{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# Install the necessary dependencies\n",
    "\n",
    "import sys\n",
    "import os\n",
    "!{sys.executable} -m pip install --quiet matplotlib numpy pandas ipython jupyterlab_myst scikit-learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sklearn\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import train_test_split\n",
    "from IPython.display import HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "---\n",
    "license:\n",
    "    code: MIT\n",
    "    content: CC-BY-4.0\n",
    "github: https://github.com/ocademy-ai/machine-learning\n",
    "venue: By Ocademy\n",
    "open_access: true\n",
    "bibliography:\n",
    "  - https://raw.githubusercontent.com/ocademy-ai/machine-learning/main/open-machine-learning-jupyter-book/references.bib\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective of this session\n",
    "\n",
    "We have already learnt how to use Linear Regression and Logistic Regression models.\n",
    "\n",
    "The code might seem quite easy and intuitive for you. And you might naturally ask:\n",
    "- What's behind the ```.fit()``` function?\n",
    "- Why sometimes it takes quite a bit for this ```.fit()``` function to finish running?\n",
    "\n",
    "In this session, you will learn that the ```.fit()``` is the training of ML models, \n",
    "i.e. tuning of parameters for ML models. And the technique behind is called \"Gradient Descent\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Video\n",
    "\n",
    "The corresponding video (in Chinese) for this notebook is [👉 available here on Bilibili](https://www.bilibili.com/video/BV1SY4y1G7o9/).\n",
    "You can (and should) watch the video before diving into the details of gradient descent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<p style=\"text-align: center;\">\n",
       "<iframe src=\"https://player.bilibili.com/player.html?aid=642485873&cid=764796592&page=1&high_quality=1&danmaku=0\" width=\"105%\" height=\"700px;\" style=\"border:none;\"></iframe>\n",
       "video. <a href=\"https://player.bilibili.com/player.html?aid=642485873&cid=764796592&page=1&high_quality=1&danmaku=0\">[source]</a>\n",
       "</p>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "\n",
    "display(\n",
    "    HTML(\n",
    "        \"\"\"\n",
    "<p style=\"text-align: center;\">\n",
    "<iframe src=\"https://player.bilibili.com/player.html?aid=642485873&cid=764796592&page=1&high_quality=1&danmaku=0\" width=\"105%\" height=\"700px;\" style=\"border:none;\"></iframe>\n",
    "video. <a href=\"https://player.bilibili.com/player.html?aid=642485873&cid=764796592&page=1&high_quality=1&danmaku=0\">[source]</a>\n",
    "</p>\n",
    "\"\"\"\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's be playful ... to gain some intuition\n",
    "\n",
    "- [Tensorflow Playground](https://playground.tensorflow.org/#activation=sigmoid&batchSize=10&dataset=circle&regDataset=reg-plane&learningRate=0.00001&regularizationRate=0&noise=0&networkShape=&seed=0.71864&showTestData=false&discretize=false&percTrainData=50&x=true&y=true&xTimesY=true&xSquared=true&ySquared=true&cosX=false&sinX=false&cosY=false&sinY=false&collectStats=false&problem=classification&initZero=false&hideText=false)\n",
    "- [Gradient Descent Visualization](https://github.com/lilipads/gradient_descent_viz)\n",
    "- [Optimization Algorithms Visualization](https://bl.ocks.org/EmilienDupont/aaf429be5705b219aaaf8d691e27ca87)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some mathematics ... to gain more insight"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Abstract\n",
    "\n",
    "The idea behind gradient descent is simple - by gradually tuning parameters, such as slope ($m$) and the intercept ($b$) in our regression function $y = mx + b$, we minimize cost. \n",
    "By cost, we usually mean some kind of a function that tells us how far off our model predicted result. For regression problems we often use `mean squared error` (MSE) cost function. If we use gradient descent for the classification problem, we will have a different set of parameters to tune.\n",
    "\n",
    "$$ MSE = \\frac{1}{n}\\sum_{i=1}^{n} (y_i - \\hat{y_i})^2 \\quad \\textrm{where} \\quad \\hat{y_i} = mx_i + b $$\n",
    "\n",
    "Now we have to figure out how to tweak parameters $m$ and $b$ to reduce MSE.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Partial derivatives\n",
    "\n",
    "We use partial derivatives to find how each individual parameter affects MSE, so that's where word _partial_ comes from. In simple words, we take the derivative with respect to $m$ and $b$ **separately**. Take a look at the formula below. It looks almost exactly the same as MSE, but this time we added f(m, b) to it. It essentially changes nothing, except now we can plug $m$ and $b$ numbers into it and calculate the result.\n",
    "\n",
    "$$𝑓(𝑚,𝑏)= \\frac{1}{n}\\sum_{i=1}^{n}(y_i - (mx_i+b))^2$$\n",
    "\n",
    "This formula (or better say function) is better representation for further calculations of partial derivatives. We can ignore sum for now and what comes before that and focus only on $y - (mx + b)^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Partial derivative with respect to $m$\n",
    "\n",
    "With respect to $m$ means we derive parameter $m$ and basically ignore what is going on with $b$, or we can say its 0. To derive with respect to $m$ we will use chain rule.\n",
    "\n",
    "$$ [f(g(x))]' = f'(g(x)) * g(x)' \\: - \\textrm{chain rule}$$\n",
    "\n",
    "Chain rule applies when one function sits inside of another. If you're new to this, you'd be surprised that $()^2$ is outside function, and $y-(\\boldsymbol{m}x+b)$ sits inside it. So, the chain rule says that we should take a derivative of outside function, keep inside function unchanged and then multiply by derivative of the inside function. Lets write these steps down:\n",
    "\n",
    "$$ (y - (mx + b))^2 $$\n",
    "\n",
    "1. Derivative of $()^2$ is $2()$, same as $x^2$ becomes $2x$\n",
    "2. We do nothing with $y - (mx + b)$, so it stays the same\n",
    "3. Derivative of $y - (mx + b)$ with respect to **_m_** is $(0 - (x + 0))$ or $-x$, because **_y_** and **_b_** are constants, they become 0, and derivative of **_mx_** is **_x_**\n",
    " \n",
    "Multiply all parts we get following: $2 * (y - (mx+b)) * -x$. \n",
    "\n",
    "Looks nicer if we move -x to the left: $-2x *(y-(mx+b))$. There we have it. The final version of our derivative is the following:\n",
    "\n",
    "$$\\frac{\\partial f}{\\partial m} = \\frac{1}{n}\\sum_{i=1}^{n}-2x_i(y_i - (mx_i+b))$$\n",
    "\n",
    "Here, $\\frac{df}{dm}$ means we find partial derivative of function f (we mentioned it earlier) with respect to m. We plug our derivative to the summation and we're done.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Partial derivative with respect to $b$\n",
    "\n",
    "Same rules apply to the derivative with respect to b.\n",
    "\n",
    "1. $()^2$ becomes $2()$, same as $x^2$ becomes $2x$\n",
    "2. $y - (mx + b)$ stays the same\n",
    "3. $y - (mx + b)$ becomes $(0 - (0 + 1))$ or $-1$, because **_y_** and **_mx_** are constants, they become 0, and derivative of **_b_** is 1\n",
    "\n",
    "Multiply all the parts together and we get $-2(y-(mx+b))$\n",
    "\n",
    "$$\\frac{\\partial f}{\\partial b} = \\frac{1}{n}\\sum_{i=1}^{n}-2(y_i - (mx_i+b))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final function\n",
    "\n",
    "Few details we should discuss before jumping into code:\n",
    "\n",
    "1. Gradient descent is an iterative process and with each iteration ($epoch$) we slightly minimizing MSE, so each time we use our derived functions to update parameters $m$ and $b$.\n",
    "2. Because it's iterative, we should choose how many iterations we take, or make algorithm stop when we approach minima of MSE. In other words when algorithm is no longer improving MSE, we know it reached minimum.\n",
    "3. Gradient descent has an additional parameter learning rate ($lr$), which helps control how fast or slow algorithm going towards minima of MSE\n",
    "\n",
    "That's about it. So you can already understand that Gradient Descent for the most part is just process of taking derivatives and using them over and over to minimize function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time to code!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear regression With gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LinearRegression:\n",
    "    def __init__(self, learning_rate=0.0003, n_iters=3000):\n",
    "        self.lr = learning_rate\n",
    "        self.n_iters = n_iters\n",
    "        self.weights = None\n",
    "        self.bias = None\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        n_samples, n_features = X.shape\n",
    "\n",
    "        # init parameters\n",
    "        self.weights = np.zeros(n_features)\n",
    "        self.bias = 0\n",
    "\n",
    "        # gradient descent\n",
    "        for _ in range(self.n_iters):\n",
    "            # approximate y with linear combination of weights and x, plus bias\n",
    "            y_predicted = np.dot(X, self.weights) + self.bias\n",
    "\n",
    "            # compute gradients\n",
    "            dw = (1 / n_samples) * np.dot(X.T, (y_predicted - y))\n",
    "            db = (1 / n_samples) * np.sum(y_predicted - y)\n",
    "            # update parameters\n",
    "            self.weights -= self.lr * dw\n",
    "            self.bias -= self.lr * db\n",
    "\n",
    "    def predict(self, X):\n",
    "        y_predicted = np.dot(X, self.weights) + self.bias\n",
    "        return y_predicted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'lr': 0.0003, 'n_iters': 3000, 'weights': array([0.36114314, 0.15172482, 0.01138062, 0.07103796, 0.10143793,\n",
      "       0.14812986, 0.09146885, 0.00270041]), 'bias': 0.014542612245156487}\n",
      "0    -1.470137\n",
      "1    -1.226722\n",
      "2    -1.633534\n",
      "3    -1.145394\n",
      "4    -1.385705\n",
      "        ...   \n",
      "92    0.985388\n",
      "93    1.125408\n",
      "94    1.936285\n",
      "95    1.776223\n",
      "96    1.680470\n",
      "Name: lpsa, Length: 97, dtype: float64\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHn0lEQVR4nO3de1xUdf4/8NcMV0EYxRsgqKh5Qbxx89Jly9TVWtMuW5oparlp2mXddltzt+zK9qu22q8rqYVillZrWm2tZVle8gKoKF7yiooKIqDDTQaYOb8/bEhwhjln5sycc2Zez8eDxyNhmHkzkud1Pp/35/PRCYIggIiIiEgGeqULICIiIu/BYEFERESyYbAgIiIi2TBYEBERkWwYLIiIiEg2DBZEREQkGwYLIiIikg2DBREREcnG39MvaLFYcP78eYSFhUGn03n65YmIiMgJgiCgsrIS0dHR0Ovtj0t4PFicP38esbGxnn5ZIiIikkFhYSFiYmLsft3jwSIsLAzA1cLCw8M9/fJERETkhIqKCsTGxjZex+3xeLCwTn+Eh4czWBAREWmMozYGNm8SERGRbBgsiIiISDYMFkRERCQbBgsiIiKSDYMFERERyYbBgoiIiGTDYEFERESyYbAgIiIi2TBYEBERkWwkBYuFCxdCp9M1+YiMjHRXbURERKQxkrf07tevH7777rvGP/v5+claEBEREWmX5GDh7+/PUQoiIiKySXKPxbFjxxAdHY24uDhMnDgRJ0+ebPHxJpMJFRUVTT6IiIi0IK/wMh56bxcqauuVLkUzJAWLIUOGYOXKlfjmm2+wbNkyFBcXY/jw4SgrK7P7Penp6TAYDI0fsbGxLhdNRETkbj/8XIJJS3di2/FS/PPbo0qXoxk6QRAEZ7+5uroaPXr0wF/+8hfMmzfP5mNMJhNMJlPjn63nuRuNRh6bTkREqvRJbiHmf5YPs0XAb3p1wOLJiQgNktw94FUqKipgMBgcXr9depdCQ0PRv39/HDt2zO5jgoKCEBQU5MrLEBEReYQgCPj3D8fxxi8jFPckdsZr9w5AgB93ZxDLpXfKZDLh8OHDiIqKkqseIiIiRZgtAp77/GBjqHjs1h548/cDGSokkjRi8fTTT2PcuHHo0qULSkpK8PLLL6OiogJpaWnuqo+IiMjtauvNeGpNHjYcLIZOBywc1w9pw7spXZYmSQoWZ8+exaRJk1BaWooOHTpg6NCh2LlzJ7p27equ+oiIiNzKWFOPmStzkX2qHIF+erw9cRDu6M+ReGdJChZr1qxxVx1EREQed/7yFUxbno2jF6oQFuyPZVOTMbR7O6XL0jTfbnElIiKfdfRCJdIys1FkrEWn8CBkzUhFn0iuVnQVgwUREfmc7IJyPJKVg4raBvTs2BpZM1LRuU0rpcvyCgwWRETkUzYcKMITa/JQ12BBcte2eC8tGW1CApUuy2swWBARkc/4YMcpPPfFQQgCMDq+E/41aTCCA3iYppwYLIiIyOsJgoA3vj2Cf/9wAgAweUgXvDg+AX56ncKVeR8GCyIi8mr1Zgue/Swfn+4+CwD406hemDuiJ3Q6hgp3YLAgIiKvVVPXgDkf7sEPRy7CT6/DKxMSMDG1i9JleTUGCyIi8kplVSbMyMrFvsLLCA7Q498PJuL2vp2ULsvrMVgQEZHXKSyvwdTMbBSUVqNtSADen5aCxC5tlS7LJzBYEBGRVzlwzohpy3NQWmVC5zatsPLhVPTo0FrpsnwGgwUREXmNbcdK8egHuaiuM6NvVDiypqegY3iw0mX5FAYLIiLyCp/nncPTn+5DvVnA8B7t8O6UJIQHByhdls9hsCAiIs1btuUkXvn6MABg3MBovPH7AQjy58ZXSmCwICIizbJYBLzy9WG8v60AAPDwTXFYcEdf6LnxlWIYLIiISJNMDWb8+dP9+GLfeQDAgjv6YuYt3RWuihgsiIhIcypr6/HoB7ux/UQZAvx0eP2+gZgwuLPSZREYLIiISGNKKmqRtjwHh4sqEBrohyVTknHTDe2VLot+wWBBRESaceJiFdIys3H20hW0bx2EFdNTkNDZoHRZdA0GCyIi0oQ9Zy7h4RU5uFRTj7j2ocianoou7UKULouaYbAgIiLV+/7wBcz5aA9q6y0YGGNA5rQUtGsdpHRZZAODBRERqdrHOWfw7LoDMFsE3Na7A/49OREhgbx8qRX/ZoiISJUEQcCiTcfx5sajAID7kmKQfk9/BPjpFa6MWsJgQUREqmO2CHju8wP4cNcZAMDc23riT6N7Qafjxldqx2BBRESqUltvxhOr9+LbQxeg0wEv3tUPU4Z1U7osEonBgoiIVONyTR0eycpF7ulLCPTX418TB2FMQpTSZZEEDBZERKQK5y5fQVpmNo6XVCE82B/vpaUgNS5C6bJIIgYLIiJS3M/FFZiWmYPiilpEGYKRNSMVvTqFKV0WOYHBgoiIFLXzZBlmrsxFZW0DbujYGlkzUhHdppXSZZGTGCyIiEgxX+cX4ak1eagzW5DaLQLLpibDEBKgdFnkAgYLIiJSRNb2U1j45UEIAvDbfp3wzsTBCA7wU7oschGDBREReZQgCHj9myNY/OMJAMCUoV2x8K5+8NNzjwpvwGBBREQeU2+24K9r87F2z1kAwNOje2HObT258ZUXYbAgIiKPqDY14LEP92Dz0Yvw0+uQfk9/3J8cq3RZJDMGCyIicrvSKhNmrMjB/rNGtArww78nD8aIPp2ULovcgMGCiIjc6nRZNdIys3GqrAYRoYHInJaCQbFtlC6L3ITBgoiI3Cb/rBHTV2SjtKoOMW1bYeWMVHTv0FrpssiNGCyIiMgtthy9iNmrdqO6zox+0eFYPj0FHcOClS6L3IzBgoiIZLdu71n8+dP9aLAIuKlne2Q8lIiwYG585QsYLIiISDaCIGDplpNI/9/PAIDxg6Lx+n0DEeivV7gy8hQGCyIikoXFIuClrw5h+U+nAAAzb47D/LF9oefGVz6FwYKIXGa2CMguKEdJZS06hgUjNS6Cuyj6GFODGfM+2Yev9hcBAP52Z188cnN3hasiJTBYEJFLNhwowgtfHkKRsbbxc1GGYDw/Lh5jEqIUrIw8paK2Hn9YmYudJ8sR4KfDG78fiPGDOitdFimEk15E5LQNB4owe9WeJqECAIqNtZi9ag82HChSqDLylAsVtbj/3R3YebIcrYP8sWJ6KkOFj2OwICKnmC0CXvjyEAQbX7N+7oUvD8FssfUI8gbHSypxz+Lt+Lm4Eh3CgvDxo0NxY8/2SpdFCmOwICKnZBeUXzdScS0BQJGxFtkF5Z4rijxm9+lLuO/dHTh3+Qq6tw/FZ7OHo1+0QemySAXYY0FETimptB8qnHkcacfGQxfw+Oo9qK23YFBsG2ROS0FEaKDSZZFKMFgQkVPE7qDInRa9y+rsM1iwLh8WARjRpyMWPTgYIYG8lNCvOBVCRE5JjYtAlCEY9haV6nB1dUhqXIQnyyI3EQQBb393FPM/uxoq7k+OwdIpSQwVdB0GCyJyip9eh+fHxQPAdeHC+ufnx8VzPwsv0GC24Nl1+Xj7u2MAgCdG9MRr9w6Avx8vIXQ9/lYQkdPGJEQh46FERBqaTndEGoKR8VAi97HwAlfqzJi1ag9WZxdCrwNenpCAeaN7Q6djYCTbOIZFRC4ZkxCFUfGR3HnTC12qrsPDWTnYc+Yygvz1+Nekwfhtv0ilyyKVY7AgIpf56XUY1qOd0mWQjM5eqkFaZjZOXKxGeLA/3p+WgpRu7JchxxgsiIioicNFFUjLzEZJpQnRhmBkzUjFDZ3ClC6LNILBgoiIGu04UYY/rMxFpakBvTuFYcWMFEQZWildFmkIgwUREQEA/rv/POZ9vA91ZgtS4yKwbGoyDK0ClC6LNIbBgoiIsPynArz430MQBGBsQiTeemAQggP8lC6LNIjBgojIh1ksAl775mcs2XwSAJA2rCueG9ePq3rIaS7tY5Geng6dToennnpKpnKIiMhT6s0WPP3pvsZQ8ZcxvbHwLoYKco3TIxY5OTlYunQpBgwYIGc9RETkAVWmBsxetRtbj5XCT6/Da/cOwH1JMUqXRV7AqRGLqqoqTJ48GcuWLUPbtm3lromIiNzoYqUJk5buxNZjpQgJ9MP7ackMFSQbp4LFnDlzcOedd2LkyJEOH2symVBRUdHkg4iIlHGqtBr3ZmxH/jkj2oUGYvXMobi1d0elyyIvInkqZM2aNdizZw9ycnJEPT49PR0vvPCC5MKIiEhe+89exvTlOSirrkOXiBCsnJGKbu1DlS6LvIykEYvCwkI8+eSTWLVqFYKDgx1/A4D58+fDaDQ2fhQWFjpVKBEROe/HIyWYuHQnyqrrkNA5HGtnD2eo8DJmi4AdJ8rwed457DhRBrNFUKQOnSAIol95/fr1uPvuu+Hn9+vaZrPZDJ1OB71eD5PJ1ORrtlRUVMBgMMBoNCI8PNz5yomISJS1u8/imbX70WARcPMN7ZHxUBJaB3G3AW+y4UARXvjyEIqMtY2fizIE4/lx8bKdMiz2+i0pWFRWVuL06dNNPjd9+nT06dMHzzzzDBISEmQrjIiIXCMIAjI2n8D/23AEAHD34M547d4BCPR3aacBUpkNB4owe9UeNL+YWxcNZzyUKEu4EHv9lhRZw8LCrgsPoaGhaNeunahQQUREnmG2CHjxy4PI2nH1ZvDR33THM7/tAz33qPAqZouAF748dF2oAAABV8PFC18ewqj4SI/tT8LYSkTkZWrrzXh89R5k7TgNnQ547nfxmD+2L0OFF8ouKG8y/dGcAKDIWIvsgnKP1eTyJNuPP/4oQxlERCQH45V6/GFlLnYVlCPQT49/PjAQvxsQrXRZ5CYllfZDhTOPkwO7d4iIvESxsRZpmdk4cqESYUH+WDI1CcN7tFe6LHKjjmHiVmiKfZwcGCyIiLzA8ZJKTH0/G+eNtegYFoQV01MRH80GeW+XGheBKEMwio21NvssdAAiDcFIjYvwWE3ssSAi0rjcU+W4N2MHzhtr0b1DKD57bDhDhY/w0+vw/Lh4AL+uArGy/vn5cfEePViOwYKISMO+OViMye/tgvFKPRK7tMHaWcMR0zZE6bLIg8YkRCHjoUREGppOd0QagmVbaioFp0KIiDTqw12n8ff1B2ARgJF9O+L/JiWiVWDLmxSSdxqTEIVR8ZHILihHSWUtOoZdnf7w5EiFFYMFEZHGCIKAtzYexb82HQcATEqNxUvjE+Dvx0FoX+an12FYj3ZKl8FgQUSkJQ1mCxasO4CPc6+eu/Tk7TfgqZE3QKfjHhWkDgwWREQacaXOjLkf7cH3P5dArwNentAfDw7ponRZRE0wWBARaUB5dR0ezsrB3jOXEeSvx6IHEzEqvpPSZZEHmS2CKnooHGGwICJSucLyGqRlZuNkaTXahATg/bRkJHX13L4EpDxPnF4qF3b6EBGp2MHzRtyTsR0nS6vRuU0r/GfWMIYKH2M9vbT5mSDFxlrMXrUHGw4UKVSZbRyxICISSY6haCnPsf14Kf7wwW5UmRrQJzIMWTNS0Snc+a2ZtTKUTr9S4+mljjBYEBGJIMdQtJTn+GLfefzpkzzUmwUM7R6BJVOSYWgVoGj95HlSTi9Vw1JTgFMhREQOyTEU3dJzzFq1B+98dxSf553DjhNlWLblJJ5YvRf1ZgF3DohC1oxUl0OFlobS6VdqPL3UEY5YEBG1QI6haEfPAQBvfXfsuq9NG94Nz/0uHnoXhri1OJROv1Lj6aWOcMSCiKgFUoainX0Oe4bEtXUpVIh5bTH1k3Ksp5fa+y3Q4eqUlidPL3WEwYLIS5ktAnacKGscXjdbbN2zqpda6pdjKNqZYWodgBf/e9jln1uLQ+n0KzWeXuoIp0KIvJDWG/XUVL8cQ9HODFPL1ZSnxaF0KXxhpYv19NLm/09EqvT/aQYLIi9jbdRrfp9rbdRT4hhlKdRWv3UouthYa7NPQYer/8C3NBTt6Dla4upIghz1q5WaAqi7qen0Ukc4FULkRcQ0Cb7w5SHVTouosX45hqJbeg5HXB1J0OJQuhi+uNLFenrp+EGdMaxHO9X+nTFYEHkRrTfqqbV+61B0pKHpRT7SEHzdCIq93hB7z2GPnE15UurXAjUGUPoVp0KIvIjWG/XUXL+YoWhHQ/NjEqJgvNKA+Z/th0UAurYLwemyGuiAJhdJd4wkaGko3REtbhrlSxgsiLyI1hv11F6/dSjaFke9IYsnD8aJi9V449ujAIB7EjvjtXsH4PvDFzzWlNdS/Vqi5gBKDBZEXkXrjXparV/M0PyfPt2PmjozAGD2rT3wl9/2hk6n86qRBE9RewD1deyxIPIiWm/Uc9TkKACYmBLr0ZrEELMBVk2d+eoul3f1wzNj+kCn+/Un1EpTnlpocdMoX8JgQeRltN6o56jJ8a3vjuGm1za5vetfygZdYofcp93YDWnDu8lUoe/SeoD2djpBEDzaNltRUQGDwQCj0Yjw8HBPvjSRT9H6xkFmi4BFm47jre+OXvc160/hrqAkdX+EHSfKMGnZTofPu3rmUK/ocVALX9rHQg3EXr8ZLIhIlcwWATe9tsnuFIO132LbMyNkDUz2mjBbCjPWWlvaACvKDbWS9gO0loi9fnMqhIhUSYk9LZzdH+HaoXlbdPDtoXl3nvvC/hT14aoQIlIlJZYUurI/wpiEKPzhljgs2VLQ5PO+PjTP6Qrfw2BBRKqkxJJCV8LMBztOYenWq6EiuWtbPJASi5i2IT49NK+2c1/IMxgsiEiVlNjTwpkwIwgC3vz2KBb9cBwA8OCQLnjxrn7w9/PtmWZH00o6XJ1WGhUf6bPBy1v59m8+EamWEksKpe6PUG+24Jm1+xtDxbxRvfDKhASfDxWAes99Iffjbz8RqZan9+SQEmZq6hrwh5W5+CT3LPQ64B/39McTt9/QZOMrX8Ztt30Xp0KISNU8veW1Ncy0dH5HWZUJM7Jysa/wMoID9Fg0KREj4zu5pR6t4rbbvovBgohUz9OHZ7UUZgrLazA1MxsFpdVoGxKA96elILFLW9HP7Sv7Lmj13BdyHYMFEZENtsLMgXNGTFueg9IqEzq3aYWVD6eiR4fWDp/LGiY2HirG+rzzKK+ua/yaty69tE4rzV61xyPHwpN6cOdNIvIIqXfqth4PQLG7/W3HSvHoB7morjOjb1Q4VkxPQfvWQQ7rsbWPw7XcvT250riPhffglt5EpBpSLy62Ht8mJAAAcLmmXtRzyOnzvHN4+tN9qDcLGN6jHd6dkoTtx0sd/kz29nFozl3bk6uFr0z/eDsGCyJSBalnb4i9GLf0HHJatuUkXvn6MADgdwOi8Ob9A/HDzyUOf6ZR8ZEtnnViCw8pIzXjWSFEpDipZ2+09HhbWjq/w/p8zp5RYbEIePm/hxpDxcM3xeFfEwfDX68X9TPtPFkmKVQAXHpJ3oHNm0TkNlLP3nD0eDHPYeXK3L6pwYw/f7ofX+w7DwBYcEdfzLylO4CrR6SL+Zl2nCiT9HMAXHqpZZzu+RWDBRG5jdRNkly5Y7/2e105o6Kyth6PfrAb20+UwV+vwxu/H4gJgzvbfJ2WiR8d4dJLbWODalOcCiEit5G6SZIrd+zW73X26HMAKKmoxQNLdmL7iTKEBvph+fSUJqFCSo3DurdvcXtwKy691DZriG0+imUNsRsOFClUmXIYLIjIbaSeveHo8WKew9kzKk5crMI9GdtxqKgC7VsH4eNHh+HmGzo4/TMN7dHO7vbg13LX9uTkfq6EWG/GYEFEbiP1ILGWHm+Lredw5oyKvWcu4b6M7Th76Qq6tQvBZ7OHI6Gzweb3SfmZ7J11EhEagIdv7IbVM4di2zMjGCo0iget2cYeCyJyKzFnb4h5vK19LGw9h9Tpl+8PX8Ccj/agtt6CgTEGZE5LQbvWQbL9TJ4+68QRNhnKhwet2cZgQURuJ/Xiau/xgOOdN6WcUfFxzhk8u+4AzBYBt/bugMWTExESKO6fRSk/k6fPOrGHTYby4kFrtnGDLCLyOtaGOsD2GRWLJw/G8ZJqvLnxKADgvqQYpN/THwF+zs8Oq30kQOpGZeSY2SLgptc2OQyx3rKjqtjrN0csiMjrtDRV8fc7+2Lb8TJ8uOsMAGDubT3xp9G9oNM5/w+/2kcCHDUZ6nC1yXBUfKRXXAA9hQet2cYRCyLyWs1HEQbEGPDHj/Pw7aEL0OmAF+7qh6nDurn0GloYCdhxogyTlu10+DhuKe4ctQdLuXDEgoh83rW9DZdr6pCWmY3c05cQ6K/HOw8Mwtj+rv2jr5WRADYZupfaGnSVxmBBRF7v3OUrSMvMxvGSKoQH++O9tBRZdrkUu9xwxU8FaB8WpNgFh02G7qeWBl01YLAgIq92pLgSaZnZKK6oRWR4MLJmpKJ3ZJgszy32Dv+lrw43/rcSQ+RSVsoQuYobZBGRZjk6vXTnyTLc9+52FFfU4oaOrfHZY8NlCxWAc3f4Smz1LHWjMiJXcMSCiDTJUcPc1/lFeGpNHurMFqR0a4v3pqbA8MsmW3JxNBJgi1K9F1I3KiNylqRVIRkZGcjIyMCpU6cAAP369cNzzz2HsWPHin5BrgohIlc5Wolxf0oMPsk9C0EAftuvE96ZOBjBAX5urQWQcp7pVUqswlD7fhukXm5ZFRITE4N//OMf6NmzJwAgKysL48ePx969e9GvXz/XKibSEP7jrBwxBz99nHMWADBlaFcsvKufW/9u7I0EiKHEKgw2GZK7SQoW48aNa/LnV155BRkZGdi5cyeDBfkMX1mzrlaOVmJY3Z8cixfH93Np4yuxmi83LK00NWnYtIerMMgbOd1jYTab8emnn6K6uhrDhg2z+ziTyQSTydT454qKCmdfkkhx9obgrQ15atgMyduJvcu/sWc7j4QKq2tHAswWAe9tK+AqDPJJkleF5Ofno3Xr1ggKCsKsWbOwbt06xMfH2318eno6DAZD40dsbKxLBRMpRcwQ/AtfHrpuZQLJS+xdfvvQlk8odSeuwiBfJjlY9O7dG3l5edi5cydmz56NtLQ0HDp0yO7j58+fD6PR2PhRWFjoUsFESpGyGRLDhftYV2I4uiT/6dN9Hl3S2Zy19yLS0DQIRRqCObJFXs3ls0JGjhyJHj16YMmSJaIez1UhpFWf553Dk2vyRD1WiZ4LX2oo3XCgCLN+WYlhj1rO6vClvxfybh47K0QQhCY9FETeSkqjnad7LnytoTQk0B9B/nqYGix2H6OWszq4CoN8jaSpkGeffRZbt27FqVOnkJ+fjwULFuDHH3/E5MmT3VUfkWqIHYIHPNtzYW0obT5No8QOj56wbu9ZzFiRA1ODBf2iWx71tE5PZReUe6Y4IpIWLC5cuIApU6agd+/euP3227Fr1y5s2LABo0aNcld9RKrRUkOeLZ64qPlSQ6kgCFiy+QT++PE+NFgEjB8UjRk3dhP1vTy1UzpH26UT2SNpKuT99993Vx1EmuDMZkjuvKiJbSjNLijX9HC8xSLgpa8OYflPpwAAM2+Ow/yxfbFLZGjjfhG22ev/8LWpNZIXzwohksi6GdKKnwpc3gTJ1cY+saFFy3fspgYz5n2yD1/tvzql87c7++KRm7sD4KmdrrAXHu4aGIWlWwq4Vws5jcGCyAl+eh2m3Rjn0iZIctwVir0T1+ode0VtPR5duRs7TpYhwE+HN34/EOMHdW78unV6avaqPdCh6Vkd3C/CPnsbvRUZa7FkS4HN71FLMyypH49NJ3KSK5sgydVw6aihVIerYUWLd+wXKmpx/7s7sONkGVoH+WPF9NQmocJKi/tFKNm/0FJfjiNshiUxOGJB5AJnjqJ21HAp5a7QW+/Yj5dUIS0zG+cuX0GHsCCsmJ6CftEGu49vflaHmveLULp/QexZKy3R8tQauR+DBZGLpF7U5G64dCbcqNnu05fwcFYOLtfUI659KFbOSEVsRIjD79PCfhFqOGtGjlCg1ak18gwGCyIZSLmouaPhUkt37C3ZeOgCHl+9B7X1FgyKbYPMaSmICA1UuixZyDlS5QpXQgGbYUkMBgsiD3NXw6UW7thbsjr7DBasy4dFAEb06YhFDw5GSKD3/BOllqXBjlbSWHnT1Bp5Fps3iTzMmxsuHbHVtCgIAt7+7ijmf3Y1VNyfHIOlU5K8KlQA6lka7KjpWAfg0VviNNUMS+riXf/nEmmAtzZcOmKraTEyPAg3dAzD1uOlAIDHR/TEvFG9oNN5188OqGtpsJi+nL+M6av5qTVShsunm0rF002JrlJ6dYAn2WtatNLpgJfGJ+ChoV0lPa+WTg41WwTc9Nomh/uebHtmhMd+Bi29f6Q8j51uSkTO8ZaGS0fE7JvQplUAJqV2kfS8Wgtmahyp0npfDqkTeyzI5ym5WZH1H/bxgzpjWI92XhcqAHH7JlyqqZe06ZJWT3TV4mZeRFJxxIJ8mtbuerVI7qZFtSzbdJavjFSR72KwIJ+lhs2KvEVLc/VyNy2qZdmmK8ROQbAHgrSIwYJ8ktbvet3FmQuZo1Gf1LgItGkVgMtX6m1+v9RNl9SybNPdOJpGWsVgQT7JG+565ebMhUzMqE+RsRbGFkIFIK1pUU3LNt2Fo2mkZWzeJJ/kK3e9YjnTDFnXYMGz6w7YHfURAPzpk32NI0O/6dUBkeFBTR7nTNOi1jYYk9oc7Gg0Dbg6mubJJmMiKThiQT7JF+56xXJmWmjDgSI8uy4f5dW2RyKsquvMAIA//7Y3Hru1BywCXO4ZUOOyTXucGQXiaBppHUcsyCdp7a7XnaRcyIBfRzcchQqrB1NjMee2ntDpdLItrxWzbFPJZcSA80tiOZpGWscRC/JJWrrrdTcpFzIxm101N25gZ+cKu4atptKWlm0q3fjoSnMwR9NI6xgsyGeJOS/BF0i5kInZ7Opacoz6OAoJzacD1ND46Mp0hqPTR3l0OakdgwX5NE9vVqTGfQmkXMj+u/+86OfVwfVRH6khQS3LiF2ZzuBoGmkdgwX5PE+dl6D08Lw9Ui5kYkc32oUG4pW7E2z+XGLDlTMhQS2Nj65OZ3A0jbSMwYLIA9QwPN8SsRcyR6MbABARGoAd829HoP+vveHWMLHxUDHW551HeXVd49fshStnQoJaGh/lmM7g1t+kVQwWRG6mluF5R8RcyKyjG7NW7bH5HDoAr97dv0mosDVScy174cqZkKCWxke5pjN4+ihpEZebErmZ1OWcYrhrKaWY5aCj4iNxW+8O130+ysZmV/aWXF7L3qZPzoSE1LgIRIbb/z5PLiPmSabkqzhiQeRmcg/PK9mrUVtvxrxP8vDDkYsAgClDuyK5W1uboxtSlqbamtZwZjph46Fi1DaYbb6GEo2PzkxnqLHBl0gKBgsiN5NzeF7JXg3jlXr8YWUudhWUI9BPjzfvH4hxA6PtPl7q0lSgabiSOp1g772xahMSgPR7+nt8pEDKdIZaG3yJpOBUCJGbybXLp5JnSBQba3H/uzuwq6AcYUH+WDEjpcVQATjXINk8XImdThAzOhLkr8eo+EjJNXmKszt1EqkNRyyI3EyuRj6lllIeL6nE1Pezcd5Yi45hQVgxPRXx0eEOv09Kg2RLqyTETCeIGR0prjCp9nwNrTT4EonBEQsiD5CjkU+JpZS5p8pxb8YOnDfWonuHUKydPVxUqACsjZRBjh/4i5bClaOmUrUsM3WWOxp8iZTCEQsiD3F1XwJPL6X89mAxHl+9F6YGCwZ3aYPMtBS0DQ0U/f1+eh0mpXbBW98dc/jYp0b2cqmHQC3LTJ2l9WBEdC0GCyIPcmVfAk+eIfHhrtP4+/oDsAjAyL4d8X+TEtEq0M/mY1taxdCtfaio1+vWPsSlerV+vobWgxHRtRgsiDTCE2dICIKAt747hn99f3WUYWJKLF6ekAB/P9uzprZWMUSEBuDuQZ0xMj4S7VuLmwpx9YKp9fM1tB6MiK6lEwRB/hbyFlRUVMBgMMBoNCI8XNxcLZE9vrjm311LEhvMFvxt/QGsySkEADx5+w14auQN0Olsv5+OlncCQGR4EGobLDDW1Ld4wdz2zAj46XUu/31qebmm9f0EbAcjbqpFShN7/WawIM3S8kXEVXIHqit1Zsz9aA++/7kEeh3w8oT+eHBIlxZf/6bXNjlciXHt6IG9kQTrBVOuv08th01f/p0m9WOwIK9m726Zd3fSlVfX4eGsHOw9cxlB/nr836TBGN2v5f0edpwow6RlO0U9vw6AISQAwf5+KK6wfcF09Pf57wcHo21okCbDglRaDkbk3cRev9ljQZrDNf/yKSyvQdrybJy8WA1DqwBkTktGUlfH8/hSVicIAC7X1OPDhxOh1+uuu2CK2fhr7uq9uHbfL2++i+fBY6R13MeCNIdr/uVx8LwR92Rsx8mL1ejcphXWzh4mKlQAzjVbfnOoGADwuwHRTfaiELO5VfPNRLkbJZF6MViQ5nDNv+u2Hy/FA0t24mKlCX0iw7B29nD07Bgm+vsdbVNuy8odpzFp2U7c9NqmJoHAmb8nd29hTkTOY7AgzeGaf9d8se880pZno8rUgJ4dQvH06N7oECZ+h0zg1+Wdzmg+2uDs3xNHpojUicGCNEeuQ7180XtbT+KJ1XtRb756l3/8YjUeWZl73SiCGNZtyqMM0oJB89EGZ0Y/rsWRKSJ1YbAgzbn2brn5xUgLmyEpwWIR8OrXh/HyV4dtft3ZnoUxCVHY9swIrJ45FDNu7IYIkVt+Xzva0NLfpxgcmSJSFwYL0iQ5DvXSArNFwI4TZfg87xx2nChzqp+grsGCeZ/kYemWk3Yf40rPgnUVw3Pj+iFnwUisnjkUU4d1FfW91tEGe3+fLeHIFJE6cbkpaZarh3qpnRybJVWZGjDrg93YdrwUfjodzC1sWyPHsevXLpVcueO0w8dfO9owJiEKv+nVEamvfofK2oYWv48jU0TqxREL0jRHx2krxdWRBuuGUc2XYUqZsiiprMXEpTuw7XgpQgL9MPOWOFGvLUfPgjN9MBsOFOHG1753GCoAICI00KtGpoi8CUcsiGTm6kiDHBuAFZRWY2rmLhSWX0G70EAsn56CapMZ7262Px1iJUfPgtRDwcScO3Ktv93Zl6GCSKU4YkFeS47+BKnsjTQUGWsxa9UefL3/vMPncHUDsLzCy7g3YzsKy6+ga7sQrJ09HANi2jSOIrQkMjxItp4FsX0wLQUpu3UaWslSIxHJjyMW5JUcHeftjl4MMRfIuav3YhF0uGOA/bttVzYA++FICR5btQdX6s3o39mA5dNTGo8u99PrcNfAKCzZUmD3OStrG7DxULFsowFi+mDE7LxpxePDidSPIxbkdeyNGpRX1+P9n07Z3P1RDmK3pn7so5Z7JJzdAOzT3EI8kpWLK/Vm3NKrA9b8YWhjqACuBp/P81r+mavrzLJvle2oD0ZqTwcbNonUjcGCvIrYYXV3nDUh5QLZ0rJOqY2PgiDg3z8cx5//sx9mi4B7Ejvj/bRkhAY1HZBctOlYk9NF7REc1Cc3sUGqHRs2iTSBwYK8ithhdXecNSGl6bGlHgkpG4CZLQKe+/wgXv/mCABg9q098ObvByLAr+n/2hsOFOGt747JUp/cxOy8GREagB3zb2eoINIABgvyKlKP85bzAiqmOfJa19bavNF0VHykw8bH2noz5ny4Bx/sPA2dDlg4Lh7PjOkDna7pJdo6iiOVp7bKdhSkdABevbs/Av35zxWRFrB5k7yKM0slfzp+UZZmTusFctaqPaIeb621peWp254ZYbPx0VhTj5krc5F9qhyBfnq89cAg3GmnIVRKc6St+jzBuoKk+fsQKXFDMCJSnk4QWtiKzw0qKipgMBhgNBoRHh7uyZcmH2C2CLjptU0oNtZKWr4odUfLlny9/zzmrt4LezMs1pUN254ZgY2Him3u32CNOLZ6CoqMV5CWmY2jF6oQFuSPpVOTW9wp8/O8c3hyTZ7o+q+tz9NNkmaL4LU7qRJpndjrN8cWyas4e5y3nM2cdwyIxqJJiTa/dm2PBIAWN8Kyfv3aHpCjFypxz+LtOHqhCp3Cg/DJrGEOt9+WMvKg9FbZat1JlYjEkxQs0tPTkZKSgrCwMHTs2BETJkzAkSNH3FUbkVOcOc5b7mbOOwZE4V0bNVzbIyF1I6ycU+W4L2M7ioy16NmxNT577Eb0jXI86iflWHJvO8SNiDxPUo/F5s2bMWfOHKSkpKChoQELFizA6NGjcejQIYSGhrqrRiLJrt2YaeOhYqzPO4/y6roWv0eOQ7js1WBraF/KRlgbDhThiTV5qGuwIKlrW7yflow2IeKOKG9pe22rh2/s5raNw4jIt7jUY3Hx4kV07NgRmzdvxi233CLqe9hjQUowWwS8tfEoFv1w3OFj35k4COMHdXbqNaT0B+w4UYZJy3Y6fN4ZN3bD8u2nIAjAqPhO+L9JgxEc4Ce5PjlOS3UX9lYQqZ/Y67dLq0KMRiMAICLC/va6JpMJJpOpSWFEnuan1+HGnu1FBQtnVkM4c9G2TlG01GjaOsgPmT+dAgBMSu2Cl8b3g7+fc61Raj1mXs2Bh4ikc3rEQhAEjB8/HpcuXcLWrVvtPm7hwoV44YUXrvs8Ryy0R+t3lY5WjDi7GsLeyZwtrexo/r2A7SkKqz+O7IUnbu953R4VWufKe0dEniV2xMLpYDFnzhx89dVX2LZtG2JiYuw+ztaIRWxsLIOFxji6q9RK6HB0IV/84GDcMSC6yeda+tmsYcVeE6aYsGLrvQ3y18PUYIFed3VzqImpXaT/sAqQ8nsgx3tHRJ7j1qmQxx9/HF988QW2bNnSYqgAgKCgIAQFBbX4GFI3e3eV1iWaf7glDl/sK9LEULa9jZisXvrqMPR6XWPdjgKVlJUd9hpCr52iOHmxCpk/FeDExWoEB+ixaFIiRsZ3cu2H9hCpUxpyvHdEpD6SJmsFQcDcuXPx2WefYdOmTYiLi3NXXaQSLR3qJfzysWRLwXUXCHcc8iWXMQlR+Pudtve6uLZue6ekXvsYV444v5afXoeYtq3w3raroaJNSAA+fGSopkKFo/eqObneOyJSF0nBYs6cOVi1ahU++ugjhIWFobi4GMXFxbhy5Yq76iOFObsdtDsO+ZKL2SLgpa9sn51hrXThFwex8IuWA9Wz6/IRIXLJp6OG0APnjLh78XYUlFajc5tW+M+s4Ujq2va6M0TU9l4CjsMnYPv3wNnj4YlI3SRNhWRkZAAAbr311iafX758OaZNmyZXTaQirtwtqnUoW8wQfHGFye7Xrcqr6/HEmr1oExIAY019iw2h1iPObdl2rBSzVu1GlakBfaPCsWJ6CjqFB2tmtYSzUxqOVsWIee+ISH0kT4XY+mCo8F5y3C2qbShbznou1dTj8i+hoqUjzgHYHHn4PO8cpq/IRpWpAcO6t8PHjw5tDBVSpxaU4uyUhpTj4YlIO3i6KbVIzF4LjqhtKFvuenQADCEBCPb3Q3HF9SdzArhu9UOUIRjDe7TD2j3nAAC/GxCFN+8fiCB/P4dTCzpcnVoYFR+piouuK1MaPNWUyPswWFCLxGwHbY9ah7LFDMF3Cg8CoMOFCseBSgBwuaYeHz6cCL1e12Sppb3TS4uMtY2hYsaNcfjbnX2h/yUkaG21hKtTGmrduIuInMPTTckh611lZLMDtaIMwXj0ljjooK2hbDFD8Avv6oeFd0k7JbW02tTkZE7A/umlVmHB/nj2jj6NoQLQ3moJOaY0eKopkffgiAWJ0tJd5eAubTU3lC12CD7joUQ8uy4f5dX1Dp+z+VC/mBU1lbUNyDl1qcnIgxZXS3BKg4isGCw0SKldLq13lc1pdShbTN1jEqIwok8nDE3/3u7pqPaG+p0dedDqagmt/h4QkbwYLDRGrUsQ7YUOtRNTd6C/Hq/enWBzK/CWhvqdHXloqa9FzVNMgGu/B1rZFp6IWubSsenO4LHpzuOBTcqSGurMFgGpr3yHMgcjHfbOwlBriHQHX/pZibTK7YeQOYvBwjk8sEkdpNxVf3/4Amat2o168/X/i4kNg75wF8/ATKQNbj2EjDxPa0sQvZXYof5Pcgoxf10+zBYB/aLDUVplwoVrdvMU29So1SkmsbS2ZwcROcZgoRFaW4LoqwRBwKJNx/HmxqMAgPuSYpB+T3/odTqvH3lwBgMzkfdhsNAILS5B9DVmi4DnPj+AD3edAQDMua0Hnh7dGzrd1QDBC+P1GJiJvA+DhUZodQmir6itN+OJ1Xvx7aEL0OmAF+7qh6nDuildluoxMBN5H+68qRE8sMk2NRwrfrmmDg+9twvfHrqAQH89Fj+YqHioUMP7IoY1MNv7rdXh6uoQBmYi7eCIhYZwd8Om1LBE8fzlK0jLzMaxkiqEBfvjvanJGNJd2SkPNbwvYml5zw4iso3LTTXIF5YgOqKGJYpHiiuRlpmN4opaRIYHI2tGKnpHhrn1NR1Rw/viDC2FISJfxX0syGupYU+PnSfLMHNlLiprG3BDx9bImpGK6Dat3PJaYqnhfXEFAzORunEfC/JaSi9R/F9+EZ78OA91DRakdGuLZVOT0SYkUPbXkUrp98VV3r5nB5GvYLAgzVFyiWLW9lNY+OVBCALw236d8M7EwQgO8JP9dZzBpZtEpAYMFqQoZ4a/lViiKAgCXv/mCBb/eAIA8NDQLnjhrgSPDtU7eq98cekmp0+I1IfBghTjbMOep/f0qDdb8Ne1+Vi75ywA4OnRvTDntp6NG195gpj3ytf2OmHDJ5E6cR8LUoR19ULznoBiYy1mr9qDDQeK7H6vJ/f0qDY1YObKXKzdcxZ+eh3+370DMHfEDR4PFWLeK1/a68SV3x8ici8GC/I4RwdPAVcPnmppUyfrnh6RhqbD+pGGYNmWVJZWmTBp2U78eOQiggP0WDY1CfenxLr8vFJIfa888b4oTY7fHyJyH06FkMfJtXphTEIURsVHumWO/UxZDaZm7sKpshq0DQlA5rQUDO7S1uXnlcqZ98qd74saaH31C5G3Y7Agj5Nz9YI7lijmnzVi+opslFbVIaZtK6yckYruHVrL+hpiOfteefPSTa5+IVI3Bgsb2GnuXmpevbD12EXM+mA3quvMiI8Kx4rpKegYrtwqCjW/V0rhe0KkbgwWzbDT3P3Uunph3d6z+POn+9FgEXBjz3Z496EkhAUHeLSG5tT6XimJ7wmRurF58xrsNPcMta1eEAQBSzafwB8/3ocGi4C7BkZj+bRUxUMFoL73Sg34nhCpG4PFL9hp7llqWb1gsQh46b+Hkf6/nwEAM2+Ow9sPDEKgv3r+11DLe6UmfE+I1IuHkP1ix4kyTFq20+HjVs8c6rVNcUpQsp/F1GDGnz7Zh//uvzoSteCOvph5S3ePvLYz2PtzPb4nRJ7DQ8gkYqe5MpRavVBRW49HV+7GjpNlCPDT4Y3fD8T4QZ09XocU3rzSw1l8T4jUh8HiF+w09x0XKmqRlpmNn4sr0TrIH0umJOHGnu2VLssu3pUTkZYwWPyCnea+4XhJFdIys3Hu8hW0bx2EFdNTkNDZoHRZdnGVEhFpjXo61BTGTnPvt/v0Jdz37nacu3wFce1Dse6x4aoPFVylRERaw2BxDXaae6/vDl3A5Pd24nJNPQbFtsHa2cMRGxGidFl2cZUSEWkVp0Ka8fZzFnzRmuwzeHZdPiwCcFvvDvj35ESEBKr7V5/nYRCRVqn7X1eFsNPcOwiCgHe+P4a3vzsGALg/OQav3t0f/n7qH6jjKiUi0ioGC/JKDWYL/v75QazOPgMAeHxET8wb1Qs6nfwjT+5YtcFVSkSkVQwW5HWu1Jnx+Oq9+O7wBeh0wIvjEzBlaFe3vJa7Vm1wlRIRaZX6x4SJJLhcU4eH3t+F7w5fQKC/HhmTk9waKty1aoOrlIhIqxgsyGucvVSDezO2Y/fpSwgP9seHjwzBmIRIt7yWJ1ZtcJUSEWkRp0JI06z9DXsLL2HZlpO4VFOPKEMwsmakolenMLe9rqdWbXCVEhFpDYMFaZat/gZ/vQ5PjOjp1lABeHbVBlcpEZGWMFiQJln7G5pPNDRYBDy77gDahga6daqAqzaIiGxjjwVpTkv9DVbu3pXSumrD3oSEDldXh3DVBhH5GgYL0pxdJ8tE9ze4C1dtEBHZxmBBmlJvtuCt746Keqy7d6Xkqg0iouuxx4I0o8rUgNmrdiPn1CVRj/dEfwNXbRARNcVgQZpwsdKEGStykH/OiFYBfggO0ONyTb0qdqXkqg0iol9xKoRU71RpNe57dzvyzxkRERqI1X8YivR7+gNgfwMRkdowWJCq7T97GfdmbMfpshrERrTC2tnDMSi2DfsbiIhUilMhpFqbj17E7FW7UVNnRr/ocCyfntKkb4L9DURE6sNgYYM7jsEmadbuPotn1u5Hg0XAzTe0R8ZDSWgddP2vK/sbiIjUhcGiGXcdg022NQ9xKd3aYtnWAry24WcAwIRB0fh/9w1EoD9n7YiItMArgoVcIwz2tom2HoPNuXt52QpxIYF+qKkzAwAevaU7nhnTB3qOFhERaYbmg4VcIwyOjsHW4eo20aPiIzktIgN7Ic4aKn6fFIP5d/T1fGG/4HQYEZFzNB0s5Bxh8NQx2CTurI9tx0thtgiKXMw5HUZE5DzJE9dbtmzBuHHjEB0dDZ1Oh/Xr17uhLMccjTAA0g6i8uQx2L7OUYgD3H/Whz3WsNq8PmtY3XCgyOM1ERFpieRgUV1djYEDB2LRokXuqEc0KSMMYvAYbM9Ra4iTO6wSEfkiyVMhY8eOxdixY91RiyRyX5xS4yLQJiQAl2vq7T6mTUgAj8GWgVpDHKfDiIhc5/YeC5PJBJPJ1PjniooKWZ5XiYuTr7fuydXQaLxiP7wBnj/rw0qtIylERFri9mCRnp6OF154QfbnTY2LQJQhGMXGWlkOosouKG9xtAIALtXU++zdqlwNjR/uOo2/rz9g9+tKnvWh1pEUIiItcfuuQ/Pnz4fRaGz8KCwslOV5/fQ6PD8uHoA8B1HxbtU+ORoaBUHAPzcexYJ1B2ARgIkpsfj3g4MRpaKzPqxh1d5vjA5XwxSnw4iI7HP7iEVQUBCCgoLc8tzWg6ia30lHOnEnzbtV2+TY36PBbMHf1h/AmpyrofKJ22/AH0feAJ1OhzEJUarZL8IaVmev2gMd0ORn5qmpRETiaHofC0C+g6jknlrxFq42NF6pM2PuR3vw/c8l0OuAlyYkYPKQro1fV9tZH3KGVSIiXyQ5WFRVVeH48eONfy4oKEBeXh4iIiLQpUsXWYsTS46LE+9WbXNliqi8ug4PZ+Vg75nLCPLX4/8mDcbofpFylyg7nppKROQ8ycEiNzcXt912W+Of582bBwBIS0vDihUrZCtMCbxbvZ6zU0SF5TVIW56NkxerYWgVgPfTkpHcTTujPWobSSEi0grJweLWW2+FIHjvBkG8W23KmSmiQ+crkLY8GxcrTejcphWyZqSgZ8cwj9VMRETK0XyPhTvwbvVXUqeIth8vxaMf7EalqQF9IsOwYnoqIg2+1fBKROTL3L7clLTPOkXUPCA0Xxr6xb7zSFuejUpTA4bEReDjR4cxVBAR+RiOWJAojqaI3t9WgJf+ewgAcEf/SPzz/kEIDvBTsmQiIlIAgwWJZmuKyGIR8I8NP2PplpMAgGnDu+Hvv/O91TNERHQVgwU5ra7Bgr/8Zx/W550HADwzpg9m/aY7dDqGCiIiX8VgQU6pMjVg1ge7se14Kfz1Orx27wDcmxSjdFlERKQwBguSrKSyFjNW5ODAuQqEBPph8eRE3Nq7o9JlERGRCjBYkCQFpdWYmrkLheVX0C40EMunp2BATBulyyIiIpVgsCDR8govY8aKHJRX16FruxBkTU9Ft/ahSpdFREQqwmBBovxwpASPrdqDK/Vm9O9sQOa0FHQIc8+ptUREpF0MFuTQp7mF+Otn+TBbBNzSqwMyJiciNIi/OkREdD1eHcguQRCw+McTeP2bIwCAewZ3xmv3DUCAHzdsJSIi2xgsyCazRcDCLw7ig52nAQCzftMDz4zpzT0qiIioRQwWdJ3aejOeWpOHDQeLodMBz/0uHtNvjFO6LCIi0gAGC2rCWFOPmR/kIrugHIF+evzzgYH43YBopcsiIiKNYLCgRkXGK0jLzMbRC1UIC/LH0qnJPD6eiIgkYbAgAMDRC5VIy8xGkbEWncKDsGJ6KvpGhStdFhERaQyDBSHnVDkeXpGDitoG9OgQiqwZqYhpG6J0WUREpEEMFj5uw4FiPLFmL+oaLEjq2hbvTU1G29BApcsiIiKNYrDwYR/sPI3nPz8AiwCM7NsJix4cjOAAP6XLIiIiDWOw8EGCIODNb49i0Q/HAQCTUrvgpfH94M+Nr4iIyEUMFj6mwWzBs+vy8UnuWQDAH0f2whO39+TGV0REJAsGCx9SU9eAuR/txaafS6DXAa/e3R8TU7soXRYREXkRBgsfUVZlwoysXOwrvIzgAD0WTUrEyPhOSpdFRERehsHCBxSW12BqZjYKSqvRJiQA76elIKlrW6XLIiIiL8Rg4eUOnDNi+oocXKw0oXObVsiakYqeHVsrXRYREXkpBgsvtu1YKWat2o0qUwP6RIYha0YqOoUHK10WERF5MQYLL/V53jk8/ek+1JsFDOveDkumJiE8OEDpsoiIyMsxWHihZVtO4pWvDwMA7hwQhX/ePxBB/tz4ioiI3I/BwotYLAJe/fow3ttWAACYcWMc/nZnX+j13KOCiIg8g8HCS5gazPjzp/vxxb7zAIBn7+iDmTd358ZXRETkUQwWXqCyth6zVu3GT8fL4K/X4fXfD8Ddg2OULouIiHwQg4XGlVTUYtryHBwqqkBooB8yHkrCLb06KF0WERH5KAYLDTt5sQpTM7Nx9tIVtG8diBXTU5HQ2aB0WURE5MMYLDRq75lLmLEiB5dq6tGtXQiyZqSia7tQpcsiIiIfx2ChQZt+voDHPtyD2noLBsYY8P60FLRvHaR0WURERAwWWvNJTiHmr8uH2SLgN706YPHkRIQG8a+RiIjUgVckjRAEAYs2HcebG48CAO5NjME/7u2PAD+9wpURERH9isFCA8wWAc99fgAf7joDAJhzWw88Pbo396ggIiLVYbBQudp6M55csxffHLwAnQ5YOK4f0oZ3U7osIiIimxgsVOxyTR0eycpF7ulLCPTX4+0HBuGO/lFKl0VERGQXg4VKnb98BWmZ2ThWUoWwYH+8NzUZQ7q3U7osIiKiFjFYqNCR4kqkZWajuKIWkeHBWDEjBX0iw5Uui4iIyCEGC5XZdbIMj6zMRWVtA27o2BpZM1IR3aaV0mURERGJwmChIv/LL8KTH+ehrsGClG5tsWxqMtqEBCpdFhERkWgMFiqRtf0UFn55EIIAjI7vhH9NGozgAD+lyyIiIpKEwUJhgiDgjW+P4N8/nAAATB7SBS+OT4CfnntUEBGR9jBYKKjebMH8z/Lxn91nAQB/GtULc0f05MZXRESkWQwWCqk2NWDOR3vw45GL8NPr8OrdCXggpYvSZREREbmEwUIBZVUmzFiRg31njQgO0GPx5ESM6NNJ6bKIiIhcxmDhYWfKajA1cxdOldWgbUgA3p+WgsQubZUui4iISBYMFh6Uf9aI6SuyUVpVh5i2rZA1IxU9OrRWuiwiIiLZMFh4yNZjFzHrg92orjMjPiocK6anoGN4sNJlERERyYrBwgPW7T2LP3+6Hw0WATf2bId3H0pCWHCA0mURERHJjsHCjQRBwLKtJ/Hq1z8DAO4aGI03fj8Qgf56hSsjIiJyDwYLN7FYBLz81WFk/lQAAHjkpjg8e0df6LnxFREReTGnbp0XL16MuLg4BAcHIykpCVu3bpW7Lk0zNZjxxJq9jaFiwR198bffxTNUEBGR15McLD7++GM89dRTWLBgAfbu3Yubb74ZY8eOxZkzZ9xRn+ZU1NZjWmYO/ru/CAF+OrwzcRBm3tJd6bKIiIg8QicIgiDlG4YMGYLExERkZGQ0fq5v376YMGEC0tPTHX5/RUUFDAYDjEYjwsPDpVesYhcqapGWmY2fiysRGuiHJVOScdMN7ZUui4iIyGVir9+Seizq6uqwe/du/PWvf23y+dGjR2P79u02v8dkMsFkMjUpzBsdL6lCWmY2zl2+gvatg7BiegoSOhuULouIiMijJE2FlJaWwmw2o1OnpttPd+rUCcXFxTa/Jz09HQaDofEjNjbW+WpVShAE/Pk/+3Du8hXEtQ/FuseGM1QQEZFPcqp5s/npm4Ig2D2Rc/78+TAajY0fhYWFzrykqul0OrzzwGCM7NsR/5k1DLERIUqXREREpAhJUyHt27eHn5/fdaMTJSUl141iWAUFBSEoKMj5CjWiS7sQvJeWonQZREREipI0YhEYGIikpCRs3Lixyec3btyI4cOHy1oYERERaY/kDbLmzZuHKVOmIDk5GcOGDcPSpUtx5swZzJo1yx31ERERkYZIDhYPPPAAysrK8OKLL6KoqAgJCQn4+uuv0bVrV3fUR0RERBoieR8LV3nzPhZERETeSuz1m6dhERERkWwYLIiIiEg2DBZEREQkGwYLIiIikg2DBREREcmGwYKIiIhkw2BBREREsmGwICIiItkwWBAREZFsJG/p7SrrRp8VFRWefmkiIiJykvW67WjDbo8Hi8rKSgBAbGysp1+aiIiIXFRZWQmDwWD36x4/K8RiseD8+fMICwuDTqfz5Eu7VUVFBWJjY1FYWMgzUBzgeyUe3yvx+F6Jx/dKPL5XvxIEAZWVlYiOjoZeb7+TwuMjFnq9HjExMZ5+WY8JDw/3+V8+sfheicf3Sjy+V+LxvRKP79VVLY1UWLF5k4iIiGTDYEFERESyYbCQSVBQEJ5//nkEBQUpXYrq8b0Sj++VeHyvxON7JR7fK+k83rxJRERE3osjFkRERCQbBgsiIiKSDYMFERERyYbBgoiIiGTDYCGTxYsXIy4uDsHBwUhKSsLWrVuVLkmVtmzZgnHjxiE6Oho6nQ7r169XuiRVSk9PR0pKCsLCwtCxY0dMmDABR44cUbosVcrIyMCAAQMaNzAaNmwY/ve//yldliakp6dDp9PhqaeeUroU1Vm4cCF0Ol2Tj8jISKXL0gQGCxl8/PHHeOqpp7BgwQLs3bsXN998M8aOHYszZ84oXZrqVFdXY+DAgVi0aJHSpaja5s2bMWfOHOzcuRMbN25EQ0MDRo8ejerqaqVLU52YmBj84x//QG5uLnJzczFixAiMHz8eBw8eVLo0VcvJycHSpUsxYMAApUtRrX79+qGoqKjxIz8/X+mSNIHLTWUwZMgQJCYmIiMjo/Fzffv2xYQJE5Cenq5gZeqm0+mwbt06TJgwQelSVO/ixYvo2LEjNm/ejFtuuUXpclQvIiICr7/+Oh5++GGlS1GlqqoqJCYmYvHixXj55ZcxaNAgvP3220qXpSoLFy7E+vXrkZeXp3QpmsMRCxfV1dVh9+7dGD16dJPPjx49Gtu3b1eoKvI2RqMRwNULJtlnNpuxZs0aVFdXY9iwYUqXo1pz5szBnXfeiZEjRypdiqodO3YM0dHRiIuLw8SJE3Hy5EmlS9IEjx9C5m1KS0thNpvRqVOnJp/v1KkTiouLFaqKvIkgCJg3bx5uuukmJCQkKF2OKuXn52PYsGGora1F69atsW7dOsTHxytdliqtWbMGe/bsQU5OjtKlqNqQIUOwcuVK9OrVCxcuXMDLL7+M4cOH4+DBg2jXrp3S5akag4VMmh8BLwiCVx0LT8qZO3cu9u/fj23btildimr17t0beXl5uHz5MtauXYu0tDRs3ryZ4aKZwsJCPPnkk/j2228RHBysdDmqNnbs2Mb/7t+/P4YNG4YePXogKysL8+bNU7Ay9WOwcFH79u3h5+d33ehESUnJdaMYRFI9/vjj+OKLL7BlyxbExMQoXY5qBQYGomfPngCA5ORk5OTk4J133sGSJUsUrkxddu/ejZKSEiQlJTV+zmw2Y8uWLVi0aBFMJhP8/PwUrFC9QkND0b9/fxw7dkzpUlSPPRYuCgwMRFJSEjZu3Njk8xs3bsTw4cMVqoq0ThAEzJ07F5999hk2bdqEuLg4pUvSFEEQYDKZlC5DdW6//Xbk5+cjLy+v8SM5ORmTJ09GXl4eQ0ULTCYTDh8+jKioKKVLUT2OWMhg3rx5mDJlCpKTkzFs2DAsXboUZ86cwaxZs5QuTXWqqqpw/Pjxxj8XFBQgLy8PERER6NKli4KVqcucOXPw0Ucf4fPPP0dYWFjjiJjBYECrVq0Urk5dnn32WYwdOxaxsbGorKzEmjVr8OOPP2LDhg1Kl6Y6YWFh1/XphIaGol27duzfaebpp5/GuHHj0KVLF5SUlODll19GRUUF0tLSlC5N9RgsZPDAAw+grKwML774IoqKipCQkICvv/4aXbt2Vbo01cnNzcVtt93W+GfrXGVaWhpWrFihUFXqY126fOuttzb5/PLlyzFt2jTPF6RiFy5cwJQpU1BUVASDwYABAwZgw4YNGDVqlNKlkYadPXsWkyZNQmlpKTp06IChQ4di586d/HddBO5jQURERLJhjwURERHJhsGCiIiIZMNgQURERLJhsCAiIiLZMFgQERGRbBgsiIiISDYMFkRERCQbBgsiIiKSDYMFERERyYbBgoiIiGTDYEFERESyYbAgIiIi2fx/HNUSsa3A9aoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prostate = pd.read_table(\n",
    "    \"https://static-1300131294.cos.ap-shanghai.myqcloud.com/data/ml-fundamental/parameter-optimization/gradient-descent/prostate.data\"\n",
    ")\n",
    "prostate.drop(prostate.columns[0], axis=1, inplace=True)\n",
    "\n",
    "X = prostate.drop([\"lpsa\", \"train\"], axis=1)\n",
    "y = prostate[\"lpsa\"]\n",
    "\n",
    "regressor = LinearRegression()\n",
    "\n",
    "regressor.fit(X, y)\n",
    "y_pred = regressor.predict(X)\n",
    "\n",
    "print(regressor.__dict__)\n",
    "print(y - y_pred)\n",
    "\n",
    "plt.scatter(y, y_pred)\n",
    "plt.plot([0, 5], [0, 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LinearRegressionWithSGD:\n",
    "    def __init__(self, learning_rate=0.0003, n_iters=5000):\n",
    "        self.lr = learning_rate\n",
    "        self.n_iters = n_iters\n",
    "        self.weights = None\n",
    "        self.bias = None\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        n_samples, n_features = X.shape\n",
    "\n",
    "        # init parameters\n",
    "        self.weights = np.zeros(n_features)\n",
    "        self.bias = 0\n",
    "\n",
    "        batch_size = 5\n",
    "        # stochastic gradient descent\n",
    "        for _ in range(self.n_iters):\n",
    "            # approximate y with linear combination of weights and x, plus bias\n",
    "            y_predicted = np.dot(X, self.weights) + self.bias\n",
    "\n",
    "            indexes = np.random.randint(0, len(X), batch_size)  # random sample\n",
    "\n",
    "            Xs = np.take(X, indexes, axis=0)\n",
    "            ys = np.take(y, indexes, axis=0)\n",
    "            y_predicted_s = np.take(y_predicted, indexes)\n",
    "\n",
    "            # compute gradients\n",
    "            dw = (1 / batch_size) * np.dot(Xs.T, (y_predicted_s - ys))\n",
    "            db = (1 / batch_size) * np.sum(y_predicted_s - ys)\n",
    "            # update parameters\n",
    "            self.weights -= self.lr * dw\n",
    "            self.bias -= self.lr * db\n",
    "\n",
    "    def predict(self, X):\n",
    "        y_predicted = np.dot(X, self.weights) + self.bias\n",
    "        return y_predicted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'lr': 0.0003, 'n_iters': 5000, 'weights': array([ 0.4438507 ,  0.21653236, -0.00196098,  0.08407097,  0.14496092,\n",
      "        0.13755223,  0.1197549 , -0.00640268]), 'bias': 0.021768646451838486}\n",
      "0    -1.108122\n",
      "1    -0.759352\n",
      "2    -0.798181\n",
      "3    -0.658291\n",
      "4    -1.016655\n",
      "        ...   \n",
      "92    1.736796\n",
      "93    1.300247\n",
      "94    2.055891\n",
      "95    2.676134\n",
      "96    2.001249\n",
      "Name: lpsa, Length: 97, dtype: float64\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABK+ElEQVR4nO3deVzUdf4H8NcM9zkKiICi4i3iyaGYuWXqZuVqt64HSlmadrltdmz3Qe22u/XbNlMTwSztMNMuS9fSTBFQ8b5PVA5BYThkgJnv7w8bFGSY7wzfme/3O/N6Ph48HgkDvAfDz2s+n/fn89EIgiCAiIiISAJauQsgIiIi18FgQURERJJhsCAiIiLJMFgQERGRZBgsiIiISDIMFkRERCQZBgsiIiKSDIMFERERScbT2d/QZDLh/PnzCAoKgkajcfa3JyIiIjsIgoCKigpERUVBq7U8L+H0YHH+/HlER0c7+9sSERGRBPLz89GxY0eLH3d6sAgKCgJwpbDg4GBnf3siIiKyg16vR3R0dMM4bonTg4V5+SM4OJjBgoiISGWstTGweZOIiIgkw2BBREREkmGwICIiIskwWBAREZFkGCyIiIhIMgwWREREJBkGCyIiIpIMgwURERFJhsGCiIiIJGNTsHj55Zeh0WgavUVERDiqNiIiIlIZm4/07tu3LzZs2NDwZw8PD0kLIiIiIvWyOVh4enpyloKIiIiaZXOPxdGjRxEVFYWYmBhMnDgRJ06caPHxBoMBer2+0RsREZEa5OWXYcpH26GvqZO7FNWwKVgMGTIEy5Ytw48//ojFixejsLAQw4YNQ2lpqcXPSUtLg06na3iLjo5uddFERESO9vOhYkxalIUtx0rwr5+OyF2OamgEQRDs/eSqqip069YNTz/9NObNm9fsYwwGAwwGQ8Ofzfe5l5eX89p0IiJSpM9z8/HsV3thNAn4Q892+GDyYAT42Nw94FL0ej10Op3V8btVP6WAgAD069cPR48etfgYHx8f+Pj4tObbEBEROYUgCPjvz8fwzu8zFHcN7oC37+4PLw+eziBWq35SBoMBBw8eRGRkpFT1EBERycJoEvDimv0NoeKRm7rhn/cOYKiwkU0zFk899RTGjRuHTp06obi4GK+//jr0ej1SUlIcVR8REZHD1dQZ8cTKPKzbXwiNBnh5XF+kDOsid1mqZFOwOHv2LCZNmoSSkhK0a9cOQ4cORVZWFjp37uyo+oiIiByqvLoOM5flIvvURXh7aPHuxIG4rR9n4u1lU7BYuXKlo+ogIiJyuvNllzF9aTaOFFUiyNcTi6clYGjXULnLUjX3bnElIiK3daSoAinp2Sgor0H7YB9kpiahdwR3K7YWgwUREbmd7JMX8WBmDvQ19egeHojM1CR0aOMnd1kugcGCiIjcyrp9BXhsZR5q601I6NwWH6UkoI2/t9xluQwGCyIichsfbzuFF9fuhyAAY2Lb4/8mDYKvFy/TlBKDBRERuTxBEPDOT4fx35+PAwAmD+mEV8fHwUOrkbky18NgQURELq3OaMJzX+3FFzvOAgD+Mron5o7sDo2GocIRGCyIiMhlVdfWY84nO/Hz4Qvw0GrwxoQ4TEzqJHdZLo3BgoiIXFJppQGpmbnYnV8GXy8t/vvnwbilT3u5y3J5DBZERORy8i9WY1p6Nk6WVKGtvxeWTE/E4E5t5S7LLTBYEBGRS9l3rhzTl+agpNKADm38sOyBJHRrFyh3WW6DwYKIiFzGlqMlePjjXFTVGtEnMhiZMxIRHuwrd1luhcGCiIhcwpq8c3jqi92oMwoY1i0UH06NR7Cvl9xluR0GCyIiUr3Fm0/gje8PAgDGDYjCO/f2h48nD76SA4MFERGplskk4I3vD2LJlpMAgAeGx+D52/pAy4OvZMNgQUREqmSoN+KvX+zB2t3nAQDP39YHM0d0lbkqYrAgIiLVqaipw8Mf78DW46Xw8tDgH/cMwIRBHeQui8BgQUREKlOsr0HK0hwcLNAjwNsDC6cmYHiPMLnLot8xWBARkWocv1CJlPRsnL10GWGBPsiYkYi4Djq5y6JrMFgQEZEq7DxzCQ9k5OBSdR1iwgKQOSMJnUL95S6LmmCwICIixfvfwSLM+XQnaupMGNBRh/TpiQgN9JG7LGoGgwURESnaZzln8NzqfTCaBNzcqx3+O3kw/L05fCkV/2aIiEiRBEHA+xuP4Z/rjwAA7onviLS7+sHLQytzZdQSBgsiIlIco0nAi2v24ZPtZwAAc2/ujr+M6QmNhgdfKR2DBRERKUpNnRGPrdiFnw4UQaMBXv1TX0xN7iJ3WSQSgwURESlGWXUtHszMRe7pS/D21OL/Jg7ErXGRcpdFNmCwICIiRThXdhkp6dk4VlyJYF9PfJSSiKSYELnLIhsxWBARkewOFeoxPT0HhfoaROp8kZmahJ7tg+Qui+zAYEFERLLKOlGKmctyUVFTjx7hgchMTUJUGz+5yyI7MVgQEZFsvt9bgCdW5qHWaEJSlxAsnpYAnb+X3GVRKzBYEBGRLDK3nsLL3+yHIAB/7Nse700cBF8vD7nLolZisCAiIqcSBAH/+PEwPvjlOABg6tDOePlPfeGh5RkVroDBgoiInKbOaMIzq/Zi1c6zAICnxvTEnJu78+ArF8JgQURETlFlqMcjn+zEpiMX4KHVIO2ufrgvIVruskhiDBZERORwJZUGpGbkYM/Zcvh5eeC/kwdhZO/2cpdFDsBgQUREDnW6tAop6dk4VVqNkABvpE9PxMDoNnKXRQ7CYEFERA6z92w5ZmRko6SyFh3b+mFZahK6tguUuyxyIAYLIiJyiM1HLmD28h2oqjWib1Qwls5IRHiQr9xlkYMxWBARkeRW7zqLv36xB/UmAcO7h2HBlMEI8uXBV+6AwYKIiCQjCAIWbT6BtB8OAQDGD4zCP+4ZAG9PrcyVkbMwWBARkSRMJgGvfXcAS387BQCYeWMMnh3bB1oefOVWGCyIiKjVDPVGzPt8N77bUwAA+NvtffDgjV1lrorkwGBBREStoq+pw0PLcpF14iK8PDR4594BGD+wg9xlkUwYLIiIyG5F+hqkpGfjUGEFAn08sXBqPG7oHiZ3WSQjBgsiIrLLseIKpKTn4FzZZbQL8kHGjET0jdLJXRbJjMGCiIhstuP0JTyQmYOy6jp0DQtAZmoSokP85S6LFIDBgoiIbLL+QBEeXbETNXUmDIxug/TpiQgJ8Ja7LFIIBgsiIhJtRfYZPL96L0wCMLJ3ON7/8yD4e3Mooav4fwMREVklCALe+99RvLvhKADgvoSOePPOfvD04MFX1BiDBRERtajeaMILa/ZhRXY+AOCxkd3x5Oie0Gh48BVdj8GCiIgsulxrxKMrdmHDwSJoNcCr4+MwZWhnucsiBWOwICKiZl2qqsUDmTnYeaYMPp5a/N+kQfhj3wi5yyKFY7AgIqLrnL1UjZT0bBy/UIVgX08smZ6IxC4hcpdFKsBgQUREjRws0CMlPRvFFQZE6XyRmZqEHu2D5C6LVILBgoiIGmw7XoqHluWiwlCPXu2DkJGaiEidn9xlkYowWBAREQDg2z3nMe+z3ag1mpAUE4LF0xKg8/OSuyxSGQYLIiLC0t9O4tVvD0AQgLFxEfj3/QPh6+Uhd1mkQgwWRERuzGQS8PaPh7Bw0wkAQEpyZ7w4ri88tDyjguzTqiPT0tLSoNFo8MQTT0hUDhEROUud0YSnvtjdECqevrUXXv4TQwW1jt0zFjk5OVi0aBH69+8vZT1EROQElYZ6zF6+A78eLYGHVoO37+6Pe+I7yl0WuQC7ZiwqKysxefJkLF68GG3btpW6JiIicqALFQZMWpSFX4+WwN/bA0tSEhgqSDJ2BYs5c+bg9ttvx6hRo6w+1mAwQK/XN3ojIiJ5nCqpwt0LtmLvuXKEBnhjxcyhuKlXuNxlkQuxeSlk5cqV2LlzJ3JyckQ9Pi0tDa+88orNhRERkbT2nC3DjKU5KK2qRacQfyxLTUKXsAC5yyIXY9OMRX5+Ph5//HEsX74cvr6+oj7n2WefRXl5ecNbfn6+XYUSEZH9fjlcjImLslBaVYu4DsFYNXsYQwU5hEYQBEHsg7/++mvceeed8PC4urfZaDRCo9FAq9XCYDA0+lhz9Ho9dDodysvLERwcbH/lREQkyqodZzF/1R7UmwTc2CMMC6bEI9CHpw2QbcSO3zb9n3XLLbdg7969jd43Y8YM9O7dG/Pnz7caKoiIyHkEQcCCTcfx93WHAQB3DuqAt+/uD2/PVp00QNQim4JFUFAQ4uLiGr0vICAAoaGh172fiIjkYzQJePWb/cjcdhoA8PAfumL+H3tDyzMqyME4F0ZE5GJq6oyY93kevt9bCI0GeOH2WKQOj5G7LHITrQ4Wv/zyiwRlEBGRFMov1+GhZbnYfvIivD20+Nf9A3BH/yi5yyI3whkLIiIXUVheg5T0bBwuqkCQjycWTovHsG5hcpdFbobBgojIBRwrrsC0Jdk4X16D8CAfZMxIQmwUd96R8zFYEBGpXO6pi3ggMxfll+vQtV0AlqUmoWNbf7nLIjfFYEFEpGI/7i/EYyt2wVBvwuBObbAkJRFtA7zlLovcGIMFEZFKfbL9NF74eh9MAjCqTzj+M2kw/Lx5nhDJi8GCiEhlBEHAv9cfwf9tPAYAmJQUjdfGx8HTgwdfkfwYLIiIVKTeaMLzq/fhs9wr9y49fksPPDGqBzQaHnxFysBgQUSkEpdrjZj76U7871AxtBrg9Qn98OchneQui6gRBgsiIhW4WFWLBzJzsOtMGXw8tXj/z4MxOra93GURXYfBgohI4fIvViMlPRsnSqrQxt8LS1ISEN85RO6yiJrFYEFEpGD7z5dj+tIcXKgwoEMbP2SmJqJ7eJDcZRFZxGBBRKRQW4+V4KGPd6DSUI/eEUHITE1C+2BfucsiahGDBRGRAq3dfR5/+TwPdUYBQ7uGYOHUBOj8vOQui8gqBgsiIoX56NcTeP27gwCA2/tH4l/3DYCPJw++InVgsCAiUgiTScBb6w5h0eYTAIDpw7rgxTtiodXyjApSDwYLIiIFqK034ekvd+PrvPMAgGfG9sbDI7ry4CtSHQYLIiKZVRrqMevjHdhyrASeWg3+fk9/3DW4o9xlkcoYTQKyT15EcUUNwoN8kRQTAg8ZZrsYLIiIZFRcUYMZS3Ow/7we/t4eWDAlHn/o2U7uskhl1u0rwCvfHEBBeU3D+yJ1vnhpXCxujYt0ai28sYaISCYnS6pw94Kt2H9ej7BAb6x8aChDBdls3b4CzF6+s1GoAIDC8hrMXr4T6/YVOLUeBgsiIhnk5Zfh7gVbkX/xMjqH+mPV7GHo37GN3GWRyhhNAl755gCEZj5mft8r3xyA0dTcIxyDwYKIyMl+PlSMSYuycLGqFv076rBq9jB0Dg2QuyxSoeyTF6+bqbiWAKCgvAbZJy86rSb2WBAROdEXufl45qu9MJoEjOjZDgsmD0aAD/8pJvsUV1gOFfY8Tgr8v5mIyAkEQcB/fz6Gd346AgC4a3AHvH13f3h5cOKY7BceJO6Id7GPkwKDBRGRgxlNAl5eux8fZ50GAMy+qRue/mMvnlFBrZYUE4JInS8Ky2ua7bPQAIjQXdl66iyMykREDlRTZ8ScT3bi46zT0GiAV/7UF/Nv7c1QQZLw0Grw0rhYAFdCxLXMf35pXKxTz7NgsCAicpDy6jpMW5KNdfsL4e2hxX//PBgpw7rIXRa5mFvjIrFgymBE6Bovd0TofLFgymCnn2PBpRAiIgcoKL+MlPRsHCmqRJCvJxZNTUByt1C5yyIXdWtcJEbHRvDkTSIiV3SkqAIp6dkoKK9B+2AfZKYmoXdEsNxlkYvz0GoUEV4ZLIiIJJR98iIezMyBvqYe3cMDkZmahA5t/OQui8hpGCyIiCSybl8BHluZh9p6E+I7t8WSlAS08feWuyyXppSLt+gqBgsiIgl8vO0UXly7H4IAjI5tj/9MGgRfLw+5y3JpSrp4i67irhAiolYQBAHv/HgYL6y5Eir+PKQTFkwezFDhYEq7eIuuYrAgIrJTndGE+av24P2fjwEA5o3uiTcmxMGTp2k6lBIv3qKr+H8/EZEdqmvr8dCyXHyeexZaDfDWXf3w2C09ePCVEyjx4i26ij0WREQ2Kq00IDUzF7vzy+DrpcX7kwZjVGx7uctyG0q8eIuuYrAgIrJB/sVqTEvPxsmSKrT198KS6YkY3Kmt3GW5FSVevEVXMVgQEYm071w5pi/NQUmlAR3a+GHZA0no1i5Q7rLcjhIv3qKr2GNBRCTClqMluH/hNpRUGtAnMhhfPTKMoUImSrx4i65isCAismJN3jnMyMhGVa0Rw7qF4rOHh6J9MKfZ5aS0i7foKi6FEBG1YPHmE3jj+4MAgDv6R+Kf9w2AjyfPqFACJV28RVcxWBARNcNkEvDm9wfx0ZaTAIAHhsfg+dv6QMtBS1GUcvEWXcVgQUTUhKHeiL9+sQdrd58HADx/Wx/MHNFV5qqI1IHBgojoGhU1dXj44x3YerwUnloN3rl3ACYM6iB3WUSqwWBBRPS7Yn0Npi/NwYECPQK8PfDh1Hjc2KOd3GURqQqDBRERgOMXKpGSno2zly4jLNAHGTMSEddBJ3dZRKrDYEFEbm/XmUtIzcjBpeo6dAn1x7LUIegU6i93WUSqxGBBRG7tfweLMOfTnaipM2FARx3SpyciNNBH7rKIVIvBgojc1mc5Z/Dc6n0wmgTc1KsdPpg8GP7e/GeRqDX4G0REbkcQBLy/8Rj+uf4IAOCe+I5Iu6sfvDx4GDFRazFYELkoo0ngiYTNMJoEvLhmHz7ZfgYAMPfm7vjLmJ7QaPizIfvx9+0qBgsiF7RuXwFe+eYACsprGt4XqfPFS+Ni3foOhZo6Ix5bsQs/HSiCRgO88qe+mJbcRe6ySOX4+9YY5/2IXMy6fQWYvXxno3/kAKCwvAazl+/Eun0FMlUmr7LqWkz5aDt+OlAEb08tPvjz4BZDhdEkYNvxUqzJO4dtx0thNDV3QTe5O/6+XY8zFkQuxGgS8Mo3B9DcECjgypXSr3xzAKNjIySdplX6NPC5sstISc/GseJKBPt64qOURCTFhFh8PF+Bkhhy/b4pHYMFkQvJPnnxuldO1xIAFJTXIPvkRckublL6IHy4sAIp6dko1NcgItgXmalJ6BURZPHx5legTQcL8ytQXslNZnL8vqkBl0KIXEhxheV/5Ox5nDVKnwbOOlGKez7cikJ9DXqEB+KrR4a1GCqsvQIFrrwC5bIIAc7/fVMLBgsiFxIe5Cvp41qi9EH4+70FmLYkGxU19Ujs0hZfzhqGqDZ+LX6OLa9AiZz5+6YmNgWLBQsWoH///ggODkZwcDCSk5Pxww8/OKo2IrJRUkwIInW+sLSaq8GVZYqW+gvEUvIgnLn1FOZ8uhO1RhP+2Lc9Pn5gCHT+XlY/T+wryx/2FaiyoZMNqdJy5u+bmtjUY9GxY0e89dZb6N69OwAgMzMT48ePx65du9C3b1+HFEhE4nloNXhpXCxmL98JDdBoNsH8j99L42IlaSRT4jSwIAj4x4+H8cEvxwEAU4d2xst/6iv6+Yp9Zbls22ks23ZaUb0k1ii9F0aNnPn7piY2zViMGzcOt912G3r27ImePXvijTfeQGBgILKyshxVHxHZ6Na4SCyYMhgRusaDZITOV9LGQ6VNA9cZTXjqiz0NoeKpMT3x6njxoQKw/gq0KaX0klij9F4YNXPW75ua2L0rxGg04osvvkBVVRWSk5MtPs5gMMBgMDT8Wa/X2/stiUikW+MiMTo2wqFbQM2DcGF5TbN9Fhpc+cfVGdPAVYZ6PPLJTmw6cgEeWg3S7uyH+xKjbf46Lb0CbY4athRyS6TjOeP3TU1sbt7cu3cvAgMD4ePjg1mzZmH16tWIjY21+Pi0tDTodLqGt+ho23/Zich2HloNkruFYvzADkjuFir5P3LmQRjAda/wnTkNXFJpwKTFWdh05AJ8vbRYPC3erlBhZukVqCVS9JI4svdByb0wrsTRv29qYvOMRa9evZCXl4eysjKsWrUKKSkp2LRpk8Vw8eyzz2LevHkNf9br9QwXRC7CPAg3XbuPcNLa/ZnSakxL345TpdVo6++F9OmJGNSpbau/7rWvQH/YV4Bl205b/Rx7e0kc3fugxF4Ycm02Bwtvb++G5s2EhATk5OTgvffew8KFC5t9vI+PD3x8fFpXJREpllzTwHvPlmNGRjZKKmvRsa0flqUmoWu7QMm+vvkVKABRwcKeXhJnHMaltF4Ycn2tPnlTEIRGPRRE5H6uHYSdYfORC5i9fAeqao2IjQxGRmqiwwZGqXtJzMefF5ZfxmvfHXR474OSemHIPdgULJ577jmMHTsW0dHRqKiowMqVK/HLL79g3bp1jqqPiKiR1bvO4q9f7EG9ScAN3UPx4ZR4BPlaP6PCXlJuKWxu2cMSqY6D5pZIcjabmjeLioowdepU9OrVC7fccgu2b9+OdevWYfTo0Y6qj4gIwJXZ0YWbjuPJz3aj3iRg/MAoLJ2e5NBQYSbFlkJLWz6tkaL3gVsiyZk0giA49eg1vV4PnU6H8vJyBAcHO/NbE5FKmUwCXvvuAJb+dgoAMPPGGDw7tg+0Tn6Vbe8trkaTgOFvb7Q5VADAiplDJVtmUvottKRsYsdv3m5KRIpmqDdi3ue78d2eK4c4/e32Pnjwxq6y1GJvL4m1LZ/NcUTvg7N7Ycg9MVgQkWLpa+rw8LId2HaiFF4eGrxz7wCMH9hB7rJsZutyBnsfqDlqmXFisCAiRSrS1yAlPRuHCisQ6OOJhVPjcUP3MLnLsoutO1acdQ4IqYea7nphsCAixTlWXImU9GycK7uMdkE+yJiRiL5ROrnLspu1LZ8AEBLghRfu6IuIYOW+EiV5OOO8EynZfKQ3EZEj7Th9Cfd8uBXnyi4jJiwAX80epupQAVg//lwD4M07++HOQTwOmhqzdtcLcOW8EymPgW8tBgsiUoz1B4ow+aMslFXXYWB0G6yaPQzRIf5ylyUJbvkke6jxrhcuhRCRIqzIPoPnV++FSQBG9g7H+38eBH9v1/onirdgkq3UeNeLa/3WEpHqCIKA9/53FO9uOAoAuC+hI968sx80Gg22HS9VzQAstmOfWz7JFmq864XBgohkU2804YU1+7Ei+wwA4NGR3TFvdE/8uL/wug74kAAv3DmwA0bFRiguZNjSsa+WLYOkDGq864UnbxKRLC7XGvHoil3YcLAIWg3w6vg4TBna2WIH/LWUtM3OUr3mqHBt/4SatgyScpj/HwOav+vFWT06YsdvNm8SkdNdqqrF5I+ysOFgEbw9tfhgcjymDO3cYgf8tczb7NbtK3BKvZbY0rFv6a4QpTwXUi61Nf5yKYSInMK8BHCoUI/Fv57A+bIaBPt6Ysn0RCR2uTKNK/boaymvFW8NsR37WcdLWwwgSngupGxqavxlsCAih2tuCUCrAZ4Y1aMhVAC2dbabB+2M304iLMhHln9oxda77USJ6C2DrtjYyb4Saail8ZfBgogcylIPgkkAXvv2IKLa+DVM5drT2f7adwcb/tvZ/Qri6xU3iDp7y2BLA75UYYB9Je6HwYKIHEZMz8S1SwBijr5uibOPOBZTb0SwD5K7heL9n49Z/XrO3DLY0oAPQJIwoLajqEkabN4ksoPRJGDb8VKsyTuHbcdLFXWcrpLYemrgtUdf28PZRxy3dFS3WU29CeXVtYjU+Vp8jAZXBm5nbRlsqZF01vKdmCVBk6kaj6ImaTBYENlo3b4CDH97IyYtzsLjK/MwaXEWhr+9kV39TZhMAj769YSox167BGDugI/U2ffq3dlHHJvr1fl7Nfvx8uo6PPLpLsR3bmvxHALAeVekixnwm2NrGFDjUdQkDQYLIhtwy6A4dUYTnvpiN/53qFjU45suAdwaF4kt80dixcyhSL2hC0ICvG2uwZn9CqNjI+Dr6dHsx8xD8Ld7mv9/w9lbBsXuvGmOLWFAjUdRkzTYY0EkkrVXetwyeEWloR6zl+/Ar0dL4KHVINDHA/rL9RZfDbfx92p2CcDcAZ/cLRTP3x7b0EhYUmFo1LBpiTP7FbJPXkSh3vYB8slRPTB3ZA9F7mRp7ddQ41HUJA3OWBCJxKld6y5UGDBpURZ+PVoCPy8PfJSSgLfv7t/iFHtZdR3WHyhs8euaQ8b4gR0w/YYYRfUrAPYN1hoAK3PypS/GCikGcjFfw9zYqqS/J3IOBgsikRwxtetKTaCnSqpwz4dbsfdcOUIDvLHyoaG4uVc4RsdGoI2F/gPg6kyP2OfeUsOks/sVzOwZrOUKotYG/JbYEgaU+PdEzsFgQSSS1FO7rtQEuudsGe5esBWnS6sRHeKHL2cPw4DoNgCuzPSUVddZ/Fx7BlilHXHcmsHa2T0GYgb8lj5mSxhQ2t8TOQd7LIhEkvKWQVfa3//L4WI88slOVNcaEdchGEunJ6FdkE/Dxx3VxKekI47Ng/Xs5TuhQcu7K5qSo8fAPOA3PasiooVzLCLsPNRKSX9P5BwMFkQitTR42PJqzpWaQFftOIv5q/ag3iTgxh5hWDAlHoE+jf9ZcWQTn5KOOLY0WFsi93XX1gZ8KcOAkv6eyPEYLIhsYO2VnphXc7Y0gSr1H2NBELBg03H8fd1hAMCEgVH4+z0D4O15/eqqlDM9Std0sD5VUo13NxwBYH8QdaSWBnyGAbIXgwWRjVo7tav2/f1Gk4DXvj2AjK2nAAAPj+iK+bf2htbC87d3pketF1c1HZB7RQRKtqxApAYMFkR2aM2rOTXv76+pM2Le53n4fu+V7aEv3BGLB4bHWP08W2d61HJxlZjwwx4DcjcMFkROptalgfLLdXhoWS62n7wIbw8t/nnfAIwbECX688UOsGppbBUbftQ680JkL40gCE7dOK/X66HT6VBeXo7g4GBnfmsixTAPnkDzSwNKGTzNCstrkJKejcNFFQjy8cTCafEY1i1M8u9jNAkY/vZGiz0o5tC1Zf5IWQdnS+Gn6d+fWmZeiMQQO37zHAsiGahpf/+x4grc9cFvOFxUgfAgH3z2cLJDQgWgjtNNxd7a+f0e3itD7olLIUQyUcPae+6pi3ggMxfll+vQtV0AMmckITrE32HfTw2NrWLDz9/W7HOJLcVEtmKwILcn5xq4krf0/bS/EI+u2AVDvQmDOrVBekoi2tpxy6gt1NDYKjbUXKyqtfgxNWwpJrIXgwW5NVdaA5cyIH2y/TRe+HofTAIwqk84/jNpMPy8m78WXErWGlsBICTAC4X6Gmw7Xor4zm2x4/Qlp4ZCKUONUrcUE7UGgwW5LbXsPhBDqoBUbzTh6VV78NXOcwCA+xKi8eadcfD0cE47lpijsS9W1eHJz/IAAFoNcO3dZc4IhWJ29bQN8MLFKsv3o5gpcUsxUWuxeZPcktgGPDXcNmoOSK1tEvxuz3n0e/mnhlABAJuPFGPDwSJJ67XGUmNrc5r+9TijMVLMJV6vj4/jleHkthgsyC2pYfeBGFIFpLV55zDn0124XGds9P4ivUGWHQy3xkViy/yRWDFzKP593wCEiOztcFYotLar57b+US2GDwHAxMRofLvnPLYdL1VFgCUSi0sh5JbUsPtADCnuHblQYcBfvtht8fPl2sFgbmzddry0xUbIppzVGGltV4+l00Z1/l4AgH9vONrwvmuXcHigFqkdgwW5JTXsPhCjtQEp/2I17l+4DXVGy6+Y5d7BYG+4K66ocfggbW1Xj6VLySz19Tw0IgZrdxe4RDMxuS8GC3JLaj1Wu6nWBKT958sxfWkOLlQYRH2N1szetGaAtzfcnSqpvu4UTzkGaXP4MJ8q2tKy1cLNJ6/7mBqbicm9sceC3JKYBjy5r7QWwxyQbG0S3HqsBPcvzMKFCoPoA6/sHeDX7SvA8Lc3YtLiLDy+Mg+TFmdh+NsbRfdtWHuOTWkAtPH3wrsbjijq1Etry1aWqK2ZmIjBgtyWmo7VtsSegLR293mkLM1GpaEeQ2JCsHbODQ7bwSDFjpWWnmNztZqHXqXt+GnNjI9amomJAAYLcnPX7j54b+JArJg5FFvmj1RFqDCzJSB99OsJPLZiF+qMAm7vF4nM1CS0DfB2yOyNlFt6LT3HpiVF6Hzx5KgeKKu2fIaEXIO0FP06Sm8mJgLYY0Gk6GO1xbK2Q8FkEvDWukNYtPkEAGD6sC548Y5YaK3sYIhoRU+CFDtWrD3H5k7e/HbPeVH1OXuQFnOqqDVKbyYmAhgsiFyGpYBUW2/C01/uxtd5VwbcZ8b2xsMjukKjafxyX+pL0Ryxpbe552j+s7lB9GhRhaiv5exBWsypopaopZmYCGCwIFIlsbssKg31mPXxDmw5VgJPrQZ/v6c/7hrc0eLXlXL2xplbeps70twSOQdpSzNDkTpf/GlAJBb9vivk2tChpmZiIoDBgkh1xN4LUlxRg9SMHOw7p4e/twcWTInHH3q2c1qdztrSa+nOl+YoYZBuaWZoUKe2ki5HEclBIwiCU1uj9Xo9dDodysvLERwc7MxvTaR6lgZR8xBpbtY8WVKFaenbkX/xMkIDvLF0RiL6d2zj5Gqv1gs0/ypc7O4bSzM05rMhxG7jVMNhUzx5k5RK7PjNGQsilbC2y8J89HZ4kC8eXJaLi1W16Bzqj8wZSegSFuDkaq+Qoim0pRkanZ+3qFAx9+buuKF7mCoG6WuXoxgySI0YLIhUQuwui0mLs2CoN6FfBx2WzkhEWKCP84psRtOp/7AAH0ADlFQasO14aYuDpbWr7VNv6CKqhh7tA1W380fskheR0jBYEEnAGa8sxe6eMNSbMKJnOyyYPBgBPo77FbflOZtfha/bV4CnvtwtarAUM0OzOu9cMx+9ntq2aVoLVGo5wI3cE4MFUSs565Wl2MHxxh5hWJKSAC8Px51/Z89ztnWwFDNDc7GqDiEB3rhUVavqO1+uJXbJy9m3zRKJxZM3iVpBiiOrxRJzZ0aAjweWTk90eKiw9Tnbcwqn2BmaCQOjAKj7zpdr2XKwmD2MJgHbjpdiTd45bDteyvtHSHIMFkR2kvLIajHE3Jnxz3sHwNOBocLe52zPYCl2hmZ0bITq73y5liMOFjNr7YVwRGJwKYTITlIfWS2GeZfFS2v3o0h/9brzNv5eeOuufg4fRO19zvYMlracg+Gh1Uh6aqicHHWwGPs2yFkYLIjs5MhXli0ZEN0Gwb5eKNIb4OulxV/H9ML0G2KcMoja+5ztGSxbOgK7uWWOpqeGmqf81RY0HHGwGPs2yJkYLIjs5Mwjq82OFFUgJT0bBeU1aB/sg4wZSegT6byD5ux9zvYOlrfGReKhETFY/OtJXHuUn0YDzLwxpsVGUbVu1bQ1UIkhx+wauS+bFmPT0tKQmJiIoKAghIeHY8KECTh8+LCjaiNSNGvNlBpcGcyk2pGQc+oi7lmwFQXlNegeHoivHrnBqaECsP85t9Qf0tJguW5fARZtPommbSomAVi0+WSzvQHObKh1FEvXxNvbNyLX7Bq5J5uCxaZNmzBnzhxkZWVh/fr1qK+vx5gxY1BVVeWo+ogUy97B0h7r9hVg8kfboa+pR3zntvhyVjI6tPFr9de1lbUGUgHAC7f3afY52zpYtjR9b9a0UdTZDbWOdGtcJLbMH4kVM4fivYkDsWLmUGyZP9KuGRc5ZtfIfbXqrpALFy4gPDwcmzZtwogRI0R9Du8KIVfj6Gn3j7NO48U1+yAIwOjY9vjPpEHw9fJo9ddtjZZuE7X23MUerLXteCkmLc6yWsuKmUMbpu/t+Rx3YL5TxdpS1Jb5I9ljQRY55a6Q8vJyAEBIiOWpXoPBAIPhave6Xq9vzbckUpyWbqtsDUEQ8M+fjuD9n48BACYldcJr4/s6dDupWLfGRcJkAh75dOd1Hysor8Gs5TvxoYUpe7FXs9szfc8p/+Y5om+DyBK7/4USBAHz5s3D8OHDERcXZ/FxaWlp0Ol0DW/R0dH2fksixTIPluMHdkByt9BW/wNdZzRh/qo9DaHiyVE98eadcYoIFcCVV8CvfXegxcf85YvdqK032f097Jm+55S/ZVL3bRBZYvdSyJw5c/Ddd99hy5Yt6Nixo8XHNTdjER0dzaUQIguqa+sx55Od+PnwBWg1wJt39sPEpE5yl9WI2CWHQB9PvHNvf7sGLXum7znlbx1vTCV7iV0Ksevlz6OPPoq1a9fi559/bjFUAICPjw+Cg4MbvRFR80orDZi0eDt+PnwBvl5aLJqaoLhQAYhfSqg01Nu9E8Oe5lhnNtSqldSza0RN2RQsBEHA3Llz8dVXX2Hjxo2IiYlxVF1Ebif/YjXu+XAbdueXoY2/Fz55cChGxbaXu6xm2bKUIMD+nRj2TN9zyp9IXjY1b86ZMweffvop1qxZg6CgIBQWFgIAdDod/Pycv/WNyFXsO1eO6UtzUFJpQIc2fshMTUL38EC5y7IoKSYEbfy8UHa5TtTjW3P4kj3NsY5qqLUHlx7I3djUY6HRNP/LsHTpUkyfPl3U1+B2U6LGthwtwazlO1BpqEefyGBkzEhE+2DlNxe+t+EI/r3hqPjHTxyI8QM7OLAi5VHzCaBETTmkx0IQhGbfxIYKImpsTd45zMjIRqWhHsldQ/HZw0NVESoAYO7IHmjj7yX68WGBPg6sRnlc4QRQInsoY+8akYOZL6Rak3cO246XKuLkxcWbT+DxlXmoMwq4o38kMlITEewrfqCWm4dWg7fu6if+E+T/kTuNK50ASmQrXkJGLk9p09Emk4A3vz+Ij7acBACk3hCDv93eB1oVrrvfGheJB27ogiW/nbL62JIqg9XHuApe+kXujMGCXJp5Orrp60LzdLSzdwkY6o346xd7sHb3eQDAc7f1xswbu1rsX3IUc0NhYfllXKyqRUigDyKC7WssHBUbISpYuNOhVDwBlNwZgwW5LGvT0RpcmY4eHRvhlC79ipo6zFq+A78dK4WnVoN37h2ACYOc28xoNAl4f+MxLP3tZLM7OuyZybH3SnQxtap1NwVPACV3xmBBLkvO6eimg2KXUH88kJmLAwV6BHh7YMGUeIzo2U7S72nNun0FeOarvSirtrxFtMCOmRxH3EOhtOUrWzkqbBGpAYMFuSy5pqObGxQ9NBoYBQFhgd7ImJGEuA46Sb+nmJpmLb/+wjBLbJ3JMR9K1fR5R4gMA9cGsVMl1Xh3wxG7l6+UMNPBS7/InTFYkMuSYzraUk+H8ffjYh6/pYdNoUKKQdK8JCSWvTM59h5K1dIV7E3rsrZ8paSZjtaGLSK1YrAgVWtp4HX2dHRLPR1mH/xyHH8e0llUOJBqkLS2JGSJPTM5Yq9EN7MUxCxpKfQorVEXUNYJoETOwmBBqmVt4HX2dLSYAVzsTICUg6S9Sz2ObiwUE8QsafqclNaoey1bwxaR2vGALFIlMacaGk0CdH7emHFDF7QN8G70OEdcSCVVT4fUhyvZGhA0uBLQHN1YaO9MCnD9c7KlUZeIHIszFqQ6Yl6dPvPVXry89gAK9VcHm5AAL9w5sANGxUY4ZDo6LEDckdXNDfTXLumUVBgk3c1ibUmoOc5oLLRnJsXS8hXPjSBSDgYLUh0xr06vbKlsvK3yUlUd0n87hUQHhIqaOiMytp1q8TGWBkWxzYtNiR0kW1oSasqZjY72zKQAzYcenhtBpBwMFiQJZ27xs/dV57Vr7SN7t8eO05ckqbesuhYPZuYi9/QleGo1qDcJons6bG1evJYtg6SlHQo6X0+Mjm2PG3q0s/vkTXuZZ1LEBqqQAG+8cWdcs6GH50YQKQeDBbWas7f4teZVp3kZYWja/3Cxqrbh/fbWe77sMlLSs3G0uBJBvp74aFoCLlXXitpiaG/zor2DpNJ2KJhnUsSer/G32/tY/PvhuRFEyqERBMGp1+uJvc+d1MHSK27zP9+O2OJnNAkY/vZGm3oGrLGn3sOFFUhJz0ahvgYRwb7ITE1Cr4ighhqtDeDbjpdi0uIsh9epdO9tOIJ/bzhq9XErZg4VtZtGKedYELkaseM3ZyzIbnJt8bOlZ0AsW+vNOlGKmctyUVFTjx7hgchMTUJUG79GNVobBO1Z0nHFw5XmjuyBFdn5jRptr2XLDI3SZmWI3BGDBdlNzrs4LJ5qGOyDmnoTyqvrbA4cYuv9YW8BHv8sD7X1JiR2aYvF0xLQxt/b4uMtEbuk88LtfRAW5OOyg6SHVoPxAyOxcPNJi4+xZRmD50YQyYvBguwm9xY/S69O1x8obNVsRkv1Lv3tJF79fZYmsUtbZMxIQoCPfb9GYhsOp98Q43Jh4lrr9hVgUQuhYlRsuEvN0BC5Oh6QRaIZTQK2HS/Fmrxz2Ha8FGGB9p/bIBXzq9PxAzsguVsoPLSahtmMCF3j7xsS4CXqazZXryAImPVxbqOln5xTlzDqX5uwbl+B3bW/NC4WwNXeCTN3aTgU08C6/kAxvt9j38+YiJyPzZskSnNNcdaWHcyvuLfMHynL4Ni0gTK+c1v84R8/W50haFpvndGEqUu2I+vE9ac2StFM6c4Nh2IbWEMDvJH9/CiXDllESsfmTZKMpZ0fRXpDw/uUuMWvubV2W7ckVhnq8cgnO5sNFYA0Taru3HAodpmstKrWIb06RCQ9BgtqkZidHzp/L/h6ejTq6lfq7gVbrrIuqTQgNSMHe86Wt/g1zU2fWcdLodVq7AoH7tpwaMsyGY/jJlIHBgtqkdjjsz95YLDdg6qziZkhOFNajWnp23GqtBoB3h6oqjVa/bpzPt2JsstXjxF3l+WM1kiKCUFIgBcuVtVZfSyP4yZSBzZvUovEvkosqTJc10SpZM01fZrtPVuOuxb8hlOl1ejY1g+vjY8T9TWvDRVA45tWqXkeWg1eF/HzdcZtq0QkDQYLapG7Xe7069ELmLhoG0oqaxEbGYyvZg/DHQOiEBJg+zkV9lxx3lTTnTj2fh0lu61/FB4eEWPx4xq4/u4YIlfCpRBqkTtd7rR611n89Ys9qDcJuKF7KD6cEo/fjpXglW8ONLpXxBatOSTMVXaLiDne/NnbYjGgY1v8bc0+Se5wISL5MFhQi9zhcidBELBo8wmk/XAIAPCnAVF4594B2HioyOrNo238vK5bAmmOrY2HlnbimJdX1HJXSHPhqI2fF2bc0AVzR/Zo9P/Nbf0j8cc499wdQ+RKuBRCVlk6cCpC56uaAc4Sk0nAa98ebAgVM2+Mwbv3D4SHVmP14KaQAC/8Z9IgUd/HlqUiaztxgNYtrziLORw1bf4tu1yHf284ivjX11/Xf9JS7wsRqQNnLEgUVzxrwVBvxF8+341vfz/V8fnb+mDmiK4Arhzc1NJuGAC4WFUHrUYj+VKRnHewSEXMiZpl1XWqmn0hInEYLFRIzJq1I7jSWQv6mjo8vGwHtp0ohZeHBu/cOwDjB3Zo+Lgtu2GkXiqS+w4WKVgLR2YCHHMDLhHJh8FCZVyloU9ORfoapKRn41BhBQJ9PLFwajxu6B7W6DG27IZJ7hYq+tAtMaHQFXbi2BJ6lD77QkS2YbBQEVdp6JPTseJKpKRn41zZZYQF+iBjRiLiOuiue5ytu2HELBWJDYWusBPH1tCj5NkXIrINmzdVwlUa+uS04/Ql3PPhVpwru4yYsACsfmRYs6ECsO/m0ZYaDy01MjZ3iJYr3HpqDkdiKXn2hYhsw2ChErY09NH1NhwowuSPslBWXYeB0W2wavYwRIf4t/g5Uu2GsScUqn0nzrXhqCUa8FRNIlfDpRCVcIWGPrmszD6D51bvhUkAbu7VDv+dPBj+3uL+15diN4y9uzzUvhPn1rhIfDhlMJ75ai/Kqq8/60Mtsy9EZBsGC5VwhYY+ZxMEAe/97yje3XAUAHBfQke8eWc/eHrYNlHX2t0wrQmFat+JYw5H7288hqW/nWx0mJhSb8AlotZhsFAJV2joc6Z6owkvrNmPFdlnAACPjuyOeaN7QqNx/itjdw+FHloNHh/VA3NHdlft7AsRicdgoRLucLS2VC7XGvHoil3YcLAIGg3w6vg4TB3aWbZ6GAqvUPvsCxGJw+ZNFVF7Q58zlFXXYsqS7dhwsAjenlosmBwva6gAXGOXBxGRWBpBEJy6P1Gv10On06G8vBzBwcHO/NYuQ66TN5Xu7KVqpKRn4/iFKgT7emLJ9EQkdlHOLAAPNyMiNRM7fnMpRIU4pXy9gwV6TF+ajSK9AZE6X2SmJqFn+yC5y2pgNAnQ+Xnj6T/2wsWqWoQE+iAimKGQiFwPgwWp3rbjpXhoWS4qDPXo1T4IGamJiNT5yV1Wg5ZmKhgqiMjVsMeCVO2b3ecxdcl2VBjq0TsiCCseGqq4UCH2xE0iIlfAYEGq9fSXu/Hoil2o//3EykOFFbj9/35VzGDNY9iJyB0xWJDqCIKAhz/Oxee5Z6/7mJJmAngMOxG5IwYLUpU6owlPfr4bP+4vavbjSpoJ4DHsROSOGCxINSoN9UjNyMHXu861+DilzAS4+4mbROSeuCuEVOFChQGpGTnYe64c3h5a1BpNVj9HzEyAI88EUfKJmzwLhYgchcGCZCVmgDtVUoWUpdk4XVqNkABvPDmqJ15Ys8/q17Y2E+DoA6uUegw7D+oiIkfiyZskGzED3J6zZZixNAelVbWIDvHDstQh6BTij+Fvb7Q6E7Bl/kiLg7Z5G2jTzzc/Wqoj0o0modmbPeUayJ31vInI9fDkTVI0SwOceVfHgimD4eftidnLd6C61oi+UcFYOiOxYRaiNTMB1raBanCl+XN0bESrZhOaC05t/Lww44YumDuyh81fu7XLF8563kTk3hgsyOnEDHDzV+1FlaEe9SYBN/YIw4Ip8Qj0ufq/q/lCtqYDd4SImQBbtoHae3S6peBUfrkO7244il4RQTbNDEixfOGM501ExGBBTidmgCv/fdlgwsAo/P2eAfD2vH4D061xkRgdG2Hzq3hHbwOVemZAzOxO03DR3OyGK25/ZRMqkfIwWJDTiR24RvYOx7/uGwhtCwOFPReyOXobqJQzA0aTgJfX2hZSLM1uTEyMFlW/Wra/sgmVSJl4jgU5ndiBa+aNXVsMFfYybwO19JU1uDJA2bsNVMqZgfc3HkWhXvzpnS3dTfLvDUfRxt/LYc/bmXgHC5FyMViQ01kb2AHHDnDmbaAArqtBim2gUs2IrNtXgH9vOCrqaxVX1IhagjFzxPN2Ft7BQqRsNgeLzZs3Y9y4cYiKioJGo8HXX3/tgLLIlV07sDdHA8cPcObmzwhd48E9Qufb6i2XUsyImAdPscKDfEUtwZRV1+GJUT0d8rydhXewECmbzT0WVVVVGDBgAGbMmIG7777bETWRG7g1LhIvjuuD1749iGtfWDpzjdze5k9rpDgYy9rgeS1zSPl2z3lRj+8S5o8t80eqtunRFZtQiVyJzcFi7NixGDt2rCNqITeSe+oi3t1wDCYBiNL54sEbY9AnUuf0Ac6e5k8xWrMdFrBtUDSHFFuWYBz1vJ2Bd7AQKZvDd4UYDAYYDIaGP+v1ekd/S1K4n/YX4tEVu2CoN2FQpzZYkpKIkABvucuSXGtmRMQOik+O6tkQUpR8N4mU3OV5EqmVw5s309LSoNPpGt6io8VteSPX9Mn205i1fAcM9Sbc0jscnz441CVDhZl5ZmD8wA5I7hYqejZGTINrRLAP5o7s3uh7ObIpVSnc5XkSqZXDg8Wzzz6L8vLyhrf8/HxHf0tSIEEQ8K/1R/D86n0wCcDExGgsnBoPP28PuUtTJGuDpwbAy3/qe93g6cimVCVxl+dJpEYOXwrx8fGBj4+Po78NKVi90YS/fb0PK3OuhMrHbumBJ0f1gEbDV5QtsbdPw1FNqUrjLs+TSG148iY51OVaI+Z+uhP/O1QMrQZ4bUIcJg/pLHdZqmHv4Knm5kxbuMvzJFITm4NFZWUljh071vDnkydPIi8vDyEhIejUqZOkxZG6XayqxQOZOdh1pgw+nlr8Z9IgjOkbIXdZqsPBk4jUxOZgkZubi5tvvrnhz/PmzQMApKSkICMjQ7LCSN3yL1YjZWk2Tlyogs7PC0tSEpDQxTW79HkRFhHRVTYHi5tuugmCwKNyybID5/VIWZqNCxUGdGjjh8zURHQPD5K7LIfgRVhERI3xrhCS1NZjJbh/4TZcqDCgd0QQVs0e5tKhghdhERE1xmBBklm7+zxSlmajwlCPITEh+Ozh5Ou2A6qF0SRg2/FSrMk7h23HS6+70IoXYRERNY+7QkgSS7acxGvfXrk067Z+EfjXfQPh66XOMyrELG/YchEWGy+JyJ1wxoJaxWQS8Ob3BxtCxfRhXfCfSYNVHSrELG/wIiwiouZxxoLsVltvwtNf7sbXeVdu1Zx/a2/M+kNXmw6+UtKOCmvLGxpcWd4YHRvBi7CIiCxgsCC7VBrqMevjHdhyrASeWg3evrs/7o7vaNPXUNqOCluWN3gRFhFR87gUQjYrrqjBxEXbsOVYCfy9PfBRSoJdoUJpOypsWd7gRVhERM1jsCCbnCypwt0LtmLfOT1CA7yx8qGhuKlXuE1fQ6k7Kmxd3uBFWERE1+NSCImWl1+G1IwcXKyqRedQf2TOSEKXsACbv46UOyqk7NGwZ3mDF2ERETXGYEGi/Hy4GI8s34nLdUb066BD+vREtAuy79ZaqXZUSN2jYV7emL18JzRAo3DR0vIG7/IgIrrKJZZCrB1mRK3zRW4+HszMxeU6I0b0bIeVDw21O1QAti85NMdRPRpc3iAiah3Vz1gobWeBKxEEAR/8chz/+PEwAOCuQR3w9j394eXRujza2h0VtmwLtWdJgssbRET2U/WMhRJ3FrgKo0nAi2v2N4SKWX/ohn/eN6DVoQJAq3dU2NKj0Zoak7uFYvzADkjuFspQQUQkkmqDhSN3Frj70kpNnRFzPtmJj7NOQ6O5Msg/M7a3TQdfWdOaJQeeeklEpFyqXQpx1F0N7r60Ul5dh5kf5yL75EV4e2jxr/sH4I7+UQ75XvYuOfDUSyIi5VJtsHDEq1bz0krT+Qnz0oqrN+8VlF9GSno2jhRVIsjHE4umJTh8t4M9Oyp46iURkXKpdilE6let1pZWBLj2NdhHiipw1wdbcaSoEu2DffD5rGTFbqHkqZdERMql2mBhftVqaejQ4MoShthXrdaWVoDWNwQqVc6pi7hnwVYUlNegW7sArJo9DH0ig697nJJ6T7gtlIhImVS7FGLvYUaWFOrFLZmIfZxarNtXiMdW7kJtvQnxndvio2kJaBvg3czjlNd7wm2hRETKo9oZC0DaV60XKw2SPk4NPs46jUc+2YHaehNG9WmPTx4cYjFUKHVbL7eFEhEpi2pnLMyketUa0syA2prHKZkgCPjnT0fw/s/HAACTkjrhtfF94dnMGRWOPoyKiIhci+qDBSDNXQ0ROj9JH6dU9UYTnlu9F5/nngUAPDmqJx67pbvFMyocta2XiIhck0sECymYm0FbGkRtaQZVouraesz9dBc2HiqGVgO8eWc/TEzq1OLn8DAqIiKyhap7LKR07RZGS9S8hbG00oBJi7dj46Fi+HppsWhqgtVQAfAwKiIisg1nLNxA/sVqTEvPxsmSKrTx98KSlETEd24r6nN5GBUREdmCMxa/MzcpWmJuUlTbAVn7zpXjrgVbcbKkCh3a+OHLWcNEhwqAh1EREZFtGCx+54wbM51ty9ESTFyUhQsVBvSOCMJXjwxD9/BAm78OD6MiIiKxuBTyO1drUlyTdw5PfbEbdUYByV1DsXBaPIJ9vez+ejyMioiIxGCw+J0rNSku3nwCb3x/EABwe/9I/Ou+AfDx9Gj115ViWy8REbk2BovfuUKToskk4M3vD+KjLScBAKk3xOBvt/eBlrMKRETkJOyx+J3amxQN9UY88VleQ6h47rbeeOEOhgoiInIuBotrqLVJsaKmDqkZOVi7+zw8tRr8+/4BeGhEN4unaRIRETkKl0KaUFuTYrG+BtOX5uBAgR4B3h5YMCUeI3q2k7ssIiJyUwwWzVBLk+KJC5WYlp6Ns5cuIyzQGxkzkhDXQSd3WURE5MYYLFRq15lLSM3IwaXqOnQJ9UdmahI6hwbIXRYREbk5BgsV2nioCI98shM1dSYM6KjDkumJCAv0kbssIiIiBgu1+TwnH8+u3gujScAferbDB5MHI8CHf41ERKQMHJFUQhAEvL/xGP65/ggA4O7BHfHW3f3g5cGNPUREpBwMFipgNAl4cc0+fLL9DABgzs3d8NSYXtxOSkREisNgoXA1dUY8vnIXftxfBI0GeHlcX6QM6yJ3WURERM1isFCwsupaPJiZi9zTl+DtqcW79w/Ebf2UeUgXERERwGChWOfLLiMlPRtHiysR5OuJj6YlYEhX5Z+tQURE7o3BQoEOF1YgJT0bhfoaRAT7IiM1Eb0jguUui4iIyCoGC4XZfqIUDy7LRUVNPXqEByIzNQlRbfzkLouIiEgUBgsF+WFvAR7/LA+19SYkdmmLxdMS0MbfW+6yiIiIRGOwUIjMrafw8jf7IQjAmNj2+L9Jg+Dr5SF3WURERDZhsJCZIAh456fD+O/PxwEAk4d0wqvj4xR7myoREVFLGCxkVGc04dmv9uLLHWcBAH8Z3RNzR3bnwVdERKRaDBYyqTLUY86nO/HL4Qvw0Grw5p1xuD+xk9xlERERtQqDhQxKKw1IzcjB7rPl8PXS4oPJgzGyd3u5yyIiImo1BgsnO1NajWnp23GqtBpt/b2wZHoiBndqK3dZREREkmCwcKK9Z8sxIyMbJZW16NjWD5mpSejWLlDusoiIiCTDYOEkvx69gFkf70BVrRGxkcHImJGI8GBfucsiIiKSFIOFE6zedRZ//WIP6k0Cbugeig+nxCPI10vusoiIiCTHYOFAgiBg8a8n8Ob3hwAAfxoQhXfuHQBvT63MlRERETkGg4WDmEwCXv/uINJ/OwkAeHB4DJ67rQ+0PPiKiIhcmF0vnT/44APExMTA19cX8fHx+PXXX6WuS9UM9UY8tnJXQ6h4/rY++NsdsQwVRETk8mwOFp999hmeeOIJPP/889i1axduvPFGjB07FmfOnHFEfaqjr6nD9PQcfLunAF4eGrw3cSBmjugqd1lEREROoREEQbDlE4YMGYLBgwdjwYIFDe/r06cPJkyYgLS0NKufr9frodPpUF5ejuDgYNsrVrAifQ1S0rNxqLACAd4eWDg1AcN7hMldFhERUauJHb9t6rGora3Fjh078MwzzzR6/5gxY7B169ZmP8dgMMBgMDQqzBUdK65ESno2zpVdRligDzJmJCKug07usoiIiJzKpqWQkpISGI1GtG/f+Pjp9u3bo7CwsNnPSUtLg06na3iLjo62v1qFEgQBf/1yN86VXUZMWABWPzKMoYKIiNySXc2bTW/fFATB4o2czz77LMrLyxve8vPz7fmWiqbRaPDe/YMwqk84vpyVjOgQf7lLIiIikoVNSyFhYWHw8PC4bnaiuLj4ulkMMx8fH/j4+NhfoUp0CvXHRymJcpdBREQkK5tmLLy9vREfH4/169c3ev/69esxbNgwSQsjIiIi9bH5gKx58+Zh6tSpSEhIQHJyMhYtWoQzZ85g1qxZjqiPiIiIVMTmYHH//fejtLQUr776KgoKChAXF4fvv/8enTt3dkR9REREpCI2n2PRWq58jgUREZGrEjt+8zYsIiIikgyDBREREUmGwYKIiIgkw2BBREREkmGwICIiIskwWBAREZFkGCyIiIhIMgwWREREJBkGCyIiIpKMzUd6t5b5oE+9Xu/sb01ERER2Mo/b1g7sdnqwqKioAABER0c7+1sTERFRK1VUVECn01n8uNPvCjGZTDh//jyCgoKg0Wic+a0dSq/XIzo6Gvn5+bwDxQr+rMTjz0o8/qzE489KPP6srhIEARUVFYiKioJWa7mTwukzFlqtFh07dnT2t3Wa4OBgt/+fTyz+rMTjz0o8/qzE489KPP6srmhppsKMzZtEREQkGQYLIiIikgyDhUR8fHzw0ksvwcfHR+5SFI8/K/H4sxKPPyvx+LMSjz8r2zm9eZOIiIhcF2csiIiISDIMFkRERCQZBgsiIiKSDIMFERERSYbBQiIffPABYmJi4Ovri/j4ePz6669yl6RImzdvxrhx4xAVFQWNRoOvv/5a7pIUKS0tDYmJiQgKCkJ4eDgmTJiAw4cPy12WIi1YsAD9+/dvOMAoOTkZP/zwg9xlqUJaWho0Gg2eeOIJuUtRnJdffhkajabRW0REhNxlqQKDhQQ+++wzPPHEE3j++eexa9cu3HjjjRg7dizOnDkjd2mKU1VVhQEDBuD999+XuxRF27RpE+bMmYOsrCysX78e9fX1GDNmDKqqquQuTXE6duyIt956C7m5ucjNzcXIkSMxfvx47N+/X+7SFC0nJweLFi1C//795S5Fsfr27YuCgoKGt71798pdkipwu6kEhgwZgsGDB2PBggUN7+vTpw8mTJiAtLQ0GStTNo1Gg9WrV2PChAlyl6J4Fy5cQHh4ODZt2oQRI0bIXY7ihYSE4B//+AceeOABuUtRpMrKSgwePBgffPABXn/9dQwcOBDvvvuu3GUpyssvv4yvv/4aeXl5cpeiOpyxaKXa2lrs2LEDY8aMafT+MWPGYOvWrTJVRa6mvLwcwJUBkywzGo1YuXIlqqqqkJycLHc5ijVnzhzcfvvtGDVqlNylKNrRo0cRFRWFmJgYTJw4ESdOnJC7JFVw+iVkrqakpARGoxHt27dv9P727dujsLBQpqrIlQiCgHnz5mH48OGIi4uTuxxF2rt3L5KTk1FTU4PAwECsXr0asbGxcpelSCtXrsTOnTuRk5MjdymKNmTIECxbtgw9e/ZEUVERXn/9dQwbNgz79+9HaGio3OUpGoOFRJpeAS8IgktdC0/ymTt3Lvbs2YMtW7bIXYpi9erVC3l5eSgrK8OqVauQkpKCTZs2MVw0kZ+fj8cffxw//fQTfH195S5H0caOHdvw3/369UNycjK6deuGzMxMzJs3T8bKlI/BopXCwsLg4eFx3exEcXHxdbMYRLZ69NFHsXbtWmzevBkdO3aUuxzF8vb2Rvfu3QEACQkJyMnJwXvvvYeFCxfKXJmy7NixA8XFxYiPj294n9FoxObNm/H+++/DYDDAw8NDxgqVKyAgAP369cPRo0flLkXx2GPRSt7e3oiPj8f69esbvX/9+vUYNmyYTFWR2gmCgLlz5+Krr77Cxo0bERMTI3dJqiIIAgwGg9xlKM4tt9yCvXv3Ii8vr+EtISEBkydPRl5eHkNFCwwGAw4ePIjIyEi5S1E8zlhIYN68eZg6dSoSEhKQnJyMRYsW4cyZM5g1a5bcpSlOZWUljh071vDnkydPIi8vDyEhIejUqZOMlSnLnDlz8Omnn2LNmjUICgpqmBHT6XTw8/OTuTplee655zB27FhER0ejoqICK1euxC+//IJ169bJXZriBAUFXdenExAQgNDQUPbvNPHUU09h3Lhx6NSpE4qLi/H6669Dr9cjJSVF7tIUj8FCAvfffz9KS0vx6quvoqCgAHFxcfj+++/RuXNnuUtTnNzcXNx8880NfzavVaakpCAjI0OmqpTHvHX5pptuavT+pUuXYvr06c4vSMGKioowdepUFBQUQKfToX///li3bh1Gjx4td2mkYmfPnsWkSZNQUlKCdu3aYejQocjKyuK/6yLwHAsiIiKSDHssiIiISDIMFkRERCQZBgsiIiKSDIMFERERSYbBgoiIiCTDYEFERESSYbAgIiIiyTBYEBERkWQYLIiIiEgyDBZEREQkGQYLIiIikgyDBREREUnm/wHlyEF5rrbhSwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prostate = pd.read_table(\n",
    "    \"https://static-1300131294.cos.ap-shanghai.myqcloud.com/data/ml-fundamental/parameter-optimization/gradient-descent/prostate.data\"\n",
    ")\n",
    "prostate.drop(prostate.columns[0], axis=1, inplace=True)\n",
    "\n",
    "X = prostate.drop([\"lpsa\", \"train\"], axis=1)\n",
    "y = prostate[\"lpsa\"]\n",
    "\n",
    "regressor = LinearRegressionWithSGD()\n",
    "\n",
    "regressor.fit(X, y)\n",
    "y_pred = regressor.predict(X)\n",
    "\n",
    "print(regressor.__dict__)\n",
    "print(y - y_pred)\n",
    "\n",
    "plt.scatter(y, y_pred)\n",
    "plt.plot([0, 5], [0, 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic regression with gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogisticRegression:\n",
    "    def __init__(self, learning_rate=0.001, n_iters=1000):\n",
    "        self.lr = learning_rate\n",
    "        self.n_iters = n_iters\n",
    "        self.weights = None\n",
    "        self.bias = None\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        n_samples, n_features = X.shape\n",
    "\n",
    "        # init parameters\n",
    "        self.weights = np.zeros(n_features)\n",
    "        self.bias = 0\n",
    "\n",
    "        # gradient descent\n",
    "        for _ in range(self.n_iters):\n",
    "            # approximate y with linear combination of weights and x, plus bias\n",
    "            linear_model = np.dot(X, self.weights) + self.bias\n",
    "            # apply sigmoid function\n",
    "            y_predicted = self._sigmoid(linear_model)\n",
    "\n",
    "            # compute gradients\n",
    "            dw = (1 / n_samples) * np.dot(X.T, (y_predicted - y))\n",
    "            db = (1 / n_samples) * np.sum(y_predicted - y)\n",
    "            # update parameters\n",
    "            self.weights -= self.lr * dw\n",
    "            self.bias -= self.lr * db\n",
    "\n",
    "    def predict(self, X):\n",
    "        linear_model = np.dot(X, self.weights) + self.bias\n",
    "        y_predicted = self._sigmoid(linear_model)\n",
    "        y_predicted_cls = [1 if i > 0.5 else 0 for i in y_predicted]\n",
    "        return np.array(y_predicted_cls)\n",
    "\n",
    "    def _sigmoid(self, x):\n",
    "        return 1 / (1 + np.exp(-x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LR classification perf:\n",
      " [[88  9]\n",
      " [40 16]]\n",
      "LR classification error rate:\n",
      " 0.3202614379084967\n"
     ]
    }
   ],
   "source": [
    "heart = pd.read_csv(\n",
    "    \"https://static-1300131294.cos.ap-shanghai.myqcloud.com/data/ml-fundamental/parameter-optimization/gradient-descent/SA_heart.csv\"\n",
    ")\n",
    "heart.famhist.replace(to_replace=[\"Present\", \"Absent\"], value=[1, 0], inplace=True)\n",
    "heart.drop([\"row.names\"], axis=1, inplace=True)\n",
    "X = heart.iloc[:, :-1]\n",
    "y = heart.iloc[:, -1]\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.33, random_state=42\n",
    ")\n",
    "\n",
    "regressor = LogisticRegression(learning_rate=0.0001, n_iters=1000)\n",
    "\n",
    "regressor.fit(X_train, y_train)\n",
    "y_pred = regressor.predict(X_test)\n",
    "perf = sklearn.metrics.confusion_matrix(y_test, y_pred)\n",
    "print(\"LR classification perf:\\n\", perf)\n",
    "\n",
    "error_rate = np.mean(y_test != y_pred)\n",
    "print(\"LR classification error rate:\\n\", error_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your turn 🚀\n",
    "\n",
    "Modify ```LogisticRegression``` so that the training will use SGD instead of GD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [optional] At the frontier of Machine Learning Research "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<p style=\"text-align: center;\">\n",
       "<iframe src=\"https://www.youtube.com/embed/mdKjMPmcWjY\" width=\"105%\" height=\"700px;\" style=\"border:none;\"></iframe>\n",
       "video. <a href=\"https://www.youtube.com/embed/mdKjMPmcWjY\">[source]</a>\n",
       "</p>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "\n",
    "display(\n",
    "    HTML(\n",
    "        \"\"\"\n",
    "<p style=\"text-align: center;\">\n",
    "<iframe src=\"https://www.youtube.com/embed/mdKjMPmcWjY\" width=\"105%\" height=\"700px;\" style=\"border:none;\"></iframe>\n",
    "video. <a href=\"https://www.youtube.com/embed/mdKjMPmcWjY\">[source]</a>\n",
    "</p>\n",
    "\"\"\"\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bibliography\n",
    "\n",
    "- [Gradient Descent, Step-by-Step - StatQuest](https://www.youtube.com/watch?v=sDv4f4s2SB8)\n",
    "- [Stochastic Gradient Descent, Clearly Explained!!! - StatQuest](https://www.youtube.com/watch?v=vMh0zPT0tLI) \n",
    "- http://43.142.12.204:12345/05-ML_04-Under-the-Hood.html\n",
    "- http://43.142.12.204:9999/GradientDescentAnimation.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myconda1",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}