1ship.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "orig_nbformat": 2,
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "colab": {
      "name": "ship.ipynb",
      "provenance": []
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EMl5OP9t3CeN"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/hikmatfarhat-ndu/CSC645/blob/master/1ship.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QzAGoawvkOyE"
      },
      "source": [
        "# Supervised Learning\n",
        "\n",
        "Machine learning tasks can be loosely grouped into three categories\n",
        "1. Supervised Learning\n",
        "1. Unsupervised Learning\n",
        "1. Reinforcement Learning\n",
        "\n",
        "In this course we will practice mostly supervised learning. In supervised learning we are given N data points.\n",
        "\n",
        "$$Data=\\{(x_1,y_1),\\ldots,(x_N,y_N)\\}$$ \n",
        "\n",
        "presumably generated (or sampled) by some __unknown__  function __y=f(x)__. Our goal basically is to __learn__ (an approximation of) __f(x)__. If we do it successfully then for any input __x__ we can compute __y=f(x)__. We type of supervised learning is __classification__. In this case __y__ belongs to a discrete set of classes _C_. For classification problem usually we assume that each pair $(x_i,y_i)$ was sampled from some __unknown__ probability distribution __p(y|x)__ and our task is to learn an approximation of __p(y|x)__.\n",
        "\n",
        "In this lecure we give the first example of **supervised learning**. We are given a set of (image,label) pairs (CIFAR10)   where each image can be in one of the __ten__ classes: ship, horse, car...etc and so each label is a value between 0 and 9 denoting the class of image. For example an image with an associated label of 8 is that of a ship. \n",
        "\n",
        "Our goal is \"learn\" a probability function __p(y|x)__ which, given an input image, it gives us the probability that that image is a ship or not (you can use others).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KWRP93VKkOyI"
      },
      "source": [
        "## Logistic Regression\n",
        "\n",
        "In this module we introduce Logistic Regression which can be regarded as the simplest neural network, a single \"neuron\". This type of network is sometime called a Perceptron, but the method used for learning is different from the way a Perceptron learns. \n",
        "\n",
        "As can be seen from the figure below the input is a vector of size _n_ and it feeds a single unit (a neuron or perceptron). To obtain the output we perform the **dot** product between the matrix **W** and the input **x** and the result is fed into some function (usually nonlinear) _f_\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "z&=\\sum_iw_i\\cdot x_i+b\\\\\n",
        "\\hat{p}(y|x)&=f(z)\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "Since $z$ depends on $w$ and $b$ so does $\\hat{p}$. The input and _f_ are known whereas _W_ and _b_ are parameters to be determined. Our goal is to find the _optimal_ _W_ and _b_ such that the output is as close as possible to the label associated with the input.\n",
        "![title](https://github.com/hikmatfarhat-ndu/CSC645/blob/master/figures/perceptron.png?raw=1)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-bYPx_7YkOyK"
      },
      "source": [
        "## Needed Packages\n",
        "We need a few packages to help us in this module.\n",
        "1. Tensorflow is needed just to retreive the data. We can always downloaded  \n",
        "    from the web and read it but it is convenient to use tensorflow.\n",
        "1.matplotlib is used for graphs. \n",
        "1. all matrix computation is done with numpy (or cupy). \n",
        "numpy is the de facto standard of matrix computation with python. cupy has almost the same interface but it uses NVIDIA cuda to perform the computation on the GPU. If you have access to a GPU ( like we have on Google colab) it is better to use cupy. If you are using Google colab go to Runtime->change runtime type and choose GPU"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4Uk5mr3L3CeS"
      },
      "source": [
        "import cupy as np \n",
        "#import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import tensorflow as tf"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mKskktTfkOyc"
      },
      "source": [
        "## The data\n",
        "\n",
        "The CIFAR10 data set contains 50000 images for trainning and 10000 for testing. Each image has dimension 32x32 pixels and each pixel has an RGB value. Also each image has an associated label from 0 to 9 describing what type of image it is. The keras function returns tensors of type numpy.ndarray.\n",
        "After reading the data we print some of its properties and plot the first ten images."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "tags": [],
        "id": "m70hk4cc3Cej",
        "outputId": "c9061792-62c7-4ce3-f4d0-95b7eff80f8e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 69
        }
      },
      "source": [
        "(img_train,label_train),(img_test,label_test)=tf.keras.datasets.cifar10.load_data()\n",
        "print(img_train.shape)\n",
        "print(label_train.shape)\n",
        "print(type(img_train))\n"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "(50000, 32, 32, 3)\n",
            "(50000, 1)\n",
            "<class 'numpy.ndarray'>\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bsDWCG2Sokl3"
      },
      "source": [
        "## Plotting the images\n",
        "\n",
        "We plot the first ten images using matplotlib where the title on top is the corresponding label. As can be seen 9--> Truck, 1-->car, 8-->ship... etc\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "StAPFRYWkOyq",
        "outputId": "3ce10781-020b-4cf8-9037-0a2a9f74f774",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 233
        }
      },
      "source": [
        "fig=plt.figure()\n",
        "fig.tight_layout()\n",
        "plt.subplots_adjust( wspace=1, hspace=1)\n",
        "\n",
        "for i in range(1,11):\n",
        "    img=img_train[i]    \n",
        "    t=fig.add_subplot(2,5,i)\n",
        "    t.set_title(str(label_train[i]))\n",
        "    t.axes.get_xaxis().set_visible(False)\n",
        "    t.axes.get_yaxis().set_visible(False)\n",
        "    plt.imshow(img)\n"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADYCAYAAADVjzngAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9WYxk2Xmg953tbrHnVplVXUuzu7lTJCWS0oys4UgjLyPZhv1i+8keY2zAGMyDAQM2MC+GMbANP3qRYQxgP42XB3lgzBiwRwsMS/JYIy5NUk32wu6uqu6qyso1tht3O5sfbnZzbbpLJLuSUn6NWx1xMyLzjz/u/c85/3ZEjJErrrjiiis+eOTTFuCKK6644i8qVwb4iiuuuOIpcWWAr7jiiiueElcG+IorrrjiKXFlgK+44oornhJXBviKK6644ilxZYCvuOKKK54SPzUDLISIQoiNEOI/fZ+v/5tCiPLifc//tOS6bFzp6f3zZ9DVf3Lx+iiE0D9t+S4LV9fU++NS6CnG+FM5gAg8/33n/iXgJaAE/gnw8ffzvj/Px5WefjxdfdfP/s2Ln/8733f+zsV5/bTlf5p6Av4e8CoQgL/xpPr983hcBj19YC4IIcQLwP8I/HvAFPhHwD/8izQzeT9c6enJEULMgL8DfPNpy3KJ+Trwt4CvPm1BLjkfqJ4+SB/wPw/8YYzxj2KMDvgvgBvAFz9AGX4WuNLTk/OfA/8VcPq0BbmsxBh/K8b4+0DztGW5zHzQevqgg3Di+x4L4JMfsAw/C1zp6X0ihPgC8Dngv3vaslxxxZPyQRrg3wO+KIT4q0KIhH7JmADFByjDzwJXenqfCCEU8N8CfzvGGJ62PFdc8aR8YAY4xvgK8G8B/w1wCOwA3wIefFAy/Cxwpacn4m8B34gx/vHTFuSKK/4sfKCBnRjjbwO/DSCEmAJ/E/jSBynDzwJXenrf/DX61cJvXDzfAj4rhPhMjPFvP0W5rrjiffGBGmAhxC8AX6O/UX4L+IcXM74rvosrPb1v/gaQfdfzf0A/cP33T0WaS8yFO0vSxxOMECIDuivXzffyQevpgw7C/ZfAgj7Pbg78ux/w3/9Z4UpP74MY4yLG+PidA+iAVYxx+bRlu4T8DlADf5k+17UG/spTlehy8sHq6aeY5NwAS+Dvvs/X/9v0RqcBPvS0k7Q/qONKTz9VXf3HF69vAPW05b/EevoLeU1dBj2Ji198xRVXXHHFB8xVM54rrrjiiqfElQG+4oorrnhKXBngK6644oqnxBOloQ1G47i1tUVwHSFEhFRokyKkQGuDEOBsR9c2eGfRSmJ0gk4MWhu01mitEVISYyCG8K4z2jnHer2maWoQEaIgxj5IGELAB4/WCikEQvTnfQhUm4aus+L/X/oPjp2dnXjnzp2fyO/yzmK7jhB8r3MhyPMCqRXfW7H8o7l37x6np6eXSk8ASoqojYQIQvxw8YQQSCGIIeC8RyqFkhIpJUJKBAJrLTEGhBRIIYlADAEhJUYrBKC1IkZw3ve/R0qkkIQYUFIgpYAIZdXSdu5S6Wo8mcbdawf9k3fCNlJd6Ewg3g3lRGK4CPDI/iOE/kYCep1B5L0jP5GIuPgb33lVHysS3/Pes+NHrFfzS6WnYZLErTxHiEiUEicVNsIsTehsS91ZnA+Ei48juNAPAkRECoExmulkfHG9OaRShBBpu46qbgCBVpLEGIiREDwxBoIPIETfO0AJhFQEJLa1zNflaYxx9/vlfSIDvLW9w7//H/5HLE8eIYViONtHpkN8sNy6fYssUTTlgqZckGjB/vYe29sHJKOMYlBQZBmZMUgkVb3BOtvfTEpRrkv+8A//iP/7D34fIR1N45D0xrppNhyfPGa2NSFNNSE4mqalrlr+33/yjZ/IF/eT5M6dO3zpSz9e3USMEUHk3it/yisvfplqtcRHGE63+cKvfJHtGzeIKKIQ9P/9aD7/+c//WPL8tDBGcutgRJamCCnx3uF9eCfqzGg8ZDIaoqXg7PiETd2yvbdDIg1aabIkJUkTDh8fEYRlMBiQ5znWWuqmZjwZMx4WNHXFbDJC64S6sZwvVhijSVJJ1zQkSpGZFGMM/+vvfOUpa+UH2dq9xt/5z/5r6qpmva7YVC0qK8iLMUYYQBBjwHtP6Bx1VZElCVFEWtfQOocUCdf2rnNtfxuj+1v/uyL8AIQQL873j9/5eQgBIUQ/sBHxMfJ3/4N//Slo4kczynL+tc9+Fo0lZDl1NmUdFb95Y8aLh6/x7eM5q01LcP3nMcbQOkcUAq0VSWqYzab8pV/6POPcUG3WbG1P6GzDvbfe5huvvEGM8NFnn+PZ688ggme1OmO5OOXs7JRsuAU6Ic0jgyJFqgFf+cbb/E//+Pfv/zB5n8gAxxgJzhMDoAQSkEQG4yGjyZAsMQyHOSLskKUaLVO8Llhbz8njEwyRnfGE/e1tnLNUdUXwIKXGWc94vEMxmFBu5oAEIYgEQrAIGUgSCQSsdXSdwzn/439jl5BIJBKI3nJ2/JhqOcc1NSFCaxKq1YLtgwMQkki/IvhZRQrJoMgZj0dEoGlbhBA0VU2eZYzHQybjAUqCtw3TrQmDwQAVNUrId1dU29MxUXmKQY7Smq5TJIlgMEiRMmASiY+e4Dq8DxityXNDkgq0MEgvkFEgfORHTA+fGrbrePPeXZq6pSwbyqojKkNejNEYogAfAy74fkZmPZlKUBpc7Oicw3vJqmqwsWVrOiPNUrTS/Uw5xt7o0q8+f6gOYuS7L7ZLqCY8gmPvIXQYqdEaTJJSLTccL0tOyw2bpl8tASTBEEP/mZVXOCJiU/PyG/e4dbDPoDAIZciUYDabsru7S72pydKMYjhiPBwy3Z7x8C2omxqRjGjkABkahA8MZSTa9zazT1YJd+E2MDpBSYWRgq3ZhN0b18inE6RSxODxtmPd1hw/OsJaiZOWB/fvIm3HL3z84+x+YRshBN4HVqsNtvM4G+ksDEczVusFUhlCCHRdQ9NVpKkmTQ11XeGcJ/j+63+vZevT5s8qVz/zhRgDbVNxcnQEIZCnCdY6XFVyfvSY688+h0guZjH098Vl1cWPQkrJZDxiNB4RBaRdSmIMy/MFe9s7GCMxiUTIyNb2tJ/ddhaDweiEGAJN2zIeFpBEtNF47zBaomQC0eJcJEkM1nXYztO2Ae8Cg4FGyohSgAdvHbbtLpbrl4umbXnt9dchSnwQWCewUaA3HSoogoAuBiweFwO5TsiEJk0UQvWG2XaB1cOapt1wcLDPbDplMBySZRnqwp0Rib07I3JhlC9mxBf/XEwN+vOXUE9RRDoVEHi0giQ1xDwhdCVN62htoCMQZMQHT+08GnkxAAeCkmjveHx2jkwznhteR2cDikQhdQrpmPt332JVNjRRcuvgBnmqqboNh6dnrG3C0ueME8VIdUQfiT55T3mfbAYcIlJqssGY7a0pz975ELvXn6G0jjcePmJd19RVRVWuWayXHJ8sGU92QXS89NUvI7uO2WDCr/yiRqmErdk2UiasliWr1Tmvv/4mWhnywRBrHW1bEy/8c8PpmBgc63KFtwLQGGNQ6s9jHDESo2exWPDw4SOa+TmZkojgiETuvvYKN57/CLP9GyAvHFlP4A++TCityLOU6D0mSxmMhqQmQUXYnk4JoaNqSqzvyIscpQWhiygJUkRssLTNBkREREVnBc45EKCUpN00FIOcEAKrdUkMAiUv4hXOUp6ViAi+DbgmEELgMubG+xDpwkVvUiHRiUKiUFKRkPRGhYANkaa1NHVDKiR5TNC691lKk+Bsx4OjEx6fnVHkBbPplN3tbWazKZPpDKM1qi8Q+M7fjpFAb5wJF+uzS6gjAIFnJ3MMUYzyFDXQPI4VeRZIZCAVEp+khETivKNpW2IQZGmKNoZiNGDvYI/96zc4LxtOy4YXBlPGowHFsOPDP/d5fqf5fb70pa+wfavho5Mdtg92SR/c48F8w8pOaDxMb++TbxlsvWZ7O3tPeZ/IAAsBXg1oleQ8pnTHc7587zHz8w2Hj48RRLQEET2dt7StpdhRLM7PUa5Dx8jifMHr9++zv7+NMZrJ1pTZ9pTZ7oST+THf/tbb7OzNeHy8wlc12EDoAkEHVKJJtaGyHiEkg8EA8ecukaOffjjb8ejhQ+4/eMTDu/cptGZrlLE1yrEu8K0Xv8qnf3nEcDJDiJ9dHQjA2w7XRZy3dL6jkQolBeVqgRAB51oa27BYrxgPhuTK0IWuXy0IgUkV1lmc7UCId/2VSgjSxECEuumQyqCMwSiNEoIkMSwXC9brNUYkGN3fhEKpp62WHyDGiI8Xs9QYicECFkFAxH4wsSKAkgyyAeW6Ym0bGu8xxpAaQ64UUho2PiKDpF3XrNcNjx+fUAxy9vcPuHPrNoM0IzUGKQTee1wURAThe6q46IN1lwwhBc/tXeO6MGzJnCAU3yzPyCcptalxxrK/d53bz95hPBpx7+5dHj44BKEgOBLpGeeaL3zu05zW8LU//qd8+/X77G9NUVLSqoxl1bGxgrfP1hyvG6Z7ko0VnC472nTArf1nuHZwDT3U1OmGL34x57f+wd/7ofI+kQGW0lA3CYenZ5T1EmINSuIaS11uEAS8710GZVVRVhsePH6dNMkYZANmeUG7WfMnf/z/cPP2be7cuc10NibPE4yRjEYJwq1oOkldtVRljd1sCLElyw2bckOWZmgtcR6stYRwWXuJvCOX+L7/8+5Srl/CveNXeyeUJi6CHg7rLGXd8Ph0QbSOxTin3ZuyFwXffvGrbO/ucecTnyQfTYhIZOwTSJDv+OfCJbxFvp94kZ0QiSJiXUfnA5O0QEuBlhqCx0qNNoKu7bBdhc4TtOkza5CS4CNZkuGsw2QZaZqijUbg2VQVzjqiNCTa9NFr5+jqjmAlWhTkwwFpmuJ9uJTXVIiRTV33A4tSaG0uLqlIEP3R2Za2CyRJJDWK4ATWdng8joCLASPku++TUvaTquhZb9as31hzPj9nkA6ZTWbMZjOGwyHaZIAgBI8P/WzYXxjly0aIgoEskFXDqjxCEPnCszfoXM1+BJXAzaHhU9MRN7d22EQ4TRKacoXraoTwFPWSmw/fIl06xqMMcXpIOD+iRfDGn77E6w8fUrU1p8cnHB8f88Jzt/nFn/8lbk5v8D/8b79PdOe8dW+D6ypCt+CTX7jznvI+mQFWiqZtOT85oVyfI2mou4qqnBOdo24bqraldQ6UYDAaIrRmsnOTnTwn6TYcPz7EaE8Eys2GZ565zsH1axR5wv7OjM987DneeHBK12qSEHGZIMQCpWE+X6O0Ic0GKK1p25bL28wpft8j8Z0n7/rWIhF/YXo18Tu5REil2bt2jSQvqNqOtunogscJSHWCcUfce+lrDGdjdm89S5qPiP4igAIE8bPilOjTv7Q2/WAuQQlJsB4hFePBiNYnhLp3wWTKYDcVxWRIFLw72x1PJ/guIFEYY0hMitYKaxuMTiiKlKYLhNDHHkTwdJ0jBklicpIkRWpNa2ucv3zB3RgjnXfECDoE+jtIIBIJ0hKEAxkQ3uNsi5YJeaJpXYcn4KPtfcQRNAqJBNG7Enzs3RdCwMnijLOw4vD4lKIYMJtO2d7ZI8vyfrATioDARfFuxsRlQiHYQbGVpJRqw6Ld0IaGna0pvzLZx8mKaRvZuX/I4PAM1zZc71pEEAipiUJBFzDffIWs7XBD08/8baSzAhMNp2VJIiJ5aKGcY+sV+/tbFGnCx778dZalp+xaqm6NXZ3y+MF7r1CfyACH4Dg6fZuTk7cpF3Nk6EAGhkPJ7v4+603HvGypfWQ4GnLtxnW29m+gkglF6Kgfb3h8eEjnGhyG4WhCtampq4Y8S0iUwK4WHIyHmGTKUXbK2WlDVTc462iajhBahExIMgWX2Bf1/TPe7zyO7xrgGD0hOpRSvQ/zu2bLSmlmW9s8+6Hneenrr7CpTlhuWgCmJqHYgUd37zK7doBJB+zcyBBSf48/ONJHtS+zIQ4hUNYdaSLQqURKhYgChWA4yBkMUmgjqdNoJKMsQxU5apCyKle01mISgxQe6ywiKqILBNGvjqKPKKkRzhKcBx+JaFKtMUVCudyQFWnvqlCCJFGXMqskAi72s3MXPTK4PqdXCrzqjSgioC6ykxQBqQSgcAKiBAj4GBEEFAoXQFwYXy/6gaxtLXgFQbKolszLc06X5+RZ0WdNaINUGiEVztmnqpMfhpGSSdtSRA9G4WRCt15h8oxndYpJAqlzpPMzJBc6lf31pgRIJNFHwnJJqj2hNoQAwUasE0SV8UmpKNMEMczZsQ2cnRALTZIqnru5z9m6Zd5YinGCWycU0/F7yvtEBrhta07PH1J3Jda16BDZ2p7wwvMH3Lx+i/nScbSoKbvIZGvKM7dvMZhtc7osYXlC11rOz+cs6gqTT/jMpz9LlhVIqYkeXN1y71svc+eFT3Dz2i7aBLrunLJybKoNMYBzAeccyjus64iXcLkIvDsTfSfZ+zuBnd742q6laRt8sIzHY5S6mL2+Y0CFJsuHfOrnPsNLX3+F9arq/ZWh4ZFcMB5kcHTK2y+/RjGcUYwmDGczAiAv5tQ+OFzwvUvjkq4UfIicrzYkiaMoMrI8RQZBYgyD3CCkBzxKQKIkuYpMZ1NWbQPOIYMnkQbnWqQALQ0iRIR1EAIg0InGtit8JxDKIDQYpfoiIRUZ5AlSy36JrtRFscLlw4VAiKHPAIoRpSSBiFa971sI+sAcEaJHxIimzycPAVyAEHufrhCyn8HGQBABL9/JdhAI+oKVaCOVXbFYn2N0QpL0Lh+jJFpKuq5+yhr5QRQRvV4R8b1bM0mwy5IQBEXXkkqH8hbpHDJAROKkJBrzboqdCBFJQA8HSCvwPuK1ROaamZCkWmO9wmU5ybokvP4mm80atbPFJE8Q2mBsIB/mxMmI8dZ77yb2RAa4rjZEAbsHN1C7grFJePZDN/nQs/vYWhBFRxM30DQonSPlAOc1be1QrSX4gPWR5bJkcnJGmuTs7e0zngzxtqErWx6+dpdUpjy/t8ON6xOC2KdqSk5Pz0iSjOEoI8kMLjiapvqRNT1Pk/j9D2J8p8YIHxyn56ccHR/hneXDH/kIw2Ff5ffOWyKSiOITn/wUv/wrDzg7PadaV3RNy1EsKYYLkmg4fO0+eT4hH09JRwNqGdARDIKqq6i6ms5aOtt98Ep4H8QIy3WJMS2dsxSdJRGa4VDStRVt52l9b0SkD4TOY0Tk/OQMH3xf3JOkdDYyGY0RXuNbh4wBSWTTehKtsF1L9BpjMozSdG2Ld44sMxS5obGOtmnw3l1KAxxixF64IPpKt4B1HU1bk6RJnxMtFUZphIh9upr3xOAASQC64OlCQCmFQ6C5qASTASc8MfZuRoRAvrMVLBCCxLY1m9qDt+AaYlf32SeXDG87muUpXkZsavBZQtZ62sUjogqEVCGFR+IwHoTMCUL3NQdS4b5r+aPG18hKhU8EdpQipinaVqTe49eB2gbs/UOq+49ww5zh83dYn83p0jHWKbqqxdYbyuS9770nMsDOOqbTGwzG1yl0wiRV7OwMmC9LHt2d0zpNEz1SBpQKeDpwDb5rkM71KWajLeZVjTTFRZiqT3yXSAbZiFv7N8iNQNEw29lle/fD7GzvUG0ajh7PmU630YnEhY48z9D6cm4UEb+TOHkRtXYXsTbBydkxL379Rf70pT+l6zpO5+f8lV/5IlnWX/XhotLIh0gxGPDFX/0iy/kc33W8/K1XWFnHtw+PmYoc0wjuvfQKLYFQaMTumGq5QgXPup6zqOZsqg2bev10FfIeaC0ZF30FmjESvEVpSZb12Qtt52hsBCHJ85QXnrvD0dExbetJiwFpluODI2LIBwW2iUgkMga8t2wWS1b1ikFR0LQgZCTQF2I479m/fkBEsKjmrDd1X558GQ1wCLRti4gKoiC4jrY5J0kis8EBWZToGBDRI5NAFJ7V4pymKpnNdsmGMzadZLVak+gc5yz9/DlA+M516iJo3MXKIxKcIAbRF/3YhrA5pV0csT4/xl1GAxw8h4tjpPeIIifbmpGPJszrNXsyIXES2Th8Y2l9hO0diuc/RBca2nJDrC0mgOw6urMS0gliWCAKjQqBuOlQH7sDyZB80dGdLaiOjti88ib+rSMGsyFubKlKz+l8RdMsuZnsvae8T5gFoUmEpKmWeJNizJDaBZoW0skA4wXaWbrgQEVaW5FmCiE6ogSZDhjM9knEEpVPITEE2RFjgwx94EkXhnyQEEPDZnFEFIJr0zG/+c/9Oi++dJ+qcXR2TWxdf7FcwpShHg8I2q5lvV6xWa8QSnBydspLr7zMN199hfv37tPWDXXb8JGPf5S93R2MSfA+UJYN6/WaW9evs3tti3/2r/8aUgfW7YJvvf4GbaV48/AYvRWpbMdZs+akXXH7c59ksSlp6pI2lthQY53FXkJ/HfRpzM/fOiBJe/+i84H58Rned+SDbZplADwSDVFQrkrOTs9oGk+SKZwT1E1LiI66rqjKllSnGCWxWBAeKQSj4ZBioBFCYoymKJK+Z0SMPDo6puks0qQopfCXMLgEEddZhlozSHKSYkgzTBHCklQdqbeMiyFZPqBLEqyz5FmBzDKS6Qwz2mE7mWDrrjfmMdLGyLquWK0WuGaNjg7jW1RwBN8Qo0PJjEBKlBqakub8iOr0Ac163qf9XTJcjBzVG8SmIfWerb1teGaPZDoiKT3mdIFvKmxVUyuJH2SYm88gFWRlSzhd4R6dEDpHIx35557DB4s7PYajYwgSHp/RxSVquk3+wi1Gz15n9a2M5d23mNSS0TVF3XUsRQuKPsD8HjyRAdbGoLSh6Tpa54naoIcS6zRRKzrXYAO9i18lOGmYDAZMxyWBBrsWOB8QUZGmGVJflBpHD1EgtQQlqdua4B0qGrx31NWcQTLiFz79ER4cHvPgUctmbak3zaWMxEIkBAsC2m7D6dljvvnSN3l8dsLZasHpcs56s0FkGuUMp+dnvPj1F7nxzAFFntN1jvn5mq5r8baiqTc41zLeSbn14X2+/eANXBU5Wm3IVMJe6BjFlu4bLSrTxK0BG98QYm90Q1B9ueUlRArBbFj0uapZX0acC6jKknv3H9C5gBCGPMvRecHx8THz8zl1p9FJAwJi6Bs3rVYrbBfoRIfRCqNhPBkihMTHCMHjoiOEjkhfkmy7jq5rUVqTGo0Q76ezxlMgRgZaM0oN3jWcrmtaJbHBIc8es1cMSZMUEQJvnZwQAqSbmkGqOV49ptErssk1Uil49NZdXDog3b1BMhqjtwrq4/uoqmRAhW02dM0Ka1uUytl0mnQ0Y2QiDkeD43IWIoPRGp0Nka1AoolOYIRh3Zb8wfEJk82aXVsyCg2di7QnJ4RX72ElxMEQnyd0d3bobMaHnn2OLp/SHR0R5zWqavG5xh4eYuctanJGO5uhhgVie8pot2B9Mme0ZbgzS7FZ4N7DE5LR1nvK+4SFGJIQY+8/ExKTZmzq3i9pm5bobB8syRKGg5TBIGNrmOFigdMW2Q3Ymo4p24CMjhjcRYDoIlFLK/LBACk6og8oITFGInyksxW5idw6GJOaA0RwrOdrvL98wSVrO05ODwkhUFUVZ+cnHJ+f8Mb9u5BovBKMtqeMZjPWxwtOz075yotf5f5bE9I0wfvIZtNibcfrr79EdB2DYUbnG6bXBjxz8xpvvXpCFRyPNytEbBEhw0fLw2+9SrM3YIXrR+sAwXma5nLOgLVSDIuMPMv6QVkq9NaEdpDz6muvYW0kMRn4SCk8gzRijILGU1clREeaGRIjKTKNGRQkaUZ2cS46x9npAh8FnQ946y96iPRLbOu6vppSaQLqYkC/hMYlRiZpCrGj9S2d9oQsI5JgV5J8Z4aaDtho2BiJQKJjSuUdi6qixpKLhFRFvG9Ylh6ZjhllQ0bjEaNrN7EPStzymGZ5wmYzp+0aiIZVI8inJdnWgNBu8F2Lc+4yDlNoo0myDNUFdJKiUCghmK9Kfvf+22zJwM9L+CgSGT2uLFm+eY+yyFhNRlR5htdQbOdc39mhbhrk2RK1LInWEqVABkG7XKE2NXa+IAwK7CglO5iQ7+8hvGUoNAejHHdjl61rO+8t7xN/Qu8RwWKUojCBYRIQWtEpBT4iRSBLBeORYndnwLW9EXFL0pUFs8KgECTZgJ2dIWmiL/xtkigiQkvyQUFwfbBKCUFiFIRA7Bw2lIzyAXpvSnSezbq8lPXoTdvw6msvQ4Q06Ze1WTFgud6we2MPpQSD8ZBEJsgmcvfum7zyysuMC0OSGqQ2IBQxWo60IzWa6XTCYJyztTviox/7EPPDNctFzbyxJLEljf3gd37vAYuFZm0EiIxgBXVV09bt01bLD0VKQZ4ZiiIjUYboA7JIGQ4L2s7RdQGiIlGKunbszWYMhtcQZw0BQZYnDIYFaaqQynOwf4vJdEaWJ0gZaMo1m/WG1vV+daMFIUqk7CcUIjF0LhJQOC/ovL2UpchSCCZ5xrpcEXGYXJPmCT4KzPY2124coNOEqqlRRYKIkA0SqBuUs+DAiQ7XbhiNcsrS09iWTVOTjAaMZnuYakV3cpdqvaStlnS27V08dV+lOBETfLXEd+1FL5bLZ4KFFCgVSVOJNH3loIueednw8tmCrUHCbp5zXSfksqPzLfb8jGUY8Ci0LNYabQw7xQ53z07IbWBclwwISKPxdQcRnLW4ck1czHHGUA9SUrfHzu41KDcEJRkmKTenQwbmvfX0ZD5gIRhmGkVCmmoOtjOevTVjVBRIJE1V0jYbYrBkheHZmxMObuygzD5tU3O6t810OOL6wXWmOwN2tickSd+GMARPlJBkKa7pkPSlpFKKPtd4nFFuNthmTaoH3LqxS5oa/pff/sc/9pf2k6Zpar7ylS+RJzkffuGjvPDhj7IoO5rmjzAqRUpPphNyk5PqlLrtODs7wW8UgyIjHRQkeUE2iExnmiLPL1K0CrQ2DH9uyOFbx3zjG6/Rto5F25HiSY3BeE9TOeKoQBWSro2cPz7FWfe01fJDeaf5EDHine1XUVIQlQShEFL2lV9SoBUcHOywtb3F1jsGOE3IixQp4fj4Ec/dvsW1/eukRYIPLavzM1f/gKMAACAASURBVB49OOR8URGiJ0kTtFKECMF7jDasyg0uKCQB27pLOagrJdmfzUgImE4QUo1ODM569p65xe3dfTZdS9V1zBKBtS07WyOkHxATha4dIVGs5xuG4wMm0hGjwzUr1mWCSvc5uPECp/deItDHVaQQEAMyeHy1YjNvka4hevduNd1lI8SAiBu0CCiVgAhUrmXdRhAGl+Ss04IzCUMCTkZidFR0rBrPRigKnbA6jnz56IQ7165xPYAsEoSHdr7ENyXBtjSbFc52eCHw6xTXbHD3HpEBQUuYThnfvo07efye8j5ZQ/Yi49f+8s+zLkuUjFzbm3H71k32dnZRUVBXJW1X46NFSsmgKBgMUpTJqLRgPYcYSj7x4RvcfuE2LjpA9BH/6PsWl0YTm77tJTEi9YVPTgk655BCEUODUQnX9sYUg/zH/Mp+8nRdx5tvvsm1nWs896xHCM3R0SmPHhwyGhYEHNXpObnOqZa9H/P27Vvc2hoxGhf4AGVZMZlK9m8MqDZNv3IIglQaRttj/pkv/hJn5wsO3z6ibjrOmhYjS3ZGI6SIbGdDiq1tFpuao4cPLm2+tLWWk5MzsmxDkSTkWhKFx3eatCgIm46iyNmaFAyywO1b+6SJQSc5USiSRJEkum+usxniqwZb1QyKlOFwCLbhxsEBxqxZV2sQHq00Eol3Dql6N4VHY31AiUBinnxh+NPGaMWv/fLnqNcli+Wcpmv6lqyt4+a1GwghWNUVI6UpXUdVbbg+2yYQUQjSTUVQkhETVAjMBgndqma1rCibCtt4Rns32HnusxAqzg+hPj+hrTcQPcNUYURLvMiOcF24jOMUgogJjpEQpDIShUM5S90Gdmcztq7vgJKctA0RSZIZhFO46BjkKUmSI3ygPDuBec3poqbKU06JbDUNYnEMtUN6SecamtBhBSA78kawOnxESp+zLYYjRnXDjGffU94nutKKPOMv/cKngIhQkiTP+jFQCKQyjPPtfpOjdzrNh4B3HpzDdg1pprnz7A0yk9NUq77sT/RNxWPwhBgJom/6YdveB2pQSCWQSMrzmkcPDrl5c5+da9tY3yffXzaCD1RNjdQKT+B8ec6Dh28zGY+IdQfOUp6uOC1rTs8WCNnxr/7L/yLKVrz67Vd445U3mD+YMx0bjl4XXLt2HWgIzFGpYTrb5ZMf/gTuNxL+93/0f/Lo/gPWrSfaDVEntDawPp8zSxMckBYpSl/ObJEQIuuqprUOOSwY6JzgHa2zXD844O7dt4kxUBQJ+9dG7OxMeyOZgpSK5KJNqZQCGkdX1yzPzgixYzDOSLQgz1NGo3hRQOQRUZCbDKETnLMM0wwHdN4jo7mUlXCKyM9/4oBxMaRrHXXncFFgXcQ3LW3b0VpLaz1N59lUDcZoVpsNEzOkrPt88DDb4/DoMXfffsT1tGAWBfP1KcJFQjaiuP5JPnP7Dm51xOmD+7z12iucHr/BQCzA1jQerFd93xFx+QZ1GSRzOyZzFYOuodjM4UThywUvPP8MO3duwWrD6viELZGAjuhEktVN/34tEUpy7/E5WxvJ7etjVqmiXC7xD89IqznCRYSXdDJgJVghcM5y7kryfETpLE3bETYd2eMT5MFPyAcspWSQpZhE94ZWKQK823kqftd2KMSIkBLHO530BcokFKMxIYIPvaWOxIvE8T6a7ZXujXpwCO+QPpBEhfaSvJVw3lGqBaMkY7w9upS90EIMNG3N4+NDvvaNF3n51Vf56pe+jG87zh49Zn1yhq8b8IEueMzOkBe/8lVCveLtRw94/PCU5syyPItMtlPmp4G2dUgBWW4IvMmX/uQVRFKgTEo6nNDGisoFHpctIYWsbAiLJWmRMt0a99/ZZUQIus6ihMI7j/MWbQTCRoJrqdcl2nWIOOHG/jZZKhhkA0bjtG8MEz0QkFKwNRtyft6y3lQcnp9hfcPB/i5n5yVHx+e0riI1kkwrjAgkSYpTEtt1yOAxylAk+hJ6Nnt3yfHDN8hv3mI43GJrewuZ5ng0m/k56/W6r1IzKXXrqOuOuq77RkQhYp2naRratmWWaHTnuH7zNulgwrzqOFq0hCRBdBaRjBhe22J48HH2P/3X6OZHrN/8Go/vvszhoweEdo40LdjyaavlB/Ah8pXSc8tEPk7gWVuTnp2Q2g0f/fgvsnvzDg9ee4ujtx6xdi1BOxIUOkRkWxI3DVlRcHM8w7k1Ojc88/EP4cuKjY+U3z6nC5GgEjocUqfobEBMMxqjiFtTOiTLsmZdrnGriuUbd99T3iduxiN1gicSvCc4h3WOrunwwWGd77tJ+X5rFO89Tdv0+8R1LV5ETJqQZDmDrK9Iiu8kg8eAJJDlCTY6XNfgXAsxQ0QB3jFIDduTEUKAbRuUGF3KpHkpJdkgpXUdjw8f4RrH0VtvIZGkQuGrltg58AGhBNvTLYzQ1A7ydMh4GvFtRVWV5OSclo62bojWoglEIejUgjYIqtbhoyCoBIunipHgIkYopEnJigGZoO83cQnRSrEz3SZL+tJjlUoiHq0imfZsjwqKVDPKU3a3RhRGMchTnMxxMWKtxXYtSZYw2pbMVy3LZcOj0wXLxYLTs5r5+ZpN1fdlLbKk763sHCL2VWNGK0TvDSNcTCQuG1IIqtUZy7OU6CNDFHma9WXUftjv+iEFaZIQYoGLkq5tefToECk1aZYhgLqquTkd8fwzNzBpQVAJZec5X1vmm8h81VAvViyJtMpgTYrJ98mf/2Vu7T7P6Og+xw9e5+jBa1SHbzxttfwAMXjeXi1YKg8xJdeG67Ijl44b+3tMp9scm8e01tI1FqM9FosSGuE8CIf1AS3AKljXa+z8DBUg5Ckbk1MVCTYraGyHTFNknlHFyCZGNIHFuqVxkRaF7zrS5fI95X0iA9x1lq9/89WLtncVzrZ456mrmnJTYi86WGmToLWCEKnKkrOzE+qmRieG8daULC+YjLaZTibkRd+pKs8UKjqKzNAl/Y0RQt+sRggBUjDeKnBumy70qUdKgVKXb74ilWQ8m+AaR1fWbM5WGGsRGELs+1hFqRBRopVERMHyfIW0AeE1QiR4bWmFZt3RD2wWhI0o6/BE1l1FZT06yftCJqXeXXEE4bE+4iMIqTGJubS7ZQghGQ9HpBf7s2kd6bqGGCJFlnB9f49US8ajIUWeYrRCK0lUF1kMUgEKk2iUESgJ58sVj08WlKuSugpUmxpCy+7uEKMNUkCgb95Td+7dPrc+RPw7uz1cMpSUxK6lXJxTbVr8oxNqB9OtXXZ2tjEyoFTExQapUpRIyE1kOkwwJiXPC4zW2DbHjXNuHuzRWE/jIpsuUKRrlF8ykprVGuabmvNywcoGjtUUlQ0YFDvkzyQcFDk6EZwc/9Btzp4qiYwc5J55Zbm/EgyFJptK8lQwLUZkss8QccHRdpZU9A2rjLzYHNh76k2HiALftKyPTqhch1EauW5o04JHSrKQAiclWkRkdLQEOh8JqyWLTSTonCAkSkTSH7FpxBMZ4LLc8H/83h9SNSXl8oxqvSTRmhgij4+OcM5TjEZMt7cYTcZIH9EucPj4IetyjTSGydYUHyJaZ2xvz9jZnXHjxj63nrnGTqHJtSKORgQpkErhY5/apnRgsDVAJBLnPMIIlBGXdgY8nI1xq4aq7JBtx14+BGlYbCpaFYhSARIRYDEveePVe2zlBVVbU7YdPoIVivN103cIQ6Jl37Gr7VpW64ogFMr0+6pJ3TdjEfS7GVjnqKqWNLcMlLmUAZN3EEohlUJIQcRjfR/gybOcbHeAlpDlOf2uvhofAlHaPpimgSjQOhI8aCU5Pz/n+Owc7yA4gbN9PwiiwgeJl31P27ZtWZY1Uup+1RY8LlzOLAhBZLNY4VvLqjrkweMzvn33ATeffY5PffpT3L6+xyiVffRfZwiZIKQkTzSpERjpUTGgdYDCkGeaxkZq6wmLktXpA+aHR+xfu8HO/oiyNhyenhIeHnNUzTmXQ4ZFxt7QMNva5bb5OC995f962mr5AXIj+NU7Q7729prTleO185qxyUiTFOPAlhW2a/DRU9uOYWLwFztrewR1jJTWEjpPumnowhwW64tuaZI2H3LfNjxoG5IIqVckscMkGuEiVVlTxxwhEySRoVLspO+dKPBEBnhTN7x67wghHMcP3+bx23e5cW2fra0ZpyfndNayLQWD7QlWBJanJ/zq53+Jj3/sw3TO4mPfDu/B4SNef+Mu9+6/QZpoPvKRFzDic+zcuY4KsD3dwUqBSN/Z9hMcAak8emAwJCADQYZLGTCJoj+UlBggj5LxYIIVgsWmQmYpJCnRQz0v6dY11aJkESSt6zDjIVu3bhCQLOcr0kHap7tEegMcLa2LCBnRUiO0RniPF6Hfssf1wZHNpgW5xtpwSctr+zS0qmn7LIYY8L7Fh0CeDfrdG+oKgiOtYF01dA4CDqP67efBE30HHmzdkmpJtSk5m88J0dBpS3Kxh2HbRjZVR6sCWmtWreVkuSJEedG9LiLl9248eVlomoY/+eMvM5lOWZQVr99/wKOjM0bTKb/7e7/LF37uY9zYGTNIFUlaoE1B01bMJlsU2QCJfLd4QiiJRYIyzNcb/vRbr/L3//7/zPnpgs989nP8lV//TW7dfobt8YDCNux4wcvLgJeesN7QJoKd2S5Jmj5ttfwAqRb89U/tsp8rvvrmitdOGr702PKpZw7YPDqkMRnNakHwjtI27MW+xN1FTxk9SylZ5QWN9oxQFOMhIkDYNIimQaUFJ13L0kemWlMYQ57n5GmKkC1dV2NEB92ciY98CMWgeu8U0CcywLPZFr/xr/wbdF3J1/7pH1AtV1y7SIHJ0hFC1ewe7PL8x59jOJ1Q727xL/z6X2U4HmCDxwZHoL/RGtdwdrbk6PCUrvPkJuXxgzPeevkNRNNweD5n++Yuz33iWfaubeGCR6YJwgQIfSoMwl/GewXvA4vFmrSTJOmAyZZk/uiMNx884tQ2DA/2SfIc23mqEPEWfN1Rd33rSCNq4vE5m01NdJHc5DQh4qUHJfGdJSkKYohY7y92NQhY75AmwaQJRV6Q5QOMSbH2cu70C/D/sfemsbJl133fb+29z1DTnd7cr9ns5ixSHiJIgjLJMi1/CGwEUBIgn+IpAxDDyGAkhpEvQaAERhJkMhIjhiUEcWAoCRzHlhMrEhDbgS1Z1ECatCg22exmj2+8c03nnD2sfNh1XzfFborPYneV5P0D6r1X99Wts2vXOf+z9lprrxVixMfIuu+BiJjE3v4ei/maENYMvcei2Hnka6+8hpB9nXfu3mQ2s2iCbhkJ/UD0HdWmo/Kbbz3gYum5e3iDo9k+Vqas1jlmoQyIMyzWHefrjqS5bq4BHHEnu22HmPj6wxPq03MWyxW9j/zwZ38f/8Gf+Q/5iZ/8n/jC53+Vb9RwY39MVTvG0xkxRvamexzMDhmPJlhrqesaEcsqRrwRPv/lr/C3fvbv8NWXXqKpan7m5/4v2sMbfN/v/X72moZR0zDrAncm4K2hXwuWhNVEXb13s8ltkTSxVyc+fXcKKni94MXHC776+Jzn12ui2XSPjrDsBxgGXONy0kDSzRw1LBi4fHaP/U98lNpVnL52j+OvvsLNfuCZg+vI4HEmEWzkvFuwXC25OZ5w++gwd2NedcwWS+7Gnrv736V6wCLw+qsvc7k449Gjh2hKDP3AarUCoK4sRiP94pIYeo7fvMff/rt/m1Xfs1wvSSijyYTpbI/RtObBw4eMmhn7s2v49cAv/sLPc/r1L5MGz1unZ1TfeJ2vvnGfZ5+7y3Q6Ym9/RtM2jFtH0xhc7VDdQQVOQF8T1bKWRGcjD8XyICSWQyJcLKmGSEqwTrlkehcCmmLedOADi9OzjRAIp5cXYCSn46lSjUZMqyrnSgNqBGcES42xhqqqMHWdS+ylhLGb3MAdxPvA8dk5hsR43HBwMGO17pkvVoSoCAZrEmeLNYvlKVYsj+4/4mBvxnQ2Y+gHHj58iBjhD372+6h1ynhUM9ubsPArzi4uUK84Z5gvGmwV6P2CPnrEVXQ+gXGoQgwDaQdr3AJUdY2Z7jGkgVTXjEYVN5+5jYpy5+YtPv8Lf49XF6fcOBjTNjXNqEVVaFzFqBnRtnkZXlc1TdWS6op56Pniiy9xfHrCj/zI7+PTn/k0/9v//n/w8te+yLW9lk9+5OPYRjh9dMqXX3kF09ZM2zGHkwnJTXb0jBIkJQ5Gjo/dGOOHSBcH7s1XrA1cu3GDyfkKbE0XEt1igfMR5xpGIqShJy6WTGOO31wGZbo3w10umSnU/Zo7cUKNoG2D1AZRxfaRG3ZEK2B8wq962n7JwVh54dnvUlPOFCNf/Nzf5d6D+1ycPcYvL4l9R0yJYRhAIn3oefTwMU3bcvfOs6zUcbZa8ejRORcXZ7mNULKcXz7izbde5c7tO3z/9/0gH7n7PL/2pS/RH5+yGgbmfsDPPS8f97S/9oBJrTRtjWtqRm3FbDbm1p1nCHH3TgNBMGrxPtF1PfP1movgCdaQjGHwgajdpstBRMTSTMY5nc9YxOVOB7ZWjLGYyiEm7wq8KmkpxpDcpouvACb/zBjJgSkxuXttisSYl9e7iKJcLpeM6goXLKt1n33/GNZ9v+mSFzDiGddwcr7ixa+9xahpEbHEEBh8T1VbXn3jAUfNIaKRvdmIVTCsTi5RUc4vLjg6mtAYwUelGyL92hMiJB2o6ho0EWJiJ701kqtXGFtRteCMY90PnJ1fMF+seHx8zvL8hOVyRV01VxUmqZ2jrRz1pn9e2zQ0zQh1hpPlnLc2rYd+/2c/yw/8wA/w+utv8MV/9GVe+vqv4YynX3ZcnpzRnT2gGrcMfsxlX7O+tHi/m9vb0YTTxLVG+fhhxelyxKL3pMEzGU+ZTPcx9YhVVMIwkELE2wHEUgFTFKfKsJxjQ8CcXTBSZSSJmkTbL9iTmtXaQ6oxtoIgmMslXegJQYn9grZWbtze5+6dg/cc6lNvRbbRE7oVoonReAzG4CqHcRZjHcbV+NgwrQ544WPfy961u6g2PHr4mMfHZzw+PqZpp9x/8Cb3H7zB2fk5+3sHzEYz1qbmwuwj04raGipnSSir1ZL52QldvyKJwVjDeDrlw8+zkwKsqnTrnr4bWK97lus1QZV6PMrtT6zFGEHF0NgK69ym0eJVMMpsChTlrr8i2d+Zmx9r3uASs7sixkhCURGMNU/SzZxzOGff3s69q/aKGExVYes67wBcr3FVTfCJwYfNjWeTFpTg7HLN2WVPnFV5A4ZrsabCJ8+9h6fIniFFT1tbptMGupbRpMV3HZfzOft2hKsbTEgkP5BU6IeBum1omjZXxdr2pLwLuXVThxXQlHfsvfnWfd66/4DX33qLi+WKbojYTpHec3m5IIRA5Sy1yx0sjDXUdUVTN0SBy/WSzkfufuh5nnvuOQ4PDnjuuQ/xuV/5ZV555WsMqznDuuPRyYLTi0uqUcOiqqnrmnbcEsJubm9P5DlqRLk5dnz8+oz7FwvSYo6EgHMO09QM5C4hEn0+uTC5ZYIIxkAdPKwSGvp8bUFOl10pYxqiBoYqp+olFaIEVBJ9TKADN0Y1BzcPmI2/S0E4YwzP3f0Qrm7ofYczuWCOWIOmhHMNYitc3XBwcMAzd1+gqSrG9Zg4BM5Oz3jw+DG3744JxiJ1w8VyyWv37nHtxrPUN+5izITpdEZdVVSVw/s1Jw9eYxVWzBcLlsNAkorFMGDaBYPfPX9dTInlckm37um6gX7wWGsYT8fYKgukSLZUxTqsy2liVeVwziHGvN3+e1PU/UqIrwT4qiO0sZak6ZsEWFVxrqJyub17tpx3VICN5BtQVaEkOu9xmouAJ9kUzXF5BRA0Ml/2RCzVaMx0f4+6rvA+sO6WnM879s0Ka4TJqGatCWYBaw2VGbFcLWjHlr3xDGsTKQU05TY8TVMxmYxYi2Fzp9spUkosuz4Xpoo5b75/9TV+9Qtf5CsvfpWL5ZoUQIZIipHTy1Uus2kMzuQOF4Jgq5zGl1CG4KnqlrsYKudYLhbUTcMwDDw8Oeb85GRz4+rojMO0Nc5Y2nbE/tHhTnaPRgRxDRoGjElMG8vzN0a8eerx8zO6+QU+DCQj9AheBZsiUSMShYTBi2FwOS+CqIhPhGTyaixGkuS6yesQWTpLV1fExiFNrjUdoqetldE0r9QlvPc8PZUAxxj5yAsf5fmPfQzr8j5853LFoZgiRgWNOZcuxIHlyX0Wq4Hzs0vuv3WP4+Nj5ssFhylBVWPrMX5Y8eD4lK98/Rs886GPc2N2nUYsravR6Lk4f8y6GxiNJ0xDT3e5oh0fUTUzlp2w2mwh3CVijKzXa4YhPBFK1zpcW2HXOd0q57Dmi11VCDEgVqhql9upi8kFajYnuWx6xcmmbX3Xddn61bTJK5YnIq2qm7/ZFCDf3TzglJT5cklVOZq6wojQ+8BkPAMTEaNYp1gTkRDo1hHX1LR7I0xjciudKlHZmmhqVt3A/t6MW9UIf7LEr9asFnM+8txHODt5gA+eEBOo0K16VCxtW9M2NdbKprnlDs6VZPfDEAb6fsD3PWcXc/7aX/8bvPXaPZbLFSkpF2uPppQDTZCrF0ZhIx8wJFTCkwJIo2RZrdacnJwwalsuLy7p+4GLywXLyxUScyfl1DhEc8pjZR1Wqp2sxS1iaKZ7BOlQ4xGFa67mozeXvNhdcHn6iMv5wHJYs9JEj1Apm56JhqiwFmGl4NSgXkk+Fx8SEUjK2uUV6DqBNyn3YZREYypMVCYEbk1rnt1rmEiku3jvbjRPaQELGgdeffUNhtBzeDjj6PAA7z0XF5ekrsPEwLh13HzmBrf3J9x7/Ih+2dO2LR967sPc1LuYespytebGjTucnx4zDMr5+Zwb15egymoIJFPl/mlDRz0aUwP+9BjFMt2/xt7BLUC5d++rv5Xv630hb8VOWGOpKrBWqds2f4mWTe6y5C3Zmq0bayy2spjK4CqLNdmStSk9sYZz25jswtjb2yOEwOBzet9VLPLKLRFjQCNoMqgmUtq9lQLkTStK3mziKoexloqc311VBrFkAZaIquCTR6yhnVYggRQTtrZcnF7i3JQ9FVxdM26E2bhCru1xks4ZjyrkaI+kkdVyyRASRmCyt8dk2iIauLw45+RssZNZENY5Dq9dy7W3+6t6vMLF+QUHh0fM9g43bawSST0xDMQQCCESg6IxnxfDMGzOp4SQC1w9fPyYX/6VX2Y0HvG1l77OYrEihJyZZ9SSxBBiIvqAhlyIZz7v8DtaYU9su7k5e/oYcZXy7K0xr90LpLgmBM+QPBfBcyKWmbE5f16EQYWLpDz2misyYjYrUJ48f6yBSxKrBGMMRyjTIWC7gRu24dO3r/HcrTG39gzNekEf3ttX/pQ+YPDDnBdf/Ie89torWBKTtsH7SNetsSoc7O/xyU9/gk/+Ux/nhbu3uLh3n7M3LhhS4OBgn2Y04ux8xWc+/mk+9onn+ZXP/SL33ryPTZZh1eXmiSES6wZXVVRty7Mfeo6Tey+DsTQjx91nn+HOMx9Go+frL753u49t4Zzj2vUbuYtz3FwUkstUjswI2bTCSUnxKZGSUpG3C4vkpp2q8Td0ZsjvE0Out2GdJYSQVxsp+4GNMU8sYGPMkzoZV/7iXcSIMJlNSCh931M5S1NXCJEq+7hoRyMqZ+mXCwY34JrEEPpNqxeL72FYJx52J+zduoXxFwyamIgynlmu793g7Pw+k6bFuYYI2BS4cec2iqPrAv0Q0GQ4nO1hzcNtT8u3oChBA+IMjZ0wFou1loODaxAiKekTF1UMHSn6XBYg5U0tKSWCD6zWK/q+JwSfb9QhElPgb/3Mz/DrL77I57/wD7lcrFBx2YImt7zXmNM/I9kSNMZtGn7uGsKw6pCUu4wLoD4yGY+5NvOcLeasHs2Zn19waoTPpcCBwBRHu2lN3xtl7pV+Ew0wIjgxVCLkDd15h1trc7Au+og3MCIxmwRsWtIvEgsaqrFF6+9SHnBKCUX5+Cc/ya07tyHm1vQpJJIIxuY95wfXD3l8ucav3+C48wTnOLm85MGDe/jBc+P6M3zmuRcQze3Ba2vQGOi6ARFLEqVPWXAmVcv1GzfoF6e88FHHSy9/g8X5McdWSCHsZF8qawzj8QhVQ4w5mLTq17jKYpxBU/bjkhSXDFFz/YzcLnKzVEzyZEfW1QW0MYJztbjOE0LIFwf5V574jVVzyUZrsMZgJLs8dhERmIxGm7rAihPFGkjJk1K24BHwRuhXK4REXdt8MxoGQhcZukhlW0Z7U6RpCF1PSgHncocEdY7e1zTWMR6NwVVcdEuMVdarJX0X8829bkEEu4sl9sgrHzGCJEHIriUU1IR3fMeKRpMzOkIgbYpjZQ+WcKgHDIPPq6KY8jzGwDB4Hjx8yPXrNxhNDuh7T9h0YY6qRFLOld0cxxjD2cmjLc/Ju6Mpf17jDKIOhiyAhyPHq/fu8fjNCxbnc7xYjhGWwTNKiRahNtnNEpzLmUUiGJM7RXuUkDY7NUlYAZMUDUqyglhFCSy7BRUBkzyOEZLeuw7LU+YBC+2o4u6zd7kRbqLB4xJYNWhd4+oKY8AYiF5ZI7SHN7kWaqaPH2OOj4l9QA0sV5fYSnGVZTJpQZVhWOF9xA9droZmLQf+iIOjGefLFfPTY9bLBWeP7pOGjrYd7WwgwBibRVUjMXlSiojkgEgkkNKmXGfKAipm4+GVtzsjX6G5XNGmfRP5mWgOYGGy+2GjGTEqyFUGBZvgy4622SGfU7XNrYCMKMZACB6RTSH2ylFVLgd7m4YQA9euHWEah2DBGRiRF4eNAiYH9cQgtcUYQ7SGvT1HhaOpRyRrqZNHxFO5CkY1Yh2KEEPY0ZkCwTzx81/FAqy1INXbgVZVsAYj4GIixIhuDKeruIBxMTdAUH3SvTvGhG78oN0QGXwghNwyTI1sVmwpF2nfrLIWtDjDVwAAIABJREFUl6dbnpF3xzibzZiN0aExksLAUVsx0h7Tr5C+o1EYgIjQSbb0OyVnO0hC0ialNOdG5GqPGlE27gkxVEaoxdAADYpNSgqBoe/ojLKuoDHvvWHlKWsURqqqZ29UM/ia9TIxLHqMaTC2oh61NE4wGvHe07Yj2nbMtPPsHcy4dv2I6XSMrYQHD+8R9fome0IYhoHF4pTFYsWwXhKGDjGWrl8x259xMV+wuJwzGY0YNzXjpuZgbx9rdrPKl2rMu7zCgPc93nsg5/HGTd2BfMLnbgeqbE6ajeGrkn1wCBiwRklRSWT/pDEG4+zmpNi4J2LeIWiNyXEkVTTFXc1AAzZ+NSVbExsXzBD7nAZU11TOYa3BiuCamhAD08kIbwWjBmdqaps7SfdxTYwBZ3NBKLuxlE1VgQMTBBELCM46EKWuc82RXLIx5ADdLqIQ41W5bcWQb+DOuXwTElARRNnkg1tsLYgPT1IaYaPPKQtzSm9btZWAsy6/51VVw6RPKsNdBXavaq/EGLF2B0ucimCranMR5XmRWCFdx6yGO3s1d2c1jy4Mug4sNbdO3NyHckBc8vVrVRGEtDGMrkrnGsCJpTaWkTFMjDARZWKhVsEkzdkSweMHg/82XtKnmkFNHuIjrG3RIXH26Jyvv/gqlopqlgX25vV99qfj3EpGW0JUFosTUlwznVRMxo6T01MenzwmpUC3XLJcLFiv5vT9wHKxJPQdkjzW1YgVHj08onYVt28/w/WjfQ729tnbO2A62+fnf+UXfytf1/uD5k4PYeMiCcFvNgyE/CVuvlBjLa6yYAwpxifuA8gXkWyCUVe+3RgCBMkBmN9wMcQhEjat5+u6yiJMdnPIlSDvIEIWhCf+603pzJxK5zY3qx5roDGGUdsgAtY4dBMrcxuLZz0kfEyINVTGYozNFa/aMX3vGfoBI5obcKpibD5m8IHFuqMfcr+vXTSBVcEPm8JU5m0rNKlulsqQo0UgmsNFlaswm6aUV1k0bGy5GOPGRZEbIcQQ8SmvoHxKJNmYgiq5xvdG7N9Z1tTsqKvGuLyp5uqk15Dz7lsCH705ZYiK1cjXHl7weB1YR+hUiBjSRnBjykFakavV6abu+cZgGBlLaxytFcYGJiZx4AytESqJOJO/qxQD/bfZr/J0tzBVJPbYIDSADAu+9PnP8ejxCbia7/meT/K9n/4kd27fYLlc8uV1x6rvee3+A9588ICzszOGfgCEZnadxWJJt5yzXl4ydOt8MRrYm465ce0m+/tH7B3e5Nr1W9y++RkOZ2Mqa7Nfc1NNzOygBaybOrVxY2UIWVTqzVjtRlTZCE6IeQtySin7r0yu+vXOz6abnN9qswTt+35jheQ0wKu/U8q+v6qpqVxNZa8SkHbxYskntTMb6wKwVcW0ck8srvVygfd9fk3TYGuHT4Z+iGjKeehhGJiMxk+0s48ePDSidN1AMB2n8yWr8wWVazBVxbxb0LQGsFxedCzXPUmVpmmedOneNVSz/5+UwAh1Xec0xyhv55cjWDUMIebiTWxEZLO54MrFZSubm79uuArU+uBzD7iUz90UE0bzMv6dBsLmnT+4D/+dYgRc9cSyz6XEA0oOGh7UymduTzmoDDdr5RuP5jxaBS4DdMkwqMGrEiQ3HZVNbr2IIKpUCazCxFhqMVSiNCYxtYl9J4yN0FaWyuaViveB7tt0DpGniY6LyGNg14qAflhVb2x7EO+kzNN3Tpmr74wyT98ZOzpP8B5z9VQCXCgUCoXvHruZm1QoFAr/BFAEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgSRYALhUJhSxQBLhQKhS1RBLhQKBS2RBHgQqFQ2BJFgAuFQmFLFAEuFAqFLVEEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgSRYALhUJhSxQBLhQKhS1RBLhQKBS2RBHgQqFQ2BJFgAuFQmFLFAEuFAqFLVEEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgSRYALhUJhSxQBLhQKhS1RBLhQKBS2RBHgQqFQ2BJFgAuFQmFLFAEuFAqFLVEEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgSRYALhUJhSxQBLhQKhS1RBLhQKBS2RBHgQqFQ2BJFgAuFQmFLFAEuFAqFLVEEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgSRYALhUJhSxQBLhQKhS1RBLhQKBS2RBHgQqFQ2BJFgAuFQmFLFAEuFAqFLVEEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgSRYALhUJhSxQBLhQKhS1RBLhQKBS2RBHgQqFQ2BJFgAuFQmFLFAEuFAqFLVEEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgSRYALhUJhSxQBLhQKhS1RBLhQKBS2RBHgQqFQ2BJFgAuFQmFLFAEuFAqFLVEEuFAoFLZEEeBCoVDYEkWAC4VCYUsUAS4UCoUtUQS4UCgUtkQR4EKhUNgS75sAi4iKyFJE/rPv8PX/uogsNr/3sfdrXLtGmafvjH+MefrRzTwlEfnR93t8u0Q5p74z/jHm6T/ZvF5FxH1XBqGq78sDUOBj73j+CeBvAI+BU+BngU/+Zr/3O/3xLvP0zwOL3/BQ4F8u8/TNnxf4LPB54BJ4Bfi33uX3XgV+dNvj3+ZcAdeBnwdOgHPgHwD/7Hcyx7+TH9/u8wJ/ZPP//8Zv+Pnzm5+778YYPkgXxAHw08AngVvAL5EFufAOVPXvqer06gH8YbII/z9bHtpOISIV8H8CfxHYB/5V4L8Wkd+z1YHtJgvgTwA3gEPgPwf+5nfNivsdhogcAv8R8OX3+1gfmACr6i+p6k+q6qmqeuC/AT4pItc+qDH8NuWPAn9VVZfbHsiOcQTsAf+LZn4Z+Arw6e0Oa/dQ1U5Vv6qqCRAgkoX4aLsj21n+HPDngeP3+0DbDML9MPBAVU+2OIadRkQmwL8C/M/bHsuuoaoPgZ8C/riIWBH5p4EPA39/uyPbXUTkS0BHXon+hKo+2vKQdg4R+UHg+4H/8YM43laWICLyLPA/AH96G8f/bcS/RL4L/3/bHsiO8lPATwD/3eb5v62qb2xxPDuNqv5uEWmBHwPqbY9n1xARC/wF4E+pahKR9/2YH7gFLCI3gJ8D/oKq/tQHffzfZvxR4C/rxvtfeBsR+RTwv5KDJTXwGeDPiMgf2urAdpyNO+KngD9b/OXfwp8EvqSqv/hBHfADFeCNc/vngJ9W1e8o9eOfVETkQ8CPAH95y0PZVb4X+Jqq/qyqJlX9KvB/A//Clsf124UK+Mi2B7Fj/AHgx0TkgYg8AP4Z4L8Skf/+/TrgB+aCEJE9curZz6vqn/2gjvvbmH8N+AVVfXnbA9lRvgB8XEQ+C/wdspj8YeC/2OqodhAR+SHytf5LgAX+HXIm0ue2Oa4d5I8B7Tue/zXgrwI/+X4d8IP0Af8Y8APAZ0Tkj73j559W1dc/wHH8duGPAP/ltgexq6jqyyLyJ8jR6g8DF8BfIfuEC99MQ56njwAe+EfAH1LVe1sd1Y6hqufvfC4iA3Cpqhfv50HfryTnjnxR/Ph3+Po/Tk4S74CPbDtJ+4N6lHl63+bpD2zmaQ38/m2Pf8fnqpxT39nr/+PN6zvAfjfGIJs3LhQKhcIHTCnGUygUCluiCHChUChsiSLAhUKhsCWeKgtidtjo9NAxdIBarK0xxmAMWHFY47DWYq0Bk4h0qOlRibl+kIIgGAP5X4IYwQiIGDQ5NBlijPghkBQEiyL4EDAGUkwkVVBFSVyerOlX/v3fsvIUGGvV2RoVMEYwzmCtIfgAGEQEERCBylnatiWmyDAMAFhrcc7ijMGIwVjJ06f5s6eYEIQQAt4HVECM5F3+KGLMkzk1mwPNzxasV91OzRPAZP9ID249S65N9XY8Qsif50mMQvMJZJ7MnaDkB4Ci+RzjN/mIV9OEsjmNvoXzR2+yujjdqbmaTCd6dHiQrzFnMcai5O9YRL5pnlIMDH1PjJGUIiFEYvCkFDfXn4DqZh5Nfpj8MNZgrMUah7F2c31bjN28RgQQRIT79+9xfna2U/PkKqdNmzf5yeZPMbK55oR8muRzRYwgmq9DyH8bI4gYYohZZ96JSD7/zGa+N6dciilfh/BkfrI+5XMspYQfhmNVvfEt432aD3d4u+VP/re/h0dvWXS4xWz8LNYo89PH9BfC0f4tjq7vc3jDMb6+oB99Fd17QB8j/dwiXcvUjhnPYDRqGdXXESNYE3BAWh1x+vAul29d8LWXfp219owOp0wPp5ycXpJCz3o5ECMMvme1vuRn/uKvPu139L5TuZo7dz7CoAk3bji4ucezdw659/p9fG9pqpa6tuzvt7zw0Vt8z+/6FBeLOW+8+SY+Bg4OJnz4uVs8c+saaECcglFUQJMSek/qI6uLJYv5iigG09SYxhCdx00srjJU1lEZh4jw43/qA9na/tQc3HqWf/PP/zSQxTZFJSWwAkYgCZAiJiWsKLapSRqprEWMBWMBIWgipbcvuvyGm4OIkDRLtbm64WsihIhPPDn2FX/p3/0XP5DP/jTcunWbf/9P/3uM2xGHB4fs7x9gqhofE8451us1KSUq52ibltD3PLx/j9XlBUPfsVpecnFxSmUds/GYt954hde/8RKX5+dIzDd9qRxV29KMRtRVQ101NE1LPRoxmY6Z7B3QTg5pJ3uMx2N+/M/9p9uelm+hbhyf+r3PoyFixFBVDVXTMp1OGY8neO/p+g7VyHg0onYVDiEMHZVNVC7fkObzJWCIYlAjIAZjDM45rAgxelKCFGG16jk+PqXrFTEOYwwhBEIIxBjp+543X33ttXcb71MJsPcRcTU3njkkdjPC4PGhIxnPydkly/XAYBaYvRbDgkHPWM5PCLFCuzF1UKIIhobJ+BCjMy4vIykEGgdOLb20LDs4fljhbYeZRPZrpR3XzM9XqCohJPrOs1p27GIWhxhoRwaXDLZyNFLRupqRq4mLAY09JIuMLVaUbrXA92v29qZUtePw2ozbd69x/fYRIawwLmFttlxEhdgH4jow7LfE/gDUINYhrcFMBB2D2pTv+hsxqmq75Vl5d5R8U7kyo7Jlm4iqG/3MVko2NxKaIhHBkS1hRDfGsT55v7wQeIdh9mTFBCllk0dUSR/MR/yu4Jxj//Aafdex6tY0TUW7sf8dBicJnwKahLoZM5sdoArztkVTIMabDH6grsfsTWYc3nqGIIbVV36N0wf3MSIYY7CuwlmLIe/YMNYizlHVjmY8pR7t07RTRqMxy/l8y7PyLohkq1Pz+iineyWcs1SVJaWIqwxJASNYZxGFuq1xJqHJs1qvsrVrHQkFI5jN9ack1v2AkJ6szq0V2rbBh56YvlWP6vq9y248nQAPynwhHOy1qHrmi7fou4696T64xPn8jHapHMQ9RmlN1ydOlxFVGDtlVAnOCU29T1sdsLisOH2shL5lNm6oRzWLtmXR7DP31xBVjLFYlzBWyOsGZfA96/WaoQ9P+/V8IFhnODhsSMGAOCpV+kWPSQbtAz504GBFx+O3hG69oPcDB0dHTG4cMW1aRm1D01a01YS6gsYZBCX5iDc9imHV90QP6iPEXGlwvD8jOiFWCbWKyibf0OzUSvGb+Bb3gyia0sa6zVqaVCF6XEy0VYPTlG9GqsTNe4jwDqv3yR+bg1xdjJslKJD07V+4Kryyizd0AEE5ODjg9OSYvl8zv4zEGBiPp4gkKiukACRFN5bZaDSm71csFxckTUxmU2azG0zGB7R7+1ysVzw+O+Wte29BjNS2oo6gEhDVLDJGSEBCUTkBcdnKE8NysXsCnF0NJrtQJIe4VBOKElNESRibf56IJLKh0jQNzip9Dz51VHULCKop34msgCG7/dJAbezGlWqoasdo1LBceXxIWJtdsSklUkq0bftew31aC1h5682e5cEp1inBrNEqMT24wXMffoaTkwWmNpzPV5y+fEmQxMH1DyGyorJKUyWms5a22We5FM5OB4Y1OKkxMuXSH/CV5V2G0SHD4TFTbTBywcXFJQ8fzolDxA89625F3+cl1xMHzg5RN47nXjgiLiOXZx1nZ2e8dnyCUcewXCOhBxNZDnBveYatHVGV/aNDLq4fMnz4DmMi12czrt89ZNwIjQH1Ae/XrAbPMO+Zv3XG4uSSfrFmGDzeJA4/8gzmcERzNMbuV8jIbizh3USBuHEPZMs1/zvGlJfFIqgIScEYaJ0waSGuO7zUJKnRjU9SN35iRZ4I8W/83PrEsuaJIH/zaHaTvutYXJxD9MSh46Jb0PUBW40QmwhJ8AkMiqT8sNYyGbfce/Nljh8/YG865ROfnDEMnpCgne7R7h0SqzGr9SUjm0hGqGz2MWNcXjXoxv2lCslD6AElxt0zgLL41tR1BSnh/UA/9KxWS4zkm3lMESxUjUNJ1HWDdTbHbJoRowMHalnMl6gkrAgqiaCQSBwcHVBj0ZTdZd5Hht4TY0A1x7CMkY3/3Hz3BFhVmJ8EwnpO0wrNxLJ3/QCtA9demNK7niF41l1DFfc5PeuZ1nvcun2AqS7xEjkPKxZnJzRmj+VacLWjrnp6Izyc3+Dn/uZA0oc8o2P2928Qe8PpwyW+V4wVhhAISbFVTV2PkB1M5JgdTPnBH/oeXvyVF1k8WkDXsVokUjK0KNPaMqqESWXYd2P2RjNwFtRizjyXqzf5wuv3+cbnv8wP/+gP8cmPfYi2rTE+4BcRcxo4ffOM/qWHrB6fsl6tWPVrjtcrXnvpDezRHqM7ezzzu5/j2vM3sLVDd3i9LQL2yndgNlbLRkONy4FInMFGIaRAt1ywfvSY0dFNRvstCSEkyHaifrP74R3o1eNJwOrdXrSbIhxi4PP/4P+FfmA6m7F3/TpDiIwmY0aj24AhqRBSZN0tEDMh+J5Xv/YVfu2Xfp77b77G0bXr3Dy6wc1nRzjneO7Z52j/uT/IC89+gkcP7nFxfsJicc56MWe9XtF1HX3fE2LcrBoEZx2Vc1RVhXWLbU/LuyCMx3vU1iEkQvRE7XBW8GFAREhkK7VpR4SuZ9V1AFRtizGGiGG97mhne3TrC9BEZR110zB4j/ceVU/lamJKxBRp2hrVhDUG1Qg4nMvxlxDe+0b1VAJsEELnGe3NIAmnj+ZcrOZgXuFTn3iBT3z/lHEzppIxw7rhlVc6lpdnnDcO2w54u8ZfXHAw8dzerxnvNzgVkgrnK8+rb73Bq784xy9fp792wepwybWjgYObFc2sBtNjbEXTGIIkok87KcCjUcX1O/vEkFAPR6MpKXhOl+fcmFU8OxszcYITqI1jf9ZSj8YkDK6qqBtHXSvLx5d8/ef+PtP7n2I43KMyhjgEUg9updRR8YtL/GpFHDzhcsn8ZM7oZIU/m3N4NMNfnyHXxtuekvfEDwP3772VMz6cBWcxRqhchUkJZwWsRayhFiAmglZUh3cYXIXvenoxWFOhoqRN0wcj2WVDyrIMVwL8jmcbN4Vs7G/VbOnsogRbawl9x+WjR1yenXCxWjI5PELE0DYN+4fXqeuavh/o+zVtW/GNV17iC7/6yzy+f5/QdZwfH/Pg3j2a2TXq0YSmaXn+Ix/lhec/ln2b3ZL1asF6ueDy/ISTRw+5f+8er77+OuPRiBs3b3Kwf8R4MqVtW+79pd0ru6EKkcSQBowmSJ4UOgZNVHWLqSqqpqYeNQiWmBxoxFpL3w2kGEhhoFutmE4n7I+niIYcf0lCjLBYdaQYmYw3wUsLyUTaSc16mcU2aYSUV27fzvZ5KgFOSYnRs+oXuMZiXKJbekLqefP1+8wvaibjmvFoymR8yLVnKvr7a04XF1jvaUYwmrWMRx1iHuNcQyVjVCfYYUwYImk4586dMz70zCVHex2zWWR2WONjQwiObnHB0A2QDE6q38p39b5hrbBarYhBqE2DdYZae0QHahFutBVH4wZrhSiCMtCvArZuEQsVhpmzHNYzWoTh4QlnlwuMuUpBq9EoWAfN2JGiJcZIYx3nyxWdv2DUOBoPNQY1VxbM7tH5wFfuPcaIYI3Bmpw6V7kKR8IZIRkwzrE/brk1GXFUO9pRk8U2DmgSVn5JHwZSjIixT/xwimKNZfD+7TQ1VXwIOSDjHJWr85Ibs0kd2sG5EmE8mXCeEuvVChZzln3P/HLOa6+9xt1nn2cy3c8Wlx+onOPXv/QFHj16yLrr8rI45iBmSjltCoSqrqjrhrpqqdqWg8MbGIHz44fUriUEuHbjDtevX2O2d0A7nuBcTj+t6mbbs/Iu5O9WnKFyhsrUxJAwtsI1LYvVms73jDUxaiZUtiaFnspZYlRiAiuJ2kCKA6OqRkPapPSlzeeuSSGACNY5jFO6oWcyGbFeXlBVNSIWTTn75tv5/57OBZETc+n9CqkrmqYmrWtWy46z48Bq4anbBePJBQfXzhk1E9wIDJaqFiZtw8F4n1kzwUpE6bJjXBOKpV/v0TTKp7634uZh4GCamM4cTWN4eF9ZrSL9omPoeipb4dpqJ912IkL0Ht8NhN5DMlgCGhWiw1HTuIqqFoLAYtmBcTg7wjlBiCSfcI2lsQ6/7liFAdFIiIEUDH0fEBNpakuqHUMfqesGjSuGYWBxvmB1vsKvByqq3RQVICrMo8k5qTFbqwZFjMeQqCXnU1orEJVrtWM0Nuy3gjEOHyKX68Cj5ZJHizWqOQvlyv9mxeCsZeh7ILs7NCmD98QQsM7S1C3GbPzN8G2XjNtCFepmhBhL9J4w9ITgOT8/5+x8wcuvfIN2NMFaBwqTUcOje6/TrxYMXYfGQNtuov4ipBTxPl/PzlU4V2184kK3XnFyes69Bw94dHzM4eERSWG1XjNExbgKRPDBb3tavoWrfF8RwVpDVRmaSrCuRm0FqzUx5DzpoeuRYDCaaIwQjSFojj/sTRuy7boJYKccq3B1RaVpYzhpvh5TIsZAXVeIJWdMiAFrSClb5O/FU5ajzJFEHxJ1Eipbo85y4YVhbUkR+qGn80t6OWPczJhNjpi216mcZdy0TN0+Y+dAPD56Uoio9PjguTirmcwMn/lde8yaS8Yjw2jkSFF4dG+gX60JfYSoefOHkScLyl1CBGzU7BYYBpypaGsDagipYQiWEMFpRDAE75nuj5lOW8RaSJG+j1hxjEeOECKaApo8oe/pVp755YrxtGbcVnhrNzciQ0zgfaI/nXP88Iz9syVHz05+0/0JW0MMuIarTLIrV3Aip44FwKREpYr3nr4fiMmSoseIRTUx+IHL1cD9y54ohitngxGoTN7QErx/2/mwyQEOISAiVNWQLeBNupsPcUuT8W0QQTZWPaqEMBCtYbWc8+D+PV5/4w2UHDxzznG0N0V9R+0s/TAgKT45B8QYwiZpOqWUNx2kSPCei8s5b77xOi9/7df5xstf4/TslGtH1zg+fYyKYGyNrbIFvFqutjsn70IOwm02jGxWVZVrsK6iDzkwmTQhSRlWa2yyNDYyrkZ4VfooVLVlMpmwDptzSTZJ6cpmQ4o8yfW9soxTSlhX5U1oXG1yAZBNts2789QuiPW6Q5wh9AZvVggdTe0xEmjHE2wjqOvphgXrVWI2ukPrjnDiMIMjmjFu2oAb8HFFIqA4usFx/55jdlTx0Y/NOBjdQtMS7z2rOaBzNPbU1uGqhrp2OCe7uZk6gQmJSoW2qpiNx1TBIacX+GRY+8S6j4gkjHNoChzt7zE73GPVrVkvO3wIhAAxBFIlpOhR3xO6nm45sLyco6mhdXtoVAYfWK0H1kPEB6UfOh49OOXw8TmH6QY7mwehiTQMXN1Kr4I9+exNJFEsWYCJgZNFjxC5v+hJGLoAyy6w6CJDUuYpomqQzUlvRYHwTVkWqgLqSMnklUHMr9ErAX6XXM5tY50jac5xjTESuw6pHc4KdWXxvifEQBSLpshyEUnDmslkgg8BUqLrB+bzOTEEeh9p25YkiXVYs5ovWC7nfP2Vl/jqi7/OG2+8yunpY/qu5979N3GVRVUQcYixGBEWO5iGlnN1lRgTKULSHMRVhaHPO03rqsobeaKiQ49UiVlTM+/WBFUmkzGz2Qy7jnRDIMSUc4Elp+TFlDb7EfIN/GqDhqa8gzWlvEsXhZi+fb75UwmwCPQ+0D0eWJwI+1PH4XUYtw3zbk6SnsnMYWrHxXmDTRPiYBlWPWKFyo3BTzk/XWPHifPlwHq1IqYxj08s9x9abt+9pJ3OcabDd4Y4zDBxxDN3BmajEcf31hjjcJUBE3Bu9xQ4xsDy7IIGy7QZU7uKLg6oS6yl59IbZrGiEYczwmw8Ym9vxHRak1LH6rJnNV9hmWCPZiDgh0DsA8krw5BYLHoWyyXGGFSgT5HT5ZLzrqcPSu8HmpNLbp7O8YPfWR9w3l8RnyQuaMqnq2yyISJQm0gtUItyulzTe4cxPT5CUqg2uapjk/A+EpJDksWQg2rJZGFPpM0WUnlibb9zE0d+1W7OVF3VuKZl2fUcn51CJdy8c/vJ7r5s8elmu3ok+sS4qVl3a9adxyC01YovfulLvHGyYDrbZzTK52boPa+99g3Ozo95441XuLg8Y/ADSdOTqUkpoon82PxsF10Q2dWWd6qt1x4NPXY6AslxkLquid4jUZmOR1gbcGlNXK3pV9l91yTHyfKYi3VA6hpXV6jJOecheLqhpxJD0zY0dZNXA6s1lWtZ7q+ZX66yz10EV1mif+8V1VO7ILSPzMZTGiyuU+LCoY1j6BdcnPcE6RjPLNNmwtH+HW4eXedw74iKGqObbXp2yXx1zP3jNzh+dMLF6R6PH1ri8HGm+8c8Pn/0bBzWAAAgAElEQVSJUFtG5hazyS1GB3d45u4e+IblWU4FidoRteOvj3awq4oqYbFm7CoGCVycX3I2LLj2oQP2xxXz0wtOuoaq3qN2jqPDfabjlsrBeORYNpbTocOtIKVrrJYdvl8ThkDoExfzjovFQNIB3BntdIS3wiIFLmOkV2FQ2E+WPgkhvUfRgx3BmCfrNXSTz5sFhY1lbIFEbZSFbZgHZdQI4pRKhMpYFn1g7AzjynK8jCySJyE4kZwpI+/MhtAnFvIVqtnft6sYMdy8/Sxf3d/Dnz9CROkH5WLeMyTB1DWVJlKIJO8JpkJtgzUVLho0RpZD/P/be88nvbL7zu9z0g1P6oyMCZjE0XBEURR3V5tUq9rafWGXX9t/pe0qb/mFtxyqTNuStUtRYhjODCdiADQ6PfGGE/3i3G5gSIIryJTxbLm/VahudAD6uX3v75zz+30DP/v0Uy7++m+pqzFGZ7FCCJGua4kxEKMffEUMV8S9lBBSPRPzXJqJv8oL8jtQ1zXeQnA93nuss0hVECMgVZZXC4nvLaWU0LWcPe6Q9YxKjWhXa1arNWuXmBzu4ElEkYgpoQqDKQz9umVU11fzAmM0SkrKwrAir1RSKYwxxPTiTeLLFeCUMBEqZRgpg/KR6CIUlsoY1quAswk6w/7N29zcvU+pRwRbIihRQpNEpGkTX3zzmOP5I2zb4VY1rAQ79Lw+E6Suw6USJVuEPkPIiC4bDqb77E5u0Hee3rWZf2e2j2IlEGghiS7StT2bpsXKyPf/9Lu89eCAL/6vv2Hz9Zx17zmcTpiOR3jvCE2P85YQHLa3bGJisVhCyvQY5wJdE1muLQGDUIrFpscXhlgWpLJgE9f0UeKRRFMginK7JbeXT/FzAomrOPCYiAJ8innRbTfAFGVqdmpNbSTTsuRoOuFoVlCXEpXgk6dLfnG85PPzjo0fxAkCQngmvhC/9kNc7oi3FTFGDg4Ome0dYI4f0TVrFss1TdvhvR8k3TEX0RBwwHKzyUfiQQnWWUdje9q+wzmPQj5b+MjvxhiRaTA8gsFM5tkJ4fLtNl8vIYdeLRqZBN4nSi0pCpWl1QiiC2zaniQ1E5M3dc18jY8JGT1Ka6pRwWRnj8Y2OGczgcaFPONREucsIUiUVhhtaJs1SgrM0FYMMSF8oDAvZmu9XAtCCuq6whiN0YbCQJCeGDuwnrRJ4CTCFEhXU8R9pJ8R3ZgkahAaFy2bTU833yc2liIEVDrCh5Kq/ZI7hWBX38AMTixdf0EbLhDdhlpLai2QpUHEQLBiK3nAAsF0dAhyQRscnfLs3Jrw4Q8e8N67t9grBZ//+Fcsvzynj5bgI8uzJT55QgqsV5bGJWxwrLoOjSeFQN97msZjnUfqghA9m94SbEQWCqtqulRgARstwShUZTCl3kbBYEbKD73gebYuXDlKydyGiHh0bHhtVHDv9dc52i0plaTQmp264u5BgZGJFCTjUqGVoPfwy/MOhqGLGGhoeScnhqlfvFLQXf3fW1hZUkwUpmBcjxmPJhAjEkmMASUELniSz/eJFJdOXAHnwhB9I4gpYJ0jpaxiu9Qf5uL7TE2YnivAv+2++XbbZrtw9dpElk6mmPvAzgcmo4qyKrHW0tgeLxJeCWI5QuhId7Fkvd6gSezszSinNbLUyGSQKRJDzFJvEmVpsM5mJ0IhidHjXY9IAiUFKZIFLMlTqxfT9V5SiAGFzt+SOy2SZDPTWHeJyikUBcaN0WGKsEfgDxFpgpAVCEnyFmVL9oqSQpyTXMT2higStfuKA3GHG6M7NGlDZy0+9LiwIoQ169FTxMSBrHBJs2k1MW4fZYgEo9EeyRhCJan3R7z1/fu8851b3L63z6h4D10oPtefsPnilJOLBa7LU/kkYeM8GxspC4lLiegj3nm6ztFYTyQ7enUhsfHQd54oeuZ9pEsanyJdSpkqOK4oq2IrJdsw9GbjM+OglL6920oxkRQIpTGm4MZI8vbRmL1phTEKkiCGSNt0bPoGYcbIBFOj2S310LxIpDhYWX5r2BdJKWQ+8WDPuI3FF3KLpNlsUFIyGY3yBiTk3VhdlSgncEAYXkNKmSkUQh4G5R1yLgpXxkSXhjXP21nClYVpnoWKZycSnqkIt/QykVIkOHdleYvItp3OeyAhZe6ThxSJShC0xCpNURakRYtzefAbBIxqQxSBJNNVC+vSGCrTHHOxlVKQQiA4Rwoqt7dSvndjiITflxKOBLHNXFSrBESJ7R3aCCoqVJQoxow4ohR3SN0tfDrC6CkIRRgGLjNzg+nRWzg6+sZy4ed05gsm+ueUcZfdyQFd2yK8yXK+5HBB4Cx08XOEEFhR0/kxPnR//9/WPxBiSPQBQimZ3J5x460D/vG//gNu351hysTRawcII1ES/ub4lE++fIzsTX5YpMCmRJ8Co8kIVWis7WianqYPWC+udj4r39NL6Lue1abneOGwKtPQLKCqgmpU5yPRtj4yw40qhuk1XLYJUi6WMiGSQqgRoRzRKsFitcYTqaoCHxKn8w3riwvas0fs3X8biWbdOXQI7BSwsICQeWAlLgtvApHyAh4j+cnKHlfbWF5C8Dz8+nNs31AYTSgLbPAoJZhOpnjn6Pse2/e4EOh7i1K5SGQ/4DhYfeaG1Lcl22lwEXuuuA799+fxW0Iqtw8pZVaNMkiZ+eX5NBDpe0tMAR8iSchs8SoFNgbGVU1VjyiKBpEcKUWECATf411P9G4wB5MQE653GJ1pgRKZ7VB9QKR8JldSEmS+hs69eFj5cgXYQzgHLz3JSHRRYvQY6SLJW1KQ1PVtdqp3kP4+Z8clRmt8lSWB3jqcC5RVidQwne0zLjR+NmG8u0FXFet+wWnXoaaK2WyC9CXejlHxFgnBfPUxUnm0uMHu7gFKb5/Noguen3z6E3xl+f6ffMA/+fMfcPfBDqgWaxucC8wOxjz48G0effqYn330v1PYEd7mHp6Rib1Jwe0beyChdY6lDaw6aL3K/FadWGuHnpSsm54nizUXXeDw/hGPHp4QvCChCT6xWW6yiGYLMThKkmIcWAnPJcYOBtoIifWCh43i44Xjs/ljZjtjhJTY3rNZrfGPP0XNv+Tf3L2P79fUwfNgt6QYjfmLr+aoBGOjGBtBZTTGSIRSWOvoup5N71lZWDi5heUX2nbDT3/6V4jksX1D121QRcloNOLmjSPadcN6taQVEEj44LKTWciGNM4GRJLfWlueGZF/m01/ycX+zxIpETpLkimrHAuFqRRq8GUQSKqyQunEcrWmntQYIajHJdPJFLfp6Nslwfds1hfMDg8g9BAsWmr0YErUtD3TgylFUWSqbLOh0AWFrrE20fUBYiCS3dlehJcqwFUx4t17PyBqzagesTOZUU4miCg4O3vIonHI6iboe3TtjM73lNWC3hZ46+m7jq7rrpQjk/GYalTRdj1P52s6VfGkOWV0njCpzA7+faSSJTv1fbQRBKeQjJjUtzjce4fCvNhp6FVBl4Y0i7z21i3e/+Eb3Ht3j5AaQuyxwYEAVUqmh1PufOdNNuZH9C6xbi0iwa29CW+8fo/X33iT5WZNe24560Oe7ruIVAGl14yPJN//0z/g6WnD4m8+Y7Pa8F/+1/+KH/1vf8lP/u9vuJivODtZ4u3+q74kvxOShBJZoeVTxA5HtstdcD4yB/qYOO8ThYpUXYsRAml7xv05fVriUyQsnnIxn7NYLagnE/7kj/+Uw7ricGy4tTNiUkoKIylMNksJMRKs5WTe8JOvzvn3H53w+RZW4L5refTVFxwe7ucZTFHQO49zPVoJpIhZSB0TTd+DkJRlRdc1eTgXfBawPPdvXhXawS3vP7Wrff7zYktbWkIICmMy1zvLedBGk6lpMbethMIUkv2DfQSBopIkPEIE6tJQmwmL9pxu3lBPxuA9hZTUZUFwPlMeTYXddCQ79NgDGKlZLxf07aCCtYEksiz+RXi5AlxN+PDDf4GoK6qqYlTUqLJAofnksxGbJ2fMG83p0zVat5S1pPAbFqHAtnmyH/BorWg3a755/JBCK6SQtF3HxjnO1he87u7Qnicu5ies5yeYoNkZn3Hj3oR1OKOuJqi4YCwv+N1WF68GSmv+q//mv0AXksluhZcNgoBAUZZjICIkqEJw6/5N3nz7DY5/eQHJ512KrGhQ/PyrbzhdLHl6uuZs2bNyiRaQsWdcOr77w+/z4Z99j88+fsLHD0/o5yvGOyV//m9+yNPHHWen53z+6de8/8c3B6+D7YMQUBrJqNAgFK0P9N4/T45ACtAiK4pUityeFLx1MCE5y3rV0NDhY+BsveE//PjH3H39TW7ceZ3JdMz+ZMTNg30Ox5qduqTQAq2gLAxSSZx1rNuOVdMQQkf0lm1sQXjncJuOtCtQVYU2gvOLU0bliLZdIUQgpUBvHb0XlOWU9XpFdHnnlgqJ75+1dr7FfuA3e70MHtLPBnPfvibbW4Al1WSM8xbvLS5a+j4gRHYE7DuHLkvqyYSdvR1cv6KzDXVh0Dogpc9eGzFBUGwWDc5mHnXqUxaCCc26aXF2Q1EWlGWFRrNeLrm4WGCKMUZpghL4oWf8IrxUAS7LmvtvvE8yEqUjRiUkJbJQhC80y1XPyXyFjWuqsoCLSFWOmNQzCl1S6gLrE9F6bGdZzZfIFJAp4kNL052yCY4m9hAgeMlmrTg/uWC229IUU1LRUds1nXVsNg63hT3gRGLv5g7IiDZZAiClIDpBShIQgwQb9m7s8E/+xQ/5n5/+n7TLSHAJryQrFUmjgsXasy4SfqohRgoNs9kO771/mw/+9A/Zu7/PdLFmtj/BfSx4enLCH//RA1577YCPP/6GRw+POT9dILeQLQLZs2FUV3hvsdbS+SzrlMMDnlsUCZEiisSD/REf3p0xNYp14+lEJAaL7RqKsubo5i0evP0eRzduUpUFVZnbYJWRFEZRFSpL2FPg/PScn//qCz473/Dp6YZfnfWsfdzK43eMER8ifdejjMmUUJ1bKdZaSm1IJFzfEX0kiEhwNg8xUzafD4NT3OWR+HfNBXLdyD3PS9XXJcSvDe22CUJKbIokFZAq5cKZctRQSALnI7btcYA0Ei09m2bF7nTC7v6IIgWeNEuSTbgIq0WLRCBTIghP0AUSR9tkP3IlNMIItNC43mWzLCHQSqIVBP+7N4gvVYCV0cx2jwjSg2yRskekhB4lvGuZn5yymDfZ1FjvMT9f0OoxYiegZruYokDGyGqxZLNZ0KyXCO+RyRHTBi8WeBmZN3OqWmMKw2i0z7kKNM5zsW4ZzQJRRKw9p91IfLD/b35f/yBIKSL0paoqHwsTkpgkKapscJ0SiIiZat793rv8xf/yU9KjFS4EYq0Z3dvl7Q/fpH6yy95FQ9MF+hQRBo6OZnzwwRvcfnAbakE5Vezsj0lR883DY/78z97jjdcP+frLY5pVy2bVvupL8kKklPAuSz6tDwz+MJdOkTk0kUQhE7Pa8Cdv3eTeboHrerxTaBHBW3Z3Zrz77h/w7nvvcO/+fXZnOxiVj9ZaSUqT8M5ytu5YrpY8PT7m448/5hffnHFaHnDiDEsrsfEltUn/HyGm3NdtmxVKZe/a2WxCaSS2b3DO46wjBoeOARUEk7qi7R29c/gQh7icZ7va38ZueO4DXEotfj0tZFt3v0BmWimJUgKtc7SS23hEyiwuQSJ6T98GGiMYjwxd1+G8ZVaXMC6yyrYhS5CFoyoLUkz01mFGJpu8dzZ7B4fMXlJKEWK+ZkIKlJDZfc4N9LUX4CV5wNmvVqRB3qkcUTiMlsRuweb0Ef26Z7y/Ty3GPL44pRMrRG9JzlGNRvjgOT99QtOsaNo1IgZEcgi1Qc86dFmw2qxBFtRmyv7hLiHUmTieLih0ltV6J2it20pXfoYhCCIf96KUg2GHIiWJH0KikQlRJPbuH1AdjGnxtNFxMNPc/uAOP/y3P+SNsxMWizVN67A+IBTs7025dXMPqxyOjqKW7O5NUbLg/GRFWQrefHCTn//0K1pvcf0WmssMiDGxXreZ6zukzj7rluQSYCTsVJr3bs/4l+/f4WzRsBGJ4Av6qqAfVdx98Boffv8fM5mMKY3BiIQWESUFJMt6seb09JRvHj3iV198wS9++TE//duf0UzvMv3wn5PKST4uiu3ki+SFqqVtEoWBqjJMdncYj0asFguOnxxjraPUCqlBaUExmeAXm+wPEvyVTBd+s4g+32a4nH2mlK4sGLODWrz6vm0uwkprikpTFgkVYbUJiBAg5GSh7F4QSc5myXoaPFeiB5FQSiKQRB8RJZSlJjiP67NvsOs93geUBusd0nUYURCEIApBGvQt0kNKnhBe3P576eU+pZz26bFE1RF1T1p7UnNBbC/QQCkm9OdntGePc1TKZs1yfo42JUJD161p2xWb9iLvcmSkGDluH2kObkypyoQQCR97tE7cu3MLHe/Qul8gisf4AMJMqeUhiU/+/r+pf0DkWBKZhQQh0ntPShrIAxGpQEiFF1BNJaMbU5bB0gWHGCn278+Y3KkxR7e4DwTrc3/PucHzNiBSwiidjaN3p4wnJYUZUY2nvPnuPW7d+xVPnp5DEFv8wCRCyj1xrYbdWQSPxIjMrxwbxRuHY/7sD+5wY6JZXgR2RxXTQrE/Knhw74i33nzAZLqHsz0508DjbcfZ+RlfffU5f/2Tv+Xnv/yEr756yOnZKU3bEoVm94c/wFIjg0DGmKPXt7ACpxhJYZNPBGmEFBGtEweH++zvHnD85ISiKNmbTTEqT+k3LuJDygGUQn6LCv7bCvBv/J/Dxy+/9rIQv+jrtwICUgiIGNFCYLSiKiqUS5ACUkmqshgYEppCSfZmO4yqETFB23V0vc3hpA7K0jCdjQb7yg5ERGqJj9n8XxCRIhIJmfevJFEOC5WKSAXp97UDTinRdR3WdfjQEtOGEBYEHO1iRRIBVRVEmVicLTg7foRNHbO9JePpIaqcEBJE29J2p3ThlKRz4qipNYf3bvDg7dd5enGC1hFdbCBu8P2Sm7sfgLwNo5YnJ+c4C7PZTbR8ceLoq0JM4BOUSoEUuBhoekvXr4dhBxhTYCiIQiBEx+zmLpQlWEFZj5jNZrjgcXh86BHJZ0MREXDJ42PmzaqkKcsJ050d9g4Kbt25TZRjZjc1t97YxyYLcWu90AAQeAQShabQhklh6BkWMR85rAWv7RQ8OJrRdhYRekpTMTs8oK5vI0WkMAXRdcwvTnj4zTd8+dWXfPbZZ/zil7/k8y++YL1eZaFLCKSYUEJS791geng7K5mCzz3PLbXjKbTkzTu7TGd7VKNdRFFjfc/p2Sn3br3BrRu38d5S1wXW9jz66GPOFmvc4GP7G0O25/CsmIpn1LRnDYirr7t0/oL8u9lGiJQoEdBGrEtEFQk+YaRCIQkip1hUVcVoXCNEZLeeUhUG33X01pJIaKXQRhCTpGkiwQmKYsTRjSOUMhwv/hZtFLLSUEi8CIymY0bjUT6pdD0hBkxRYLvf0xAuq2dAqwpFxPUtdtFz4ZaUsyl/+MM/5MK2nC9XnCzPOXgjF5gYHFE1SF1RmDHny2M6e8r9dyfoiaQLHdJEZgcVCEPfiHxxyBHSu7MphwcSKfdZ22N2xgEndlCU+P9Ek/tVICXofMIlB0TCoNYzhSGGPGTyPuCjxflIqQTj6QhdanQPRTHBqB1sG/GyJ9oOFXPwccxqSKL3tF2LEZp27lmuVvS+ox5VnC9WWfmlDVpXNE14Ts+/XRAISmO4uTvlzv4Ot/em3JwWLJqG1aZFhZ5RaNmlx/YBawPjSU1laozSaKEojGS5POVHP/oLfvzjn/DZZ19x/PSM5WpFZ1uc9yQ/mGsLgTQGVYzQe3cRxQhJAKFJQm6tKU9ZaG7fPMTUUxoLp2cnLJqGrv0Vp3cvmE6nJBF5cnrK4ydPefT0jIgcXtOl18NvDs8u/55tHMmybbhqxMcUyN3ToX2xzSs5QAwI20MIOCJeCaQ21FVNSJHG2eznqyQxRFK0zDcrJpMxWgjKsmB3dwfnEnYFjQtsGotREmU0i9UiL2iFRlWGoCHncmd/7zQo7spSQxJsbE4eeRFeegdMkgOhWyOSRhuDrivqoxFlbLk1K7l1z7B4IAlxgrPw+S/nzC8MSIMpKrq+oR5L3nzvBi0Llm2DLhWHNw+YTg94/f6UpDe4oEmxJSI4W3/BdDLC2ZpCSqQeD8KF7XtYYoLVpkcp8MHlnUlp2JlNiSr32kKMWOvpOktUcRjcCWxwzJcNJ8cNe7csomxJwZKCwLlI5xzWO2LMRuSFNMwvGk5OL9i0DVJLNk2DFAbrYLnqefT4fCuvE8CoNPzT91/j9u6Yg7pgpiU7RaIzhrb2xN4Se5iVl74FicpIjILQbdisOy76NX/90Sf8d//uf+Kbr79hvWnprcP5YeGJObheSDUsSgZdjjD7N0FXg3TXXynAthFSSUw9JaoShyOQlahSSNabJTZalqsVJ6enXFwssb+2Q03Pub99uwhfTj0FSpJTWwbDojgM7HwIhBSzr+7ld23pdYKEZsh/UzKbORHorUVpRVGYK0N5BBTaMK0q6jp7RHgsxmRf36Oq4uHTE5RWSJ0d5TrbZUWvgCSHtyINeYQBKQCZUEmiDWgdaZvfkxIuJfA25B0DHqkiVV0Q0oa2nzM/f8x4XDAbV4T9Hh8Dhdjl+OGKvqsRcopQJdZ5DmcjRhND11mKyjPbrSmM4JuvT0AmzKgFmXO8lssNS87o/A1UrNBygi4rfMoBjtuGXBw9RmhCzD3gEBNSGmJ0hBDy1Np5YggEIiEFhIY+9hyfnvCzn37K3o2avaPca4ohezE3NotZwiApxUhW64aLxQrrPcoorMs8WusCm6bnyfHZs6/fMkyqgn/2nbvIFGnnc5bH59SjPK+W0bE6P0aREKMcEdOtzvFdR9v0LOZzlosLFps1X55uON84uiBIMsdVyYE5kGyDEkW+VxIgFbIaoWdZoBJjXgCvmr9bGCEtpcIlxWLZsG4trXWZUuUi8+US0TZsmoZV09A5hxiMaC7d364K57eUhly1xPLAU4JSV6rJ4cP4mH2Bc9G9pDNu3zWC/Pq0FiiR6XaXXiMhJUQUJCkJMZsSSSEoqoJxYXJuYAKSyNaSukCmRBQ+256WEqIm+UyhLUuDNhqtFegcdd/3HTujcTbGDxEhElorUnoxU+ul7SiD7fP0XrUk0RBiw+rsKc3FGYunp7hCQajodlc01nJzd0ylI+O6IFHhoiYimO4XJOFxfYeSkXGpaZYLPv3sa0Z7cHg3MaoNMmqiTaDWtNZRyyMOZ1OmownJiK0swJeDC6kkUhmUBKM1MSZs764KaMqbjazGIqJKRTkpaG3Lz37+Ebv7Ne+9f4NqLJE695Kdza5oIeSBQvKO+XzNfL7M5iqKPCDwnpgiLjjcyhHDdu6AK6N448YuT87mfHF6ysOPf0mzW2CKAmTk4vyEyXTK/tENuq7l/Okxp0+OefL0lOOTU07nSzY2MDl6jQfvfofm5k2SD7gQ6Lxj03XYZk3oGvrNEts1BKlI4z1kPSGGfjD2GgxXBpO0bYOQkuXG8uR0zrrrCQl8H7C9Q/YOqRXW2exzkLIXQUwipzNcWX2Kb7UhUgIxMHVyLqMEqfOgKeW5QYyXbnE5NTrbRgxJ0tv36AE5ky3FkFMxBvaGlHqIX8oqXO89WgqMGuUlZXA6u2R6VXWBbRtMKVBFwlQSLUr6BuxqQ1kUaGNQRoJK+JDofUtZ7ZE6h/eeGNO3Bpe/DS+ZiJGIoUEJQQxrGnvM2elHfPrjjzByysgbmrML2pNj1A3BRdcxe3Ofgx1NSolV44gOTKk5uhNRpocgkakAG3n89Tf84qePuf2uREwEnR0h7JjQaHb3NY8Xj+lXK77zjmZ3d8p4JkFs30oshaAqM0leSYWoCozS9L1lvW7o++4qCJEY6GyPkFDWhjfffoOiMpydPeXf/ff/nrPHf8S7373L3q0RupSkpFCygJRdllZtw9OTM07PLxBKEpLHh4RzeeUuR2q4t7awqpCPuSdN4OFFx6ePz/j08y84LSM7dclOEZjNZozHY6SpOJ+veXrR8ZNPvuSr46esWk9Ikno84w/ffod33n2dUgpKrWl7y+l8ztOLBYt1yye/+Dmf//wndG0LRY3au00UitAtQahs3qLUs2K8ZUgJzi4WnJ5dsOp6fASiILhAPRqhQxZqpDBwp4ed7RW9bLBNvBRhxPgst0wMg8cYAkpmDpUSAjWsSJd0tNzG8llKv7WG7IIkc+hoTJDIoaxGG3pn8T7iU8RLSaEUUkpiCMgk0EqTTELESFkZUhc4vDEjCUGhJVpKgld4bxlVBYXKmyyh8+dt0WGMIUQz2H76bOKufk9eECE5VvYxYhNp+1POVh9xMv8rzp4uuKHfY+IEeIvr1mhGFF3gcfqMN/7RHfzY0j6aE5qSw5uR9/9IUI0LQjzi/HzNfN5R1SPeeec20zstxvREr3ENXDxq2JwrnO1ZNUtmI83eniPuV7iwnSIDGR0FAiUViJhpRBGKosSYXHyFkKTU4AOUdUmUib2jHd58cI+2u8X/+j/8Jf/jf/uXrNff5YN/dJe9mxMiOcHgckeyWKyYz9ckJLdu32K92fD09CwPqIrA/tEIKPhV/PxVX5Lfio3z/MeHCy5OVqy6Hl1EzpsVXxyfcGsy5s8efMCb736Xm3fvoouag9v3OXj7LTrrmNQVs3HFtNRMq5qiLBmVJUYpfAw01jFvHceLnv/jYB9RlHz56BEr64gi0Fw8IQoo6xGKEYnMqNnGtSpHybv8s0WILgDZDlGlhEqDqZGQUBjcFfUpXb15FhSZcZk1mTNHIlIkRMjcaaMUhZIopRBC4Qcxx2VSMCTOt7ACR0DIEm2GU2CIGFNQVSXrbkOSAi3J/r59RwgWpbNjnFG5rbDpNjQJJgcjtDXEkJWAzllQkf2jPVxvB2LTM1MAAApySURBVOc1n0OCdUKIiNYav05EL4DBgU/9nlgQPjpOV58Q2oauO2PTfkPsW2Z1olt9QTGWTGtNSFO0qZn5HeSNmoOjMVYrTk48TWMBycVTiY9ztN5jMlFcLDxN25C05bDMjCyFwLfQNIHHnzSMTMn9d+5QKMNidYE0BWI7SZt0zQaRAlVdURQF2hQoGTFD/HfwAdtboo3IaPB9wDvLer0EGXnnO2+wOnP81Y8+4uGjY3YeKnp2Kesx4/EU7x2rVcNy3TDemfCdD3eYTveZ7WmWqxVKKkY7I14bVQSvrrLWtg0hBFarBt/3qARaGLwsubmXuPPGB7z+3T/m5mv3mIxy9tZsIjjafxtDNvG5TLuQ5MwuksM7SwSUEIwLw42Z5gff+x6TnX3+489+xk8//oin56fE2OJNBUKhk0IqrmhW24YUA75twXtkjEjiIIxQJNejBChAKAm6wnd9TtMe7DUTudUlhp5whiAOCjEJ1CpRGxiXmqoyFCafCpRWSDEe+uS5aCslOV1+/YquxouRUkKoAi0EIXakYBEp4UPAlAVJ5LpipKQuTGZ/GEnbrIhJIwU4b3lyfMburX1cH7FdT4yeRLb0VEqCyeIqHyPRBWSR6AeXutzSyQu5MuqKJfHb8HJDuOhJnCHKjsK0RCOwTUW17/CzJUKPmFV7IEouLpY0X2948+gtJlPFvgjcvNmzPl/z5Gvo1gpZNZSmZ2wMUU65OFvgohvacZGiFBR7Y/ydKWefH+NdSbeKzC82JL2mtfp3Wr29KsSYWG16rAProKoSZZmTWom5nxd9JNiA7zxtGzh/umBcj9iZzSiKgqgi99474HS5R11rnBc0m0gQHl0Egk9IU7Kzd8Bs52BIfhUoLXBB0LQNyqhsNJ2Gp2YLoaXkYFLgxJjRckpVTjl1AjPZ4Xvff5/bNw4pjMGFRAqeZAVGCWQ2uCIBKgniYI6dqVMx53+JZ45qs+mE1+7e5nx+wfHTJ5yenhBiuPKcyOPsNCxU27eop5SY1JKUNHWZ8CEP2bTSFGpI+U0CmyJNiBRK44vMnAkhEbMP+zMZPCAlaJkwGsZVyeFszLTKfhlaq3y8FgKlsnMcCIRMw8ckxjx+tRflt0BISeMsQQi0FOiyAJkHiaaoBkP6iFAKU1ckskFP3gEPM5kUGNUTUtDYYLExoqVgVFaoqqBZtay9w4Y8CEUmjCyYTWdYawe3R4e1HlC5rfMCvFwsfbB4f0oUjqQjZmKY3dpDOIfXgeRKpDO4DfTW0c17Th7P6W8pvFHsHgre+cAwm07Z24tUu1DNCqQZ0/f7nD/dZbF5RFHMkUkg0ShTUI81aElRQIw9XRthGQgDx3HbEGKiaRzeS5wXdJ1HiQbneqRQWc4tJMFH2t5hbUQazesP7nFwcMRsb4LQktnhmLfev0epC3b3ZpSjAlGIrGqSiUKWkHK/K6VI8AFlDLoYIaREqtzPe15Cum0wSnJ7p6bTiWI6Jkx30KMxNx7c5979+8MCMvhppKzqEkiiGmTLYjjqXb2+ZxP8Z3LaRKEV01HN0f4etw8PefT4MeebHlSBEENx+bUj+nYhMRtr6krgQx7oChRaDsbj5NfqQszcU5lPWs5GvEuEkDcGlwVYCIlWAqMjpUmMq4Kd6ZjSKIySSCXyPSQkUmqkHEqFuHRJk88Wry2CGBZgP9wBkoQRiqIscTERhnivgMCHhI+J4LNMWw49XRUUWiqICutDFkspgRSSyhT0tEiZrjL0Usr8/KqqsL0b/h7w3pFi+J1OhC83hEsBSZsJKFKhqpLJTNK3K2JI2HOJ7RxuJaBVpFZy9vSCMC+RowKRDAd7ip3RDvXYYmqNqUFXgulswt70BhcrRa9yHpqIFaAxlWbvaMROmTC1IyRP3wd0p7dSEpkS+CARHkL0tK2D4OjaFjmocMqyHKhFiaQ1s4MprxnDeDJlPKlAB4QquX3vJrWpqasaZMIli0xZMBBjzHFFwRNTyN6nQF1VKK2J5OKzraolyDvgo0nFkojd2WF66w6HdcFbH7zD3u4+Sum8ywWQAgUomVsOclB3yWG6LwcGw5XuLyXCwADQUlIVhsP9XV6/d5fT+Zzuq8d0aCTZElUMPbtthADqSjISEoTJu1AMUqrhuBsJIeJ9pLKSsi4RIuJ9InhBjFyJTLIoIxdgrRJGJcpCX0Wsy0Hxlv8opMwNnyuXJJHf2cZ2jRjuC0R2f0spYoRGGYN3EakVSYbsL+0D2nqidzmgdBhACimJISE9eB8QkquVWZLbQVVlQEKIIvPyQ0IYQQgWkhqCOwVRZBe2F+GlCrAWgklV00eFE4JERFSR1I+xZ5HVw4a4CKi+ovYzbA2uAFJBtwK76imT5M79HWIq2JwHRN+h6hXVuOT20W2mu5KVg4v5kt4KYlKYMos2ZJiT6LAh5KPi9rLBEaIgoXEuZJcq72gbi9IShEapAlUI0AZVKGZHO5hJMbAjcoKvFopRPaI2dW7uB0v0ETUYhayWK3rX58JkNEorgs/0tEIa2s7hXcx0pN+hR3+VUFIwNobOFJT7h+zVE27vTXjz/i4GidQKJQRagFbZSEUgUCIN7lbkXZmQOWwzj76JpEzVimkYG0m0Vuzv7fLgrTfpEjxdWY4vmmEYJXhmgbmd95VSEqUkWufJvlaXQ8M0FN+Ad5FCF4wnFSH1w6KiBhZEGmYm4oolMRDKnuMEyyvJsRR6KMIaITRy8AdGZKHyNrb/BPm5kWLodQsQKve6k8gUM4jEZLHeIjuLCI5RYbCDWjJE8C6gVY6rSsTBRjbPLGKKTKcTjBOE2BAagKEnHCxalmgjqUVJSorV8vckxNBJcTC9x2LRcnG+YNWuCJVDuQniscecB6KXkCpSmjB6vaZ8E5Sa4Oc9/aMl8/MFftlw9MYMkQSF2cH5xHy1RoYn7B8dcbDzgORPODk9p2layqpitl8Q+hKlBe0GutZi19uZTZUSdL1H9AHnHM76nB+lR9npS2h8Gkx6QsDZHojoWiJExIUe7yKhj9jG4qRFK0NvO7q+ZTadkID56QVd31ONx4xm2SN2s1oipUAqzWLRsNn0xJiGCfYWIkWs61AicrQ7obqxx9GsoBIRqQTqkm6aEvJStHXpU3ApMpAyC15iTj3wIWFDyjLSvqeziZAUfYh0AXQ14ujWHXZnX3C6aHJhSjH/uWpebB+kVGitKQqd1Voqx52n4ZQjBSiR+apa5qEbIiJlQoosWhLDrvVykbl8fMRAU5NDER7covIOWKj8uSGI8rLts40tCImgSDJTQlTmNEeRsN4Tks6vW+ZNkDRgtEKGHN3aNo6YstlOVVYgwKSBshfzycGmiCg0ZVWwbjYA1NUIUxpCaNBKMZg2PPcTvRjiZQqYEOIE+PLvdWX+4fB6SunoVf8Qz+P6Ov3dcX2t/m64vk5/N2zpdYIXXKuXKsDXuMY1rnGN3x+2r4lzjWtc4xr/P8F1Ab7GNa5xjVeE6wJ8jWtc4xqvCNcF+BrXuMY1XhGuC/A1rnGNa7wiXBfga1zjGtd4RbguwNe4xjWu8YpwXYCvcY1rXOMV4boAX+Ma17jGK8L/A/nHMQMW4ZU6AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 10 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u91W-0xKo1wS"
      },
      "source": [
        "## Preprocessing the data\n",
        "\n",
        "img_train is a tensor of shape (50000,32,32,3). This means we have 50000 images with of size 32x32 with RGB encoding. Since our logistic regression assumes the input is a vector we flatten the data. Also the return value from keras is of type numpy.ndarray it would cause a problem with cupy so we convert it to cupy array.\n",
        "Finally, since the data is a row vector we convert it to a column vector by taking the transpose:  X.T is the transpose of X\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "tags": [],
        "id": "QRAYR6DU3Cev",
        "outputId": "d48f51e8-9015-4839-c663-613d61598e69",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 69
        }
      },
      "source": [
        "img_train=np.array(img_train.reshape(50000,32*32*3)/255.)\n",
        "img_test=np.array(img_test.reshape(10000,32*32*3)/255.)\n",
        "label_train=np.array(label_train.reshape(50000,1))\n",
        "label_test=np.array(label_test.reshape(10000,1))\n",
        "\n",
        "print(img_train.shape)\n",
        "print(type(img_train))\n",
        "print(label_train.shape)\n"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "(50000, 3072)\n",
            "<class 'cupy.core.core.ndarray'>\n",
            "(50000, 1)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wtlXsJXtp3Ua"
      },
      "source": [
        "## Binary classification"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nsu_YVsm3Ce-"
      },
      "source": [
        "In this exercise we would like to recognize if an image contains a ship (label=8) or not. We convert anything that is 8 to 1 and everything else to 0"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zB5AbuWa3CfA"
      },
      "source": [
        "\n",
        "feature=8\n",
        "for i in range(label_train.shape[0]):\n",
        "    if label_train[i]==feature:\n",
        "        label_train[i]=1\n",
        "    else:\n",
        "        label_train[i]=0\n",
        "\n",
        "for i in range(label_test.shape[0]):\n",
        "    if label_test[i]==feature:\n",
        "        label_test[i]=1\n",
        "    else:\n",
        "        label_test[i]=0        \n"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Q5Lj6AbgqGFG"
      },
      "source": [
        "## Sigmoid\n",
        "\n",
        "So far we have not specified the function _f_ that our  model depends on $\\hat{p}=f(z)$. In this example we use the **sigmoid** function. Given an input _z_ it has the form\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\sigma=\\frac{1}{1+e^{-z}}\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "We can plot the sigmoid function using matplotlib. As you can see below the values of $\\sigma$ go from 0 to 1 which we interpret as a probability. For example, if $\\sigma=0.65$ then the probability of the image being a **ship** is 0.65 and of **not** being one is 0.35 so decide it is a ship.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nITpqcXBavpA",
        "outputId": "e3fd0e2b-e7ab-4c81-8d46-4e0f98fcd1c6",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        }
      },
      "source": [
        "z=[1/(1+np.exp(-x)) for x in range(-10,11)]\n",
        "plt.plot([x for x in range(-10,11)],z)\n",
        "plt.xticks([t for t in range(-10,11,2)])\n",
        "plt.show()"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAcZklEQVR4nO3de3zU9Z3v8dcnmVwkhIskoJJgQMFCxaKk1NparYhV6BG7bRV3tbW1sqdbe2z1dNdue9xdu499bNvTantqT0u1XtpVH9BalyN0vVRt7UUkeOGqELklgBBAbsGQTOZz/pgfOIRABjIzv5nfvJ+PRx7zu3yT72eSmTc/PjO/+Zm7IyIiha8k7AJERCQzFOgiIhGhQBcRiQgFuohIRCjQRUQiIhbWxDU1Nd7Q0BDW9CIiBWnJkiXb3b22t32hBXpDQwNNTU1hTS8iUpDMbMPR9qnlIiISEQp0EZGIUKCLiESEAl1EJCIU6CIiEdFnoJvZz81sm5ktP8p+M7MfmlmzmS01s/MyX6aIiPQlnSP0B4DLj7H/CmBs8DUb+L/9L0tERI5Xn+9Dd/c/mFnDMYbMBB7y5OfwvmhmQ8zsVHffkqEaRSSi4t0JDsQTxLudrkRw252gqztBPJFcjnc78USCznjy9uCYeMKJJxx3xx0S7iSCW09ZTjhHjEnuT247+AHiyWU/tJy89ZR9R44jZdth6z3vaI8BU8eP4H31Q/r3y+tFJk4sGgm0pKy3BtuOCHQzm03yKJ5Ro0ZlYGoRCVP7gTgbd+5n9ztd7OuIs+9AnL0H4sFyF3s7ksvvbgvGBPs7uhJh34WcMXt3efigyrwN9LS5+xxgDkBjY6OurCFSADq6utm4cz/rtrezbns767e3sza43bb3wFG/r7TEqK6MMbAi+VVdGaNmYDkNNVWH1qvKY1SWlRArLaGs1IiVlBArNcpLk7exkmB7sL+stIRYSXBbapSaUVJilJhRYlBihgW3B7fZYduS64e2kwza5NLhoXtw+8FtybGWsnz4/bWeG0KQiUDfBNSnrNcF20SkQMS7E7S+/c6h0F63vZ31O9pZ29bO5t3vHNYxqBlYTsOwKi4aV0tDTRUNw6oYWlVGdUUZAyvfDe+KWElehFwxyUSgzwduNrNHgQ8Au9U/F8l/Xd0Jnnt9G3ObWvn96m10db+b2tWVMcbUVNHYMJTRNXWMrqlidE0VDTVVDKosC7FqOZY+A93MHgEuBmrMrBX4J6AMwN1/AiwEpgPNwH7gc9kqVkT6b83Wvcxb0spjL7eyfV8ntdUVXH9+A+85tZoxQXCfXFWuo+sClM67XK7tY78DX8pYRSKScXs6unjitS3MbWrh1ZZdxEqMqeOHc3VjPReNqyVWqnMMoyC0j88VkexKJJxF63Yyr6mFhcu30NGVYNyIgXxzxniuOnckNQMrwi5RMkyBLhIxm3e9w6+XtDJvSSsbd+6nuiLGJ8+r4+rGes6pG6xWSoQp0EUioDOe4KmVbzG3qZUX1rThDhecMYxbp43jY+89hZPKS8MuUXJAgS5S4N5u7+SGBxbzWssuThtcyZcvGcunJ9dRf/KAsEuTHFOgixSwrXs6uP6+RazfsZ8fzJrEx885jdIStVSKlQJdpEBt2NHOdfctYue+Th743Pu54IyasEuSkCnQRQrQ62/t4fr7XqKrO8HDN52flc8FkcKjQBcpMC9vfJvP3b+YyrIS5v3tBxk7ojrskiRPKNBFCsgf12xn9i+aqK2u4Jc3fkAvfMphFOgiBeK/lm/hfzzyKmNqq3joxikMr64MuyTJMwp0kQIwt6mF23+9lEn1Q7j/hikMHqAPyJIjKdBF8ty9L6zlXxes4sKxNfz0+skMKNfTVnqnR4ZInnJ3vv/0av7Ps81Mn3gKd10ziYqYzviUo1Ogi+ShRML5l/+3ggf/soFrGuv5t7+aqBOGpE8KdJE809Wd4GvzXuPxVzdz04Wj+cfp4/WBWpIWBbpIHuno6ubmh1/mmVXb+NrHzuLvLj5DYS5pU6CL5Im9HV184cEmXlq/k29ddTbXn3962CVJgVGgi+QBd+emh5pYsuFt7r5mEjMnjQy7JClAuu6USB547o1tvLh2J//03yYozOWEKdBFQubufO+p1Yw6eQCzpowKuxwpYAp0kZA9ueItVmzewy1Tx1KmizVLP+jRIxKi7kTy5KEzaqu46ly1WqR/FOgiIXpi6WZWb93HV6eN04lD0m8KdJGQxLsT3P3MGt5zSjXTzz417HIkAhToIiF57JVNrNvezq3TxlGio3PJAAW6SAg64wl++Ls1nFM3mGkTRoRdjkSEAl0kBHObWmh9+x1unTZOp/ZLxijQRXKso6ubHz3bTOPpQ7loXG3Y5UiEKNBFcuzhRRt5a08Ht16mo3PJLAW6SA7t74zz4+ebueCMYVxwRk3Y5UjEKNBFcuihv2xg+75ObrtsXNilSASlFehmdrmZvWFmzWZ2ey/7R5nZc2b2ipktNbPpmS9VpLDt7ejiJ79/k4vPqmXy6SeHXY5EUJ+BbmalwD3AFcAE4Fozm9Bj2DeBue5+LjAL+HGmCxUpdD//43p27e/itmlnhV2KRFQ6R+hTgGZ3X+vuncCjwMweYxwYFCwPBjZnrkSRwrdrfyf3vrCWyyaMYGLd4LDLkYhKJ9BHAi0p663BtlT/DFxnZq3AQuDLvf0gM5ttZk1m1tTW1nYC5YoUpp+9sJZ9nXFuVe9csihTL4peCzzg7nXAdOAXZnbEz3b3Oe7e6O6NtbV6/60Uhx37DnD/n9YzY+KpvOeUQX1/g8gJSifQNwH1Ket1wbZUNwJzAdz9L0AloPdkiQA/+f2bdHR185VLdXQu2ZVOoC8GxprZaDMrJ/mi5/weYzYCUwHMbDzJQFdPRYre1j0dPPSXDXzi3DrOHD4w7HIk4voMdHePAzcDTwKrSL6bZYWZ3WlmVwbDbgNuMrPXgEeAG9zds1W0SKH48XPNdCecW6aODbsUKQKxdAa5+0KSL3ambrsjZXkl8KHMliZS2DbteodHXmrh0431jBo2IOxypAjoTFGRLPnRs2sA+PIlZ4ZciRQLBbpIFqzf3s7cplb++gOjOG3ISWGXI0VCgS6SBT/83RrKSo2/u/iMsEuRIqJAF8mw5m17efzVTXzmgw0MH1QZdjlSRBToIhl21zNrOKmslL/9yJiwS5Eio0AXyaCVm/ewYOkWPveh0QwbWBF2OVJkFOgiGXTXM6uproxx04U6OpfcU6CLZMiGHe08vXIrN354NIMHlIVdjhQhBbpIhixYtgWAT02uC7kSKVYKdJEMWbhsC++rH0LdUJ0VKuFQoItkwIYd7SzftIcZE08JuxQpYgp0kQw42G6ZPvHUkCuRYqZAF8kAtVskHyjQRfpJ7RbJFwp0kX5Su0XyhQJdpJ/UbpF8oUAX6YeNO/ar3SJ5Q4Eu0g8H2y1XnK12i4RPgS7SDwuWbeZ99UOoP1ntFgmfAl3kBKndIvlGgS5ygtRukXyjQBc5QWq3SL5RoIucALVbJB8p0EVOgNotko8U6CInYOGyLbyvbrDaLZJXFOgix2njjv0s27SbGefo6FzyiwJd5Dip3SL5SoEucpzUbpF8pUAXOQ5qt0g+U6CLHAe1WySfpRXoZna5mb1hZs1mdvtRxlxtZivNbIWZPZzZMkXyg9otks/6DHQzKwXuAa4AJgDXmtmEHmPGAl8HPuTu7wW+koVaRUJ1sN2iC1lIvkrnCH0K0Ozua929E3gUmNljzE3APe7+NoC7b8tsmSLh05WJJN+lE+gjgZaU9dZgW6pxwDgz+5OZvWhml/f2g8xstpk1mVlTW1vbiVUsEhK1WyTfZepF0RgwFrgYuBb4mZkN6TnI3ee4e6O7N9bW1mZoapHsU7tFCkE6gb4JqE9Zrwu2pWoF5rt7l7uvA1aTDHiRSFC7RQpBOoG+GBhrZqPNrByYBczvMeZxkkfnmFkNyRbM2gzWKRIqtVukEPQZ6O4eB24GngRWAXPdfYWZ3WlmVwbDngR2mNlK4Dnga+6+I1tFi+SS2i1SKGLpDHL3hcDCHtvuSFl24NbgSyRSFi5Xu0UKg84UFenDgqVqt0hhUKCLHIPaLVJIFOgix6B2ixQSBbrIMajdIoVEgS5yFGq3SKFRoIschdotUmgU6CJHsWDpFs5Ru0UKiAJdpBeHrkyko3MpIAp0kV6o3SKFSIEu0ouFy9RukcKjQBfpYeOO/SxtVbtFCo8CXaQHtVukUCnQRXpQu0UKlQJdJEXLTrVbpHAp0EVS6MpEUsgU6CIp1G6RQqZAFwkcbLfo6FwKlQJdJHCw3aL+uRQqBbpIQO0WKXQKdBHUbpFoUKCLoHaLRIMCXQS1WyQaFOhS9NRukahQoEvRW6h2i0SEAl2K3gK1WyQiFOhS1NRukShRoEtRU7tFokSBLkVN7RaJEgW6FC21WyRqFOhStNRukahRoEvRWrhsCxNHqt0i0ZFWoJvZ5Wb2hpk1m9ntxxj3STNzM2vMXIkimdeycz+vte5mxjk6Opfo6DPQzawUuAe4ApgAXGtmE3oZVw3cAizKdJEimaZ2i0RROkfoU4Bmd1/r7p3Ao8DMXsZ9C/g20JHB+kSyQu0WiaJ0An0k0JKy3hpsO8TMzgPq3X3BsX6Qmc02syYza2prazvuYkUyQe0Wiap+vyhqZiXA94Hb+hrr7nPcvdHdG2tra/s7tcgJUbtFoiqdQN8E1Kes1wXbDqoGzgaeN7P1wPnAfL0wKvlK7RaJqnQCfTEw1sxGm1k5MAuYf3Cnu+929xp3b3D3BuBF4Ep3b8pKxSL9cLDdopOJJIr6DHR3jwM3A08Cq4C57r7CzO40syuzXaBIJqndIlEWS2eQuy8EFvbYdsdRxl7c/7JEsuNgu2XUMLVbJHp0pqgUDbVbJOoU6FI01G6RqFOgS9FQu0WiToEuRUHtFikGCnQpCmq3SDFQoEtRULtFioECXSJP7RYpFgp0ibzfLle7RYqDAl0ib8FStVukOCjQJdLUbpFiokCXSFO7RYqJAl0iTe0WKSYKdIkstVuk2CjQJbLUbpFio0CXyFqw7C3OHjlI7RYpGgp0iaSWnft5rWUXMyaeFnYpIjmjQJdIUrtFipECXSJJ7RYpRgp0iRy1W6RYKdAlcua/thmA6RNPCbkSkdxSoEuk7O3o4t4X1nLh2BpOH1YVdjkiOaVAl0i5/0/reXt/F7dddlbYpYjknAJdImP3/i5+9sJaLh0/gkn1Q8IuRyTnFOgSGT97YS17O+LcOm1c2KWIhEKBLpGwY98B7v/TOmZMPJUJpw0KuxyRUCjQJRJ++oe1vNPVzVenjQ27FJHQKNCl4G3b08GDf17PzEkjOXN4ddjliIRGgS4F78fPv0k84dwyVUfnUtwU6FLQNu16h4cXbeTTk+toqNH7zqW4KdCloP3o2WYc5+ZLzgy7FJHQpRXoZna5mb1hZs1mdnsv+281s5VmttTMfmdmp2e+VJHDbdjRzrymFq6dMoq6ofoQLpE+A93MSoF7gCuACcC1Zjahx7BXgEZ3Pwf4FfCdTBcq0tMPfreG0hLjSx/V0bkIpHeEPgVodve17t4JPArMTB3g7s+5+/5g9UWgLrNlihyueds+Hn9lE9effzojBlWGXY5IXkgn0EcCLSnrrcG2o7kR+G1vO8xstpk1mVlTW1tb+lWK9HD3M6upLCvlv198RtiliOSNjL4oambXAY3Ad3vb7+5z3L3R3Rtra2szObUUkVVb9vDE0i3ccEEDNQMrwi5HJG/E0hizCahPWa8Lth3GzC4FvgFc5O4HMlOeyJHueno11RUxZn9kTNiliOSVdI7QFwNjzWy0mZUDs4D5qQPM7Fzgp8CV7r4t82WKJC1t3cVTK7fyhQvHMGRAedjliOSVPgPd3ePAzcCTwCpgrruvMLM7zezKYNh3gYHAPDN71czmH+XHifTL959ezZABZXz+ww1hlyKSd9JpueDuC4GFPbbdkbJ8aYbrEjnCkg07ef6NNv7h8vdQXVkWdjkieUdnikrB+N5Tq6kZWM5nL9B5ayK9UaBLQfhz83b+/OYOvnjxmQwoT+s/liJFR4Euec/d+d7TqzllUCV/84FRYZcjkrcU6JL3fr+6jSUb3uZLl5xJZVlp2OWI5C0FuuQ1d+f7T69m5JCTuKaxvu9vECliCnTJa0+v3MrS1t3cMnUs5TE9XEWORc8QyVuJRPLovGHYAP7qvGN9fJCIgAJd8tiCZVt4/a29fHXaOGKleqiK9EXPEslL8e4Edz2zmrHDB/Lxc04LuxyRgqBAl7z0+KubWdvWzq3TxlFaYmGXI1IQFOiSd/64Zjt3/OdyJo4czMfee0rY5YgUDAW65JX/Wv4Wn39gMaNOHsB9n22kREfnImnTOdSSN+Y1tfAPv17K++qHcP8N79fH44ocJwW65IX7/riObz2xkgvH1vCT6yZTVaGHpsjx0rNGQuXu3PX0an74bDNXnH0Kd8+aREVMp/eLnAgFuoQmkXDufGIlD/x5PVc31vFvn5io95uL9IMCXULR1Z3g73+1lN+8sokvfHg035gxHjO9ACrSHwp0ybmOrm5ufvgVnlm1lf952Ti+9NEzFeYiGaBAl5zadyDOFx5czItrd3LnzPfymQ82hF2SSGQo0CVndrZ3csP9L7Fi8x7uvmYSV52rD9wSySQFuuTElt3vcP19L9Gycz9zrp/M1PEjwi5JJHIU6JJ167e38zf3LmL3O108+PkpnD9mWNgliUSSAl2yauXmPXzm5y/RnUjwyE3nM7FucNgliUSWAl2yZvH6ndz4wGKqKmI8OvuDnDm8OuySRCJNgS4Ztb8zzm+XvcXcphYWrdvJ6JoqfnHjFOqGDgi7NJHIU6BLv7k7L2/cxbymFp5YuoV9B+I0DBvA1z52Fn89ZRRDq/QhWyK5oECXE7Ztbwe/eXkTc5taeLOtnZPKSplxzqlc3VjP+xuG6mQhkRxToMtx6epO8Nzr25jb1Mpzb2yjO+FMPn0o3/7kGGaccxoD9SmJIqHRs0/SsmbrXuYtaeWxl1vZvq+T2uoKbrpwDJ+aXMeZwweGXZ6IoECXXrg7O9s7Wbe9nZVb9vCbVzbxysZdxEqMqeOHc3VjPReNq9UnI4rkGQV6EdvT0cX67e2sC74OLq/d3s7ejvihcWOHD+SbM8Zz1bkjqRlYEWLFInIsaQW6mV0O/AAoBe5193/vsb8CeAiYDOwArnH39ZktVY6Hu7O/s5t9B+LsbO9kw45kUK9ra2f9jmRwb9/XeWi8GZw2+CRG11Rx1aSRNNRUMaamitE1VZw+bIBe4BQpAH0GupmVAvcA04BWYLGZzXf3lSnDbgTedvczzWwW8G3gmmwUXKjcnXjCiXc7XYlE8rY7QVd3cjmeSNB1aJsT704QT7y7fiDezb6OOPsOxNkb3B5aPxBnX0fXoW17D8RpPxAn4UfWUVtdweiaKqa+ZwSja6toGFbFmNoqRp08gMoyXSlIpJClc4Q+BWh297UAZvYoMBNIDfSZwD8Hy78CfmRm5u69REr/zF3cwpwX1h5a7znFERP6kasHvye5fHC74/7u+sGf7YeW3x2T8OS+hDsJh0Ry52Hrfmjdew3W/qoqL2VgZYyBFTEGVpZRXRFjeHXloW3Vh/bFGHJSOacPG0BDTZXehSISYek8u0cCLSnrrcAHjjbG3eNmthsYBmxPHWRms4HZAKNGjTqhgodWlXPWiB6nkNsxV49oFxjJFsO7y/bu9xlY8BPMeo41SkqS40sMSswoCXaWHNxWYpilrJthlvyJsRIjVlpCWakRKzHKYiWUlZQQKw22lxhlpcn1stKSw8ZXxEqTIV0Zo6o8RmmJWiAicricHq65+xxgDkBjY+MJHbdOmzCCaRP00asiIj2l876zTUB9ynpdsK3XMWYWAwaTfHFURERyJJ1AXwyMNbPRZlYOzALm9xgzH/hssPwp4Nls9M9FROTo+my5BD3xm4EnSb5t8efuvsLM7gSa3H0+cB/wCzNrBnaSDH0REcmhtHro7r4QWNhj2x0pyx3ApzNbmoiIHA+duy0iEhEKdBGRiFCgi4hEhAJdRCQiLKx3F5pZG7DhBL+9hh5noeZQWHPrPkd/3jDn1n0unLlPd/fa3naEFuj9YWZN7t5YTHPrPkd/3jDn1n2OxtxquYiIRIQCXUQkIgo10OcU4dy6z9GfN8y5dZ8jMHdB9tBFRORIhXqELiIiPSjQRUQioqAC3cw+bWYrzCxhZo099n3dzJrN7A0z+1iW65hkZi+a2atm1mRmU7I5X4+5v2xmrwe/h+/kat5g7tvMzM2sJodzfje4v0vN7DdmNiTL810ePIaazez2bM6VMme9mT1nZiuDv+stuZi3Rw2lZvaKmT2RwzmHmNmvgr/vKjP7YA7n/mrwu15uZo+YWWWW5vm5mW0zs+Up2042s6fNbE1wOzRjE7p7wXwB44GzgOeBxpTtE4DXgApgNPAmUJrFOp4CrgiWpwPP5+j+fxR4BqgI1ofn8HdfT/IjlDcANTmc9zIgFix/G/h2FucqDR47Y4Dy4DE1IQf38VTgvGC5Glidi3l71HAr8DDwRA7nfBD4QrBcDgzJ0bwjgXXAScH6XOCGLM31EeA8YHnKtu8AtwfLt2fyMV1QR+juvsrd3+hl10zgUXc/4O7rgGaSF7fOWinAoGB5MLA5i3Ol+iLw7+5+AMDdt+VoXoC7gL+nl+twZ5O7P+Xu8WD1RZJXzMqWQxdEd/dO4OAF0bPK3be4+8vB8l5gFcnQyQkzqwNmAPfmcM7BJMPuPgB373T3Xbman+RHh58UXGFtAFl6Drv7H0heIyLVTJL/mBHcXpWp+Qoq0I+htwtZZ/MJ8RXgu2bWAvxv4OtZnCvVOOBCM1tkZr83s/fnYlIzmwlscvfXcjHfMXwe+G0Wf36uH0dHMLMG4FxgUQ6nvZvkP9aJHM45GmgD7g9aPfeaWVUuJnb3TSSftxuBLcBud38qF3MHRrj7lmD5LSBjF0nO6UWi02FmzwCn9LLrG+7+n/lQBzAV+Kq7/9rMriZ5lHFpDuaNAScD5wPvB+aa2RgP/u+WxXn/kWTrIyvS+Zub2TeAOPAf2aojbGY2EPg18BV335OjOT8ObHP3JWZ2cS7mDMRItiK+7O6LzOwHJNsP/yvbEwc965kk/1HZBcwzs+vc/ZfZnrsnd3czy9j/evMu0N39RIIxnQtZZ6wOM3sIOPjC1Twy+F/VPub9IvBYEOAvmVmC5If8tGVrXjObSPKB/5qZQfJ3+7KZTXH3t/o777HmTqnhBuDjwNRM/ON1DBl/HKXLzMpIhvl/uPtjuZgz8CHgSjObDlQCg8zsl+5+XZbnbQVa3f3g/0R+RTLQc+FSYJ27twGY2WPABUCuAn2rmZ3q7lvM7FQgY63TqLRc5gOzzKzCzEYDY4GXsjjfZuCiYPkSYE0W50r1OMkXRjGzcSRfSMrqp8W5+zJ3H+7uDe7eQPKJeF6mwrwvZnY5yXbAle6+P8vTpXNB9Iyz5L+U9wGr3P372Z4vlbt/3d3rgr/tLJIXeM92mBM8flrM7Kxg01RgZbbnDWwEzjezAcHvfirJ1y1yZT7w2WD5s0DmOg+5eFU5g68Yf4JkoBwAtgJPpuz7Bsl3KLxB8A6ULNbxYWAJyXdBLAIm5+j+l5M8ilgOvAxcEsLfYD25fZdLM8m+9qvB10+yPN90ku8yeZNkyycX9/HDJF9sXppyP6eH8Le9mNy+y2US0BTc78eBoTmc+1+A14Pn0i8I3jmWhXkeIdmn7wqy60ZgGPA7kgeCzwAnZ2o+nfovIhIRUWm5iIgUPQW6iEhEKNBFRCJCgS4iEhEKdBGRiFCgi4hEhAJdRCQi/j/8aTN5ql43XwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "e-VCy7S_3CfL"
      },
      "source": [
        "def sigmoid(z):\n",
        "    s = 1/(1+np.exp(-z))\n",
        "    return s\n"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d_uOQhhO3uzf"
      },
      "source": [
        "## Learning\n",
        "\n",
        "So how do we learn the __unknown__ conditional probability __p(y|x)__ of a given input image being a ship or not?\n",
        "Note that we have already modeled our conditional probability as a sigmoid: $$\\hat{p}_{w,b}(y=1|x)=\\sigma_{w,b}(x)$$\n",
        "and clearly \n",
        "$$\\hat{p}_{w,b}(y=0|x)=1-\\sigma_{w,b}(x)$$\n",
        "Where we have explicitly showed the dependence of $\\sigma$ on the weights $w$ and the bias $b$\n",
        "\n",
        "We can combine both outcomes, and drop the subscripts for simplicity keeping in mind that $\\sigma$ depends on both $w$ and $b$, to obtain our distribution\n",
        "$$\\hat{p}_{w,b}(y|x)=\\sigma(x)^y(1-\\sigma(x))^{1-y}$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EfiBT94z5d8G"
      },
      "source": [
        "### Kullback-Leibler\n",
        "\n",
        "To obtain our best model we use the Kullback-Leibler divergence formula (or relative entropy) to measure how close our distribution $\\hat{p}$ is to the __true__ but __unknown__ distribution $p$\n",
        "\n",
        "$$KL(p||\\hat{p})=\\sum_{x,y}p(y|x)\\log\\frac{p(y|x)}{\\hat{p}(y|x)}=-\\sum_{x,y}p(y|x)\\log \\hat{p}(y|x)+\\sum_{x,y}p(y|x)\\log p(y|x) >=0$$\n",
        "\n",
        "Our goal is to obtain $\\hat{p}$ as close as possible to $p$ so we minimize the KL. Since the second term does not depend on $w$ and $b$ we minimize the first term, or maximize its negative, i.e. __maximize the log-likelihood__\n",
        "$$\\mathbb{E}_{x,y\\sim p}[\\log \\hat{p}(y|x)]$$\n",
        "\n",
        "But since we don't have access to the __true__ __unknown__ distribution we replace the expectation with an expectation over the data. Let $D=\\{(x_1,y_1),\\ldots,(x_n,y_n)\\}$ be our data set then\n",
        "$$\\mathbb{E}_{x,y\\sim p}[\\log \\hat{p}(y|x)]\\approx \\frac{1}{|D|}\\sum_{x,y\\in D}\\log \\hat{p}(y|x)$$\n",
        "\n",
        "In the current setting recall that $\\hat{p}(y|x)=\\sigma_w(w)^y(1-\\sigma_w(x))^{1-y}$\n",
        "\n",
        "Therefore we obtain the expression to be minimized\n",
        "$$E=\\frac{1}{|D|}\\sum_{x,y\\in D}y\\cdot\\sigma_w(x)+(1-y)\\cdot(1-\\sigma_w(x))$$\n",
        "\n",
        "In what follows we will refer to the above expression for $E$ as the __error__\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KsWE14PLqTMe"
      },
      "source": [
        "## Gradient Descent\n",
        "\n",
        "Now that we have an expression to optimize we need a method to find the optimal parameters. Typically, one computes the gradient and the optimal value corresponds to the value  of the parameters when the gradient vanishes. Unfortunately, for logistic regression there is __no closed form solution__ so we seek a numerical method to find the optimal parameters.\n",
        "\n",
        "Our goal is to find the **optimal** values for _W_ and _b_. To do so we give them some _arbitrary_ values and then using the expression for $E$\n",
        "In the figure below we show an arbitrary function _E(w)_. For a given value of _w_ we compute the derivative (slope) of _E_ with respect to _w_ (two different values are shown). The point on the left side has a negative slope so we need to **increase** the value of _w_ to move toward the minimum whereas the point on the right side the slope is positive so we have to **decrease** the value of _w_. \n",
        "\n",
        "In general we \"update\" the values of _w_ and _b_ as follows\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "  w=w-\\alpha\\cdot \\frac{\\partial E}{\\partial w}\\\\\n",
        "  b=b-\\alpha\\cdot \\frac{\\partial E}{\\partial b}\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "where $\\alpha$ is a parameter chosen by us, called the __learning rate__."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "weMlF8yqScvo"
      },
      "source": [
        "\n",
        "<center>\n",
        "<img src=\"https://github.com/hikmatfarhat-ndu/CSC645/blob/master/figures/gradient-descent.png?raw=1\" width=\"350\">\n",
        "</center>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7SZ1EkLBkQYn"
      },
      "source": [
        "## Gradient Descent for Logistic Regression\n",
        "\n",
        "The first thing we need to do is compute the __derivatives__ of $E$ with respect to the parameters. \n",
        "\n",
        "Let $X$ and $Y$ be a sample image and its corresponding label. And let $\\hat{p}$ be the probability that the outcome is $Y$ given input $X$ our __error__ function can be written as \n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "E=-Y\\log\\hat{p}-(1-Y)\\log (1-\\hat{p})\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "Taking the derivative with respect to _w_ and _b_ we get\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "dw &=\\frac{\\partial E}{\\partial w}=X^T\\cdot (\\hat{p}-Y)/m\\\\\n",
        "db &=\\frac{\\partial E}{\\partial b}=\\sum_{i=1}^{m}(\\hat{p}_i-Y_i)/m\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "Where $T$ is the transpose, $\\cdot$ is the dot product and $m$ is the number of samples. For details on the above **please refer to the pdf lecture notes**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AaogpCm63Cfg"
      },
      "source": [
        "# Perform forward propagation\n",
        "# Then compute the cost and gradient\n",
        "def propagate(w, b, X, Y):\n",
        "# m is the number of samples which is the number of rows in X\n",
        "    m = X.shape[0]\n",
        "    p_hat= sigmoid(np.dot(X,w)+b)   # compute activation\n",
        "    cost=-np.sum(Y*np.log(p_hat)+(1-Y)*np.log(1-p_hat))/m\n",
        "    # Compute Derivatives\n",
        "    dw = np.dot(X.T,(p_hat-Y))/m\n",
        "    db = np.sum(p_hat-Y)/m\n",
        "    cost = np.squeeze(cost)\n",
        "    \n",
        "    return dw, db, cost\n"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZZtUVolZ3CfW"
      },
      "source": [
        "def initialize(dim):\n",
        "  # this sets the initial p_hat to 0.5\n",
        "    w = np.zeros((dim,1))\n",
        "    b = 0\n",
        "    \n",
        "    return w, b"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZoR7_2Rfu-uJ"
      },
      "source": [
        "learning_rate = 0.026\n",
        "print_cost = True\n",
        "\n",
        "    \n",
        "# initialize parameters with zeros \n",
        "w, b = initialize(img_train.shape[1])\n",
        "dw,db,cost=propagate(w,b,img_train,label_train)\n"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NramD97ZqIra"
      },
      "source": [
        "## Parameter initialization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Vy2onH1IDqfh"
      },
      "source": [
        "### Gradient Descent\n",
        "\n",
        "A single step of gradient descent. Update parameters _w_ and _b_ "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i-2mL8fC3Cfr"
      },
      "source": [
        "def learn(w, b, X, Y, num_iterations,learning_rate, print_cost = False):\n",
        "    \n",
        "    for i in range(num_iterations):\n",
        "        # Cost and gradient calculation \n",
        "        dw,db, cost = propagate(w,b,X,Y)\n",
        "        \n",
        "        # update rule\n",
        "        w = w-learning_rate*dw\n",
        "        b = b-learning_rate*db\n",
        "\n",
        "        if print_cost and i % 100 == 0:\n",
        "            print (\"Cost after iteration %i: %f\" %(i, cost))\n",
        "    \n",
        "    \n",
        "    return w,b,dw,db\n"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IX4KD3f53jDu",
        "outputId": "d2b63df3-a5e3-4966-f795-024bbd828c87",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 69
        }
      },
      "source": [
        "w,b,dw,db= learn(w, b, img_train, label_train, 1,0.01, True)\n",
        "print(w.shape)\n",
        "print(dw.shape)"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Cost after iteration 0: 0.693147\n",
            "(3072, 1)\n",
            "(3072, 1)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JTfzufFAEIyo"
      },
      "source": [
        "## Prediction\n",
        "\n",
        "After we finish learning convert output >=0.5 to 1  (ship) and <0.5 to 0 (not ship)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "24YgreCj3Cf1"
      },
      "source": [
        "def predict(w, b, X):\n",
        "    \n",
        "    m = X.shape[0]\n",
        "    Y_prediction = np.zeros((m,1))\n",
        "    \n",
        "    # Compute vector \"Y_hat\" predicting\n",
        "    #    the probabilities of a ship being present in the picture\n",
        "    p_hat= sigmoid(np.dot(X,w)+b)   # compute activation\n",
        "\n",
        "\n",
        "    for i in range(p_hat.shape[0]):\n",
        "        \n",
        "        # Convert probabilities Y_hat[0,i] to actual predictions p[0,i]\n",
        "        if p_hat[i]>=0.5:\n",
        "            Y_prediction[i]=1\n",
        "        else:\n",
        "            Y_prediction[i]=0\n",
        "    \n",
        "\n",
        "    return Y_prediction\n"
      ],
      "execution_count": 18,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xMZ07QQnEZDN"
      },
      "source": [
        "## Main loop\n",
        "Define the number of iterations and learning_rate then call gradient descent. After performing num_iterations of gradient descent we predict the type of all images including the ones in the test set and compare with the given labels to gauge the accuracy of our prediction."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "tags": [],
        "id": "iCNQrbHd3Cf8",
        "outputId": "cfae10e7-bc08-437a-aa76-7386f62f6121",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 225
        }
      },
      "source": [
        "num_iterations = 1000\n",
        "learning_rate = 0.026\n",
        "print_cost = True\n",
        "# initialize parameters with zeros \n",
        "w, b = initialize(img_train.shape[1])\n",
        "\n",
        "# Gradient descent \n",
        "w,b,dw,db= learn(w, b, img_train, label_train, num_iterations, learning_rate, print_cost)\n"
      ],
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Cost after iteration 0: 0.693147\n",
            "Cost after iteration 100: 0.278183\n",
            "Cost after iteration 200: 0.274052\n",
            "Cost after iteration 300: 0.271784\n",
            "Cost after iteration 400: 0.269959\n",
            "Cost after iteration 500: 0.268390\n",
            "Cost after iteration 600: 0.266999\n",
            "Cost after iteration 700: 0.265740\n",
            "Cost after iteration 800: 0.264583\n",
            "Cost after iteration 900: 0.263509\n",
            "train accuracy:89.992\n",
            "test accuracy:89.82\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N0DJBgdZOtlq",
        "outputId": "2b36c40d-25b0-4321-c83b-a1f140df22d5",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 51
        }
      },
      "source": [
        "\n",
        "\n",
        "# Predict test/train set examples\n",
        "\n",
        "\n",
        "Y_prediction_test = predict(w, b, img_test)\n",
        "Y_prediction_train = predict(w, b, img_train)\n",
        "\n",
        "# Print train/test Errors\n",
        "\n",
        "\n",
        "print(\"train accuracy:\"+str((100 - np.mean(np.abs(Y_prediction_train - label_train)) * 100)))\n",
        "print(\"test accuracy:\"+str((100 - np.mean(np.abs(Y_prediction_test - label_test)) * 100)))\n"
      ],
      "execution_count": 20,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "train accuracy:89.992\n",
            "test accuracy:89.82\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sUCWt5S33CgC"
      },
      "source": [
        ""
      ],
      "execution_count": 14,
      "outputs": []
    }
  ]
}