unit07_optim/demo2_grad_descent.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-Cjj4lAgzVqW"
      },
      "source": [
        "# Demo:  Gradient Descent Optimization\n",
        "\n",
        "In the [breast cancer demo](../logistic/breast_cancer.ipynb), we used the `sklearn` built-in `LogisticRegression` class to find the weights for the logistic regression problem.   The `fit` routine in that class has an *optimizer* to select the weights to best match the data.  To understand how that optimizer works, in this demo, we will build a very simple gradient descent optimizer from scratch.  You will learn to:\n",
        "* Compute the gradients of a simple loss function and implement the gradient calculations in python\n",
        "* Implement a simple gradient descent optimizer\n",
        "* Visualize the effect of the learning rate in gradient descent\n",
        "* Implement an adaptive learning rate algorithm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-OSBjKmozVqX"
      },
      "source": [
        "## Loading the Breast Cancer Data\n",
        "\n",
        "We first load the standard packages."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "HmgjaBgkzVqY"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "%matplotlib inline\n",
        "from sklearn import preprocessing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IT8A_xegzVqY"
      },
      "source": [
        "We next load the data from the [breast cancer demo](../logistic/breast_cancer.ipynb).  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "V2qsHqvCzVqZ",
        "outputId": "7a3ddae2-f854-4155-b685-9829990f3c60",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        }
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "        id  thick  size_unif  shape_unif  marg  cell_size  bare  chrom  \\\n",
              "0  1000025      5          1           1     1          2   1.0      3   \n",
              "1  1002945      5          4           4     5          7  10.0      3   \n",
              "2  1015425      3          1           1     1          2   2.0      3   \n",
              "3  1016277      6          8           8     1          3   4.0      3   \n",
              "4  1017023      4          1           1     3          2   1.0      3   \n",
              "\n",
              "   normal  mit  class  \n",
              "0       1    1      2  \n",
              "1       2    1      2  \n",
              "2       1    1      2  \n",
              "3       7    1      2  \n",
              "4       1    1      2  "
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-3f695593-81bb-4e80-9e76-39544745d028\" class=\"colab-df-container\">\n",
              "    <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>id</th>\n",
              "      <th>thick</th>\n",
              "      <th>size_unif</th>\n",
              "      <th>shape_unif</th>\n",
              "      <th>marg</th>\n",
              "      <th>cell_size</th>\n",
              "      <th>bare</th>\n",
              "      <th>chrom</th>\n",
              "      <th>normal</th>\n",
              "      <th>mit</th>\n",
              "      <th>class</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>1000025</td>\n",
              "      <td>5</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "      <td>1.0</td>\n",
              "      <td>3</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>1002945</td>\n",
              "      <td>5</td>\n",
              "      <td>4</td>\n",
              "      <td>4</td>\n",
              "      <td>5</td>\n",
              "      <td>7</td>\n",
              "      <td>10.0</td>\n",
              "      <td>3</td>\n",
              "      <td>2</td>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>1015425</td>\n",
              "      <td>3</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "      <td>2.0</td>\n",
              "      <td>3</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>1016277</td>\n",
              "      <td>6</td>\n",
              "      <td>8</td>\n",
              "      <td>8</td>\n",
              "      <td>1</td>\n",
              "      <td>3</td>\n",
              "      <td>4.0</td>\n",
              "      <td>3</td>\n",
              "      <td>7</td>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>1017023</td>\n",
              "      <td>4</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>3</td>\n",
              "      <td>2</td>\n",
              "      <td>1.0</td>\n",
              "      <td>3</td>\n",
              "      <td>1</td>\n",
              "      <td>1</td>\n",
              "      <td>2</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "    <div class=\"colab-df-buttons\">\n",
              "\n",
              "  <div class=\"colab-df-container\">\n",
              "    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-3f695593-81bb-4e80-9e76-39544745d028')\"\n",
              "            title=\"Convert this dataframe to an interactive table.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n",
              "    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "\n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    .colab-df-buttons div {\n",
              "      margin-bottom: 4px;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "    <script>\n",
              "      const buttonEl =\n",
              "        document.querySelector('#df-3f695593-81bb-4e80-9e76-39544745d028 button.colab-df-convert');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      async function convertToInteractive(key) {\n",
              "        const element = document.querySelector('#df-3f695593-81bb-4e80-9e76-39544745d028');\n",
              "        const dataTable =\n",
              "          await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                    [key], {});\n",
              "        if (!dataTable) return;\n",
              "\n",
              "        const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "          + ' to learn more about interactive tables.';\n",
              "        element.innerHTML = '';\n",
              "        dataTable['output_type'] = 'display_data';\n",
              "        await google.colab.output.renderOutput(dataTable, element);\n",
              "        const docLink = document.createElement('div');\n",
              "        docLink.innerHTML = docLinkHtml;\n",
              "        element.appendChild(docLink);\n",
              "      }\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "\n",
              "<div id=\"df-096d675e-cdd9-43cb-a5a7-c335d7204b4a\">\n",
              "  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-096d675e-cdd9-43cb-a5a7-c335d7204b4a')\"\n",
              "            title=\"Suggest charts\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "     width=\"24px\">\n",
              "    <g>\n",
              "        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n",
              "    </g>\n",
              "</svg>\n",
              "  </button>\n",
              "\n",
              "<style>\n",
              "  .colab-df-quickchart {\n",
              "      --bg-color: #E8F0FE;\n",
              "      --fill-color: #1967D2;\n",
              "      --hover-bg-color: #E2EBFA;\n",
              "      --hover-fill-color: #174EA6;\n",
              "      --disabled-fill-color: #AAA;\n",
              "      --disabled-bg-color: #DDD;\n",
              "  }\n",
              "\n",
              "  [theme=dark] .colab-df-quickchart {\n",
              "      --bg-color: #3B4455;\n",
              "      --fill-color: #D2E3FC;\n",
              "      --hover-bg-color: #434B5C;\n",
              "      --hover-fill-color: #FFFFFF;\n",
              "      --disabled-bg-color: #3B4455;\n",
              "      --disabled-fill-color: #666;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart {\n",
              "    background-color: var(--bg-color);\n",
              "    border: none;\n",
              "    border-radius: 50%;\n",
              "    cursor: pointer;\n",
              "    display: none;\n",
              "    fill: var(--fill-color);\n",
              "    height: 32px;\n",
              "    padding: 0;\n",
              "    width: 32px;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart:hover {\n",
              "    background-color: var(--hover-bg-color);\n",
              "    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "    fill: var(--button-hover-fill-color);\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart-complete:disabled,\n",
              "  .colab-df-quickchart-complete:disabled:hover {\n",
              "    background-color: var(--disabled-bg-color);\n",
              "    fill: var(--disabled-fill-color);\n",
              "    box-shadow: none;\n",
              "  }\n",
              "\n",
              "  .colab-df-spinner {\n",
              "    border: 2px solid var(--fill-color);\n",
              "    border-color: transparent;\n",
              "    border-bottom-color: var(--fill-color);\n",
              "    animation:\n",
              "      spin 1s steps(1) infinite;\n",
              "  }\n",
              "\n",
              "  @keyframes spin {\n",
              "    0% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "      border-left-color: var(--fill-color);\n",
              "    }\n",
              "    20% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    30% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    40% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    60% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    80% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "    90% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "  }\n",
              "</style>\n",
              "\n",
              "  <script>\n",
              "    async function quickchart(key) {\n",
              "      const quickchartButtonEl =\n",
              "        document.querySelector('#' + key + ' button');\n",
              "      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n",
              "      quickchartButtonEl.classList.add('colab-df-spinner');\n",
              "      try {\n",
              "        const charts = await google.colab.kernel.invokeFunction(\n",
              "            'suggestCharts', [key], {});\n",
              "      } catch (error) {\n",
              "        console.error('Error during call to suggestCharts:', error);\n",
              "      }\n",
              "      quickchartButtonEl.classList.remove('colab-df-spinner');\n",
              "      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n",
              "    }\n",
              "    (() => {\n",
              "      let quickchartButtonEl =\n",
              "        document.querySelector('#df-096d675e-cdd9-43cb-a5a7-c335d7204b4a button');\n",
              "      quickchartButtonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "    })();\n",
              "  </script>\n",
              "</div>\n",
              "\n",
              "    </div>\n",
              "  </div>\n"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "dataframe",
              "variable_name": "df",
              "summary": "{\n  \"name\": \"df\",\n  \"rows\": 683,\n  \"fields\": [\n    {\n      \"column\": \"id\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 620644,\n        \"min\": 63375,\n        \"max\": 13454352,\n        \"num_unique_values\": 630,\n        \"samples\": [\n          1287971,\n          385103,\n          636375\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"thick\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 2,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          10,\n          3,\n          1\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"size_unif\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 3,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          6,\n          4,\n          3\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"shape_unif\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 2,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          7,\n          4,\n          3\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"marg\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 2,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          9,\n          5,\n          4\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"cell_size\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 2,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          10,\n          7,\n          4\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"bare\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 3.6438571604929035,\n        \"min\": 1.0,\n        \"max\": 10.0,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          8.0,\n          10.0,\n          9.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"chrom\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 2,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          6,\n          9,\n          5\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"normal\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 3,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 10,\n        \"samples\": [\n          9,\n          2,\n          3\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"mit\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1,\n        \"min\": 1,\n        \"max\": 10,\n        \"num_unique_values\": 9,\n        \"samples\": [\n          8,\n          5,\n          7\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"class\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 0,\n        \"min\": 2,\n        \"max\": 4,\n        \"num_unique_values\": 2,\n        \"samples\": [\n          4,\n          2\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    }\n  ]\n}"
            }
          },
          "metadata": {},
          "execution_count": 2
        }
      ],
      "source": [
        "names = ['id','thick','size_unif','shape_unif','marg','cell_size','bare',\n",
        "         'chrom','normal','mit','class']\n",
        "df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/' +\n",
        "                 'breast-cancer-wisconsin/breast-cancer-wisconsin.data',\n",
        "                names=names,na_values='?',header=None)\n",
        "df = df.dropna()\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hy5gdnkSzVqZ"
      },
      "source": [
        "As in the breast cancer demo, we create a data matrix `X` of various features of the breast cancer sample.  The response vector `y` is a binary indicating if each sample is benign or malignant.  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZBlhGjkUzVqZ"
      },
      "outputs": [],
      "source": [
        "# Get the predictors\n",
        "xnames = names[1:-1]\n",
        "X = np.array(df[xnames])\n",
        "\n",
        "# Get the response.  Convert to a zero-one indicator\n",
        "yraw = np.array(df['class'])\n",
        "BEN_VAL = 2   # value in the 'class' label for benign samples\n",
        "MAL_VAL = 4   # value in the 'class' label for malignant samples\n",
        "y = (yraw == MAL_VAL).astype(int)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pe7f9pQfzVqZ"
      },
      "source": [
        "We want to learn the classification rule to predict `y` from `X`.  To this end, we first split the data into training and test and then scale."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jobKO_hpzVqZ"
      },
      "outputs": [],
      "source": [
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# Split into training and test\n",
        "Xtr, Xts, ytr, yts = train_test_split(X,y, test_size=0.30)\n",
        "\n",
        "# Scale data\n",
        "scal = StandardScaler()\n",
        "Xtr1 = scal.fit_transform(Xtr)\n",
        "Xts1 = scal.transform(Xts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wHQQ0BQ-zVqa"
      },
      "source": [
        "Just as a reference, let's see the test accuracy on this split using sklearn's built-in optimizer."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "sLDfb0DezVqa",
        "outputId": "33a9f080-0176-469f-e649-685ba8250efb"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Accuracy on test data = 0.960976\n"
          ]
        }
      ],
      "source": [
        "# Fit the model\n",
        "from sklearn import linear_model\n",
        "reg = linear_model.LogisticRegression(C=1e5)\n",
        "reg.fit(Xtr1, ytr)\n",
        "\n",
        "# Measure the test error\n",
        "yhat = reg.predict(Xts1)\n",
        "acc = np.mean(yhat == yts)\n",
        "print(\"Accuracy on test data = %f\" % acc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wj0-H2e3zVqa"
      },
      "source": [
        "## Computing the Gradient and Loss Function\n",
        "\n",
        "Recall that training a logistic function means finding a weight vector `w` for the classification rule:\n",
        "\n",
        "    P(y=1|x,w) = 1/(1+\\exp(-z)), z = w[0] + w[1]*x[0] + ... + w[k]x[k-1]\n",
        "    \n",
        "From the class notes in logistic regression, we saw that the weight vector can be found by minimizing  the negative log likelihood.  The negative log likelihood is called the *loss* function.  For the logistic regression problem, the loss function simplifies to:\n",
        "\n",
        "\n",
        "    f(w) = \\sum_i -y[i]z[i] + ln(1+exp(z[i])),  z[i] = w[0] + w[1]X[i,0] + ... + w[k]X[i,k-1].\n",
        "    \n",
        "We first write a function to compute `f` and its gradient `feval`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xT3OtbHfzVqa"
      },
      "outputs": [],
      "source": [
        "def feval(w,X,y):\n",
        "    \"\"\"\n",
        "    Compute the loss and gradient given w,X,y\n",
        "    \"\"\"\n",
        "    # Construct transform matrix\n",
        "    n = X.shape[0]\n",
        "    A = np.column_stack((np.ones(n,), X))\n",
        "\n",
        "    # The loss is the binary cross entropy\n",
        "    z = A.dot(w)\n",
        "    py = 1/(1+np.exp(-z))\n",
        "    f = np.sum((1-y)*z - np.log(py))\n",
        "\n",
        "    # Gradient\n",
        "    df_dz = py-y\n",
        "    fgrad = A.T.dot(df_dz)\n",
        "    return f, fgrad"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kSCCbB5jzVqa"
      },
      "source": [
        "One issue is that the function `f(w,X,y)` depends on the parameters `X` and `y`.  Most numerical optimizers expect a function that only depends on `w`.  So, we need an alternate way to pass those parameters to the function.  There are two approaches."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "otHsJTeRzVqa"
      },
      "source": [
        "### Method 1.  Use a class method\n",
        "\n",
        "One approach to pass the parameters is to implement the function as a method in a class that contains the parameters.  In the code below, we create a class `LogisticFun` where the constructor takes the parameters `X` and `y`.  Then, the class method `feval` takes a variable `w` and computes the loss and its gradient.  \n",
        "\n",
        "You can see that one benefit of using a class method is that processing that is independent of `w` can be performed in the contructor.  For example, in the code below, constructs `self.A` in the constructor.  Then this step does not need to be repeated every time the function is called."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5UnOh_oyzVqa"
      },
      "outputs": [],
      "source": [
        "class LogisticFun(object):\n",
        "    def __init__(self,X,y):\n",
        "        \"\"\"\n",
        "        Class to computes the loss and gradient for a logistic regression problem.\n",
        "\n",
        "        The constructor takes the data matrix `X` and response vector y for training.\n",
        "        \"\"\"\n",
        "        self.X = X\n",
        "        self.y = y\n",
        "        n = X.shape[0]\n",
        "        self.A = np.column_stack((np.ones(n,), X))\n",
        "\n",
        "    def feval(self,w):\n",
        "        \"\"\"\n",
        "        Compute the loss and gradient for a given weight vector\n",
        "        \"\"\"\n",
        "        # The loss is the binary cross entropy\n",
        "        z = self.A.dot(w)\n",
        "        py = 1/(1+np.exp(-z))\n",
        "        f = np.sum((1-self.y)*z - np.log(py))\n",
        "\n",
        "        # Gradient\n",
        "        df_dz = py-self.y\n",
        "        fgrad = self.A.T.dot(df_dz)\n",
        "        return f, fgrad"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QnPNGXq2zVqa"
      },
      "source": [
        "Once we have defined the class we can create an *instance* of the class with the training data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Iu_JWRduzVqa"
      },
      "outputs": [],
      "source": [
        "log_fun = LogisticFun(Xtr1,ytr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lecCFQAXzVqa"
      },
      "source": [
        "Then, we can call the method with only the value `w`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GV52mBrPzVqa"
      },
      "outputs": [],
      "source": [
        "# Some random point\n",
        "p = X.shape[1]+1\n",
        "w0 = np.random.randn(p)\n",
        "\n",
        "# Call the function\n",
        "f, fgrad = log_fun.feval(w0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0Cf_aKyHzVqa"
      },
      "source": [
        "### Method 2.  Using Lambda functions\n",
        "\n",
        "An alternative method is to use python `lambda` functions.  The code below creates uses the `lambda` construction to fix the parameters `X` and `y`.  A slight disadvantage of this method is that the pre-processing is less easy to separate."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Qg0Wc0cqzVqb"
      },
      "outputs": [],
      "source": [
        "# Create a function with X,y fixed\n",
        "feval_param = lambda w: feval(w,Xtr1,ytr)\n",
        "\n",
        "# You can now pass a parameter like w0\n",
        "f0, fgrad0 = feval_param(w0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zZ_a0NKnzVqb"
      },
      "source": [
        "### Testing the gradient\n",
        "\n",
        "Either way you create the function, it is very important to test if the gradient is correct.  This is the number one reason people's code does not work with numerical optimizers.  The simplest method is to take two points `w0` and `w1` that are close to one another and then verify that\n",
        "\n",
        "    f(w1) - f(w0) \\approx fgrad(w0).dot(w1 - w0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "AlQC_mj6zVqb",
        "outputId": "d6470703-3ce1-4815-9544-b1722ceec44a"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Actual f1-f0    =   1.2462e-04\n",
            "Predicted f1-f0 =   1.2462e-04\n"
          ]
        }
      ],
      "source": [
        "# Take a random initial point\n",
        "p = X.shape[1]+1\n",
        "w0 = np.random.randn(p)\n",
        "\n",
        "# Perturb the point\n",
        "step = 1e-6\n",
        "w1 = w0 + step*np.random.randn(p)\n",
        "\n",
        "# Measure the function and gradient at w0 and w1\n",
        "f0, fgrad0 = log_fun.feval(w0)\n",
        "f1, fgrad1 = log_fun.feval(w1)\n",
        "\n",
        "# Predict the amount the function should have changed based on the gradient\n",
        "df_est = fgrad0.dot(w1-w0)\n",
        "\n",
        "# Print the two values to see if they are close\n",
        "print(\"Actual f1-f0    = %12.4e\" % (f1-f0))\n",
        "print(\"Predicted f1-f0 = %12.4e\" % df_est)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tpTSLqwczVqc"
      },
      "source": [
        "You can see that the two agree well.  You would have got the same answer had you used the `feval` constructed in method 2."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l4QEqUjnzVqc"
      },
      "source": [
        "## A Simple Gradient Descent Optimizer\n",
        "\n",
        "Now, we build a \"barebones\" gradient descent optimizer function with a fixed learning rate."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SJnwDPM_zVqc"
      },
      "outputs": [],
      "source": [
        "def grad_opt_sim(feval, winit, lr=1e-3,nit=1000):\n",
        "    \"\"\"\n",
        "    Simple gradient descent optimization\n",
        "\n",
        "    feval:  A function that returns f, fgrad, the objective\n",
        "            function and its gradient\n",
        "    winit:  Initial estimate\n",
        "    lr:     learning rate\n",
        "    nit:    Number of iterations\n",
        "    \"\"\"\n",
        "    # Initialize\n",
        "    w0 = winit\n",
        "\n",
        "    # Loop over iterations\n",
        "    for it in range(nit):\n",
        "\n",
        "        # Evaluate the function and gradient\n",
        "        f0, fgrad0 = feval(w0)\n",
        "\n",
        "        # Take a gradient step\n",
        "        w0 = w0 - lr*fgrad0\n",
        "\n",
        "    return w0, f0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RBf-yPdIzVqd"
      },
      "source": [
        "The above code will work, but to get a bit more insight, we can modify this code to add a history dictionary that keeps track of the loss function and values over iteration."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HYvpl8ezzVqd"
      },
      "outputs": [],
      "source": [
        "def grad_opt_simp(feval, winit, lr=1e-3,nit=1000):\n",
        "    \"\"\"\n",
        "    Simple gradient descent optimization\n",
        "\n",
        "    feval:  A function that returns f, fgrad, the objective\n",
        "            function and its gradient\n",
        "    winit:  Initial estimate\n",
        "    lr:     learning rate\n",
        "    nit:    Number of iterations\n",
        "    \"\"\"\n",
        "    # Initialize\n",
        "    w0 = winit\n",
        "\n",
        "    # Create history dictionary for tracking progress per iteration.\n",
        "    # This isn't necessary if you just want the final answer, but it\n",
        "    # is useful for debugging\n",
        "    hist = {'w': [], 'f': []}\n",
        "\n",
        "    # Loop over iterations\n",
        "    for it in range(nit):\n",
        "\n",
        "        # Evaluate the function and gradient\n",
        "        f0, fgrad0 = feval(w0)\n",
        "\n",
        "        # Take a gradient step\n",
        "        w0 = w0 - lr*fgrad0\n",
        "\n",
        "         # Save history\n",
        "        hist['f'].append(f0)\n",
        "        hist['w'].append(w0)\n",
        "\n",
        "    # Convert to numpy arrays\n",
        "    for elem in ('f', 'w'):\n",
        "        hist[elem] = np.array(hist[elem])\n",
        "    return w0, f0, hist"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0gBrgqzczVqd"
      },
      "source": [
        "We now run the gradient descent starting from a random initial condition."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CoF0PXkrzVqd",
        "outputId": "10868283-bbf6-4b19-8455-278a68be4111"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Text(0, 0.5, 'Training loss')"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3iV5f3H8fc3e0JCCAFCWLJkj4AoSoOodaN1gFVx0FK3Vlur1c6fbbXWjYu6sFpxVcGtRSKKimxEZgCRsDeEsBLu3x95SCOGcDJOnpNzPq/rOtc5537G+Yabk0+edT/mnENERAQgyu8CREQkdCgURESknEJBRETKKRRERKScQkFERMopFEREpFyM3wXURtOmTV3btm1rtOyuXbtITk6u24Kk1tQvoUd9Eppq0y8zZ87c5JzLrGxagw6Ftm3bMmPGjBotm5+fT15eXt0WJLWmfgk96pPQVJt+MbOVh5um3UciIlJOoSAiIuUUCiIiUk6hICIi5RQKIiJSTqEgIiLlIjIUVm0pZt7GEkoPaNhwEZGKIjIUXpmxivtn7mXw3yfz4H+XsHrbbr9LEhEJCQ364rWauv7EjhzYsop5u5J5aNJSHpq0lMEdM7loQA4ndskiLiYis1JEJDJDIS4miv7NY/h13jGs2lLMqzMLeXXGKq56YRZNU+I4r28rLuyfw1GZKX6XKiJSryIyFCrKaZLEzSd34sahHZmyZCMvffUdT322gienLGdAuyaM6J/D6T1akBAb7XepIiJBF/GhcFB0lDGkSzOGdGnGhp17eH3mal6e/h03vzKXP0z8hnP7ZDO8fw7dWjb2u1QRkaBRKFSiWWoCV+cdxVU/as+Xy7fw8vTvGD99Fc9/sZIe2Y0ZMSCHs3u1JDUh1u9SRUTqlEKhCmbGsUdlcOxRGfyxeB9vzl7N+OmruOON+dz19kLO6NmCiwbk0Ld1Ombmd7kiIrWmUAhQWlIclw9qx2XHtWVu4XZenv4dE+es4bWZhXRolsKI/jmc2yebjJR4v0sVEakxhUI1mRm9c9LonZPGnWd05e15axg/fRV3vbOQe95fxCndmjOifw6DjmpKVJS2HkSkYVEo1EJyfAzD+7dmeP/WLFq3g5enr+KN2at5Z95aWqUnMjw3h/NzW9GicaLfpYqIBERXadWRLs0b8YezuvHl7UN5aERvWjdJ4r6PljDo7o8Z9dx0PvxmHftLD/hdpohIlbSlUMcSYqMZ1jubYb2zWbl5Fy9PX8WrMwuZtGgDmanxnN+vFZcd25bmjRP8LlVE5Ae0pRBEbTKSufXULnxx24n8c2QuPbMb8+Qnyxh6Xz5Pf7aCEm05iEiIUSjUg5joKE7umsXTl/dn8q/y6N+uCf/39gKGPTqVOau2+V2eiEi5oIaCmaWZ2WtmtsjMFprZsWbWxMw+MrOl3nO6N6+Z2cNmVmBm88ysbzBr80ubjGSevbw/j13cl01Fezn3san87s35bN+93+/SRESCvqXwEPC+c64L0AtYCNwGTHLOdQQmee8BTgM6eo/RwONBrs03ZsbpPVrw35t/xOXHteXFaSsZet8nTJizGud0jwcR8U/QQsHMGgGDgacBnHP7nHPbgGHAOG+2ccA53uthwPOuzJdAmpm1CFZ9oSA1IZY/nNWNidcdT8u0BG4cP4eRz3zFik27/C5NRCJUMLcU2gMbgWfNbLaZPWVmyUCWc24tgPfczJs/G1hVYflCry3sdc9uzBvXDOLPw7ox57tt/PjBKTz036XsLSn1uzQRiTAWrN0VZpYLfAkMcs5NM7OHgB3A9c65tArzbXXOpZvZO8DfnHOfee2TgFudczMPWe9oynYvkZWV1W/8+PE1qq+oqIiUlNC7X8K2PQd4adE+pq0rpXmSMbJbPF0zImfY7lDtl0imPglNtemXIUOGzHTO5VY2LZih0Bz40jnX1nt/AmXHDzoAec65td7uoXznXGcze9J7/ZI3/+KD8x3uM3Jzc92MGTNqVF9+fj55eXk1WrY+TFmykd9NmM/KzcWc07sld5zRlczU8B9XKdT7JRKpT0JTbfrFzA4bCkHbfeScWwesMrPOXtNQYAEwEbjMa7sMmOC9ngiM9M5CGghsryoQwt3gTpl8cNNgbhjakXe/XseJ9+XzwpcrOXBAB6JFJHiCfUXz9cCLZhYHLAeuoCyIXjGzUcB3wAXevO8CpwMFQLE3b0RLiI3m5pM7Max3S3735nzufHM+r80s5K/n9qBry0Z+lyciYSiooeCcmwNUtokytJJ5HXBtMOtpqI7KTOHFnx3DhDlruOudBZw15jOuOK4tN5/SiaQ4jVQiInVHVzQ3EGbGOX2ymXRzHsP75/DUZys499HPWblZp6+KSN1RKDQwjZNi+eu5PfjXqAGs27GHs8dM5ZMlG/0uS0TChEKhgTqhYyZvXXc8LRoncMWzX/HkJ8t0NbSI1JpCoQFrnZHEf645jtO6t+Bv7y3ihvFzKN5X4ndZItKAKRQauKS4GMb8tA+/ObULb89bw3mPf8GqLcV+lyUiDZRCIQyYGVfnHcUzl/dn9dZizhrzGVMLNvldlog0QAqFMDKkczMmXnc8zVLjufTpaTz16XIdZxCRalEohJm2TZP5zzWDOKVrc+56ZyG/fHkOu/dpYD0RCYxCIQylxMfw2MV9+dUpnZgwdw3nP/E5hVt1nEFEjkyhEKaioozrTuzIUyNz+W5zMWePmcoXyzb7XZaIhDiFQpgbenQWb143iPSkWC55ehoT5qz2uyQRCWEKhQhwVGYKb147iH5t0rnp5Tm8Mn3VkRcSkYikUIgQqQmxjLtiACd0zOTW1+fx7NQVfpckIiFIoRBBEuOi+efIfpzSNYs/vbWARycX+F2SiIQYhUKEiY+J5tGL+zKsd0vu/WAx936wSNcyiEg5DcYfgWKjo7j/wt4kxkbz6ORlFO8r5fdndsXM/C5NRHymUIhQ0VHG337Sg8S4aJ6d+i2795Xyl3N7EB2lYBCJZAqFCGZm/P7MriTHxTBmcgG795fyjwt6ERutvYoikUqhEOHMjF/9uDNJ8dH8/f3FFO8rZcxP+xAfE+13aSLiA/1JKABck9eBP57VlY8WrOfnz89kb4nGSxKJRAoFKXf5oHbc/ZMeTFmykRtfmkNJ6QG/SxKReqZQkO8ZMaA1vzuzK+9/s4473piv01VFIoyOKcgPjDq+HduL9/HwxwU0Torl9tO66HRVkQgR1FAws2+BnUApUOKcyzWzJsDLQFvgW+BC59xWK/ut8xBwOlAMXO6cmxXM+uTwfnlyJ7bt3s/YKctJS4rlmrwOfpckIvWgPnYfDXHO9XbO5XrvbwMmOec6ApO89wCnAR29x2jg8XqoTQ7DzPjjWd0Y1rslf39/MS9OW+l3SSJSD/w4pjAMGOe9HgecU6H9eVfmSyDNzFr4UJ94oqKMf1zQixO7NOPON+fz1tw1fpckIkEW7GMKDvjQzBzwpHNuLJDlnFsL4Jxba2bNvHmzgYpjOhd6bWsrrtDMRlO2JUFWVhb5+fk1KqyoqKjGy0aa4TmOwvVR3DR+NiuWLKBnZvD+26hfQo/6JDQFq1+CHQqDnHNrvF/8H5nZoirmrexI5g9OffGCZSxAbm6uy8vLq1Fh+fn51HTZSHTsoP2MePJLHptXxAuj+pLbtklQPkf9EnrUJ6EpWP0S1N1Hzrk13vMG4A1gALD+4G4h73mDN3shkFNh8VaA9leEiEYJsTw/agAtGidy5XPTKdhQ5HdJIhIEQQsFM0s2s9SDr4FTgPnAROAyb7bLgAne64nASCszENh+cDeThIamKfE8f+UA4mKiuPK56WzZtc/vkkSkjgVzSyEL+MzM5gJfAe84594H7gZONrOlwMnee4B3geVAAfBP4Jog1iY1lNMkibEjc1m3Yw+/+NcMDYchEmaCdkzBObcc6FVJ+2ZgaCXtDrg2WPVI3enbOp37LujF9S/N5rbXv+b+C3vp4jaRMKErmqVGzurVkpWbd/GPD5fQrmkyNwzt6HdJIlIHFApSY9cO6cDyTbu4/6MltG2azNm9WvpdkojUkgbEkxozK7t724C2TfjVq3OZuXKr3yWJSC0pFKRW4mOiefLSfrRsnMDo52ewakux3yWJSC0oFKTW0pPjePry/pQccIwaN52ivSV+lyQiNaRQkDpxVGYKj13cl4INRfz61bm6D4NIA6VQkDozqENTbj/taN6bv47HP1nmdzkiUgMKBalTPzuhHWf3asm9Hywmf/GGIy8gIiFFoSB1ysy457yedGneiBtems3Kzbv8LklEqkGhIHUuMS6asZf2IyrKGP38THbpwLNIg6FQkKDIaZLEIxf1YemGndz62jwdeBZpIBQKEjQndMzkN6d24Z2v1/LklOV+lyMiAVAoSFCNHtyeM3q24N4PFjNt+Wa/yxGRI1AoSFAdPPDcpkkS1780m4079/pdkohUoVqh4N0AJzlYxUh4SomP4bFL+rJjz35uHD+b0gM6viASqo4YCmb2vJk1MrMk4BtghZndHPzSJJx0ad6I/xvWnc+XbebB/y7xuxwROYxAthR6OOd2AOcAH1J27+TLg1mUhKcLcnO4MLcVj3xcwGRd2CYSkgIJhTgziwGGAW865/YBB4JbloSrPw/rTpfmqfzy5Tms3rbb73JE5BCBhMJTwHdAOvCJmbUGioJalYSthNhoHr+kHyWljuv+PYv9pfr7QiSUHDEUnHMPOOdaOudO8e6jvAo4MfilSbhq1zSZu8/rwezvtvHARzq+IBJKAjnQfJ2ZNfJePwlMA04IdmES3s7s2ZIR/XN4/JNlTC3Y5Hc5IuIJZPfRaOfcDjM7BcgGrgb+HtyyJBL8/qyutG+azC9fnsPmIl2/IBIKAgmFgyeVnwY865ybGeByIlVKiovhkYv6sm33fn6t8ZFEQkIgv9znmtm7wFnAe2aWwv+C4ojMLNrMZpvZ2977dmY2zcyWmtnLZhbntcd77wu86W2r/+NIQ9O1ZSPuOP1oPl60gWenfut3OSIRL5BQuAL4IzDAOVcMJACjqvEZNwILK7y/B3jAOdcR2FphXaOArc65DsAD3nwSAUYe24aTjs7i7vcWsXJHqd/liES0QM4+KgWaArea2d1Af+fc7EBWbmatgDMoO60VMzPKzlx6zZtlHGUXxUHZdRDjvNevAUO9+SXMmRn3nt+TJslxPD5nr+6/IOKjmCPNYGZ/AQYB//aafm1mg5xzdwaw/geBW4FU730GsM05d/BbX0jZwWu851UAzrkSM9vuzf+9U1PMbDQwGiArK4v8/PwAyvihoqKiGi8rwXFZZ/j79AP8YuwkftYj3u9yxKPvSmgKVr8cMRQoO5bQ9+AvcjN7BpgFVBkKZnYmsME5N9PM8g42VzKrC2Da/xqcGwuMBcjNzXV5eXmHzhKQ/Px8arqsBEcesGDLh7y1bD8XnNCdYb2zj7SI1AN9V0JTsPol0LOIUg/zuiqDgLPN7FtgPGW7jR4E0rxhM6BsHKU13utCIAfAm94Y2BLgZ0mYOOeoWPq1SeeON+bz3eZiv8sRiTiBhMLfgVlm9pSZPQ3MIICDwM65251zrZxzbYERwMfOuYuBycD53myXARO81xO993jTP3Y6RzHiREcZD43ojRlcP362hsEQqWeBHGh+ATgeeNd7DHbOvViLz/wNcLOZFVB2zOBpr/1pIMNrvxm4rRafIQ1Yq/Qk7jmvJ3NXbeO+DzUMhkh9OuwxBTPreUhTgfecYWYZzrl5gX6Icy4fyPdeLwcGVDLPHuCCQNcp4e30Hi24aEAOT05ZRl7nTAa2z/C7JJGIUNWB5kermOaAwXVci8j33HlGV75YtplbXpnLezedQKOEWL9LEgl7hw0F55wGvRNfJcfH8MDw3pz/xBf8YcI3PDC8t98liYQ9jWEkIa1P63SuP7EDb8xezVtz1xx5ARGpFYWChLzrhnSgd04ad7zxNWu3625tIsGkUJCQFxMdxQPDe1NywHHLK3M5cEBnKosESyA32elZyaONmSlQpN60a5rM78/syufLNvPM1BV+lyMStgL5xf40MBN4HvgXZRevvQEsNbOhQaxN5HuG98/hpKOz+Pv7i1m0boff5YiEpUBCYSnQzznX2znXC+gHzAF+DNwXzOJEKjIz7jmvB40SY7lp/Bz27Ncw2yJ1LZBQOLrihWrOua8pGyCvoIplRIIiIyWev5/fg0XrdnLfh4v9Lkck7AQSCsvM7BEzG+Q9HgYKzCwe0MD3Uu9O7JLFJQNb89RnK/i8YNORFxCRgAUSCiMpG8H0NuB2ykY1vYyyQNAxBfHFHad3pV1GMre8Opftxfv9LkckbAQyIF6xc+4e59xZzrkznXN3O+d2OedKnXPb66NIkUMlxkXz4IjebNy5l99NmO93OSJhI5BTUgea2XtmtsDMlhx81EdxIlXp2SqNG4d2ZOLcNUyYs9rvckTCQiB3XnuWsltqzgR0uoeElKvzjiJ/yUbufHM+uW2bkJ2W6HdJIg1aIMcUdjjn3nLOrXHOrT/4CHplIgGIiY7igQt7c+CA4+aX51Cqq51FaiWQUPjYzP5mZv0rXtUc9MpEAtQ6I4k/nN2NaSu28NSny/0uR6RBC2T30fGHPIPupyAh5oJ+rfh44Qb+8eFiTuiYSdeWjfwuSaRBCuTsoxMqeSgQJKSYGX/9SQ/SkuK46eXZutpZpIYOGwpmdpH3fENlj/orUSQwTZLj+McFvViyvoh73l/kdzkiDVJVWwrp3nPmYR4iIedHnTK5/Li2PDv1W6Ys2eh3OSINTlW343zMe/5d/ZUjUnu3ndaFz5dt4pZX5/LBTYNpkhznd0kiDUYgF681NbNbzewxMxt78FEfxYnUREJsNA8O78P24v385vV5OKfTVEUCFcgpqROALOAzYFKFh0jI6tqyEbee2pmPFqxn/PRVfpcj0mAEckpqsnPuluqu2MwSgClAvPc5rznn/mBm7YDxQBNgFnCpc26fN+rq85Tdr2EzMNw59211P1fkoCsHtSN/8Ub+/NYCjmnXhPaZKX6XJBLyAtlSeM/MTqnBuvcCJ3o35ukNnGpmA4F7gAeccx2BrcAob/5RwFbnXAfgAW8+kRqLijL+cUEv4mOj+OXLc9hfesDvkkRCXiChcBXwvpkVmdkWM9tqZluOtJArU+S9jfUeDjgReM1rHwec470e5r3Hmz7UzCzAn0OkUs0bJ/C3c3swt3A7D/13qd/liIS8QHYfNa3pys0smrKB9DoAjwLLgG3OuYM35ykEsr3X2cAqAOdciZltBzKATYesczQwGiArK4v8/Pwa1VZUVFTjZSV4gtEvicAJ2TE8OrmARsWFdEqPrtP1hzt9V0JTsPrlsKFgZh2dc0uBboeZZd5h2ss550qB3maWBrwBHF3ZbAc/soppFdc5FhgLkJub6/Ly8o5URqXy8/Op6bISPMHql9xjSzjj4U8Zt9jx3k2DaJQQW+efEa70XQlNweqXqnYf3eY9P1rJY0x1PsQ5tw3IBwYCaWZ2MIxaUXYnNyjbasgB8KY3Bo64m0okECnxMTwwvDfrduzhDxO+8bsckZB12FBwzo3ynms09pGZZXpbCJhZInASsBCYDJzvzXYZZae8Akz03uNN/9jpBHOpQ31bp3P9iR14Y/ZqJs5dc+QFRCJQIMcUMLMuQFcg4WCbc+7fR1isBTDOO64QBbzinHvbzBYA483sLmA28LQ3/9PAv8ysgLIthBHV+klEAnDdkA5MWbKRO974mn5t0nVTHpFDHDEUzOxO4BSgC/AB8GPKLmSrMhScc/OAPpW0LwcGVNK+B7ggoKpFaigmOooHh/fhtIemcPPLc/j3zwcSHaWT3EQOCuSU1OHAEGCtc+5SoBcBbmGIhKLWGUn8aVh3pq3YwpNTlvldjkhICSQUdntnEZWYWSqwDmgf3LJEguu8vtmc0bMF93+4hFnfbfW7HJGQEUgozPYOGD8DzAC+omx4CpEGy8z467k9aN44gRtems323fv9LkkkJFQZCt4VxX90zm1zzj0KnAH8wjk3sl6qEwmixomxPHJRH9Zt38Pt/9FoqiJwhFDwTgl9u8L7AuecthIkbPRpnc6vf9yZd79ex4vTvvO7HBHfBbL76Csz6xv0SkR88vMT2jO4UyZ/fnsBC9fu8LscEV9VdY/mg2cYHU9ZMCw2s1lmNtvMtLUgYSMqyrj/wl6kJcZy3b9nUbyv5MgLiYSpqrYUvvKezwE6A6dTdh3B+eh6AgkzTVPieXB4b5Zv2sUfJ2oYDIlcVV1vYADOOZ3ILRHhuA5NuW5IBx75uIBBHZoyrHf2kRcSCTNVhUKmmd18uInOufuDUI+Ir24c2pEvl2/mt//5mh7ZjXW3Nok4Ve0+igZSgNTDPETCTkx0FA+N6ENcTBRXv6DjCxJ5qtpSWOuc+3O9VSISIlqmJfLwRX0Y+cxX3PHGfO6/sBe6CaBEiqq2FPQtkIh1QsdMfnlSJ96YvVrXL0hEqSoUhtZbFSIh6LohHcjrnMmf31rAnFXb/C5HpF5UdZMd3fVMIlpUlPHg8N5kpsZz7Yuz2LJrn98liQRdIFc0i0SstKQ4nrikHxt37uXG8bMpPaDxkSS8KRREjqBHq8b8aVg3Pl26iYcmLfW7HJGgUiiIBGBE/xzO79eKhyct5f356/wuRyRoFAoiATAz7jqnO71y0rj5lTksWqeB8yQ8KRREApQQG83YS/uREh/Dz5+foQPPEpYUCiLVkNUogScv7cf6HXu59sVZ7C894HdJInVKoSBSTX1ap/O3c3vwxfLN3PX2Ar/LEalTQQsFM8sxs8lmttDMvjGzG732Jmb2kZkt9Z7TvXYzs4fNrMDM5unGPhLKzuvXip8d345xX6zk37riWcJIMLcUSoBbnHNHAwOBa82sK3AbMMk51xGY5L0HOA3o6D1GA48HsTaRWrvttC4M7pTJ7yfM57Olm/wuR6ROBC0UnHNrD97P2Tm3E1gIZAPDgHHebOMou4kPXvvzrsyXQJqZtQhWfSK1FRMdxZif9uGozBSufmEmi9ft9LskkVqrl2MKZtYW6ANMA7Kcc2uhLDiAZt5s2cCqCosVem0iIatRQizPXNGfxLhornxuOht27PG7JJFaqWro7DphZinA68BNzrkdVQxBXNmEH4wpYGajKdu9RFZWFvn5+TWqq6ioqMbLSvA01H65prvxt692c8GYydw+IIGEmPAZZLih9km4C1a/BDUUzCyWskB40Tn3H695vZm1cM6t9XYPbfDaC4GcCou3AtYcuk7n3FhgLEBubq7Ly8urUW35+fnUdFkJnobcL607r+dn42bwamEKY0fmEh0VHsHQkPsknAWrX4J59pEBTwMLD7l150TgMu/1ZcCECu0jvbOQBgLbD+5mEmkITuySxZ/O7sakRRu4882vcU6D50nDE8wthUHApcDXZjbHa/stcDfwipmNAr4DLvCmvQucDhQAxcAVQaxNJCguPbYt63fsZczkAtKT4rj11C5+lyRSLUELBefcZxz+7m0/uIGPK/uz6tpg1SNSX245pRNbivfxWP4y0pPi+Png9n6XJBKwoB9oFok0Zsb/DevO9uL9/OXdhTROiuXC3JwjLygSAhQKIkEQHWXcP7wXO/bs57bX55EaH8NpPXTZjYQ+jX0kEiTxMdE8cUk/euekcf1Ls/ngG92HQUKfQkEkiJLjYxh35QB6tGrMtS/O4kMFg4Q4hYJIkKUmxDLuygF0y27Mtf+exX8XrPe7JJHDUiiI1INGCbE8f+UAurZoxDUvzuLjRQoGCU0KBZF60jgxludHHUPn5qlc9a9ZfKQtBglBCgWRetQ4MZYXRh3D0S0bcdULM3l9ZqHfJYl8j0JBpJ41TorlxZ8dw8D2Tbjl1bk889kKv0sSKadQEPFBSnwMz1zen1O7NefPby/g/o+WaKwkCQkKBRGfxMdEM+anfbgwtxUPT1rK7ybMp6T0gN9lSYTTFc0iPoqJjuKe83qSnhzHk58sp3Drbsb8tC8p8fpqij+0pSDiMzPj9tOO5i/ndufTpZs4//HPWbt9t99lSYRSKIiEiIuPacMzl/encOtuho2ZyteF2/0uSSKQQkEkhPyoUyavX30csdFRnP/E5zplVeqdQkEkxHRunsqE6wbRp3Uat7w6l9+9OZ99JToALfVDoSASgpqmxPPCqGMYPbg9//pyJSPGfsG67Xv8LksigEJBJETFREfx29OP5tGf9mXRup2c+cinfLp0o99lSZhTKIiEuDN6tmDCtYNIS4rj0qe/4i/vLGBvSanfZUmYUiiINAAds1J567rjuWRga/756Qp+8tjnFGwo8rssCUMKBZEGIjEumrvO6cE/R+ayZttuznzkU56buoIDBzQ8htQdhYJIA3Ny1yw+uGkwx7TL4I9vLeDCJ7/QVoPUGYWCSAPUrFECz13Rn/su6MXSDUWc/vCnPDq5gP0aO0lqKWihYGbPmNkGM5tfoa2JmX1kZku953Sv3czsYTMrMLN5ZtY3WHWJhAsz47x+rfjvzT/ipKObce8Hizl7zFS+WrHF79KkAQvmlsJzwKmHtN0GTHLOdQQmee8BTgM6eo/RwONBrEskrGSmxvPYxf144pJ+bC/ex4VPfsENL83W+ElSI0ELBefcFODQP1mGAeO81+OAcyq0P+/KfAmkmVmLYNUmEo5O7d6cSbfkccPQjrz/zTpO/McnPDq5gD37dfqqBM6CeWMPM2sLvO2c6+693+acS6swfatzLt3M3gbuds595rVPAn7jnJtRyTpHU7Y1QVZWVr/x48fXqLaioiJSUlJqtKwEj/qlbmwsPsD4xfuYub6UjATjnA6xHNcyhugoq/a61CehqTb9MmTIkJnOudzKpoXKoO2V/U+tNK2cc2OBsQC5ubkuLy+vRh+Yn59PTZeV4FG/1J0LTofPCzZx9/uLeHr+dvLXx/LrH3fmx92aYxZ4OKhPQlOw+qW+zz5af3C3kPe8wWsvBHIqzNcKWFPPtYmEneM6NGXCtYN44pK+mBlXvTCLs8dM5b2v1+r6BqlUfYfCROAy7/VlwIQK7SO9s5AGAtudc2vruTaRsGRmnNq9BR/cNJh7z+/Jzj37ufrFWZz0wOVU6okAAApDSURBVCe8MmOVRmCV7wnmKakvAV8Anc2s0MxGAXcDJ5vZUuBk7z3Au8ByoAD4J3BNsOoSiVTRUcYFuTlMuiWPRy7qQ3xMNLe+No+8eyfz1KfL2b57v98lSggI2jEF59xFh5k0tJJ5HXBtsGoRkf+JjjLO6tWSM3u2IH/JRh6fvIy73lnIfR8u4Zw+2Yw8tg1Ht2jkd5nik1A50Cwi9czMGNK5GUM6N2P+6u3864uV/GdWIS999R0D2jbh4oGtOaVrc7/LlHqmUBARumc35p7ze3L76V14dUYh//pyJTeOn0NqfAx9MiG57RZy26RX66wlaZgUCiJSLi0pjp8Pbs+o49sxbcUWXptZyNtzC5nyxBe0yUjinN7ZnNGzBR2bpSggwpRCQUR+ICrKOPaoDI49KoOTm2xhV3pHXp9VyMMfL+WhSUtpn5nM6d1bcGr35nRr2UgBEUYUCiJSpYQY49R+rTivXys27NjDB9+s473563gsv4Axkwto3SSJE7s0I69zJgPbZ5AQG+13yVILCgURCVizRglcemxbLj22LZuL9vLRgvW8/806XvrqO577/FviY6IY2D6DvM6Z/KhTJu2aJmsrooFRKIhIjWSkxDNiQGtGDGjNnv2lfLl8M/mLN/LJko386a0FAGQ1iueYdhkMaNeEge2bcFSmjkWEOoWCiNRaQmw0eZ2bkde5GQArN+9iytJNTFu+mS+Wb2bi3LJRa5qmxDGgXRP6tk6nV04a3Vo2IilOv4ZCiXpDROpcm4xkLs1I5tKBbXDO8e3mYqYt38y0FVv4asUW3v16HQBRBp2yUumdk0bPVmn0bNWYDs1SdFzCRwoFEQkqM6Nd02TaNU1mxIDWAGzYuYd5q7Yzr3Abcwq38/436xg/fRVQdsV1u6bJdG6eSpesVLq0aESX5qlkpyUSVYOhv6V6FAoiUu+apSZwUtcETuqaBYBzju+2FPP16u0sXreThWt3Mq9wG+/M+9+4mMlx0bTLTKZd0xQvZJLKXzdOjPXrRwk7CgUR8Z2Z0SYjmTYZyZzZ83/tRXtLWLxuJ4vX7WTJ+p0s37SLOau28s68NVQc+btJchxtM5LITk8iOy2R7LQEstMTyU5LomVaAqkJCo1AKRREJGSlxMfQr006/dqkf699b0kpq7YUs2JTMSs2FbFiUzErN+9iXuE2Ppi/jn2l3x8OvFFCDNnpSTRvFE+z1AQyU+Np1iiezJSy54NtOpahUBCRBig+JpoOzVLp0CwVyPretAMHHJuK9lK4bTert+5m9bbdrPFer9uxh2/W7GBT0V4qu8dQakIMmSnxpCfHkZ4US1pS2XPZ+4ptZa8bJ8USHxNeQaJQEJGwEhVlNGuUQLNGCfRtnV7pPKUHHFt27WPDzj1s2LmXjYc8thbvY/W2sgDZsmsfe6u4EVFcTBSp8TGkJMSQEh9DakIMKfGxpCYcfF02LTUhltT4GBLjokmMjSYpLpoE7zkxLpqk2BgS4qKIi47y9VoOhYKIRJzoKCMzNZ7M1Hi6BTD/7n2lbC3ex9bifWwr3u+93s/24n0U7S1l5579FO0toWhPCTv3lLB62+7ytp17Siitxq1Po6OMxNjo8vCo+DohNor4mLLnTrGl5NX4X+DwFAoiIkeQGBdNYlwiLdMSq72sc449+w+wc+9+du4pYfe+UnbvL/3Bc/G+UvbsL6V4Xwm79x1g9/6yeYsrzFe0t4S9JaXsLTlAVnZw7rGtUBARCSIz80Ilmmapdbfe/Pz8ultZBUG7R7OIiDQ8CgURESmnUBARkXIKBRERKadQEBGRcgoFEREpp1AQEZFyCgURESlnzgXnqrj6YGYbgZWHNDcGtlcy+6HtTYFNQSrtSA5XY7DXE+j8R5qvqumB/vsfrs2vfvGrT6qzTE37pbbt+q7UfL5Q/a60cc5lVjrFORdWD2BsIO3AjFCrMdjrCXT+I81X1fRA//2raPOlX/zqk/rol9q267tS931S3X6pz+9KOO4+equa7X6oq1qqu55A5z/SfFVNr86/v/qkesvUtF/qqt0P+q4E9jl1pkHvPqoNM5vhnMv1uw75PvVL6FGfhKZg9Us4bikEaqzfBUil1C+hR30SmoLSLxG7pSAiIj8UyVsKIiJyCIWCiIiUUyiIiEg5hYLHzJLNbJyZ/dPMLva7HgEza29mT5vZa37XIv9jZud435MJZnaK3/UImNnRZvaEmb1mZlfXZl1hHQpm9oyZbTCz+Ye0n2pmi82swMxu85p/ArzmnPs5cHa9FxshqtMnzrnlzrlR/lQaWarZL29635PLgeE+lBsRqtknC51zVwEXArU6TTWsQwF4Dji1YoOZRQOPAqcBXYGLzKwr0ApY5c1WWo81RprnCLxPpP48R/X75U5vugTHc1SjT8zsbOAzYFJtPjSsQ8E5NwXYckjzAKDA+yt0HzAeGAYUUhYMEOb/Ln6qZp9IPalOv1iZe4D3nHOz6rvWSFHd74pzbqJz7jigVru/I/GXXzb/2yKAsjDIBv4DnGdmjxNal/lHgkr7xMwyzOwJoI+Z3e5PaRHtcN+V64GTgPPN7Co/Cotgh/uu5JnZw2b2JPBubT4gpjYLN1BWSZtzzu0CrqjvYgQ4fJ9sBvRLxz+H65eHgYfruxgBDt8n+UB+XXxAJG4pFAI5Fd63Atb4VIuUUZ+EJvVL6Al6n0RiKEwHOppZOzOLA0YAE32uKdKpT0KT+iX0BL1PwjoUzOwl4Augs5kVmtko51wJcB3wAbAQeMU5942fdUYS9UloUr+EHr/6RAPiiYhIubDeUhARkepRKIiISDmFgoiIlFMoiIhIOYWCiIiUUyiIiEg5hYIIYGZF3nNbM/tpHa/7t4e8/7wu1y9SlxQKIt/XFqhWKHjDGVfle6HgjWQpEpIUCiLfdzdwgpnNMbNfmlm0md1rZtPNbJ6Z/QLAG5Vyspn9G/jaa3vTzGaa2TdmNtpruxtI9Nb3otd2cKvEvHXPN7OvzWx4hXXne3fRWmRmL5pZZQOhidS5SBwlVaQqtwG/cs6dCeD9ct/unOtvZvHAVDP70Jt3ANDdObfCe3+lc26LmSUC083sdefcbWZ2nXOudyWf9ROgN9ALaOotM8Wb1gfoRtlgZ1OBQZTdQEUkqLSlIFK1U4CRZjYHmAZkAB29aV9VCASAG8xsLvAlZSNZdqRqxwMvOedKnXPrgU+A/hXWXeicOwDMoWy3lkjQaUtBpGoGXO+c++B7jWZ5wK5D3p8EHOucKzazfCAhgHUfzt4Kr0vRd1XqibYURL5vJ5Ba4f0HwNVmFgtgZp3MLLmS5RoDW71A6AIMrDBt/8HlDzEFGO4dt8gEBgNf1clPIVJD+utD5PvmASXebqDngIco23UzyzvYuxE4p5Ll3geuMrN5wGLKdiEdNBaYZ2aznHMV75/7BnAsMBdwwK3OuXVeqIj4QkNni4hIOe0+EhGRcgoFEREpp1AQEZFyCgURESmnUBARkXIKBRERKadQEBGRcgoFEREp9/9EBjpml6tiHgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Initial condition\n",
        "winit = np.random.randn(p)\n",
        "\n",
        "# Parameters\n",
        "feval = log_fun.feval\n",
        "nit = 1000\n",
        "lr = 1e-4\n",
        "\n",
        "# Run the gradient descent\n",
        "w, f0, hist = grad_opt_simp(feval, winit, lr=lr, nit=nit)\n",
        "\n",
        "# Plot the training loss\n",
        "t = np.arange(nit)+1\n",
        "plt.semilogx(t, hist['f'])\n",
        "plt.grid()\n",
        "plt.xlabel('Iteration')\n",
        "plt.ylabel('Training loss')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z0t2UoDhzVqd"
      },
      "source": [
        "We can measure the accuracy of the final estimate by creating a predict method."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DeUgwJegzVqd",
        "outputId": "aa6ce0ab-6e7c-4be8-dca1-8f5ffd8792b3"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Test accuracy = 0.965854\n"
          ]
        }
      ],
      "source": [
        "def predict(X,w):\n",
        "    z = X.dot(w[1:]) + w[0]\n",
        "    yhat = (z > 0)\n",
        "    return yhat\n",
        "\n",
        "yhat = predict(Xts1,w)\n",
        "acc = np.mean(yhat == yts)\n",
        "print(\"Test accuracy = %f\" % acc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EAlYgYBmzVqe"
      },
      "source": [
        "You will likely not get a test accuracy quite as good as the test accuracy for the `sklearn` method, although there is some variation with the train-test split.  The reason is that the learning rate was somewhat slow and we didn't yet fully converge.\n",
        "\n",
        "To see the effect of the learning rate, the code below tries different learning rates."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dsro8B1BzVqe",
        "outputId": "d4f3b814-92bd-4f43-ba2e-8e6400d67539"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "lr=    1.00e-05  Test accuracy = 0.951220\n",
            "lr=    1.00e-04  Test accuracy = 0.965854\n",
            "lr=    1.00e-03  Test accuracy = 0.960976\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "Text(0, 0.5, 'Training loss')"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3gU1f7H8ffZkl5ISAOSECD00EPokEJVpKhYEEHFiwqoFxu2i2K5gD/LBcVCUakCogjSREqo0msglFCTAKEmJJCe8/tjFkQIISTZTJI9r+eZZ3dnZ2Y/cWS/OzNnzhFSShRFURQFwKB3AEVRFKXsUEVBURRFuUEVBUVRFOUGVRQURVGUG1RRUBRFUW5QRUFRFEW5waR3gOLw8vKSQUFBRVr36tWrODs7l2wgpdjUfil71D4pm4qzX3bs2HFBSumd33vluigEBQWxffv2Iq0bHR1NeHh4yQZSik3tl7JH7ZOyqTj7RQhx8k7vqdNHiqIoyg2qKCiKoig3qKKgKIqi3FCurykoilL6srOzSUhIICMjQ+8oNs3d3Z3Y2NgCl3FwcMDf3x+z2Vzo7aqioCjKPUlISMDV1ZWgoCCEEHrHsVmpqam4urre8X0pJRcvXiQhIYEaNWoUervq9JGiKPckIyODypUrq4JQxgkhqFy58j0f0dnkkcKeg3+x/tg8zFXr42Bvh6PZiKOd8R+PZqNQ/9Mryh2ofxvlQ1H2k00WhV//+oRfjXH0ntOGNVmd+Dm3Ewm33MdhNAgczUYczEYc7QxawbileDhcn3d9vt0/Xzvks46j2YiLgwlXe5P6h6UoReDi4kJaWlqpfV737t3ZvHkz7du3Z/Hixfkuk5mZycCBA9mxYweVK1dm7ty5XL+xdsyYMUydOhWj0ciECRPo1q3bPX3+ndYPCQnBzc0No9GIyWQq8j1bt7LJovBA59H8uvoJdlQN4KWEBbxkWsA57zYcqdaXI56duJZr5FpWDulZeaRn55KRnUt6Vi7p2dp0IS1Le55lec8y3ct4RQYBbo5mKjmacXc042Z5dHc0U8np7+faZKc9WuY72xlVQVGUm+Tm5mI0Gq2y7ddff51r167x3Xff3XGZqVOn4uHhQVxcHHPmzGHkyJHMnTuXAwcOMGfOHPbv38/p06fp3Lkzhw8fLnTWu62/Zs0avLy8SuTvvM4mi0IL/0Z4m7z5q0owTz88E3bPwnfXTHx3v057R09o8hg0exJ8GxR6m1JKMnPy/lE8/lE0sv5+TM3IISU9m5T0bJItjynp2SRcTr/xPDfvzhXGZBC4O5qp7GKHt6s93i722uP1ycUBHzdtvrujGYNBFRCl4omOjmb06NFUqVKF3bt3c+DAAat8TlRUFNHR0QUus3DhQt5//30AHn74YYYPH46UkoULF/LYY49hb29PjRo1CA4OZuvWrbRp04aZM2cyYcIEsrKyaNWqFV9//fVtxaKg9a3FJouCEIImTk2IPhvNFSd33MLfhI6vw7Fo2Dkdtk6GzV9DtVBo/iSEPAT2d77Kf32bDpbTTR7FzCelJC3z78KRkp5NyrXsf7xOTs/mYlom51Mz2XHqMueuZJKZk3fbtsxGgdf1ouFiT9VKjlTzcKTaTY/eLvaqcChFMvr3/Rw4faVEt9mgqhvvPdCwUMtu3bqVmJiYG61rOnToQGpq6m3Lffrpp3Tu3LlEc94sMTGRgIAAAEwmE+7u7ly8eJHExERat259Yzl/f38SExOJjY1l7ty5bNy4EbPZzNChQ5k1axYDBw68bbv5rQ/ad07Xrl0RQvDcc88xZMiQEvlbbLIoADR2aszKKytZl7COnjV7gsEIwVHadPUC7J2rFYjfX4blb0NIX2g2EALCwMqnboQQuDqYcXUw41/ICnO9kJxP1QrFeUvBuPl1YnI6205c4kpGzj/WtTMaqFLJgWqVHKlayZEADyeCvJyo4eVMkJczbg6Fb+OsKKUpLCzsH80t169fr0uO/Ma6F0Lccf6qVavYsWMHLVu2BCA9PR0fH59CbxdgxYoV1KlTh3PnztGlSxfq1atHx44di/un2G5RqG5XHS9HL1afWq0VhZs5e0GbYdB6KCRsh53TIOZX2DUTvOpqRw9NHteWKyNuLiQ1vV0KXDY1I5vTyRkkJl8j8XI6CcnpJF5O53RyOuuPnCfpSuY/lvdysSOosvONIlHL25navq5U93TCZFStmm1ZYX/RW8utvYTqdaTg7+9PfHw8/v7+5OTkkJKSgqen54351yUkJFC1alVOnz7NoEGDGDNmzD+2s2DBAkaPHg3AlClT7rg+QJUqVQDw8fGhb9++bN26VRWF4jAIAxEBESw+tpiMnAwcTA63LyQEBLTUpu5jYP8C7ehhxbuwcjTUu087eqgVoR1plBOuDmbq+pmp65f/KbGM7FxOXrzG8QtXOX7hKicuXOX4xausPXyen3ck3FjOzmQg2NuFun6u1PZ1oa6vK/WquFHV3UFdCFd0odeRQq9evZg2bRpt2rRh/vz5REZGIoSgV69e9O/fn1deeYXTp09z5MgRwsLCcHd3p3fv3owYMQIfHx8uXbpEamoqffv2pW/fvje26+jomO/6V69evXHz2tWrV1mxYgWjRo0qkb/FZosCQFRgFD8f/pktZ7bQKaBTwQvbu0Lzgdp0LhZ2zoA9P8GBheDmD82egKZPgEf10glvRQ5mI3X9XPMtGmmZORw7n8bhpDQOJ6Vy6GwqW45dZMGuxBvLeDrb0bCqGyHV3Amp6k5INTcCPZ1UoVDKpQ4dOnDw4EHS0tLw9/dn6tSpdOvWjVGjRhEaGkqvXr0YPHgwTz75JMHBwXh6ejJnzhwAGjZsyCOPPEKDBg0wmUxMnDgRo9FIgwYN+Oijj+jatSt5eXmYzWYmTpxI9er//P640/pJSUn07t0bg8FATk4O/fv3p3v37iXy94r8zlmVF6GhobI44ym069COjnM70qV6Fz5o98G9byQnEw4t1Y4ejq7R5tUM1wpHvfvBZF+kbOVRSno2R5JSiT1zhZjEK8ScTuFwUirZudr/X57OdjQP9CA0yIPQ6h408nfH3nT70ZXqu7/suXWfxMbGUr9+ff0CKcDdu7m4Lr/9JYTYIaUMzW95mz5SMBvNdKjWgej4aHLzcjHe6ykgkz007KtNyadg1yzYPQvmPw1FbNpaXrk7mgkN8iQ0yPPGvMycXA6fTWNfYgq7Tl1mx8nLrIxNArRTT42ruWvrVPegRXUPPJzt9IqvKIqFTRcFgMjqkSw7sYxd53YR6pdv4SycSoEQ8RZ0egOOrdFOL/2jaetACHnwrk1bKxJ7k5FG/u408nenf6tAAC6kZbLjpFYgtp24xNQNx/h2rXY0UcvbmSDHTHJ8kmhTqzLO9jb/v6eilDqb/1fXvmp7zAYzq+NXF68oXGcwQnBnbbp6AfbMgV0z4PeXYPlbWtPWNsPBxzYPv71c7OnW0I9uDf0A7aL23oQUtp24xLYTl1gfd5VV07djNgpaBnnSsY43nep4U8/PVV2TUJRSYPNFwcXOhdZVWrP61GpeD329ZL94nL2g7XCteWvCNu3aQ8yvWqFoPRTC3wQ72x4Q3cFsJKyGJ2E1tNNOf65eg3NgI9YePs/aw+cZu+wgY5cdxMfVno51vOlc34dOdXxwtCs/rb0UpTyx+aIAEBkYyfq/1nP48mHqetYt+Q8QQrvpLSAMOo+GlaNg0wStQNz3iXZRWgHAbBC0DfaibbAXb91Xn6QrGayzFIiVsUnM35GAo9lIeF1vuof4EVnPB1d1c52ilBir3nkkhKgkhJgvhDgohIgVQrQRQngKIf4UQhyxPHpYlhVCiAlCiDghxF4hRHNrZrtZeEA4AsHqU6ut/2HOlaH3RHh6OTi4wZz+MPsxuHzS+p9dDvm6OdAvNICv+jdn+zudmf1sKx5u4c/2k5d5ec5uWny0ksE/bmP+jgRSM7L1jqso5Z61b0cdDyyXUtYDmgCxwJvAKillbWCV5TVAD6C2ZRoCfGPlbDd4OXrR1Kcpq06tKq2PhOpt4Ll10OVDOL4OJraC9Z9BTlbpZShnTEYDbYO9+LBPCFveiuLn59swoFV1Ys9c4bWf9xD60Upe+mkXaw+fL7BDQaV8c3Ep+I79kta9e3cqVapEz54977hMZmYmjz76KMHBwbRq1YoTJ07ceG/MmDEEBwdTt25d/vjjj3v+/ILWz83NpVmzZgVmu1dWKwpCCDegIzAVQEqZJaVMBnoD0yyLTQP6WJ73BqZLzWagkhCiirXy3SoqMIpDlw+RkJpw94VLitEM7V6C4Vu1PpdWfQDftocTG0ovQzllMGgXokc90ICNb0by69C29Av1Z+3h8wz6fittxqxizNJYDifd3uWBUvHk5uZabduvv/46M2bMKHCZm7vOHjFiBCNHjgT+2fX18uXLGTp06D1lvdv648ePL/F7Rqx5pFATOA/8IITYJYSYIoRwBnyllGcALI/Xe4GqBsTftH6CZV6piAyIBCidU0i3cveHx2ZB/3mQkw4/3g8Lnoe086WfpRwSQtA80IOP+jRi6ztRfPNEcxr7uzN1w3G6frGOh7/ZxG+7EsnMsd4Xh1L6oqOjiYiIoH///jRq1MhqnxMVFXXXm8QWLlzIoEGDAK3r7FWrVhXYdTbAzJkzCQsLo2nTpjz33HP5FouC1k9MTGTJkiU8++yzJfr3WvNCswloDrwopdwihBjP36eK8pNfs5/bzgEIIYagnV7C19f3rv2c30laWtpt61Y1V2XBvgUEng8s0jaLzx5Do0+pfnIeAXt/Jnf/Io7VHMiZKl1B2EbHc/ntl3vlCAyoDr38HNl4Oofo+GT+Pfcy/1kAHauZCQ8w4e1kG/89S8Kt+8Td3f1Gp3P2a97DcG5/iX5enk9DMiNGF7hMamoq165dY+vWrWzevJmgoCBSU1Pp1q1bvqOyffTRR0RERBQ507Vr18jJycm3sz2A+Ph4PDw8brzv6urKyZMnOX78OC1btrwx39fXl7i4OMxmM7NmzWL58uWYzWZGjBjBlClT6N+//z+2e6f1Q0JCGDlyJO+99x5paWkFZsvIyLinf1PWLAoJQIKUcovl9Xy0opAkhKgipTxjOT107qblA25a3x84fetGpZSTgEmgdXNR1C4R8utOIWZXDJP3TaZx68Z4Onjmv2Kp6Abn38Cw5FXqHv6Gule3Qc/PoUoTHTOVjpLu5qIXkJcn2Xj0AjP+Osmy2CSWnsimawNfhnSsSYvqeu7n8iG/bi5u/HI224GxhL9GzHbY3eWXuaurK05OToSFhf3jKGHTpk0lm8XCyckJk8l0xyMGIQQuLi433jcYDLi6umI2m3F0dLwx32w24+TkxObNm9mzZw+RkdoZivT0dPz9/W/b/p3WX7t2Ld7e3nTs2JHo6OgCszk4ONCsWbNC/61WKwpSyrNCiHghRF0p5SEgCjhgmQYBYy2PCy2rLAKGCyHmAK2AlOunmUpLVGAU3+39jrXxa+lbu+/dV7Am77ow6HfYOw9WvAOTwiHsOYh4W2u1pBSawSDoUNubDrW9OZ2czqwtJ5m15RR/7E+ieWAlhnSsSZcGfhjVQEP3rsdYXT/eVrvOXrRoEcuWLSMoKIiMjAyuXLnCgAEDmDlzZvH/GCml1SagKbAd2Av8BngAldFaHR2xPHpalhXAROAosA8Ivdv2W7RoIYtqzZo1t83Ly8uTXX/uKoetHFbk7VrFtUtS/j5Cyvfcpfy/OlLu+0XKvDy9U1lFfvvFGq5mZssfNx6X7cetktVHLpadPlkt5207JbNzckvl88uTW/fJgQMH9Ali4ezsLKXUct1///2l8pl3+6yvvvpKPvfcc1JKKX/66SfZr18/KaWUMTExsnHjxjIjI0MeO3ZM1qhRQ+bk5Mj9+/fL4OBgmZSUJKWU8uLFi/LEiRO3bfdO60sp5ZUrVwqVLb/9BWyXd/heterNa1LK3UB+fUdE5bOsBIZZM8/dCCGIDIxk3qF5XMu+hpPZSc84f3P00E4fNX0ClozQOtzbNQMeGK/1uaTcMyc7E4PaBjGgdXWWx5zl6+g4Xp+/l4lr4ngxsja9m1ZVAwgpQNnsOtuabLrr7PzOXW87u41n/niGzzp9RtegrsVMaAV5ubBtqtZ81WiGh7/XBvmpIPTqOltKyYoDSfxv5RFiz1yhhpczL0UF06tJNZs/raS6zi6brNV1tvopdItmPs2oZF+J1fE6NE0tDIMRWg2B59aCiw/MfBA2fQnluLiXBUIIujX0Y8mL7fl2QAvsTQZGzN1Dly/WsjzmbL5j5SpKRaSKwi1MBhOd/DuxLn4d2XlluNuEyrXg2ZVQr6c2POgvgyHrqt6pyj2DQdA9xI+lL3Xg2wHNMQjB8zN38OikzexNSNY7nqJYnSoK+YgMjCQ1O5VtZ7fpHaVg9q7wyHSIek/rXG9qV7h0XO9UFYJWHKqw/OUOfNQnhKPn0uj11UZGzN3N6eR0veMpitWoopCPtlXb4mhy1Ofu5nslBHR4BZ6YDynxWtPVuFLsw6mCMxkNDGhdnejXw3khvBZL9p0h4tNovvjzMBnZ6g5ppeJRRSEfDiYH2lZty5pTa8iTeXrHKZzanWFINLhVg1kPw4Yv1HWGEuTqYGZk93qsfrUTXRv6MX7VEbr9bx1rD6uuSJSKRRWFO4gKjOJc+jn2XyjZW/ityrMmPPsnNOgDK9+HnwdB5u23/CtF5+/hxJePN2PWs60wCsGg77cybNZOzqZk6B1NUUqEKgp30NG/I0ZhLN3utEuCnbPWTLXLhxD7O0zpDBeP6p2qwmkX7MWyf3fg1S51WBmbRNRn0cz46wR5qstuq1NdZ2v9GYWHh9OkSRMaNmzIe++9d8/bvRNVFO7A3d6dUL/Qsts0tSBCaF1yD/gV0s7C5Ag4vELvVBWOvcnIi1G1+XNEJ5pX9+A/C/fTf8pmTl28pnc0m2NrXWfb29uzePFi9uzZw+7du1m+fDmbN28u1t95nSoKBYgMiOR4ynGOpRzTO0rR1IqAIWu1u55nPwLr/k9dZ7CCwMpOTH8mjHEPNWJ/4hW6/W8dP248ro4arMyWu86+3gEfQHZ2NtnZ2SU2vrwao7kAkYGRjNk6htWnVlOzUU294xSNR3V4ZgX8/jKs/giunIb7PgOD+j1QkoQQPNoykA61vXnr1328//sBlu8/yxePNqWKu6Pe8axm3NZxHLx0sES3Wc+zHiPDRhZq2a1btxITE0ONGjUA/TrES0xMJCBA6+TZZDLh7u7OxYsXSUxMpHXr1jeW8/f3JzExkdjYWObOncvGjRsxm80MHTqUWbNmMXDgwNu2m9/6oB0dNW3alLi4OIYNG0arVq1K5G9RRaEAfs5+NKzckNWnVvNso5IdyKJU2TnBg5PArSps/J92k1vvr0u+y2OFqpUc+fHplszbHs/o3w/Q/X/rGfdQI7qHlNoggjYlLCzsRkEAWL9+vS458rvjXQhxx/mrVq1ix44dtGzZEtC6zvbx8blt2TutD2A0Gtm9ezfJycn07duXmJgYQkJCivunqKJwN1GBUUzYNYGkq0n4OvvqHafohIAuo7Ub3lZ/qBWGh78Hk73eySqc60cNYTUq8/KcXTw/cyePhwXyn571cbKrWP/kCvuL3lpstevsm1WqVInw8HCWL19eIkVBnUO4i8hAbRCMNfFrdE5SQjq+Bt3HwcHF8NNjqmsMK6rh5cz859vyXKea/LT1FA98uYH9p1P0jlWhrV+/nt27d982WbMgAPTq1Ytp07Sh5+fPn09kZCRCCHr16sWcOXPIzMzk+PHjHDlyhLCwMKKiopg/fz7nzmljjF26dImTJ0/St2/fG5mv98Ca3/rnz58nOVnrdiU9PZ2VK1dSr169EvlbVFG4i5ruNQlyCyofdzcXVuvnofdEOBYNMx+CDPVFZS12JgNv9ajPzMGtSM3Ioe/Xm5i77ZTesZR70KFDB/r168eqVavw9/e/0Sx01KhRLFq0CIDBgwdz8eJFgoOD+fzzzxk7Vht86Oaur7t3755v19mNGzemS5cunDlz+5hid1r/zJkz9OzZk8aNG9OyZUu6dOlSYJPZe3KngRbKw1TSg+zcyWfbP5NNpzWVyRnJRf68MinmVylHe0r5bUcp0y7onUZKWXqD7OjhQmqGfGLyZll95GI5cv4emZ6Vo3ekQilrg+womuuD7NzNvQ6yo44UCiEyIJIcmcP6RH0uYllNw77w2E9w/iD8eB9cKdXRT21OZRd7pj0TxvCIYOZsi6fft38Rf0nd06CULaooFEJj78Z4OXpVrFNI19XpaulMLwF+6A6XT+qdqEIzGgSvdavLlIGhnLh4lQe+2sD6I6r/JKXsUEWhEAzCQERABBsSN5CRUwH7uKnRAQYuhPRk+L47nD+sd6IKr3MDX34f3h4/Nwee+mEb0/86oXckRQFUUSi0qMAo0nPS2XJmi95RrMM/FJ5aAnnZ2qmkC3F6J6rwgrycmf9CWyLqejNq4X7eWxhDTm756JVXqjvjy4Wi7CdVFAopzC8MF7NL+esg7174hcDTy7SuMKb31k4pKVblYm/iuydDGdKxJtP+Oskz07ZzJaMMj/gHODg4cPHiRVUYyjgpJRcvXsTBweGe1qtYd9JYkdlopoN/B6Ljo8nNy8VoMOodyTq8asOTv8KPPWF6H61IuHjrnapCMxoEb99Xn1rezryzIIYHv97E1EGhVK/sfPeVdeDv709CQgLnz6trIXrKyMi46xe+g4MD/v7+97RdVRTuQWRgJMuOL2PXuV2E+oXqHcd6qjSB/vNgRl+Y2RcGLQbHSnqnqvAebRlIoKczL8zawYNfb+KHp1vS2L/s/Xc3m83/6FpC0Ud0dDTNmjUr8e1a9fSREOKEEGKfEGK3EGK7ZZ6nEOJPIcQRy6OHZb4QQkwQQsQJIfYKIZpbM1tRdKjWAbPBXD67075X1dvAozPh3EHLnc+q6WRpaFOrMr++0BYHs5HHJm1WLZOUUlca1xQipJRNpZTXf1q/CaySUtYGVlleA/QAalumIcA3pZDtnjibnWldpTWrT622jfOptTtrHemd2gzznoScLL0T2YSa3i78OrQtgZ5OPPPjNhbuTtQ7kmJD9LjQ3BuYZnk+Dehz0/zplhvuNgOVhBBlrmvJqMAoEtMSOXzZRppthjwID4yHuJXw678gTw1WXxp83RyY93wbmgd68PKc3UxZX07H9FDKHWtfU5DACiGEBL6TUk4CfKWUZwCklGeEENf7i60GxN+0boJl3j9usxVCDEE7ksDX15fo6OgiBUtLSyvSuuZcMwLB1HVTua/SfUX67PKnOgE1n6LWgR85MymNQ3WHa72uWkFR90tFNThYknvNyEdLYtkTG8eDtc0lNphKYal9UjZZa79Yuyi0k1Ketnzx/ymEKGg0jvz+T7/tHI2lsEwCCA0NleHh4UUKFh0dTVHX/WXZLxzLPlbk9cuncFhVmSrrP6NKzfrQ9SOrfEpx9ktFFRUheWfBPuZsi8e3qj/v3F+/VAuD2idlk7X2i1VPH0kpT1sezwELgDAg6fppIcvjOcviCUDATav7A6etma+oIgMjOXT5EAmpNtaOP/I/0PJZ2PQlbPlO7zQ2w2gQ/LdvI55qG8SUDcf5z8IYNdSnYjVWKwpCCGchhOv150BXIAZYBAyyLDYIWGh5vggYaGmF1BpIuX6aqayJDNDGWKiQfSEVRAjo8QnU6QHL34TDf+idyGYYDIL3HmjAc51qMnPzKUb+spdcVRgUK7DmkYIvsEEIsQfYCiyRUi4HxgJdhBBHgC6W1wBLgWNAHDAZGGrFbMUS4BZAbY/attE09VYGIzw0Bfwawc9Pw5m9eieyGUII3uxej5ejavPzjgRembe73HSLoZQfVrumIKU8BjTJZ/5FICqf+RIYZq08JS0yIJLJ+yZzKeMSng6eescpXfYu8PhcmBIFsx+Ff63Sxn9WrE4IwYgudXAwGxm3/CC5eZLxjzXDaCjdi89KxaX6PiqiqMAo8mQea+PX6h1FH25VoP9cyLyiFYbMNL0T2ZQXwmvxVo96LN57htfn71HXGJQSo4pCEdXzrEdV56oVu4O8u/FrBP1+hKQY+OVZdQ9DKXuuUy1e7VKHX3cm8vaCfaowKCVCFYUiEkIQGRjJX6f/4lq2DXcBUbuLdvH58DL44x2909icF6Nq3xjJ7f3f99vGnfaKVamiUAyRgZFk5WWxIXGD3lH0FfYvaD0MtnwDWybpncbmvNq1DkM61mT6Xyf5eEmsKgxKsaiiUAzNfJpRyb6SbbZCulXXD6HufVpT1ePr9E5jU4QQvNWj3o37GD7/00a6YFGsQhWFYjAZTHTy78S6+HVk55btgVGszmDUOs/zqg3zBqmxnkuZENp9DI+1DODL1XF8v+G43pGUckoVhWKKCowiNTuVbUnb9I6iP3tXeGy2dsF57hOqu+1SJoTgoz4hdGvoyweLD/DbLtW7qnLvVFEopjZV2+BocrS9u5vvpHIteHgqnI2BRS9qQ3sqpcZkNDD+sWa0runJaz/vIfrQubuvpCg3UUWhmBxMDrSr2o41p9aQJ9XdpYDWIinqPxAzX+snSSlVDmYjkweGUtfPlRdm7mTnqct6R1LKEVUUSkBkYCTn0s8RcyFG7yhlR/tXoEFvWPkexNnwvRw6cXUw8+PTYfi42fPMj9s4kpSqdySlnFBFoQR09O+IURjVKaSbCQG9vwbv+jD/GbikBokpbd6u9sx4phVmo4GB328l6UqG3pGUckAVhRLgbu9OqF+obd/dnB97F3hslvZ87pPqwrMOAis78ePTLbmSns3TP2wjLTNH70hKGaeKQgmJCozixJUTHEtRv4j/wbOG1qtq0n5Y+pq68KyDhlXd+eqJ5hxKSuXF2TtVz6pKgVRRKCERARGADY6xUBi1u0DH12H3LNg1Q+80Nimirg8f9G7ImkPnVXcYSoHuqShYBsBxtlaY8szP2Y+GlRuqonAn4W9CzXBY8hqc2aN3Gpv0RKvqNwbpmbJe3dym5O+uRUEIMV0I4SaEcAL2A8eFEK9YP1r5ExUYxb4L+0i6mqR3lLLHYISHpoJTZZg3ENKT9U5kk0Z2q8f9jarw8dJYlh+VyRYAACAASURBVO4rkwMbKjorzJFCIynlFaAPsAJt7OSnrBmqvIoM1IbpXBO/RuckZZSzl9bVdkoC/PaCur6gA4NB8NkjTWhR3YMRc3ezO14VZ+WfClMU7IQQJqA38JuUMgtQV6ryUdO9JkFuQeoUUkECW0GXD+HQUtg4Xu80NsnBbGTSky3wdrVnyPTtnE1RTVWVvxWmKEwBTgEewFohRCCghtnKhxCCiMAItp3dRkpmit5xyq7WL2g3tq36AOK36p3GJlV2sWfKoFCuZuYwZMZ2MrLVAEmK5q5FQUr5hZSyqpSyq2Uc5Xgg0vrRyqeowChyZA7rE9frHaXsEgJ6fQnu1WD+YEhX3TDooZ6fG1882pR9iSm8MX+vapGkAIW70DxcCOFmef4dsAXoYO1g5VUjr0Z4O3qrU0h34+AOD30Pqadh0Uvq+oJOujb047WudVm05zRfRx/VO45SBhTm9NEQKeUVIURXoBrwAvCJdWOVXwZhICIggg2JG8jIUedqCxTQEiL/A7GLYMcPeqexWUPDa9G7aVU+XXGIFfvP6h1H0VlhisL1n3A9gB+klDsKuR4AQgijEGKXEGKx5XUNIcQWIcQRIcRcIYSdZb695XWc5f2ge/tTyo7IwEjSc9LZfGaz3lHKvrYvQa1IWP6WdtezUuqEEIx7qDGNq7nz77m7OXj2it6RFB0V5st9jxBiKfAAsEwI4cLfhaIwXgZib3o9DvhCSlkbuAwMtswfDFyWUgYDX1iWK5fC/MJwMbuoU0iFYTBA3+/A3g1+fhpDbqbeiWySg9nId0+G4mJv4tlp27l0NUvvSIpOClMUngbeB8KklNcAB/7+Ii+QEMIfuB+tBRNCCIF2kXq+ZZFpaPc/gNbkdZrl+XwgyrJ8uWM2mung34Ho+Ghy8lQHZHfl4qMN5XnhMMFxU/ROY7P83B2YPDCUc6mZvPiT6iPJVhWm9VEu4AW8IYQYC7SUUu4q5Pb/B7zB3/c1VAaSpZTXvykT0K5TYHmMt3xmDpBiWb5cigqM4nLmZXaf2613lPKhVgS0H0HVMysg5he909isJgGV+KhPCBvjLvLpisN6x1F0YLrbAkKIj4F2wGzLrNeFEO2klO/eZb2ewDkp5Q4hRPj12fksKgvx3s3bHQIMAfD19SU6Ovpuf0K+0tLSirxuYcg8iQkT0zZNI81T3dZRGMLQjkbOv+O2YDjb47PJcPTTO5JN8gEiAkx8u/YoxpQE6rtkWPXfilI01voOE3drmyyE2As0v/7r3nJ3804pZeO7rDcGeBLIQTvl5AYsALoBflLKHCFEG+B9KWU3IcQflud/WT7jLOAtCwgYGhoqt2/fXti/9R+io6MJDw8v0rqFNWzVMI4mH2XZg8sop2fCSt3mZXNpved1qBwMz/wBRrPekWxSZk4uj363mSNJqbwTZkf/nurWpLKmON9hQogdUsrQ/N4rbCsi1zs8vyMp5VtSSn8pZRDwGLBaSvkEsAZ42LLYIGCh5fkiy2ss768uqCCUB5EBkSSmJXL4sjoML6wMR1/txrbEHbD6Q73j2Cx7k5FvBjTH0c7IhF0ZpGZk6x1JKSWFKQqfADuFEFOEEFOB7RSvZdBI4BUhRBzaNYOplvlTgcqW+a8AbxbjM8qETgGdEAg1Itu9atAbQp/R+kY6slLvNDarirsjX/Vvzrlrklfn7SEvr1z/RlMKqTAXmmcC7YGllqmjlHLWvXyIlDJaStnT8vyYlDJMShkspewnpcy0zM+wvA62vF/uhzDzcvSimU8z1TS1KLr9F3waaL2ppp3XO43Nal2zMo/WtWPFgSS+WavueLYFdywKQojG1ye0X/RxwBG0X/MFXk9Q/hYZGMmhy4dISE3QO0r5YnbUhvHMSIFFL6puMHTUtbqJXk20O57XHVYFuqIr6EhhYgHTV9aPVjFEBmgX6NTRQhH4NoTO78PhZaobDB0JIRj7UCPq+rry0pxdJFy+pnckxYruWBSklB0KmDqWZsjyLMAtgNoetdV1haJq9TzUjIDlb8OFI3qnsVlOdia+GdCCnFzJ8Nm7yMpRN7ZVVPc0RrNSNFGBUew+v5uL6Rf1jlL+GAzQ5xvtdNIvz0KuagWjlxpeznzycGN2xyczZlns3VdQyiVVFEpBZEAkeTKPtQlr9Y5SPrlVgV4T4MxuiB6jdxqbdl+jKjzVNogfNp5QYzxXUKoolIJ6nvWo6lxVXVcojvoPQLMnYf3ncHKT3mls2tv31adpQCXemL+XExeu6h1HKWGFGWSncT5TdSGEKiiFJIQgMjCSv07/xdVs9Y+oyLqPBc8a8OtzWqskRRd2JgMTn2iOySh4YdZONZRnBVOYL/apwA5gOjAD7ea1BcARIUSUFbNVKJGBkWTlZbExcaPeUcovexd4cDJcSYQlr+mdxqZVq+TI5480IfbMFd5fpMbBqEgKUxSOAC2klE2llE2AFsButD6MPrNmuIqkmU8zKtlXUq2Qiss/FDqNhH3zYN/8uy+vWE1kPV+GhtdizrZ4ftmh7sOpKApTFOpLKfdefyGl3IfWQV6c9WJVPCaDifCAcNYnrCdbtaApng6vQkArWPwKJJ/SO41Ne6VLHVrV8OTd32I4nJSqdxylBBSmKBwVQnwphGhnmSYAcUIIe7QeUJVCigyIJDU7lW1J2/SOUr4ZTdpobTIPFjwPeeqctl5MRgNfPt4MZ3sTL8zcwdVM9ZVQ3hWmKAxEGwznTeAt4DRab6Y5gLqmcA/aVG2Do8mRpceW6h2l/POsAfd9Aic3wqYJeqexaT5uDkx4vCnHL1zl7QX7KOedG9u8wnSId01KOU5K+YCUsqeUcqyU8qqUMldKqZqA3AMHkwO9a/VmyfElJF1N0jtO+dfkcWjQB1Z/DKfVCHd6alvLi1e61GHh7tPM2qJO6ZVnhWmS2loIsUwIcUAIcfj6VBrhKqKnQp5CSsmMAzP0jlL+CQE9vwBnb+1u5yzVJ4+ehoYH06mONx/8foCYRPV7sbwqzOmjH4Cvgc5Ah5smpQiquVSjR40ezDs8j5RM9Q+n2Jw8oe+3cPEIrChwhFjFygwGwRePNsXT2Y6hs3ZyRQ3MUy4VpihckVL+LqU8LaVMuj5ZPVkF9nTI06TnpPPTwZ/0jlIx1OwEbYbD9qlwaLneaWyap7MdE59oxunkdEbO36uuL5RDhSkKq4UQY4QQLW8ZY0Epojoedejk34lZsbO4lq1OeZSIqFHg2wgWDoO0c3qnsWktqnvyRve6LIs5y4+bTugdR7lHhSkK7S3T56jxFErM4EaDSc5MZkHcAr2jVAwme3hoMmSlaYVB/ULV1b861KRzfR/+uzSW3fHJesdR7kFhWh+p8RSsoJlPM5r7NGfa/mlk56lzryXCpz50+QCOrIBtU/ROY9OEEHzWryk+rg4Mm7WT5GtZekdSCqmg4Tgftzy+lN9UehErrsGNBnPm6hmWH1fnwUtM2BAI7qxddD6n+vzXk7uTmYlPNOdcagav/bxHXV8oJwo6UvCwPHrfYVKKqUO1DgRXCub7mO/Jk2okqxIhhDYoj52L1kw1J1PvRDataUAl3r6vPitjzzF5/TG94yiFUNBwnF9bHv+T31R6ESsuIQSDGw0mLjmOdQnr9I5Tcbj4QJ+vISkGVo7WO43Ne6ptED1C/Bi3/BDbT1zSO45yF4W5ec1LCPGGEOJrIcSk61NphLMF3YO6U82lGlP2TVGH1yWpTjdo+S/YPBHiVM+0ehJCMO7hxvh7ODJ89i4uXVXXF8qywrQ+Wgj4AhuAVTdNBRJCOAghtgoh9ggh9gshRlvm1xBCbBFCHBFCzBVC2Fnm21tex1neDyrqH1WemAwmBjUcxJ7ze9h5bqfecSqWrh+Cdz347QW4ekHvNDbNzcHMxP7NuXQtixFzd5OXp34AlVWFKQrOUspXpZSzpZRzr0+FWC8TiLSMwdAU6C6EaA2MA76QUtYGLgODLcsPBi5LKYOBLyzL2YQ+wX3wdPBk6r6pekepWMyO8NBUSL8Mi15UzVR1FlLNnVE9G7D28Hm+WXtU7zjKHRSmKCwTQnS91w1LTZrlpdkySSASuD46yjSgj+V5b8trLO9HCSHEvX5ueeRocuSJ+k+wPnE9hy4d0jtOxeIXAp1Hw6GlsP17vdPYvCdaBdKrSVU+W3GIzccu6h1HyYe423lsIcRlwB24BmQBAu073/OuGxfCiDaUZzDaTW//B2y2HA0ghAgAlkkpQ4QQMUB3KWWC5b2jQCsp5YVbtjkEGALg6+vbYs6cOffw5/4tLS0NFxeXIq1rDdfyrvFewns0dGzIU95P6R1HN1bZLzKPxntH455ygB0tPueac0DJbr+CK+l9kp4jGb0pnfRc+KCtI+72NvHbr8QVZ79ERETskFKG5vumlLLACTDmN91tvVu2UQlYg9aRXtxN8wOAfZbn+wH/m947ClQuaLstWrSQRbVmzZoir2stn237TDae1lieunJK7yi6sdp+uXJGynE1pPymnZTZGdb5jArKGvsk9kyKrPPOUtl/8l8yJzevxLdvC4qzX4Dt8g7fqwXdvFbb8rThHaZCk1ImA9FAa6CSEMJkecsfbdAe0AbyCbB8tgnt6MSm2q8NaDAAozAybf+0uy+s3BtXP+j1FZzdB6s/1DuNzavn58aHvUPYGHeRCauO6B1HuUlB1xTetDxOzGe6a99HQghvIUQly3NHtK63Y9GOGB62LDYIrXUTwCLLayzvr7ZUNJvh4+RDr1q9WHBkARfSVWuZElfvPgh9BjZ9CUfX6J3G5vUL9eeh5v5MWH2EDUfU/+9lRUE3rw22PBa176MqwBohxF5gG/CnlHIxMBJ4RQgRB1QGrje5mQpUtsx/hb+Lkk15OuRpsvOymRU7S+8oFVPXj8Grjja28zWbOhAtc4QQfNinIcHeLrw8ZxdJVzL0jqRQuNZHCCHqCSEeFEL0vz7dbR0p5V4pZTMpZWMpZYiU8gPL/GNSyjApZbCUsp+UMtMyP8PyOtjyvk3eE1/drTpdqndhzsE5pGal6h2n4rFzgoemwLWLqplqGeBkZ+KbAc25lpXLiz/tIidXdfeit8Lc0fwuMAn4FugB/I+/T/8oVvBMo2dIy07j58M/6x2lYqrSBDq/BwcXw44f9U5j84J9XPnvgyFsPX6Jz/9UI/3qrTBHCo8CEcAZKeWTQBPAVPAqSnE0rNyQNlXaMOPADDJzVYduVtF6GNQMh+VvQdIBvdPYvL7N/Hk8LICvo4+y5pAaJElPhSkK6VLKXCBHCOEKnAVqWjeWMrjRYC6kX2DR0UV6R6mYDAboOwnsXWD+05B1Ve9ENu+9BxpSv4obI+bu5nRyut5xbFZhisIuSyui74HtwFZAddJjZWF+YYRUDuGHmB/IzcvVO07F5OoLD06G84dg2Rt6p7F5DmYjE/s3IydXMnz2TrLV9QVdFFgULN1MvC+lTJZSTgTuB56TUg4slXQ27Hq32vGp8fx56k+941RctSKgw6uwaybsKUyXXoo11fR2YexDjdh5Kpn/+0N1+aKHAouC5T6BxTe9jpNSqqOEUhIZGEmQWxBT901V3WpbU/hbENgGFo+AC3F6p7F5PRtXZWCb6kxad4w/DyTpHcfmFOb00VYhRHOrJ1FuYxAGngl5hoOXDrLp9Ca941RcRpPWm6rJHn5+CrJVe3m9vXN/fRpVc+fVebuJv3RN7zg2paBuLq63MGqPVhgOCSF2CiF2CSHU0UIp6VmzJz5OPkyNUd1qW5V7Nej7LSTtgxXv6J3G5tmbjEzs3xwJDJu9k8wcdV2ttBR0pLDV8tgHqAvcB/RDu0ehn5VzKRZmo5lBDQax7ew29pzfo3eciq1ON2j7ImybAvt/0zuNzQus7MRn/ZqwNyGFUb/tV6dQS0lBRUEASCmP5jeVUj4FeLjOw7jZufH9PjUegNVFjoJqodrdzpds8qb6MqVrQz9ejAxm7vZ4Zm45pXccm1DQTWjeQohX7vSmlPJzK+RR8uFkdqJ//f58u+dbjiUfo2YldZuI1Zjs4OHv4buOMHcgDF6hdY2h6GZE5zrEJKYwetF+6vm50jLorkO5KMVQ0JGCEXABXO8wKaWof73+OBgd+D5GHS1YnUd1rX+kpBhY8qrqH0lnBoPgf481I8DTiRdm7uRsimoIYE0FFYUzUsoPpJSj85tKLaECgIeDBw/VeYglx5Zw9upZveNUfLW7QKeRsGc27PhB7zQ2z93RzKQnW5CelcPzM3eoC89WdNdrCkrZMbCBds+gGoSnlHQaCcGdYdlISNihdxqbV9vXlc8eacLu+GR14dmKCioKUaWWQimUqi5Vua/mffxy5BeSM5L1jlPxGQxaNxiufjBvIFxVA8HorXtIFYZHqAvP1lTQIDtqBJIy6OmGT5Oek85PB3/SO4ptcPKER2bA1fPwy2BQ/VDpbkSXOkTU9Wb0ov1sPa6+pkpaoQbZUcqOYI9gwgPCmXVwFtey1Z2epaJqU7j/MzgWDas/0juNzTNaLjwHejrx3IztnLyoergtSaoolEODQwaTkpnCtAPq2kKpaf4kNB8EGz6HmF/0TmPz3B3NTH2qJXkSBk/bTkp6tt6RKgxVFMqhpj5N6R7UnUl7JhF7MVbvOLbjvv+DgNbw2zA4vVvvNDavhpcz3w5owYkLVxk+e6cayrOEqKJQTr3b+l08HDx4e8PbanS20mKyh0dnglNlmNMfUlUPnnprU6syH/cNYf2RC3ywWI2gVxJUUSin3O3d+aDdB8Qlx/Hlzi/1jmM7XLzh8Z8g/TLMfUL1qFoGPNoykCEdazL9r5NM23RC7zjlnioK5Vj7au15pM4jTD8wnW1nt+kdx3ZUaaz1qJqwTRuDQbWX193I7vXoXN+X0b/vJ1qN8VwsVisKQogAIcQaIUSsEGK/EOJly3xPIcSfQogjlkcPy3whhJgghIgTQuxVYzgUzquhrxLgGsC7G94lLStN7zi2o0FvbXCePbNhkzpS05vRIBj/WFPq+rkxfPYu9p9O0TtSuWXNI4Uc4FUpZX2gNTBMCNEAeBNYJaWsDayyvAboAdS2TEOAb6yYrcJwMjvxcfuPOXvtLOO2jdM7jm3p+AY06AN/joIDi/ROY/Oc7U18/1Qorg4mnvphmxqcp4isVhSklGeuD90ppUwFYoFqQG/gelvKaWjjNWCZP11qNgOVhBBVrJWvImnq05TBIYP5Le43Vp9arXcc22EwaKeR/EPh139BvDqFp7cq7o5MeyaMzOxcBv2wlctXs/SOVO6I0ug/RAgRBKwDQoBTUspKN713WUrpIYRYDIyVUm6wzF8FjJRSbr9lW0PQjiTw9fVtMWfOnCJlSktLw8XFpUjrlkU5MofPznxGcm4yb1d9G1dj+ezItjzuF3NWCs13voEx9xq7mn1CulPF+i1THvfJ4cu5fLItgyA3A6+3dMDeWPG6civOfomIiNghpQzN900ppVUntO63dwAPWl4n3/L+ZcvjEqD9TfNXAS0K2naLFi1kUa1Zs6bI65ZVhy8dls2mN5MvrnpR5uXl6R2nSMrtfrkQJ+XYICnHN5Uy7YLeaUpUed0ny/adlkFvLpaDf9wms3Ny9Y5T4oqzX4Dt8g7fq1ZtfSSEMAO/ALOklL9aZiddPy1kebzeVCABCLhpdX/gtDXzVTS1PWrzcvOXWRO/hoVHF+odx7ZUrgX958KV0/DTY5Cdrncim9c9pAqjezVkZWwS/1moelUtLGu2PhLAVCBW/nOUtkXAIMvzQcDCm+YPtLRCag2kSCnPWCtfRTWg/gBa+LZg7NaxnE5TNbVUBYTBg5O0pqo/Pw25qusFvQ1sE8SwiFr8tPUUY5YdVIWhEKx5pNAOeBKIFELstkz3AWOBLkKII0AXy2uApcAxIA6YDAy1YrYKy2gw8nH7jwF4d+O75El163+patAb7v8UDi+DhcMgT/3319trXesysE11Jq07xoRVcXrHKfMKGqO5WKR2wfhOV3duG6vBcp5rmLXy2JJqLtUY2XIkozaNYuaBmQxsOFDvSLal5bPaHc+rPwKHStBjHIiKd6GzvBBC8P4DDbmamcsXKw/jZGfkXx3VOOd3YrWioOirT3AfVsevZvzO8bSt2pZgj2C9I9mWDq9BejL89ZU2JkP4m3dfR7Eag0Ew7qFGZGTn8vHSWBztjAxoXV3vWGWS6uaighJC8H6b93Gxc+HtDW+Trc5vly4hoOtH0HQARI9Rdz2XASajgS8ebUpkPR/+szCG+TsS9I5UJqmiUIFVdqzMqDajiL0Uy7d7v9U7ju0RAh4Yr931vOJd2PSV3olsnp3JwNdPNKdtrcq8Pn8P87bF6x2pzFFFoYKLCoyid63eTNk3hT3n9+gdx/YYTfDQFEtheEcVhjLAwWxk6qCWtA/24o1f9jJbjfX8D6oo2ICRYSPxdfLl7fVvcyXrit5xbI/RrApDGeNgNjJ5YCgRdb15e8E+pv91Qu9IZYYqCjbA1c6V/7b/L6evnmboyqFqbGc9qMJQ5jiYjXz7ZAu6NPBl1ML9fL/huN6RygRVFGxEqF8on3b8lJgLMby0+iU1Wpsebi0M6z5VYzHozN5kZGL/5vQI8eODxQf44s/DNn+DmyoKNiSqehQftvuQrWe38mr0q2TnqRZJpe56YWj0CKz+EP54R93gpjM7k4EvH29Gvxb+jF91hHd/iyE3z3YLg7pPwcY8UOsB0nPS+XDzh7y1/i3GdRiH0WDUO5ZtMZqh73fg6AGbJ0L6Jej1pTZf0YXJaOCThxtT2cWeb9ce5dLVLL54tCkOZtv7t6GKgg16pO4jpOek8+n2T3E0OTK67WgMQh00liqDQbvT2dkL1nys3ejW7wcwO+qdzGYJIXizRz28XOz4aEksl69tZdLAUNwcbKtYq28CGzWo4SBeaPICv8X9xrit42z+PKouhIBOb8D9n8Hh5TC9D1y9qHcqm/dsh5r879GmbD9xmYe/2WRzI7ipomDDXmjyAoMaDGL2wdl8uUvdcaubls9qRwlndsOUSDh/WO9ENq9Ps2pMfyaMsykZ9J64ke0nLukdqdSoomDDhBC8Gvoq/er0Y/K+yUzZN0XvSLarYV94aglkXYWpneFYtN6JbF7bYC9+G9YOd0cz/Sdv4dedttEthioKNk4Iwbut36VnzZ6M3zmeWbGz9I5ku/xD4dlV4FoVZj4E23/QO5HNq+ntwoKhbWlR3YNX5u1h7LKDFb5lkioKCgZh4MN2HxIVGMXYrWNZcGSB3pFsl0d1GPwH1AyHxf/WxmRQo7jpqpKTHdMHh/F4WCDfrj3KwO+3cCGt4t7no4qCAoDJYOKTjp/Qrmo7Rm0axVe7vlID9OjFwR36z4OOr8OumfB9N7h8Qu9UNs1sNDDmwUZ88nBjtp+4TM8JG9hxsmJeZ1BFQbnBzmjHhMgJ9A3uy3d7v+Pl1S+TlpWmdyzbZDBC5Lvw+Fy4dAK+6wRH/tQ7lc17JDSAX4e2xc5k4NHvNvP9huMVruWeKgrKP9gZ7RjddjRvhb3F+sT19F/anxMpJ/SOZbvqdofnosE9AGb1gz/fg5wsvVPZtIZV3fn9xfaE1/Xhg8UHGDxtO+dTK87pJFUUlNsIIehfvz+Tu04mOSOZ/kv6sy5hnd6xbJdnTRi8ApoPhI3/01onXTiidyqb5u5oZvLAFrz3QAM2xF2gx/h1rD6YpHesEqGKgnJHLf1aMqfnHKq5VmP4quFM2Telwh0qlxt2TtBrAjw6C5Lj4buOsONH1aGejoQQPN2uBr8Pb4+Xiz3P/Lidd3/bR3pWrt7RikUVBaVAVV2qMr3HdLoFdWP8zvG8se4N1fW2nur3hBc2QUAr+P1lmPUwJKtBYvRU18+VhcPb8Wz7GszcfIru49exKe6C3rGKTBUF5a4cTY580vET/t383/xx4g8GLBtA7MVYvWPZLrcqMOBX6PEJnPwLJraGLd+p3lZ1ZG8y8m7PBswZ0hoB9J+yhZHz95Jyrfz1RKyKglIoQggGNxrM152/5lL6JR5f8jjjd45X4zLoxWCAVs/BsM0Q2BqWvQE/dIfzh/ROZtNa16zM8n935PlOtZi/M4HOX6xlyd4z5eq0q9WKghDieyHEOSFEzE3zPIUQfwohjlgePSzzhRBighAiTgixVwjR3Fq5lOJpX609C/sspGfNnkzZN4V+v/dj97ndeseyXZUCYcAvWlfcFw7DN221MRoyUvROZrMczEbe7FGPhcPa4etmz7DZO+k/eQsHz5aPoXCteaTwI9D9lnlvAquklLWBVZbXAD2A2pZpCPCNFXMpxeRu785H7T/i287fkpGTwcBlAxm7day61qAXIaDJYzBsGzTtD39NhC9bwM7p6pSSjkKqufPb0HZ82CeE2LNXuG/8et5bGEPytbLdpNhqRUFKuQ649Za/3sA0y/NpQJ+b5k+Xms1AJSFEFWtlU0pGu2rtWNB7AY/Ve4xZsbN4cNGDrE9Yr3cs2+XirQ3WM2SN1ox10YswOQJObNQ7mc0yGQ082bo60a+FM6B1dWZsPkn4p9FMWX+MjOyy2UpJWPNclxAiCFgspQyxvE6WUla66f3LUkoPIcRiYKyUcoNl/ipgpJRyez7bHIJ2NIGvr2+LOXPmFClbWloaLi4uRVpXud3RjKPMvjibcznnqOtQl96VehNgH3DP21H7pYRIic+5ddQ6Og37rItc8mjG8RpPkOpW+543pfZJyYlPzWPOwUz2X8zD00HQu5aZ9tVMGA3inrdVnP0SERGxQ0oZmt97ZaUoLAHG3FIU3pBS7iho+6GhoXL79tvqRqFER0cTHh5epHWV/GXlZjH30Fwm7Z1EcmYy99W4jxebvYi/q3+ht6H2SwnLugbbpsCGzyH9MtTrqXWf4VO/0JtQ+6TkbTp6gU+WH2J3fDI1vZz5d5c63N+oyj0Vh+LsFyHEHYtCabc+Srp+WsjyeM4yPwG4+WelP3C6lLMpB45EzQAAEChJREFUxWRntOPJBk+y9MGl/KvRv1h9ajUP/PYA47aO41JGxew8rMyzc4J2L8HLeyH8LTi2Fr5uA/MGQeJOvdPZrLa1vFgwtC2TnmyBySh46addRH0WzU9bT5GZo+9ppdIuCouAQZbng4CFN80faGmF1BpIkVKeKeVsSglxtXPlpeYvsbjvYnrX6s3sg7PpNr8bY7aMITEtUe94tsnBDcLfhH/vhfYj4Oga7XrDtAcgbpW6M1oHQgi6NvRj2csd+eaJ5rg6mHnr1310/GQNk9cdIy0zR5dc1myS+hPwF1BXCJEghBgMjAW6CCGOAF0srwGWAseAOGAyMNRauZTS4+vsy/tt32dB7wV0C+rGvMPzuP/X+xm5biSHLqn29Lpw8oTO78GIGOjyodaH0swH4bsOWjfdWaoFWWkzGgQ9GlVh0fB2zBgcRk0vFz5eGkvr/67ivYUxxJ0r3Z6KrXpNwdrUNYXy5ezVs8w4MIP5h+dzLecabau25ZG6j9DJvxMmgwlQ+6XU5WTBvnmw6Ss4H6uN5dB0AIQ+A17BgNonetgdn8z0TSdYvPcMWbl5tA/2YmCb6kTW88Fk1H7LW+uagioKSqlLyUxh7qG5zD00l3PXzuHj5MNDtR/iwdoPcnDbQbVf9CAlnNykXZSOXQR5OVCjEzR9gnX/3969B0d13Qcc//7uvvTYlZAQktCDhwGLR1TM2w5EBZvmMUmd1OnEretxH566ZJrU02mn43r8Z/8I0z869UwnCZ14SNo4DWNI4mSSGjtGIbWDeRmwESDbGJDAIEBCaNFrd+/pH3d1tSuthFbSalfs7zNz5t577rlnD/ubqx93996zN0ppemTkI0dqJtwID/DjI23896GLfNLdT1VJgK+ureNr6+u58P4RTQojaVKY3aJ2lIPtB9nTuoe3L7+NiLA8sJynNjzFwwseptBbmO0h5qeea/DuD+DYD6D7EjGrAE/jH8HvPQ6Lm5wfAFIzKhqzeeNMB3uOttF8rgPbwJMr/PzLn//BpPrTpJCCJoXc0tbTxr4P9rG3ZS9dsS6KvEVsX7idzy36HJvmbyLgCWR7iPnHtuHS77jy2r9R0/kODNyGUA2s/DKs+ENnziVNEDPu2u1+XjnWTtmdSzzxpYcn1cd4ScE7pdEpNU3qQ/U8u/ZZGrsbKVlRwi/O/4L9F/bz6kevUuQtYnPtZrbVb6OpronSQGm2h5sfLAsWbaa1IULN5v+Cc7+CU3vg6EvwzrehqAKWfxFWPgqLmsDrz/aI80JVSQF/u20pzc3tGelfk4LKKZZYrK9ez/rq9Ty/6XkOXz3Mm5fepLmtmdcvvo5HPKyrWse2+m1srd+a1oNxagp8hfCpx5wy0OP8XvSZn8P7e+H498EfdD5aWvoILN0OZYuyPWI1SZoUVM7ye/xsqd3CltotvPDgC5y+cZo3297kwKUD7Dyyk51HdlIXrGPT/E1smr+JDdUbqCisyPaw732B0HCCiPTD+QPwwX748A0490unzdylsOQRuO/3YcFDzq2walbQpKBmBUssGuc10jivkWfXPsvF2xf5bftvOXz1MPsv7GfvB3sBWDpnKRurN7Kmag2rK1ZTXVyNSPrzyqgJ8hVAwxecYgzc/MhJDh/92pml9fB3nXbzVsDCT8fLZueHglRO0qSgZqWFJQtZuHIhT658kpgd42znWQ59cojDVw+z74N9vHz2ZQAqCytZXbma1fNW01jRSEN5A8W+4iyP/h4l4jzbULEUHtwB0QG48i5cfMu53fXUj+Ho95y2pQugdg3UrIWaNVDzgPOMhMo6TQpq1vNYHlZVrGJVxSqebnyaiB2htbOVk9dPuuX1i6+77ReEFrC8fDnLy5fTUN7AivIVVBRW6BXFdPMGnDuUFjwIn/kHiEXh2nvOVN6XjzpzL7X8bLj93GVQuxaqVkHlKqhaCaH5TrJRM0aTgrrn+CyfmySeWPEEADf6btBys4UzN89wruscLTdb2H9xv3tMeUE5i0sXO6VksbteE6zBEv3V2mnh8cavCtYM19256VxNXDnuJInzv3GuKIYUzIkniRVQuRIq7oe5SzRZZJAmBZUXKgoraKproqmuya3rGeyhtauVs51nae1q5UL3Bd64+Aa3Bm65bQKeAAtLFrKoZBG1oVrqgnXUBmupCdZQE6zR5yemqnguLNvulCG9ndDRAtdanGVHi3Mr7EDCz1n6ipwfEiq/z0kS5Uuc5ZwFTsLQ5ycmTZOCylshf4h1VetYV7Uuqb6rv4uPuz8eLrc/5lzXOQ60HSBiR5LaVhZWUhOsoTZUS1VRFZVFlVQWVTKvcJ679Hl8M/nPmv2KymHRFqcMMQa62+Hmh07pPO98qd3R4tzxZCfMKCoeKKmF0roRpR5KayFYDYVlznMYahRNCkqNUFZQRllBGWur1ibV28amo7eDK+ErXA5fTionOk5wrfcaUXv0dMdlgTLmFTlJoqKwwuk/UOYu5xTMoTxQzpyCOQR9Qf1uIxURmFPvlCXbkvfFotB9CW6eh+62eGl3StshOH0lOWkAWF4oroRgJYSqnWWwKl4qnQfzisqd5FFYnlcP5mlSUGqCLLGoLq6murh6VMIAMMZwa+AWHb0dXO+7Tkdvh7Pee52OPmfZ2tVKV3/XqCuOIV7L6yaKUn8pIX+IkD9E0Bd014e2g/4gJf4Sdz3kD+G3/PmXVDze4Y+SUrFjEL42nCjCHc52uAPCV+H2Zed7jTvXwdip+/AH4wmiLDlZFJY5d00VlDjPbwSGlgnFH3LGOEvMnpEqleNExL3KaKBhzHbGGPqifXT2d3Jr4Ja77Orvoqu/y63rHujmcvgy4cEwPYM9hCNhDOPPVeYRD4XewlGlyFc0bn2Bp4CAJ0DAE8Dv8eP3+N31TwY/oa2nLXm/5cczWz63tzxQUuOU+o1jt7Nj0HsTeq46y75O5/uNvlsJ613O+q02Z73/1tiJJJGvKCFJBMFf7NT5ClMsU9Ul7isEbwBPtG/63qMEmhSUmmEiQpGviCJfUVrTdNjGpjfSS89gDz2RHjdZDK2HI2F6I730Rfvoi/bRGx1eDw+G6ejtcLeHyoTtG13lFW9S8hha+iwfXsvrLhPX063zWb6U+y2x8IoXy7LwiMcplsepEwuP5Rmuj++zxMJrefHI8Lrbj1hYYiHB+EdKEw6KDYNhZ+oPt9wesR2vS2wX6XOSS3cvRHqd7UgfRO5MLMkAVct2AF+Y+FgnSJOCUrOEJRZBv/NR0Xym/kSwbWz6o/30Rfvoj/UzGBtkMDbIQGyAgdiAu/7ue++ybPmypLrE9ZHLqB11S8SO0B/tJ2JHiNgRty5xmbh+tyuhTBuVVKzkxDKUUCyxEMRdB5LqRQSL+LY47cQriE+wQiEsSkFIauP0F18aG8sYxI4hxsYyNmJsd1uMTYMd4v5MvAcZ6FMpNQtYYrlXLOO2O2+xdcnWGRlTzI6NmzSGljETwzY2UTvqLI2zjNmxpPWYiRd7xPJux9v26GMT1o0xGAy2sYeXxmATXxobGxsM7nriMUOvbzBJ7RPbDPWZ6nUA5hYUZCQGmhSUUjnDY3lmz3cVWdbc3JyRfvVGXaWUUi5NCkoppVyaFJRSSrlyKimIyOdF5JyIfCgiz2V7PEoplW9yJimIiAf4D5wbb1cCfyoiK7M7KqWUyi85kxSAjcCHxpjzxphB4H+AL2d5TEoplVdy6ZbUWqAtYbsd2DSykYg8AzwDUFVVNenbssLhcMZu6VKTp3HJPRqT3JSpuORSUkg1i9eoxxuNMbuAXQDr1683W7dundSLNTc3M9ljVeZoXHKPxiQ3ZSouuZQU2oH6hO064Mp4Bxw7duyGiFwcUV0KdKdoPrK+ArgxiXFOh7HGmOl+Jtr+bu3G2z/R93+sumzFJVsxSeeYycZlqvV6rky+Xa6eKwvH3GOMyYmCk6DOA4sBP3ASWDWJfnZNpB44msV/a8oxZrqfiba/W7vx9k/0/R+nLitxyVZMZiIuU63Xc2X6Y5JuXGbyXMmZKwVjTFREvgG8BniAl4wxpyfR1c/TrM+G6RpLuv1MtP3d2o23P533X2OS3jGTjct01WeDnisTe51pI/GMk3dE5KgxZn22x6GSaVxyj8YkN2UqLrl0S+pM25XtAaiUNC65R2OSmzISl7y9UlBKKTVaPl8pKKWUGkGTglJKKZcmBaWUUi5NCnEiUiwi3xeR/xSRP8v2eBSIyH0i8j0ReSXbY1HDROQr8fPkZyLy2WyPR4GIrBCR74jIKyLy9an0dU8nBRF5SUQ6ROT9EfWppuh+DHjFGPPXwKMzPtg8kU5MjDM54tPZGWl+STMuP42fJ38BPJ6F4eaFNGNyxhizA/gaMKXbVO/ppADsBj6fWDHOFN11DE/IF5vBMeab3Uw8Jmrm7Cb9uLwQ368yYzdpxEREHgX+D/j1VF70nk4KxpiDQOeI6rGm6G7HSQxwj78v2ZRmTNQMSScu4tgJ/MoYc3ymx5ov0j1XjDGvGmM+DUzp4+98/OOXaoruWmAf8FUR+Ta59Zh/PkgZExGZKyLfAdaIyD9nZ2h5baxz5ZvAduCPRWRHNgaWx8Y6V7aKyIsi8l3gl1N5gZyZ+2gGpZyi2xhzB/jLmR6MAsaOyU1A/+hkz1hxeRF4caYHo4CxY9IMNE/HC+TjlULaU3SrjNOY5CaNS+7JeEzyMSkcAZaJyGIR8QN/Arya5THlO41JbtK45J6Mx+SeTgoi8iPgd0CDiLSLyNPGmCgwNEX3GWDPJKfoVpOgMclNGpfck62Y6IR4SimlXPf0lYJSSqn0aFJQSinl0qSglFLKpUlBKaWUS5OCUkoplyYFpZRSLk0KSgEiEo4vF4nIE9Pc9/Mjtt+ezv6Vmk6aFJRKtghIKynEpzMeT1JSiM9kqVRO0qSgVLJvAZ8RkRMi8vci4hGRfxWRIyJySkT+BiA+K+UBEXkZeC9e91MROSYip0XkmXjdt4DCeH8/jNcNXZVIvO/3ReQ9EXk8oe/m+K9onRWRH4pIqonQlJp2+ThLqlLjeQ74R2PMlwDif9y7jTEbRCQAvCUi++NtNwKfMsZ8HN/+K2NMp4gUAkdEZK8x5jkR+YYx5oEUr/UY8ACwGqiIH3Mwvm8NsApnsrO3gM04P6CiVEbplYJS4/ss8JSInADeAeYCy+L7DickBIC/E5GTwCGcmSyXMb4twI+MMTFjzDXgN8CGhL7bjTE2cALnYy2lMk6vFJQanwDfNMa8llQpshW4M2J7O/CQMaZXRJqBggn0PZaBhPUYeq6qGaJXCkol6wFCCduvAV8XER+AiNwvIsUpjisFuuIJYTnwYMK+yNDxIxwEHo9/bzEPaAIOT8u/QqlJ0v99KJXsFBCNfwy0G/h3nI9ujse/7L0OfCXFcf8L7BCRU8A5nI+QhuwCTonIcWNM4u/n/gR4CDgJGOCfjDFX40lFqazQqbOVUkq59OMjpZRSLk0KSimlXJoUlFJKuTQpKKWUcmlSUEop5dKkoJRSyqVJQSmllEuTglJKKdf/A5N/IhwO0aiJAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Learning rate values to test\n",
        "lr_test = [1e-5,1e-4,1e-3]\n",
        "ntest = len(lr_test)\n",
        "\n",
        "# Strings for the legend\n",
        "leg_str = []\n",
        "losshist = []\n",
        "\n",
        "for i, lr in enumerate(lr_test):\n",
        "    # Run the optimizer\n",
        "    w, f0, hist = grad_opt_simp(feval, winit, lr=lr, nit=nit)\n",
        "    losshist.append(hist['f'])\n",
        "\n",
        "    # Plot the results\n",
        "    t = np.arange(nit)+1\n",
        "    plt.semilogx(t, hist['f'])\n",
        "    leg_str.append(\"lr=%12.2e\" % lr)\n",
        "\n",
        "    # Measure the test accuracy\n",
        "    yhat = predict(Xts1,w)\n",
        "    acc = np.mean(yhat == yts)\n",
        "    print(\"lr=%12.2e  Test accuracy = %f\" % (lr, acc))\n",
        "\n",
        "plt.grid()\n",
        "plt.legend(leg_str, loc='upper right')\n",
        "plt.xlabel('Iteration')\n",
        "plt.ylabel('Training loss')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tLTs5wo4zVqe"
      },
      "source": [
        "You can see that increase the learning rate more to speed up the convergence time.  Again, although it depends on the train-test split, generally you will also get higher test accuracy with the faster learning rate.  But, you cannot increase the learning rate too high as the optimization may go unstable."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nDAaKpSxzVqe"
      },
      "source": [
        "## Adaptive Step Size\n",
        "\n",
        "The above example shows that gradient descent is sensitive to the step size.  We now consider a slight variant of gradient descent with an adaptive step-size using the Armijo rule in class."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hJHuXKc7zVqe"
      },
      "outputs": [],
      "source": [
        "def grad_opt_adapt(feval, winit, nit=1000, lr_init=1e-3, lr_min=1e-6):\n",
        "    \"\"\"\n",
        "    Gradient descent optimization with adaptive step size\n",
        "\n",
        "    feval:  A function that returns f, fgrad, the objective\n",
        "            function and its gradient\n",
        "    winit:  Initial estimate\n",
        "    nit:    Number of iterations\n",
        "    lr_init:     Initial learning rate\n",
        "    lr_min:      Minimum learning rate\n",
        "    \"\"\"\n",
        "\n",
        "    # Set initial point\n",
        "    w0 = winit\n",
        "    f0, fgrad0 = feval(w0)\n",
        "    lr = lr_init\n",
        "\n",
        "    # Create history dictionary for tracking progress per iteration.\n",
        "    # This isn't necessary if you just want the final answer, but it\n",
        "    # is useful for debugging\n",
        "    hist = {'lr': [], 'w': [], 'f': []}\n",
        "\n",
        "    for it in range(nit):\n",
        "\n",
        "        # Take a gradient step\n",
        "        w1 = w0 - lr*fgrad0\n",
        "\n",
        "        # Evaluate the test point by computing the objective function, f1,\n",
        "        # at the test point and the predicted decrease, df_est\n",
        "        f1, fgrad1 = feval(w1)\n",
        "        df_est = fgrad0.dot(w1-w0)\n",
        "\n",
        "        # Check if test point passes the Armijo condition\n",
        "        alpha = 0.5\n",
        "        if (f1-f0 < alpha*df_est) and (f1 < f0):\n",
        "            # If descent is sufficient, accept the point and increase the\n",
        "            # learning rate\n",
        "            lr = lr*2\n",
        "            f0 = f1\n",
        "            fgrad0 = fgrad1\n",
        "            w0 = w1\n",
        "        else:\n",
        "            # Otherwise, decrease the learning rate\n",
        "            lr = lr/2\n",
        "            lr = np.maximum(lr, lr_min)\n",
        "\n",
        "        # Save history\n",
        "        hist['f'].append(f0)\n",
        "        hist['lr'].append(lr)\n",
        "        hist['w'].append(w0)\n",
        "\n",
        "    # Convert to numpy arrays\n",
        "    for elem in ('f', 'lr', 'w'):\n",
        "        hist[elem] = np.array(hist[elem])\n",
        "    return w0, f0, hist\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fC8y2FpTzVqf"
      },
      "source": [
        "We can now run the new optimizer."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Z4rXWI_MzVqj",
        "outputId": "c63f2050-e76c-4ae5-d9f6-213afebf7974"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUdb7/8dcnvQNJSIAk9CYgIoQiohus6Lpi7+yuumJZ27rlunvXq78tV+92XQtiWWwr61pZ7KIREFAS6b0GAoFACCEFUj+/P2YSAiRhkpnJnEk+z8cjJvOdUz7k68w758z3fI+oKsYYY4zThAS6AGOMMaYpFlDGGGMcyQLKGGOMI1lAGWOMcSQLKGOMMY4UFugC/CE5OVn79u3b5vXLy8uJjY31XUGm3VjfBS/ru+Dki37Lzc3dr6rdj2/vkAHVt29fcnJy2rx+dnY2WVlZvivItBvru+BlfRecfNFvIpLXVLud4jPGGONIFlDGGGMcyQLKGGOMI1lAGWOMcSQLKGOMMY5kAWWMMcaRLKCMMcY4kgWUMcYYR2r3gBKRF0WkUERWN2pLFJFPRWST+3s3d7uIyBMisllEVorI6Pau1xhjTGAE4ghqFjDluLYHgXmqOgiY534McBEwyP01HXimnWo0xhgTYO0eUKo6HzhwXPNU4CX3zy8BlzVqf1ldlgBdRaRn+1RqjDEmkJwyF1+qqhYAqGqBiKS429OAnY2Wy3e3FRy/ARGZjusoi9TUVLKzs9tcTFlZmVfrm8Cxvgte1nfByZ/95pSAao400aZNLaiqM4GZAJmZmerN5IU2aWXwsr4LXtZ3wcmf/eaUUXx760/dub8XutvzgYxGy6UDu9u5NmOMMQHglICaA/zA/fMPgPcatX/fPZpvAlBSfyrQGGNMx9bup/hE5HUgC0gWkXzgYeAx4A0RuRXYAVztXvwD4GJgM1AB3Nze9RpjjAmMdg8oVb2+mafObWJZBX7s34qMMcY4kVNO8RljjDHHsIAyxhjjSBZQxhhjHMkCyhhjjCNZQBljjHEkCyhjjDGOZAFljDHGkSygjDHGOJIFlDHGGEdy+mzmbXKgvIpXl+S1ef1dBTVMqq0jLNTy2xhjAqVDBtSug4f59burT75gC5YeXMLfrhtFercYH1VljDGmNTpkQJ3SM4FP/vuEqf089vx/FvL6hlIuenwBj10xku+OtJv4GmNMe+uQARUWIqTER7V5/Ym9wpg2ZSL3zl7Gj//5LfM3ZvDwpcOIieiQvy5jjHEk+5ClGRmJMbxx+xncPXkgb+Tu5JK/L2T1rpJAl2WMMZ2GBVQLwkND+NmFQ3jtR+Mpr6zhiqcX8cLCbbjuAmKMMcafLKA8MHFAMh/ddzbfGdKd385dy82zlrK/rDLQZRljTIdmAeWhbrERzJw2ht9OHc7iLUVM+dsC5m/cF+iyjDGmw7KAagURYdoZfZlz9yQSY8P5/ovf8Pv311JVUxfo0owxpsOxgGqDIT3imXP3JKZN6MNzC7Zx5TOL2La/PNBlGWNMh+KYgBKRISKyvNHXIRG5X0QeEZFdjdovDnStAFHhofz2shHMnDaGncUVfPeJBfw7Z6cNoDDGGB9xTECp6gZVHaWqo4AxQAXwjvvpv9Y/p6ofBK7KE10wvAcf3ncWI9O78PM3V3Lv7OUcOlId6LKMMSboOSagjnMusEVV2z6hXjvq2SWa1340gZ9fOIQPVhVw8eML+HZHcaDLMsaYoCZOPCUlIi8C36rqkyLyCPBD4BCQA/xUVU949xeR6cB0gNTU1DGzZ89u8/7LysqIi4tr07qbD9by7IpKio4olw0M55L+4YSItLkW0zre9J0JLOu74OSLfps8eXKuqmYe3+5VQInIACBfVStFJAsYCbysqge92GYEsBsYrqp7RSQV2A8o8Fugp6re0tI2MjMzNScnp60lkJ2dTVZWVpvXP3Skml+/s5o5K3YzoX8if712FD27RLd5e8Zz3vadCRzru+Dki34TkSYDyttTfG8BtSIyEHgB6Af808ttXoTr6GkvgKruVdVaVa0DngPGebl9v0uICufx60bxp6tPY2V+CRc9voCP1+wJdFnGGBNUvA2oOlWtAS4H/qaqPwG8nfr7euD1+gci0nh7lwPe3UejnYgIV41JZ+49k8joFsPtr+Ty63dXcaS6NtClGWNMUPA2oKpF5HrgB8Bcd1t4WzcmIjHA+cDbjZr/ICKrRGQlMBn4SVu3Hwj9u8fx1p0TmX52f15dsoNLn1zIhj2lgS7LGGMcz9uAuhk4A/i9qm4TkX7Aq23dmKpWqGqSqpY0apumqqeq6khVvVRVC7ysud1FhIXwq4tP4eVbxnGgvJrvPbmQlxdvt2umjDGmBV4FlKquVdV7VfV1EekGxKvqYz6qrcM5e3B3Prr/LM4ckMT/vLeG217O5UB5VaDLMsYYR/IqoEQkW0QSRCQRWAH8Q0T+4pvSOqbkuEhe/OFY/ueSYczfuI+LHp/Poi37A12WMcY4jren+Lqo6iHgCuAfqjoGOM/7sjo2EeGWSf1458cTiY0M48bnv+YPH62nutYmnTXGmHreBlSYe5TdNRwdJGE8NLxXF+beM4lrMzN4OnsLV89YzI6iikCXZYwxjuBtQP0G+BjXtERLRaQ/sMn7sjqPmIgwHrtyJE/dMJot+8q4+IkFvLd8V6DLMsaYgPN2kMS/3aPr7nQ/3qqqV/qmtM7luyN78uF9ZzG0Rzz3zV7OA28sp7yyJtBlGWNMwHg7SCJdRN4RkUIR2Ssib4lIuq+K62zSu8Uwe/oE7jt3EO8u28UNzy2xUX7GmE7L21N8/wDmAL2ANOA/7jbTRmGhIfzk/MHMnJbJ+j2lXDVjEbsOHg50WcYY0+68DajuqvoPVa1xf80Cuvugrk7vvGGpvHLrePaVVnLVM4vYXGizTxhjOhdvA2q/iNwkIqHur5uAIl8UZmBcv0T+Nf0MauqUq2YsZpndY8oY04l4G1C34BpivgcoAK7CNf2R8ZFhvRJ4646JJESFc+PzXzN/475Al2SMMe3C21F8O9zz43VX1RRVvQzXRbvGh3onxfDmnWfQJymWW19ayn9W7A50ScYY43f+uOX7A37YZqeXEh/F7OkTOD2jG/fOXsYri7cHuiRjjPErfwSU3d/cT7pEh/PyreM4d2gqD723hr9+utFmRDfGdFj+CCh7x/SjqPBQZtw0mqvGpPP4vE08PGcNdXX2KzfGdDxhbVlJREppOogEiPaqInNSYaEh/PGqkSTGRjBz/lYOlFfxl2tGERHmj783jDEmMNoUUKoa7+tCTOuICL+6+BSSYiN49MP1lByu5tlpY4iJaFOXGmOM49if3EHu9u8M4A9XjWTRliJueO5rim1qJGNMB2EB1QFck5nBjJvGsLbgEFc/u5iCEpsayRgT/BwVUCKyXURWichyEclxtyWKyKcissn9vVug63Si84el8sot49hbcoSrnlnMln1lgS7JGGO84qiAcpusqqNUNdP9+EFgnqoOAua5H5smjO+fxOzbJ1BZU8fVMxazMv9goEsyxpg2c2JAHW8q8JL755eAywJYi+MN79WFt+48g9jIUK6fuYSFm/YHuiRjjGkTcdKFniKyDSjGNYT9WVWdKSIHVbVro2WKVfWE03wiMh2YDpCamjpm9uzZba6jrKyMuLi4Nq/vBAeP1PHn3EoKyuq4/bRIxvboHKP7OkLfdVbWd8HJF/02efLk3EZnzRo4LaB6qepuEUkBPgXuAeZ4ElCNZWZmak5OTpvryM7OJisrq83rO0XJ4Wp+9NJScvKK+c3UEUyb0CfQJfldR+m7zsj6Ljj5ot9EpMmActQpPlXd7f5eCLwDjAP2ikhPAPf3wsBVGFy6RIfz8i3jOXdoCg+9u5r/enMlR6prA12WMcZ4xDEBJSKxIhJf/zNwAbAa1x17f+Be7AfAe4GpMDhFR4Ty7LRM7j1nIP/K2cmVzyxiR1FFoMsyxpiTckxAAanAQhFZAXwDvK+qHwGPAeeLyCbgfPdj0wqhIcIDFwzhxR9mkl98mEv+voDP1u4NdFnGGNMixwSUqm5V1dPcX8NV9ffu9iJVPVdVB7m/Hwh0rcHqnKGpzL1nEr2TYvjRyzn88eP11NpEs8YYh3JMQJn2kZEYw5t3TOT6cb156ostTHvha/aXVQa6LGOMOYEFVCcUFR7Ko1ecyh+vGkluXjHffWIBuXl2YGqMcRYLqE7s6swM3r5rIlHhoVz77BJeXLjNboBojHEMC6hObnivLsy5exJZQ1L4zdy13P36MsoqawJdljHGWEAZ1/VSM6eN4b+mDOXDVQVMfXIhm/aWBrosY0wnZwFlAAgJEe7MGsCrPxpPyeFqpj71Fe8t3xXosowxnZgFlDnGxAHJzL3nLIb1TOC+2ct54I3l7DxgF/YaY9qfBZQ5QY8uUbw+fQJ3Zg1g7soCJv8pmwffWmlBZYxpVxZQpknhoSH815ShzP/5ZG4c35u3v93F5D9l88u3V7HroN2x1xjjfxZQpkU9ukTx/6aO4MtfZHH9uN68lZtP1h+/4L/fWcVuCypjjB9ZQBmP9OwSzW8vG0H2z7O4JjODN3J2kvXHbB56dzUFJRZUxhjfs4AyrdKrazS/v/xUvvhZFleOSef1b3bwnT9k8z/vrWbrvrJAl2eM6UA6x21Wjc+ld4vh0StO5a6sATydvZl/fr2DlxfncVpGV644PY1LRvYkKS4y0GUaY4KYHUEZr2QkxvDoFSP56sFz+O+LT6Gqpo6H56xh/P/O49ZZS5m7crfdJNEY0yZ2BGV8IjUhitvO7s9tZ/dn/Z5DvLNsF+8t28289YXER4Zx0ak9uPz0dMb3SyQkRAJdrjEmCFhAGZ8b2iOBX16UwC8uHMqSrUW8/e0u3l9ZwBs5+aR1jWbqqF5cfnoag1LjA12qMcbBLKCM34SGCGcOTObMgcn87rIRfLJ2D+8s28Wz87fydPYWRqQlcPnp6XzvtJ6kxEcFulxjjMNYQJl2ER0RytRRaUwdlca+0kr+s2I37yzbxW/nruV/P1jHpIHJfPfUnmQN7W5hZYwBLKBMAHSPj+SWSf24ZVI/NheW8s6yXby7bDe/eGslAKemdWHy0BTOGZrCyLQu9pmVMZ2UYwJKRDKAl4EeQB0wU1UfF5FHgNuAfe5Ff6WqHwSmSuNrA1Pi+fmFQ/nZBUNYV1DKFxsK+Xx9IU9+vokn5m0iOS6C03t3IzzUs5A6sP8InxSvIi4yjNiIMGIjQwmRY9cNCxUiQkMIDw0hPCyEiNAQIsKE8FDXz/Vt4aEhhIaAiCBAiAghIoi4Zn8PEY4+dj8XIq7lG3+v33+IHLuOiAWvMS1xTEABNcBPVfVbEYkHckXkU/dzf1XVPwWwNuNnIsKwXgkM65XAjycPpLi8ivmb9vH5+kLWFRzC0xv9lpTVsXnNHkqP1FBZU+ffon3g+MASoD63BHc7R8NM3P85pu24x+3BH3uqqq4iYsGnx7S1b4a3z87a89/UHrtKi64hK8s/23ZMQKlqAVDg/rlURNYBaYGtygRKt9iIhs+sWiM7O5ss96ulpraO8qpaaBRuilJTp1TX1lFdo1TV1lJVo1TV1rnb6qh0f6+uVerU9aVKk9/rGh4f/bm2zrXDo8sdXeaYx+Be7+i69bXWP9eoyf2zHhPW2rCdtvyW20bxz85279pNr7QeR/fTrv+mdtpPO/6b2utfVXNwj9+2Ldq+vzGPiEhfYD4wAngA+CFwCMjBdZRV3MQ604HpAKmpqWNmz57d5v2XlZURFxfX5vVN4FjfBS/ru+Dki36bPHlyrqpmHt/uuIASkTjgS+D3qvq2iKQC+3H9OfBboKeq3tLSNjIzMzUnJ6fNNTT+K9wEF+u74GV9F5x80W8i0mRAOWqqIxEJB94CXlPVtwFUda+q1qpqHfAcMC6QNRpjjGkfjjmCEtcnvC8BB1T1/kbtPd2fTyEiPwHGq+p1J9nWPiCviae6ACUetCXjOmoLtKZqa+9ttWa9ky3b1udb0+6EvvNlv3mzPU/X82S5lpZp7XPNLd/R+s4Jr7mTLeOU11wfVe1+Qqs2fHgb2C9gEq7TeCuB5e6vi4FXgFXu9jm4TvG1dR8zPWzLCfTvo7na2ntbrVnvZMu29fnWtDuh73zZb+3Rd54s19IyrX2uhf7sUH3nhNdcW/vOKa85J43iW0jToyJ9ec3Tfzxscwpf1tbWbbVmvZMt29bnW9seaL6uy99958lyLS3T2uec2m/Q8V5zJ1vG0a85x5zicxIRydEmPrAzzmd9F7ys74KTP/vNUYMkHGRmoAswbWZ9F7ys74KT3/rNjqCMMcY4kh1BGWOMcSQLKGOMMY5kAWWMMcaRLKCMMcY4kgWUB0QkVkReEpHnROTGQNdjPCci/UXkBRF5M9C1GM+JyGXu19t7InJBoOsxnhORU0Rkhoi8KSJ3erOtThtQIvKiiBSKyOrj2qeIyAYR2SwiD7qbrwDeVNXbgEvbvVhzjNb0napuVdVbA1OpaayV/fau+/X2Q+DaAJRrGmll361T1TuAawCvro/qtAEFzAKmNG4QkVDgKeAiYBhwvYgMA9KBne7FatuxRtO0WXjed8Y5ZtH6fvu1+3kTWLNoRd+JyKXAQmCeNzvttAGlqvOBA8c1jwM2u//qrgJmA1OBfFwhBZ34d+YUrew74xCt6Tdx+T/gQ1X9tr1rNcdq7WtOVeeo6kTAq49E7M32WGkcPVICVzClAW8DV4rIMzh7HrHOrMm+E5EkEZkBnC4ivwxMaaYFzb3m7gHOA64SkTsCUZg5qeZec1ki8oSIPIuXc6k6ZrJYh2hqslpV1XLg5vYuxrRKc31XBNgbnHM1129PAE+0dzGmVZrru2wg2xc7sCOoY+UDGY0epwO7A1SLaR3ru+Bk/Ra8/N53FlDHWgoMEpF+IhIBXIfrHlTG+azvgpP1W/Dye9912oASkdeBxcAQEckXkVtVtQa4G/gYWAe8oaprAlmnOZH1XXCyfgtegeo7m83cGGOMI3XaIyhjjDHOZgFljDHGkSygjDHGOJIFlDHGGEeygDLGGONIFlDGGGMcyQLKGGOMI1lAGWOMcSQLKGOMMY5kAWWMMcaRLKCMMcY4kgWUMcYYR+qQNyxMTk7Wvn37tnn98vJyYmNjfVeQaTfWd8HL+i44+aLfcnNz96tq9+PbO2RA9e3bl5ycnDavn52dTVZWlu8KMu3G+i54Wd8FJ1/0m4jkNdVup/iMMcY4kgWUMcZ4SVXZV1oZ6DI6HAsoY4zx0r9z8xn7+89Yvask0KV0KBZQxhjjpYWb9gOwubAswJV0LBZQxhjjI4oGuoQOxQLKGGOMI3XIYeam4/h8/V5ythdz9zkDiYnw7/+uX2woZOeBCr5/Rl+Pls/NK2Zl/kFuPrNfi8vN+HILmX26kdk38aTbLCg5zKxF23ng/MFEhoUe89xz87eS2bcbp/fu1tBWXF7FM19u4f7zBjX8fsoqa3hi3ibuPmcgCVHhACzYtI89JUe4OjMDgE17S/lk7V7uyhqAiLBoy37yDxzmmrEZ1NUpj8/bxFVj0slIjOGDVQWEh4Zw/rBUFm7az76yI1x+ejofrd5DiMD5w1KZ8eVWzjslheiIUGZ/s5OaOmX62f35ZlsRYSEh9E2O5bN1e7l+XG+eyd7CjycP4KkvtnDbWf14fuE2pk3owxsbqkgfVspHq/cwunc3CksrOVhRxcLN+7nlzH6s3l1CXGQ4ibERrC04RGpCJNHhoewrraSypo7E2Ahqaut489t8LhnZi+qaOsqragkPFQSoqVNq3V81dcpn6/aSNaQ7tXVQV6fUqeu56jqluqaO6to6lmwtYnSfblTV1BEWKoSHhnC4qpbiiio27i1jYEociTER7C09AoBrT8ZXLKCMIx2uquV376/lta93APDFhn3MnDaGjMQYv+xv9a4Sbn8ll6qaOpJiI/nuyJ4tLl9WWcNdr+Wy91AlA1PiOGvQCdcYArB0+wEe+3A9Q3vE8+F9ZyHS8hvYjOwtvLQ4j5FpXY+pIa+onN9/sI7vDO7OS7eMa2j/d+5OZs7fyvBeCUwdlQbAR6v3MHP+Vvonx3LduN4A/PmTjWwuLOOK0emEhgizFm3nta93cOlpvchIjOHJzzezYudBrhqTzsbCUh6ftwlV5YELhvDHjzcQFR7K+cNSefKLTWwuLOfy09P58ycbCAsNYUyfbvzfR+vZX1bJ8p0Hyc0rBlzhuWxnMdHhoUzon8Sz87cSFRbCjC+3EBYizPhyCyWHq3n9mx1s31/Oh9uqKX53DYu3FjGsZwLFFVUUlLje+D9bVwhAWIhQU3fy02irdx066TIAry7ZcdJlFrg/X2rK8Z851XpQm/GcneIzjrN6VwmX/H0Br329g+ln9+e572eSX1zBpU8uZNHm5t8s2upgRRV3vJpLUmwEI9O78ODbK8kvrmhxncc/28jeQ5Ukx0Xyu7nrqKmta3K5v322kRCB9XtKmd/CGx3Akepa3lm2C4B3luUf89zclQUALNqyn5LD1Q3t9W/c2Rv2NbQt3LTPvWwRAOWVNazaVUJZZQ3r97jeuOtDZPnOg9TVKavySyivqmXr/nLWuN/c1xYcoqKqhu1F5WwpLKOmto5Ne8vYX1bJ/rJKtheVs21/GVv2lQOwfX85+8uODrUuOVzNzgOHyS92fQEs3e7ab07eAQCW7XA9XrHzIEDDKLhNhaUN4dSYJ+EUSEdqagNdQodiAWUco65OefbLLVz+9FeUVdbw2o/G86uLT+H8YanMuXsSSXGRTHvxG15cuA1V37xR1dUp981eTuGhSp65aQxP3TAaFO6bvbzZ0Nmwp5QXv9rOdWMz+N1lw9mwt5TZS3eesNw32w7w1eYifnbhEFITInn2yy0t1vLR6j0cOlLD6N5dyd6wj6JGb/ZzVxaQGBtBda3yxXpXKB2sqCI3r5jQEGH+xn2u01R1ykJ3iC/aUoSqkptX3PCXfc72YkqPVLNxbyngCqhtReWUVtYAroBYs9sVUOsKStm4twxVqKqtIzevmKLyKgC+WF9Ida1ypLqORVtc+9tWVH7Mv2dfWSWHq2spKq9ig3t/32x3BdOKna4gqq9jtzuM6uuornV2EDXncJUFlC9ZQBlHKCg5zE0vfM2jH67nvFNS+ei+szlzYHLD8/2SY3n3x2dy7tAUfjN3LT/99wqOVJ/8zaCqpumQqfe3eZv4cuM+Hr50GKMyupKRGMPvLh9Bbl4xj8/bdMLyqspD760mPiqMX0wZyoXDezC+XyJ/+XTjMUc2AH/9dCPJcZHcPLEft5zZj0VbiliVf/Q6mZraumNOCc1euoPeiTH87xWnUlOn/GfFbgC27CtjXcEh7soaQGpCJB+t3gO4jppq65RpE/pQVF7Fql0lrN9Tyv6yKsb1TWR/WSWbC8v4elsRoSFCclwkS7cfYMXOEuoUIsNCWL7zICvzDzbUsGpXCWt2u2rcdfAw32wraniu/igO4OM1exp+/mzdXgB2Hqig9EjN0W01+rfWnwqrv5j1sLvvHH5A1Gqe/D9pPGcBZQKqtk55ZfF2LvjrfJbvPMgfrhzJ0zeOpltsxAnLxkWGMeOmMTxw/mDe/nYX1z+3hGL3X/RNmf3NDkY8/DGvLmlymi8+XrOHJ9yDAW5wf1YDMHVUGlePSefJLzaz1P0Xf723v93FN9sO8IsLh5IYG4GI8NAlwyiuqOKpLzY3LLdkaxGLtxZxx3f6Ex0RyvXjexMfGcaz811HUarKdTOXcN/sZQBs21/Okq0HuHZsBkN7JDCsZ0LD6b65KwoQgUtG9uLC4T34cuM+DlfV8tm6vSTHRXL3OQMRcQXWws2u03s/nzIEgMVbi/h66wFGpHVh4oAklm4/QG5eMSIwdVQvVu8qITfP9TnRyPQurMovYW3BITISowF4Z9luosNDCRH4cLUroEKEY05X1n/eU12rHGjUH1XNHIF2ZH/6ZGOgS+hQLKBMwKzKL+GKp7/ioffWMDK9Cx/cexbXjM1ocSBBSIhw77mDmHHTaNbsPsSVMxY1+XnRP7/ewYNvryIqPISH3lvN3JW7j3l+9a4S7p+9nFEZXfndZSNO2Ocjlw4no1sMP31jBeXu0077yyr57ftrGd27K9eNzWhYdkRaF644PZ1Zi7az66Drs5a/fbaR7vGR3DShDwAJUeHcML43H6wqYOeBCnLyisnJK+Z99+M3cnYSGiJcNSYdgCtGp7Eiv4Qt+8qYu3I3Y/sk0qNLFFOG9+BwdS2fry/ky437OGdod5LjIhmZ3pXsjYUs2LSfQSlxjO2bSHq3aD5fX8iK/INM6JfI2L7d2Huokv+s3M3glHjOGtSdypo65izfzYi0BE5L70rujmJKj9Rw1WjXv29dwSGG9oynb3Is+8uqiI8K45SeCVTV1JEcF0l8pGucVdeY8Ibfx0/OG8yN448GvrdSEyJ9ti0TXCygTLs7dKSaR+asYepTC9l18AiPXzeKV28dT99kz6fsnzKiJ6/eOp79pZVc+cyihg//AV77Oo9fvbOKyUO6s+C/ziGzTzd+8q/lLHAPHthTcoRbX1pKYmwEM78/hqjw0BO2HxsZxp+uPo2dxRU8+uE6AH47dy3llTX835UjCQk5NtAeuGAwAH/5ZCOLtxSxZOsB7vzOgGO2ffOZ/QgNEZ5fsJVZX20nPjIMAV5ZksebuflMHpJCakIUAJee1osQgcc+XM+mwjIuOc01om9cv0S6xoTzh4/XU3qkhnNPSQUga3B3lu88yNfbDjBpkOvU6MQBSWRv2Ed1rTK+f2LDMPfNhWWM7tOVURld3f1Rw8j0rpya1qXhlGPWkO4kx7mOYof2iGdIajwAg1PjGZQSB8DAlFj6d49t2H+9c09JYZL79GxyXCRR4a63mfgoV5iFh3o+FDs6PLRhdKLpfCygTLtRdX2uct6fv+Slxdu5aUIf5v30O0wdlXbS4ddNGdcvkX/fMRFBuHrGYpZsLeLzHdX89zurOWdoCjOmjaFLdDjP/2AsA7rHcfsruSzeUsSPXl5K2ZEanv9BJinxUS1u/0eT+vHqkh08+sE63lu+m7uyBjLI/WbdWFrXaH44sS9vL7RSlxUAABeESURBVMvnofdWkxIfyQ3HHUX06BLF1FFp/CtnJx+t2cMN43tz7impPL9gK/tKK485KktJiGLSoO58unYvIQIXjXAFVFhoCOedkkpeUQURYSGc5Q6jrCHdXYMZauo42z3kfeIA13MhApl9ExmcGt8QEqN7dyO9W3RDCI1M78KItC4AhIYIQ3rEc0rPBACG9khgcENAxTGwIaDi6N89zv27SiI2whXGfZJi6J3kuhygd2I06d1cP49zB2T9fvq5/yBJdJ/OjXC/GzX+X6F7fGRDjabzcXxAiUh/EXlBRN4MdC2m7TbsKeX7L37DPa8vIyUhknfvOpPfTB1Bl+jwk6/cgiE94nnrromkJkQx7YWveXltFecOTeGZm0Y3XOjaJTqcl28ZR3JcJNc/t4S1uw/x9xtOb3gDbslPLxjCwJQ4np2/lYEpcdw1eUCzy96VNYD4yDA2F5ZxZ9aAJo/Mpp/dnyPVdagqN03ow7QJfahTSImPJGvIsddSXTnadeRwxoAkuscfPc01ZXgPAM4ckNRwce7I9K50iwknPFQY3z+xYT2AYb0SSIgKJzREGNPHdZHv6D7dEJGGo6jT0rsyKDWOiLAQBnaPIyo8tFFAxTO0hyugBqXEMzDF9fOA7nH0d4dM/+6x9EmKJTkugvio8Ibr1TISY0jvFk14qHB6b9e+xrqDqv7x6e4aeieENGwXXL8TV0C5/u1t+BsmIOxaKN856YW6IhID/BToraq3icggYIiqzvVg3ReBS4BCVR3RqH0K8DgQCjyvqo81tw1V3QrcagEVnLbvL+evn21kzordxEWE8cj3hjHtjL6Ehvju3SatazRv3nEG97y+jMrSYp5uFE71UhKieOXWcdz+Si43ju/NOUNTPdp2VHgof71mFD/79woevfLUE7bbWNeYCH4xZShv5Ozk+nFNfwYzODWeazMzCA8TMhJjSOsazZkDkzhnaCphocf+vXjBsB4M75VwwswWkwYlc2paF64de3QfoSHCDeN7s6+0siG0UhOiOH9YKhP6JzUsN3VULyqqahuC5ZKRvTh0pIY+STGICJee1os+7nA5d2gK8zfuY3haFyqqahiYEsekQcl0i4lgYEocZw5MpqKqlvdXFXBKzwQuGJ7aMEovISqcswd35zuDu1N6pIbYiDAmDkzmlFV7uGpMOp+vL2TahD6syi/h5jP7kXeggim9qvlsTwQ3jO/NvHWF9EmKITE2glEZXV37G5BEyeFq1u8ppXt8JPnFh9lVfJiUhEiOVNdysKIaBX5/2QgenrOGypOM4Mwa0v2Y68d85Uh1LbGRNgeCL8jJricRkX8BucD3VXWEiEQDi1V11Ek3LnI2UAa8XB9QIhIKbATOB/KBpcD1uMLq0eM2cYuqFrrXe1NVr/LkH5WZmal2R93A2n3wMH//fBNv5OQTHir8cGI/bj+7f5Oj83zJ+i54BWvf/eXTjTzR6JKEnF+f13DU1xn46I66uaqaeXy7JzE/QFWvFZHrAVT1sHj4gYGqzheRvsc1jwM2u4+MEJHZwFRVfRTX0VabiMh0YDpAamoq2dnZbd0UZWVlXq3fmZVUKnO3VvHFjhoUyMoI43v9w+katYcVS/ecdH1vWd8Fr2Dtu4KdrqH1AiiQveArkqMd/+mJz/iz3zwJqCr3UZMCiMgAwJtbR6YBjS+7zwfGN7ewiCQBvwdOF5FfuoPsBKo6E5gJriMobxI9WP+SC6SSimqenb+Ff3y1naraOq4ck8695w5q+IC8vVjfBa9g7bvt4dv498a1xESEUl5Vy6gxYxs+p+sM/NlvngTUI8BHQIaIvAacCdzsxT6bOvpq9jyjqhYBd3ixP+NH2/aX8/Li7byZk09pZQ3fO60XPzlvUMPoLmM6uvAw19FStDugDld1vguU/eWkAaWqn4hILjABV7jcp6rezNiZD2Q0epwO7G5mWeNAdXXKlxv3MWvRdr7cuI/wUOGiET25M2uARyPjjOlI6m+xUT+AxiaM9R1PRvHNU9VzgfebaGuLpcAgEekH7AKuA25o47ZMOzp0pJp/5+TzyuLtbC+qoHt8JPefN4gbxvUmJaH564mM6QzqL0i2CWN9p9mAEpEoIAZIFpFuHD01lwD08mTjIvI6kOXeRj7wsKq+ICJ3Ax/jGrn3oqquafs/wfjbpr2lzFq0nXeW7aKiqpYxfbrxk/MHc9GInkSEdZ4Pg41pSf11bzZhrO+0dAR1O3A/rjDK5WhAHQKe8mTjqnp9M+0fAB94XqZpbxVVNXy+vpB/fr2DRVuKiAgL4dLTevHDiX0bZgIwxhwV7Q6otQWHuMB9MbXxTrMBpaqPA4+LyD2q+vd2rMkESHmlK5Q+WFXAFxsKOVJdR68uUfz8wiFcNzaDpE50bYcxnhrXzzUzxtWZ6eTkFfO3zzZx/3mDA1xVx+DJIIm/i8gIYBgQ1aj9ZX8WZtpH6ZFqPl9fyPsrC/hy4z4q3bNUXz0mg4tO7cH4fkk+nfXBmI5mYEoc2x/7rvsGk6sCXU6H4skgiYdxfY40DNdpuYuAhYAFVJAqOVzNvHV7+WBVAfM37qeqto7UhEiuH9ebi0b0ILNvooWSMa3U1NyLxjueXAd1FXAasExVbxaRVOB5/5ZlfK2koppP1u7hw9V7WLDJdQuGnl2iuGlCHy4+tQeje3c74RYSxhjPWUD5nicBdVhV60SkRkQSgEKgv5/rMl4qqajm2x3FLN1+gJztxXy7o5iaOm24LcRFp/ZkVHpXCyVjfMTOOvieJwGVIyJdgedwjeYrA77xa1WmVVSVXQcPk7O9uOGW3hv2lqIKYSHCiLQu/Ois/lw0ogcj07u06d5LxpiTCwsRaux2Gz7TYkC5J4V9VFUPAjNE5CMgQVVXtkt1pkm1dcqGPaXk5B1g6fZicrYfoKDkCABxkWGM7tON757ak8y+iYzK6Ep0hJ16MKY9TDujD2/m5Ae6jA6jxYBSVRWRd4Ex7sfb26Moc6zDVbUs33mQnO0HyMkr5tu8YkorawDokRDF2H6JZPbpRmbfbgztkWCnGowJkKjwUA7bhbo+48kpviUiMlZVl/q9GkN5ZQ15RRVsLypn2Y5ilm4vZvWukobTBkNS47l0VC/G9k0ks2830rpG2yk7YxwiOjyUmjqlprbuhBtQmtbzJKAmA7eLSB5Qjvu2J6o60q+VdVCqSlF5FXlFFew4UO76XlRB3oEK8ooq2F929E4mEWEhjErvym1n92ds326M6Z1IlxjvbpFujPGf+vn4jtTUEWcB5TVPAuoiv1fRwdTU1lFQcoS8ogryDpS7AsgdQjuKyilvNJmkCPRMiKJ3UgznDk2hd1IMfZJi6JMYy+AecS3eYtwY4yz10x0drqolzm777jVPZpLIa49Cgs3hqlp2HKggr6jc/f1oAOUXHz5mJE9EaAgZidH0SYplfL9EVwAlxdA7MZb0btF2/YQxHUSkTRjrUxbxzVBVDlZUu0+9uU7FNT4tV1h67E2FE6LC6JMUy/C0Llx8as+GAOqTFEOPhCi73siYTsBmNPetDhlQW/aVcfnTX7V5/aLiwxRnf0LpkZpj2lMTIumTGMvZg7vTJzHGfToulr5JMXSNifC2bGNMkKs/xXflM4tY+ciFAa4m+HXIgAoR8er8b12UkDUijd6JrgDqkxRDRrcYu57IGNOi+kESh47749a0jSeTxZYCx18aXQLkAD9V1a3+KMwb/ZJjeeXW8W1ePzs7m6ysET6syBjTGUTYyD2f8uQw4y/AbuCfuIaYXwf0ADYAL+Ka6dwYY4zxKU/ifoqqPquqpap6SFVnAher6r+Abn6uzxhjgpKqzcnnLU8Cqk5ErhGREPfXNY2esx4wxhi3Hl0a7ulKZU1dACvpGDwJqBuBabhus7HX/fNNIhIN3O3H2owxJqj0SYrlkpE9ARtq7gueXKi7FfheM08v9G05JxKRU4D7gGRgnqo+4+99GmNMW00ckMzclQUcqbYjKG95MoqvO3Ab0Lfx8qp6iwfrvghcAhSq6ohG7VOAx4FQ4HlVfay5bajqOuAOEQnBdU8qY4xxrIb5+OwIymuejOJ7D1gAfAa09jc+C3gSeLm+QURCgaeA84F8YKmIzMEVVo8et/4tqlooIpcCD7q3ZYwxjtUwm0SNBZS35GQjTURkuaqOavMORPoCc+uPoETkDOARVb3Q/fiXAKp6fDg1ta33VfW7zTw3HZgOkJqaOmb27NltLZmysjLi4uLavL4JHOu74NVR+m55YQ1/+7aShyZEMaBrx7+43xf9Nnny5FxVzTy+3ZMjqLkicrGqfuBVBUelATsbPc4Hmr2qVkSygCuASKDZGtzD32cCZGZmalZWVpsLdF2o2/b1TeBY3wWvjtJ3EZv3w7dfM+zUUZwxICnQ5fidP/vNk4C6D/iViFQC1Ry9H1RCG/fZ1KypzR7GqWo2kN3GfRljTLuKtFN8PuPJKL54H+8zH8ho9Dgd10wVxhgT9OoHSVTaIAmvNRtQIjJUVdeLyOimnlfVb9u4z6XAIBHpB+zCNXXSDW3cljHGOMrRW27YMHNvtXQE9QCuQQd/buI5Bc452cZF5HVcc/Uli0g+8LCqviAidwMf4xq596Kqrmlt4cYY40TRdk8on2k2oFR1uvv75LZuXFWvb6b9A1oY8GCMMcHKblroOx7dNElEJnLihbovN7uCMcZ0Ug0X6tpcfF7zZCaJV4ABwHKOXqirNLr41hhjjEtUmB1B+YonR1CZwDC1ueONMeakQkKEiNAQGyThA57MZr4a1w0KjTHGeCAyPMSOoHzAkyOoZGCtiHwDVNY3quqlfqvKGGOCWFR4KJV2oa7XPAmoR/xdhDHGdCRR4XaKzxdaDCj3zOMPqep57VSPMcYEvaiwUDvF5wMtfgalqrVAhYh0aad6jDEm6EWFW0D5gien+I4Aq0TkU6C8vlFV7/VbVcYYE8TsFJ9veBJQ77u/jDHGeCAqPJSyyppAlxH0PJnN/KX2KMQYYzqKyLBQ9pVWnnxB0yJPZpIYhOtW7MOAqPp2Ve3vx7qMMSZoRYWHUGlTHXnNkwt1/wE8A9QAk3FNcfSKP4syxphgZoMkfMOTgIpW1XmAqGqeqj6CB7faMMaYziraAsonPBrFJyIhwCb3fZx2ASn+LcsYY4KXjeLzDU+OoO4HYoB7gTHATcAP/FmUMcYEs6jwUI7U1GJzbHvHk1F8SwFERFX1Zv+XZIwxwS0qPBRVqKqtI9J9+w3Teic9ghKRM0RkLbDO/fg0EXna75UZY0yQigxz37TQTvN5xZNTfH8DLgSKAFR1BXC2P4syxphgVn/b90obKOEVTwIKVd15XJP91o0xphn1AWVHUN7xJKB2ishEQEUkQkR+hvt0X3sQkSwRWSAiM0Qkq732a4wxbRUV7j7FZ/eE8oonAXUH8GMgDcgHRgF3ebJxEXlRRApFZPVx7VNEZIOIbBaRB0+yGQXKcM1ike/Jfo0xJpCiwuqPoCygvOHJKL79wI2N20TkflyfTZ3MLOBJXLNP1K8bCjwFnI8rcJaKyBwgFNeUSo3dAixQ1S9FJBX4y/G1GGOM09gpPt/w5ELdpjyABwGlqvNFpO9xzeOAzaq6FUBEZgNTVfVR4JIWNlcMRDb3pIhMB6YDpKamkp2dfbLymlVWVubV+iZwrO+CV0fqu03FriOnb3K/pSKvrW+zwcGf/dbW35x4sc80oPGgi3xgfLM7ErkC1yjCrriOxpqkqjOBmQCZmZmalZXV5gKzs7PxZn0TONZ3wasj9V3yrhL4eiGDTxlB1vAegS7Hr/zZb20NKG8uj24q3Jrdnqq+Dbztxf6MMaZdHR0kYaf4vNFsQIlIKU0HhwDRXuwzH8ho9Dgd2O3F9owxxlEibZCETzQbUKoa76d9LgUGiUg/XBPPXgfc4Kd9GWNMuzs6SMICyhseXajbViLyOrAYGCIi+SJyq6rWAHcDH+O6nuoNVV3jzzqMMaY9RUdYQPmCX4eXqOr1zbR/AHzgz30bY0ygRNlcfD7h1yMoY4zpjMJCQwgLETuC8pIFlDHG+IHrtu92BOUNCyhjjPGDqPAQm4vPSxZQxhjjB5FhoXaKz0sWUMYY4wdR4SFU2ik+r1hAGWOMH7g+g7IjKG9YQBljjB9EhYfaZ1Be6tjT7BpjTIBEhYew91Al89btDXQpfrW1uJYsP23bAsoYY/wgOS6SrzYXcetLOYEuxa8Gdg3hNj9t2wLKGGP84NErTuXWSf0CXYbfrV7+rd+2bQFljDF+EBMRxsj0roEuw+8ObPbfUAYbJGGMMcaRLKCMMcY4kqh6c3NcZxKRfUBeE091AUo8aEsG9vuhtNZqqrb23lZr1jvZsm19vjXtTug7X/abN9vzdD1PlmtpmdY+19zyHa3vnPCaO9kyTnnN9VHV7ie0qmqn+QJmetiWE+ham6utvbfVmvVOtmxbn29NuxP6zpf91h5958lyLS3T2uda6M8O1XdOeM21te+c8prrbKf4/uNhm1P4sra2bqs1651s2bY+39r2QPN1Xf7uO0+Wa2mZ1j7n1H6DjveaO9kyjn7NdchTfN4SkRxVzQx0Hab1rO+Cl/VdcPJnv3W2IyhPzQx0AabNrO+Cl/VdcPJbv9kRlDHGGEeyIyhjjDGOZAFljDHGkSygjDHGOJIFlDHGGEeygPKAiMSKyEsi8pyI3BjoeoznRKS/iLwgIm8GuhbjORG5zP16e09ELgh0PcZzInKKiMwQkTdF5E5vttVpA0pEXhSRQhFZfVz7FBHZICKbReRBd/MVwJuqehtwabsXa47Rmr5T1a2qemtgKjWNtbLf3nW/3n4IXBuAck0jrey7dap6B3AN4NX1UZ02oIBZwJTGDSISCjwFXAQMA64XkWFAOrDTvZjdwznwZuF53xnnmEXr++3X7udNYM2iFX0nIpcCC4F53uy00waUqs4HDhzXPA7Y7P6ruwqYDUwF8nGFFHTi35lTtLLvjEO0pt/E5f+AD1XVf3fEMx5p7WtOVeeo6kTAq49E7M32WGkcPVICVzClAW8DV4rIMzh7HrHOrMm+E5EkEZkBnC4ivwxMaaYFzb3m7gHOA64SkTsCUZg5qeZec1ki8oSIPAt84M0O7I66x5Im2lRVy4Gb27sY0yrN9V0RYG9wztVcvz0BPNHexZhWaa7vsoFsX+zAjqCOlQ9kNHqcDuwOUC2mdazvgpP1W/Dye99ZQB1rKTBIRPqJSARwHTAnwDUZz1jfBSfrt+Dl977rtAElIq8Di4EhIpIvIreqag1wN/AxsA54Q1XXBLJOcyLru+Bk/Ra8AtV3Npu5McYYR+q0R1DGGGOczQLKGGOMI1lAGWOMcSQLKGOMMY5kAWWMMcaRLKCMMcY4kgWUMX4iImXu731F5AYfb/tXxz1e5MvtG+MEFlDG+F9foFUB5b6VQUuOCSj3zNHGdCgWUMb432PAWSKyXER+IiKhIvJHEVkqIitF5HYA9yzQX4jIP4FV7rZ3RSRXRNaIyHR322NAtHt7r7nb6o/WxL3t1SKySkSubbTtbPddTteLyGsi0tRkn8Y4hs1mboz/PQj8TFUvAXAHTYmqjhWRSOArEfnEvew4YISqbnM/vkVVD4hINLBURN5S1QdF5G5VHdXEvq4ARgGnAcnudea7nzsdGI5rQs+vgDNx3VTOGEeyIyhj2t8FwPdFZDnwNZAEDHI/902jcAK4V0RWAEtwzRw9iJZNAl5X1VpV3Qt8CYxttO18Va0DluM69WiMY9kRlDHtT4B7VPXjYxpFsoDy4x6fB5yhqhUikg1EebDt5lQ2+rkWe/0bh7MjKGP8rxSIb/T4Y+BOEQkHEJHBIhLbxHpdgGJ3OA0FJjR6rrp+/ePMB651f87VHTgb+MYn/wpj2pn9BWWM/60Eatyn6mYBj+M6vfate6DCPuCyJtb7CLhDRFYCG3Cd5qs3E1gpIt+q6o2N2t8BzgBWAAr8QlX3uAPOmKBit9swxhjjSHaKzxhjjCNZQBljjHEkCyhjjDGOZAFljDHGkSygjDHGOJIFlDHGGEeygDLGGONI/x8lR9eLeFUZzgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "winit = np.random.randn(p)\n",
        "nit = 1000\n",
        "w, f0, hist = grad_opt_adapt(feval, winit, nit=nit)\n",
        "\n",
        "t = np.arange(nit)\n",
        "plt.subplot(2,1,1)\n",
        "plt.semilogx(t, hist['f'])\n",
        "plt.grid()\n",
        "plt.ylabel('Loss')\n",
        "\n",
        "plt.subplot(2,1,2)\n",
        "plt.loglog(t, hist['lr'])\n",
        "plt.grid()\n",
        "plt.ylabel('Learning rate')\n",
        "plt.xlabel('Iteration')\n",
        "\n",
        "plt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CQXpHvbnzVqj"
      },
      "source": [
        "Finally we measure the accuracy."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jdnqXJ41zVqk",
        "outputId": "710642a0-5953-4c1c-d02a-a88c32deac3a"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Test accuracy = 0.960976\n"
          ]
        }
      ],
      "source": [
        "yhat = predict(Xts1,w)\n",
        "acc = np.mean(yhat == yts)\n",
        "print(\"Test accuracy = %f\" % acc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DWeQ9iB0zVqk"
      },
      "source": [
        "## In-Class Exercise\n",
        "\n",
        "Try to a build a simple optimizer to minimize:\n",
        "\n",
        "    f(w) = a[0] + a[1]*w + a[2]*w^2 + ... + a[d]*w^d\n",
        "    \n",
        "for the coefficients `a = [0,0.5,-2,0,1]`.\n",
        "\n",
        "* Plot the function f(w)\n",
        "* Can you see where the minima is?\n",
        "* Write a function that outputs f(w) and its gradient.\n",
        "* Run the optimizer on the function to see if it finds the minima.\n",
        "* Print the funciton value and number of iterations.\n",
        "* Bonus:  Instead of writing the function for a specific coefficient vector `a`, create a class that works for an arbitrary vector `a`.\n",
        "\n",
        "You may wish to use the `poly.polyval(w,a)` method to evaluate the polynomial."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "pK4YS1tGzVqk"
      },
      "outputs": [],
      "source": [
        "import numpy.polynomial.polynomial as poly\n",
        "\n",
        "a = np.array([0,0.5,-2,0,1])\n",
        "def feval(w,a):\n",
        "  f =  poly.polyval(w,a)\n",
        "  fgrad = 0\n",
        "  #fgrad = poly.polyder(f)\n",
        "  return f, fgrad\n",
        "\n",
        "w = np.linspace(-2,2,100)\n",
        "f = poly.polyval(w,a)\n",
        "plt.plot(w,f)\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9LKVGREdzVqk"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.4"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}