monopolists_problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monopolist's Optimization Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bob runs a widget factory. He is the sole producer of his particular type of widget.\n",
    "The price, $p$, he receives in the market for a pound of his widgets is determined the quantity, $q$, of pounds of widgets he puts on the market.\n",
    "\n",
    "He faces the quadratic demand curve $p(q) = 100 - q^{2}$.\n",
    "\n",
    "Which means, for example:\n",
    "\n",
    "* if he sells 1 pound of widget he receives \\\\$99 per pound: $p(1) = 100 - 1^{2} = 99$\n",
    "* if he sells 2 pounds of widgets he receives \\\\$96 per pound: $p(2) = 100 - 2^{2} = 96$\n",
    "* if he sells 3 pounds of widgets he receives \\\\$91 per pound: $p(3) = 100 - 3^{2} = 91$\n",
    "\n",
    "Bob owns the factory and equipment in it; he cannot sell it. Bob's cost per pound of widgets is 10. This gives $C(q)=10q$ as his total costs.\n",
    "\n",
    "Profit is given by $price\\ per\\ unit * units\\ sold - total\\ costs$. In our notation, Bob's profit, $\\pi$ can be written as $\\pi(q) = 90q - q^{3}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization Techniques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll consider minimizing a function that is concave over its relevant range. In order to maximize a convex function, we can simply take its negative (which is concave) and minimize it.\n",
    "\n",
    "For simplicity, we'll only consider function of one real-valued variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gradient Descent is a simple method of finding a minimum that can be used on a very wide variety of functions. Gradient Descent along with methods based on it is commonly used to maximize functions that have no second derivative or where the second derivative of the function is difficult to compute or even approximate. Gradient Descent based methods are used in order to minimize the loss function in many machine learning models, particularly deep learning models. \n",
    "\n",
    "The idea is that once we pick an initial guess, we iteratively can follow the slope (i.e. go the opposite direction of the slope) until we reach a minimum. This is often compared to a ball rolling down a hill, but without the details to make it physically realistic. (There are other optimization techniques that draw inspiration for the physical metaphor by adding in analogs to friction and gravity.)\n",
    "\n",
    "The iterative equation for Gradient Descent, for a real-valued function of one variable, is given by:\n",
    "$$x_{n+1} = x_{n} - \\gamma \\ f'(x_{n})$$\n",
    "Where $\\gamma$, called the learning rate, controls the rate at which descent happens. The learning rate is an example of a hyperparameter. A hyperparameter is not estimated by the model and must be chosen before the model is run. In cases where the function is very steep, a low value should be chosen for the learning rate. Conversely, a high value for the learning rate should be chosen for very gradual functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton-Raphson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Newton-Raphson is a root finding algorithm. However, since finding a critical point of a function (minima, maxima, point of inflection) is equivalent to finding a root of it derivative, we can use Newton-Raphson for minimization. To use Newton-Raphson for minimization, we require that the function has a second derivative over its domain. Newton-Raphson and methods based on it are commonly used for scientific and statistical applications, because in cases where the function is well understood, Newton-Raphson has better properties and converges more rapidly than gradient descent.\n",
    "\n",
    "For a real-valued function of one variable, $f$, Newton-Raphson find roots with the iterative equation:\n",
    "$$x_{n+1} = x_{n} - \\frac{f(x_{n})}{f'(x_{n})}$$\n",
    "\n",
    "Since minimizing a well-behaved function is the same as finding the roots of its derivative, we can use Newton-Raphson to minimize f with the iterative equation:\n",
    "$$x_{n+1} = x_{n} - \\frac{f'(x_{n})}{f''(x_{n})}$$\n",
    "\n",
    "Rewriting slightly, we get:  \n",
    "$$x_{n+1} = x_{n} - \\gamma(x_{n}){f'(x_{n})}$$\n",
    "where $\\gamma(x_{n}) = \\frac{1}{f''(x_{n})}$.\n",
    "\n",
    "From this last equation, we can directly see the relationship to Gradient Descent: Newton-Raphson minimization is like Gradient Descent, but with a learning rate that is determined by the curvature of the function. Returning to the metaphor of a ball rolling down a hill, Gradient Descent has a constant acceleration, while Newton-Raphson minimization has acceleration determined by the curvature of the hill.\n",
    "\n",
    "Another important difference is that Newton-Raphson has no hyperparameters. Which contributes to the simplicity of using the method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instructions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Show the expression for profit, $\\pi(q) = 90q - q^{3}$, is correct.\n",
    "2. Mathematically solve for the optimal quantity of widgets Bob should produce to maximize his profit and his maximum profit.  \n",
    "(If you need a refresher, see: https://mjo.osborne.economics.utoronto.ca/index.php/tutorial/index/1/nen/t)\n",
    "3. Use Gradient Descent to numerically solve Bob's problem. Write the algorithm yourself; don't use a package for it!\n",
    "4. Use Newton-Raphson to numerically solve Bob's problem. Write the algorithm yourself; don't use a package for it!\n",
    "4. Compare the convergence speed of Gradient Descent and Newton's method.\n",
    "5. Use Scipy's building optimization function to solve Bob's problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Profits are given by $price\\ per\\ unit * units\\ sold - total\\ costs$.  \n",
    "So:  \n",
    "$\n",
    "\\quad \\pi(q) = p(q) q - C(q)\\\\\n",
    "\\quad \\pi(q) = (100 - q^{2}) q - 10q\\\\\n",
    "\\quad \\pi(q) = 90q - q^{3}\\\\\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Profits are given by $\\pi(q) = 90q - q^{3}$.  \n",
    "The derivative is $\\pi'(q) = 90 - 3q^{2}$, by the power rule.  \n",
    "We know that a critical point occurs at $\\pi'(q) = 0$.  \n",
    "Plugging in the expression for $\\pi'(q)$ , we get $90 - 3q^{2} = 0$.  \n",
    "Solving, we get $q = \\sqrt{30}$ and $q = -\\sqrt{30}$.  \n",
    "However, the range of q is implicitly nonnegative numbers, so the negative root doesn't make sense.  \n",
    "So, the optimal quantity is $q = \\sqrt{30}$ and the associated profit is quantity is $pi(\\sqrt{30}) = 60\\sqrt{30}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to verify that we have the global maximum (as opposed to a local maximum or a minamum).  \n",
    "We can get the intuition for this by graphing.\n",
    "\n",
    "We find the roots of $\\pi(q)$, so that we can get a reasonable range.  \n",
    "The roots occur at $\\pi(q) = 0$.  \n",
    "The relevant roots are $q=0$ and $q=\\sqrt{30}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pi(q):\n",
    "    \"\"\"Profit as a function of quantity.\"\"\"\n",
    "    return 90 * q - q ** 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEuCAYAAABLSP/KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3xV9f3H8dcniySsEAiBQNgr7BFGRVGW4sS9Z7Vo1aq11qrVOlqttdZaf8W96x6oKMMJKiggIHuGsBIggxESIPvz++Me9BITSEjuPXd8no/HfeTeM9+5F+4n53vO+X5FVTHGGGPqI8LtAMYYY4KfFRNjjDH1ZsXEGGNMvVkxMcYYU29WTIwxxtSbFRNjjDH1ZsXEGGNMvVkxMcYYU29WTEzQEJGeIrJERApF5CYRWSkiJzTQtv8uIrfUctkFItKnIfYbzBry/Xe2V+vPwFnePocAYsXE+JSIbBKRAyJSJCI5IvKyiDQ5ys3dDsxS1aaq+oSq9lHV2VX2Ne4oMiYBlwPP1HKVR4EH6rqfhiAiV4rIchHZLyI7RORJEWnup30f8v421PvvrFvXzwBc/BzML1kxMf5wuqo2AQYD6cDdVRcQkahabKcjsLKBswFcCUxX1QO1XH4qMFpE2vggS41E5A/AP4A/As2BEUAn4DMRifZnFh+4krp9BuDS52CqZ8XE+I2qZgMzgL7w01+yfxKRZcA+EYkSkTQRmS0ie5xmlDOcZb8CRgP/dY5yenj/JSwi/wM6AB878293tp3tNIutFZGxNUQ7Gfj64AsRaSQie0WkwtlWkfO8WETGqWoxsAg4qbqNOft9r8q0/4jIE4dZ/rA5RaQZcD/wO1WdqaplqroJOB/oAlzsLKci0s1rvZdF5G9er+8QkQ3OvlaJyFle8zaJyG0iskxECkTkbRGJPcz7e7j3/08i8n6V3+EJEflPde8BVT4DZ/kIEblTRHJFJE9ErheRUucohiN9DsbPVNUe9vDZA9gEjHOep+I5svir17wlzvQ4IBrIAO4CYoAxQCHQ01l+NnBNdduuZl89ga1AivO6E9C1hox5wNAq064CPvN6vQ4Y5fX6CeCxGrbXEdgPNHVeRwLbgRHVLFurnMAEoByIqmbeK8DrznMFunnNexn4m9fr84AUPH9IXgDsA9p6vX8LnPmJwGrgutq839W8/22dbSc4r6OAXGBIHT6D+4C5QBs8R2LzgB1Vlqnxc7CHfx92ZGL84UMR2QPMwfPX50Ne855Q1a3qad4YATQBHlbVUlX9CvgEuOgo9lkBNAJ6i0i0qm5S1Q01LJuAp2h56wcsBxCRxkDXg68dhc56v6Cqm4HFwMG/+scA+1V1Xj1ytgLyVbW8mnnbgaTqslST7V1V3aaqlar6NrAeGOa1yBPO/F3Ax8DA2my3mv1sB77BU7zAUwzzVXVRDasc8hk4Rx+3Apeo6g5VLQCmcehnAIf5HIx/WTEx/nCmqiaoakdVvV4PbRff6vU8BdiqqpVe0zYD7eq6Q1XNAG7B89dtroi8JSIpNSy+G2haZdpPxcR5vl1Vd3vNbwrsOUyEN/i5CF7svK5PznygVQ3nlto6849IRC53rojb4xT4vngK1UE7vJ7vx1Pcj9YrwKXO80uB/x1m2aqfwVhgtXqa8g5qyS+LyZE+B+MnVkyM27wH1NkGpIqI97/LDkD2UWwLVX1DVY/F0+ykeE5eV2cZ0KPKNO9iMoBffomlAUsPk+Vd4AQRaY/nCKXaYlKHnN8DJcDZ3hOdK+NOxtMECJ4CEO+1SBuvZTsCzwE3Ai1VNQFYAchhfo9DotZx/odAfxHpC5wGvH6Ydat+Bq3wNH0BP12gcYaT19uRPgfjJ1ZMTCCZj+fL8HYRiXbuYTgdeKuW6+fgORl98J6UMSLSCCgGDgCVNaw3HTj+4AuniSUJWOVM6oXnnMnB+bHAEODzmoKoah6eL/iXgI2qurq65Wqb02nmuR/4PxGZ4Lw/nYB38ByVHPyiXgJcLCKRIjLB+/cCGuP5ws9z9n0VzsUQtfTT+1ub+eo5Qf4enkK6QFW3HGbdQz4DYA1wjIh0FpEWwFNAZ7yKem0+B+M/VkxMwFDVUjzF42Q8X5BPAper6ppabuLvwN1O880FwMPOdnYArYE7a1jvVeAUEYlzXvcDNng1x2UD54vIcOf16cBsVd12hDxvAOOoclQiIjNE5C7nZaPa5lTVR/BcnPAonnMFG/EchYxT1X3OYjc7+fYAl+A5Oji4/irgX3iOcnKc33PuEX4Hbz+9vyJyWy3nv+Ls53BNXFDlM1DVL/Ac3S3F80fGcjxF1vvS8Np+DsYPRNWG7TVGRB4CclX18VosOx+4WlWrNrn4lXNk8QAw8gh/9btGRDrgOcpoo6p7j7BsjZ+BiJwITFbV7l7TAuJzMB5WTIwJYiJyGVCmqrVtCvQb59zXY0AzVf11Pbd1M3C8qp59xIWNK2pz17ExJkCp6pGaj1zhXE6dg+dqvAkNsMk0fnny3QQQOzIxxhhTb3YC3hhjTL1ZMTHGGFNvVkyMMcbUmxUTY4wx9WbFxBhjTL3ZpcHGmKMiIsnAB0AZnt6PL3F6CzZhyC4NNsYcFRGJBFRVK0XkSqC9qv7tCKuZEGXNXMZ4EZG/i8gttVx2gYj08XWmQKWqFV7DBTTFN0MqmyBhxcQYh9Nb8OXAM7Vc5VE8fWM1xL7TROQrZ7jcDO/hdJ35iSLygYjsE5HNInJxQ+y3vkRkoNNH1o14BgQzYcqKiTE/uxKYXmXwrsOZCowWkTZHXPIwnLE6PsIzqmQiMAl4TUS8x/eYDJQCyXh6A34qEI6KVHWJqg4H7qHmXplNGLBiYsKCiDQSkb0iUiEiRc6jQkSKRWScs9jJeIYVPrhOhIjcKSK5IpInIteLSKlzBHNwvI5FwEn1jNcLzyiT/3aajr7C0zX8ZU6OxsA5wD2qWqSqc/AUsstqm/VoHW67IhLjtWgBnrFoTJiyq7lMWFDVEqfn2YtU9UQAEVkHXKOq3ziL9QPWeq32F2A80B/PoFWfArucga8OWo1nJMZDiMgnwLE1xJmjqqcdIbLw88BVPYByVV3nNX8phw4mVZusR6PG7YrIMBF5FM+VXMVAvXoGNsHNiokJJz8Nxev8td+VQ4fjTcAz6NTB8ye3Av1VdYczbRowqso2C/GMwX6IWhQLb2uBXOCPIvJvYDSeQjHLmd8EqDoWSAHOmOl1yFonR9quqi6o7z5M6LBmLhNOvMd17wdsV9XdXvN343xBA2OB1aq6yWt+S345FnxTPKMaHjVVLQPOBE7FM9riH/AMx5vlLFIENKuyWjOcwlfbrCIyW0S0hsecaqLV9j0wxoqJCSvexWQAv/xSXIanSQmgFc5Y6fDTSfIz+OWYGml4mpwO4QzNW1TDY0bV5VV1maoer6otVfUkPGOpL3BmrwOiRKS71yoD+PlS3FplVdUTVFVqeFTXJFfb98AYKyYmPDhNNknAKmdSLzxf0t6m8/N5iDXAMSLSWURaAE8BnfEqQCISCwwBPq+6P1U9WVWb1PA4uZp8/UUkVkTinfHT2wIvO9vaB0wBHhCRxiIyEpjIz+OqHzHrUfLVdk0IsmJiwkU/YIPXZb/ZwPkiMtxrmVeBU0QkTlW/AN7Fc9QxH88XaCWH3ph3OjBbVbc1QL7LgO14zp2MBcaraonX/OuBOGf+m8BvVXUlQC2z1pmvtmtCk3WnYowXEXkIyFXVx6tMPxGYrKrdvabNB65W1YBq9qkuayBv14QGu5rLGC+qelcNs9Ko0rzj3KwXiH6RNcC3a0KANXMZUztpBM+JZ19lDab3wPiZNXMZY4ypNzsyMcYYU29WTIwxxtRb2J6Ab9WqlXbq1MntGMYYE1QWLVqUr6q/6EA0bItJp06dWLhwodsxjDEmqIjI5uqmWzOXMcaYerNiYowxpt6smBhjjKk3KybGGGPqzYqJMcaYerNiYowxpt6smBhjjKm3sL3PxJhQVVmp7NhbzPaCA+zaV8bu/aUU7C+jrLISVc/8iAihWWwUTWOjaRYXRXKzWNq3iKd5XLTb8U2QsmJiTBDbWVTCsuwClmcVsDy7gI35+9iyaz+l5ZVHtb1msVF0SWpCn5Rm9ElpTv/2zUlr24zICGng5CbUWDExJogUHCjj+w07mZuRz9yMfDLz9wEgAp1bNaZ76yaM6dWaji3jaZcQR2LjGBLiYmgeH02jqAgiRIgQKK9UCovLKSwuo+BAGTsKisnafYCtu/ezLqeQqUu38fr8LQA0j4tmeOdERnZrxfjeyaQkxLn5FpgAFbZd0Kenp6t1p2KCwc6iEj5dmcOMFdv5bsNOKiqV+JhIRnRpyfDOiQxITaBPSjOaxjZcE5WqkrX7AIu37Oa7jJ18l5nP1l2eEY8HpCZwct82TByYQtvmVljCjYgsUtX0X0y3YmJM4Ckpr+CzlTm8s3DrTwWkU8t4Tu7XltE9WzMwNYGYKP9eP5OZV8TMlTuYuWIHy7IKiBA4oWdrLhiaypherYmOtOt5woEVkyqsmJhAlJlXxJsLtvD+4mx27SulXUIcZw5K4ZR+bendthkigXHuYsvO/byzcCvvLtpKzt4S2iXEcfWxnblgaCqNG1nreSizYlKFFRMTSBZt3sUzX2fy+eocIkUY3zuZC4d14LhurYgI4JPf5RWVfLUml+e/3ciCTbtoHhfNlcd04prjOjdos5sJHFZMqrBiYtymqsxem8fkWRks3Lyb5nHRXPGrjlz6q460bhrrdrw6W7R5N898vYHPVuXQIj6aG0Z349IRHYmNjnQ7mmlAVkyqsGJi3DQvcyf//HQtizbvpl1CHNcc15nz00OjiWhZ1h7++elavl2fT7uEOO45rTcn9UkOmCY6Uz9WTKqwYmLcsCK7gH/MXMO36/Np0yyW343txvnpqSF58vq7jHwe+GQVa3YUckLPJO4/ow8dWzZ2O5apJysmVVgxMf60s6iERz9bx1s/bCEhLnyagMorKnn5u008/sV6SisquXlsd64d1YWoECye4cKKSRVWTIw/lFdU8tq8zTz2+Tr2lVZw5TGduGls97DrtiRnbzH3f7yS6ct3MKhDAo+dP5DOrewoJRhZManCionxteVZBdz+/jJWb9/Lsd1ace/pveme3NTtWK6aunQb93y4gpLyCu46JY3LRnS0cylBpqZiEvxn+4wJMMVlFfz7i3U8900mrZo04ulLB3NSnzb2pQmcMSCF4Z0Tuf29Zfzlo5V8v2Enj5zb3y4jDgF2ZGJMA1qwcRd/en8ZG/P3ceHQVO48JS3smrRqQ1V5/tuNPDxzDR0S43nyksGktW3mdixTCzUdmQTkWTARiRWRBSKyVERWisj9zvTOIjJfRDJE5G0RiXGmN3JeZzjzO7mZ34Sf0vJKHp6xhgue/Z7yykpev2Y4D5/T3wpJDUSE34zqwpu/GcG+knLOnDyXD37McjuWqYeALCZACTBGVQcAA4EJIjIC+Afwb1XtBuwGrnaWvxrY7Uz/t7OcMX6RmVfEOU99x9Nfb+CC9FRm3jyKkd1auR0rKAzrnMi0m45jUIcEfv/2Uh79dC2VleHZWhLsArKYqEeR8zLaeSgwBnjPmf4KcKbzfKLzGmf+WLEGauNjqspbC7Zw6hNz2Lp7P09fOpiHz+kfEjce+lNS00b87+rhXJCeyn9nZfC7t36kuKzC7VimjgL2X72IRAKLgG7AZGADsEdVy51FsoB2zvN2wFYAVS0XkQKgJZBfZZuTgEkAHTp08PWvYELYvpJy7pyynKlLt3FM15Y8dv5A2jQPvi5QAkV0ZAQPn9OPLkmNeXjmGrJ2H+DFK9Jp2aSR29FMLQXkkQmAqlao6kCgPTAM6NUA23xWVdNVNT0pKaneGU142pBXxJmT5/LJsm3cdmIPXrt6uBWSBiAiXHt8V566ZAhrtu/l/Ge+Z3vBAbdjmVoK2GJykKruAWYBvwISROTg0VR7INt5ng2kAjjzmwM7/RzVhIEZy7cz8b9z2bmvlFd/PZwbx3QP6F59g9GEvm149dfDyNlbwrlPfc9GZzRJE9gCspiISJKIJDjP44DxwGo8ReVcZ7ErgI+c51Od1zjzv9JwvebZ+ERFpfL3Gav57euL6dq6CZ/87liO7W4n2X1leJeWvDVpBAfKKjjv6e9YtW2v25HMEQRkMQHaArNEZBnwA/C5qn4C/Am4VUQy8JwTecFZ/gWgpTP9VuAOFzKbEFVUUs6kVxfyzNeZXDy8A+9cO8LGQfeDvu2a8861vyI6MoKLn59nBSXA2U2LxhxG9p4DXP3yD6zPLeLe03tz+a86uR0p7GzeuY8Ln51HSXklb/5mBD3bhHeXNG4LqpsWjQkEP27ZzcT/ziV79wFevHKoFRKXdGzZmDd+M4LoSOGS5+eRkVvodiRTDSsmxlTjk2XbuODZecTFRDDl+mM4vodd/eemzq08BUVEuOi5+Wyyk/IBx4qJMVW8NHcjN77xI/3bNefD60eGfU+/gaJrUhPe/M1wyisqufzFBeQWFrsdyXixYmKMQ1V5ZOYa7v94FSf2Tua1a4bbTXMBplvrprx01TDyCku48sUfKCwuczuScVgxMQbPIFZ/en8ZT87ewEXDOvDUpUNCfhTEYDUwNYGnLh3MupxCJr26iJJy63olEFgxMWHvQGkF1722iHcWZnHT2O48dFZfIu1GxIB2Qs/W/PO8/nyfuZM/vLPUOocMAAHbN5cx/lBUUs7VL//Agk27+OvEPlxmV2wFjbMGtSdnbwkPz1hDl6Qm3Dq+h9uRwpoVExO29haXceWLC1iaVcB/LhzEGQNS3I5k6ujaUV3YkFvEE1+up2tSYyYObHfklYxPWDExYWnP/lIuf3EBq7fvZfLFg5nQt43bkcxREBEePKsfm3fu54/vLaNDYjyDOrRwO1ZYsnMmJuzsLCrhoufms2ZHIc9cNsQKSZCLiYrg6cuG0KZZLL95dRHb9lhPw26wYmLCSm5hMRc+O4+N+UU8f3k6Y3olux3JNIDExjG8cEU6B0rLuf71xXaFlwusmJiwkVdYwkXPziN7zwFeunIYo+yu9pDSPbkp/zxvAEu27uFvn6x2O07YsWJiwsKufaVc+vx8tu0p5qUrh/Krri3djmR84JR+bZk0qgv/m7eZD37McjtOWLFiYkJewYEyLnthPht37uP5K9IZ3sUKSSi7/aSeDO+cyJ1TlrN6u3Vb7y9WTExIKywu44oXF7Aux3OyfWQ3G9Aq1EVFRvB/Fw+iWWw01722yLpc8RMrJiZk7S8t59cv/8CK7AImXzyY0T1bux3J+EnrprFMvmQwWbsPcM+HK9yOExasmJiQVFJewW9eXciizbt5/MKBnNjHLv8NN0M7JXLz2O58uGQbUxbb+RNfs2JiQk5FpfL7t5cwN2Mn/zx3AKf1tzvbw9UNo7sxrHMi93y4wsZA8TErJiakqCr3fLSC6ct3cPepaZwzpL3bkYyLIiOExy8YSFRkBDe/9SOl5ZVuRwpZVkxMSPn3F+t5Y/4Wrju+K9cc18XtOCYApCTE8Y9z+rE0q4DHPl/ndpyQZcXEhIxXv9/EE1+u5/z09vxpQk+345gAMqFvWy4a1oFnvtnAgo273I4TkgKymIhIqojMEpFVIrJSRG52pt8nItkissR5nOK1zp0ikiEia0XkJPfSGzd8vHQb905dybi0ZB46qx8iNh6JOdTdp6aR2iKeP763lP2l5W7HCTkBWUyAcuAPqtobGAHcICK9nXn/VtWBzmM6gDPvQqAPMAF4UkRsmLww8V1GPre+s4ShHRP578WDiIoM1H/Wxk2NG0Xxz3P7s3nnfv4xY43bcUJOQP6vU9XtqrrYeV4IrAYON1DBROAtVS1R1Y1ABjDM90mN29bnFHLta4vo3Koxz12RbkPtmsMa3qUlV43sxCvfb+a7jHy344SUgCwm3kSkEzAImO9MulFElonIiyJycOCCdsBWr9WyOHzxMSEgt7CYK1/6gdjoSF68cijN46LdjmSCwO0n9aJzq8b88b1lFJVYc1dDCehiIiJNgPeBW1R1L/AU0BUYCGwH/lXH7U0SkYUisjAvL6/B8xr/2V9azjWvLGTXvlJevGIo7VvEux3JBIm4mEgePa8/2woO8OA06124oQRsMRGRaDyF5HVVnQKgqjmqWqGqlcBz/NyUlQ2keq3e3pl2CFV9VlXTVTU9Kcm6Hw9WFZXKzW8tYUV2Af930SD6tW/udiQTZIZ0TOSaYzvz5oItdnVXAwnIYiKeS3FeAFar6mNe09t6LXYWcLDTnanAhSLSSEQ6A92BBf7Ka/zrwWmr+XxVDn85rTfjetvgVubo/H58D9olxHHnlGU2mFYDCMhiAowELgPGVLkM+BERWS4iy4DRwO8BVHUl8A6wCpgJ3KCq9q8jBL08dyMvzt3Ir0d25sqRnd2OY4JYfEwUfzurLxvy9vH07Ey34wS9KLcDVEdV5wDV3Sgw/TDrPAg86LNQxnWz1ubywCerGN87mT+fmuZ2HBMCRvdszekDUpg8K4NT+7elW+smbkcKWoF6ZGLMITJyi7jpjR/p2aYZ/7lwIJERdlOiaRh/Oa03sdER3PXBcior1e04QcuKiQl4BfvLmPTqQmKiInju8iHExwTkAbUJUklNG3HXKWks2LiLdxdtPfIKplpWTExAK6+o5MY3F7N1936evmyIXQJsfOL89FSGdUrk4Rlr2LO/1O04QcmKiQlof5+xhm/X5/O3M/sytFOi23FMiIqIEO6f2IeCA2XWs/BRsmJiAtY7C7fywpyNXHlMJy4Y2sHtOCbEpbVtxmUjOvLavM2s2rbX7ThBx4qJCUiLNu/i7g9WcGy3VtxtV24ZP7l1fE8S4mO4b+pKVO1kfF1YMTEBZ3vBAa7932LaJsRaL8DGr5rHR3P7ST1ZsGkXU5ducztOULH/pSaglJZXcv3rizlQWs7zl6eTEB/jdiQTZs5PT6V/++Y8OG21dQRZB1ZMTED527RV/LhlD4+cO4DuyU3djmPCUESEcP8ZfcgtLGHyrAy34wQNKyYmYExZnMWr329m0qgunNq/7ZFXMMZHBnVowdmD2vHCnI1k7zngdpygYMXEBIRV2/Zy1wfLGdElkdtPsvHbjftuO6knAvxzpo3KWBtWTIzrCvaXcd1ri2geF83/XTTYTribgJCSEMc1x3XmwyXbWJa1x+04Ac/+1xpXVVYqt76zhO0FB3jykiEkNW3kdiRjfnLd8V1p2TiGB6ettkuFj8CKiXHV5FkZfLkml3tO682Qji2OvIIxftQ0Nppbxvdg/sZdfL4qx+04Ac2KiXHNN+vyeOyLdZw1qB2XjejodhxjqnXR0FS6JjXm4RlrKKuodDtOwLJiYlyxo6CYW95eQo/WTXnorH54Btc0JvBERUZw1ylpZObv4435W9yOE7CsmBi/K6+o5KY3f6S4rILJlwwmLibS7UjGHNaYXq0Z0SWR//tqPftL7UbG6lgxMX73+BfrWbBpFw+e1ddGtjNBQUT440m9yC8q5aW5m9yOE5CsmBi/+nZ9HpNnZ3B+envOGtTe7TjG1NqQji0Yl9aaZ77eQMH+MrfjBBwrJsZvcvYWc8tbS+jeugn3n9HX7TjG1NkfTuzJ3uJynv12g9tRAo4VE+MXFZXKzW/9yP7SCiZfbOdJTHBKa9uM0wek8NLcTeQVlrgdJ6AEZDERkVQRmSUiq0RkpYjc7ExPFJHPRWS987OFM11E5AkRyRCRZSIy2N3fwFT1ny/XMy9zF389s6914GiC2u/HdaekvJInZ1snkN4CspgA5cAfVLU3MAK4QUR6A3cAX6pqd+BL5zXAyUB35zEJeMr/kU1N5mbk839freecwe05d4idJzHBrUtSE84b0p7X522xTiC9BGQxUdXtqrrYeV4IrAbaAROBV5zFXgHOdJ5PBF5Vj3lAgohYt7MBILewmJvfWkLXpCb89cw+bscxpkHcNLY7AE98sd7lJIGjzsVERBqLiN8avEWkEzAImA8kq+p2Z9YOINl53g7Y6rValjPNuKiyUrnt3WUUFpcx+eLBxMdEuR3JmAaRkhDHxcM78N7iLLbu2u92nIBwxGIiIhEicrGITBORXGANsN05n/FPEenmq3Ai0gR4H7hFVfd6z1NPr2t16nlNRCaJyEIRWZiXl9eASU11XvpuE9+sy+Pu03rTs42dJzGh5bcndCUyQuzciaM2RyazgK7AnUAbVU1V1dbAscA84B8icmlDBxORaDyF5HVVneJMzjnYfOX8zHWmZwOpXqu3d6YdQlWfVdV0VU1PSkpq6MjGy6pte/nHjDWMS0vm0uEd3I5jTINLbhbLhUNTeW9RFlm77eikNsVknKr+FchU1Z96OVPVXar6vqqeA7zdkKHE01HTC8BqVX3Ma9ZU4Arn+RXAR17TL3eu6hoBFHg1hxk/O1BawU1v/UhCfDSPnNvf+t0yIeu647sC8PTXdt/JEYuJqh681fNHEektIj81fItIjyrLNJSRwGXAGBFZ4jxOAR4GxovIemCc8xpgOpAJZADPAdc3cB5TBw9OX0VGbhH/On8AiY1j3I5jjM+kJMRxXnoq7/yQxfaC8L6yqy5nRBOByUCaiOThucJqMNDg50xUdQ5Q05+zY6tZXoEbGjqHqbvPV+Xw2rwtTBrVheO6W1OiCX2/Pb4r7/ywlWe+zuS+M8L3isW6FJMtqjoaQETaA2lAeJdic4icvcXc/t5S+qQ047YTbRx3Ex5SE+M5e3A73liwhetP6ErrZrFuR3JFXS4Nbi4ix4hIc1XNUtXPnSMIY6isVP7wzlIOlFXwnwsHERMVkLcwGeMTN4zuRkWl8sw3mW5HcU1d/sc3AW4DFojIJhGZISL/9FEuE2RemLORORn53Ht6H+tW3oSdji0bM3FgCq/P38zOovDss6suxWSEqp6tqj2BXsCfgWW+iWWCyYrsAh75dA0n9UnmwqGpR17BmBB0/QldKS6r5JXvNrkdxRW1uWlRAFT1p+M3VS1W1cWq+j/vZUz4KS6r4NZ3lpDYOIaHz7bLgE346ta6KeN7J/PK95vZVxJ+ozHW6qZFEfmdiBxy55mIxIjIGBF5hZ/v/TBh5l+frWVdThGPnDuAFr3hf/0AABrbSURBVHYZsAlz1x3flYIDZbz9w9YjLxxialNMJgAVwJsiss3pRiUTWA9cBDyuqi/7MKMJUPMyd/L8nI1cOqIDx/ewy4CNGdKxBcM6JfLCnI2UVVQeeYUQUpubFotV9UlVHQl0xHOfx2BV7aiqv1HVH32e0gScwuIy/vDOUjomxnPXKWluxzEmYFx7fBey9xzg46Xb3I7iV3W6flNVy5zu4ff4KpAJDn/9ZBXbCw7w2AUDrTdgY7yM7tmaHslNeObrTDz3U4cHuxnA1NlnK3fwzsIsrj+hG4M7tHA7jjEBJSJCuHZUV9bmFDJrbe6RVwgRVkxMneQXlXDnlOX0SWn20wBBxphDnTEwhZTmsTw9O3xuYqx1MRGR80SkqfP8bhGZYmOthxdV5a4pyyksKeffFwy0u9yNqUF0ZARXH9eFBZt28eOW3W7H8Yu6fBvco6qFInIsnh57X8DGWg8r7y/O5rNVOfzxxJ70SLbBrow5nAuHptI0NooX5mx0O4pf1KWYVDg/TwWeVdVpgN1YECaydu/nvqkrGd45kauP7ex2HGMCXuNGUVw4NJUZK3awbU/o94lbl2KSLSLPABcA00WkUR3XN0GqslL547uennMePW8AERF2l7sxtXHFMZ1QVV75fpPbUXyuLsXgfOBT4CTn0uBE4I8+SWUCyusLtvB95k7uPjWN1MR4t+MYEzTat4hnQt82vDl/C/tLQ7uLlboUk3tVdYqqrgdwhsX9xUBVJrRs3bWfv09fzageSVxgnTgaU2dXH9uZvcXlvL8oy+0oPlWXYjK+mmknN1QQE3gqK5Xb31tGhAgPn93POnE05igM7tCCAakJvDR3E5WVoXsTY216Df6tiCwHeorIMq/HRqwL+pDm3byVkhDndhxjgpKI8OuRncjM38fsdaF7E2NtjkzeAE4Hpjo/Dz6GqOqlPsxmXHSweeu47q2secuYejqlX1vaNIsN6cuEa9PRY4GqblLVi1R1s9djlz8CGv+rrFT+9L7TvHWOjVFiTH1FR0Zw+TEdmZuxkzU79rodxydq08w1x/lZKCJ7vR6FIuKTd0VEXhSRXBFZ4TXtPhHJFpElzuMUr3l3ikiGiKwVkZN8kSmcvLFgC99t2MmfT02jnTVvGdMgLh7WgdjoCF75brPbUXyiNs1cB4/L7lHVZl6PpqrazEe5XsYzjkpV/1bVgc5jOoCI9AYuBPo46zwpIpE+yhXyDjZvHdutlQ3Ba0wDSoiP4YwBKXz4YzZ7i8vcjtPgalNMBotICnCViLQQkUTvhy9Cqeo3QG2b0SYCb6lqiapuBDKAYb7IFepUlTumLENEePgcu3rLmIZ22YhOHCirCMnLhGtTTJ4BvgR6AYuBRV6Phb6LVq0bnSvJXhSRg32ftwO8x8jMcqaZOnpjwRbmZuzkrlPSaN/Cbk40pqH1a9+cgakJ/G/e5pAb66Q2J+CfUNU04EVV7Vzl0cUPGQ96CugKDAS2A/+q6wZEZJKILBSRhXl5eQ2dL6ht3bWfh6Z5mrcuGmbNW8b4ymUjOpKZt4/vNux0O0qDqvVNi6r6WxEZICI3Oo/+vgxWzf5zVLVCVSuB5/i5KSsb8P72a+9Mq24bz6pquqqmJyXZmOUHqSp3TlkOYM1bxvjYqf3b0iI+mv99H1on4usynslNwOtAa+fxuoj8zlfBqtl/W6+XZwEHr/SaClwoIo1EpDPQHVjgr1yh4N1FWczJyOcOa94yxudioyM5f2gqn6/OYXtB6PQmXJfuVK4BhqvqX1T1L8AI4De+CCUibwLf47nrPktErgYeEZHlIrIMGA38HkBVVwLvAKuAmcANqlpRw6ZNFbmFxfztk1UM65TIJcM6uB3HmLBw6fCOVKry5vwtbkdpMFF1WFb4eUwTnOc+aQ9R1YuqmfzCYZZ/EHjQF1lC3X1TV1JcXsnfz+lnXcsb4yepifGM7tmaNxZs5cYx3UNi1NK6/AYvAfOdmwfvA+ZxmC94E/g+XbmD6ct3cPPY7nRNauJ2HGPCymUjOpJfVMKnK3e4HaVB1KqYiOeM7LvAVXju/9gFXKWqj/swm/GhvcVl/OWjFfRq05RJo/x5UZ4xBuD4HkmkJsaFzIn4WjVzqaqKyHRV7YfnXhMT5P4+fQ15hSU8d3k60ZHBf4htTLCJiBAuHtaRf8xcQ0ZuId1aN3U7Ur3U5VtksYgM9VkS4zfzMnfy5oItXH1sZ/q3T3A7jjFh69wh7YmKEN7+YeuRFw5wdSkmw4F5IrLBuQv94JVVJogUl1Vwx/vL6JAYz63je7odx5iwltS0EePSknl/cTYl5cF9EWpdruay3nhDwONfrGfTzv28fs1w4mKsP0xj3HbhsFRmrtzB56tyOK1/ittxjlptuqCPFZFbgD/i6ZU323tcE58nNA1mRXYBz32byfnp7RnZrZXbcYwxwHHdk2iXEMdbC4K7qas2zVyvAOnAcjxjvte5TyzjvvKKSv70/jISG8fw51N6ux3HGOOIjBDOT09lTkY+W3budzvOUatNMemtqpeq6jPAucBxPs5kfOC5bzeyctteHjijD83jo92OY4zxcl56eyIE3l4YvHfE16aY/DSKi6qW+zCL8ZGN+ft4/It1nNQnmZP7tT3yCsYYv0pJiOP4Hkm8uzCL8opKt+McldoUkwHeQ/UC/X09bK9pOJ4egZcRExXBAxP7uh3HGFODC4d1ILewhFlrg3N4jNqMZxJZZajeKD8M22sayHuLspiXuYs7T04juVms23GMMTUY06s1SU0b8daC4GzqslufQ9iufaU8NH016R1b2HjuxgS46MgIzhvSnllrc4Oya3orJiHswWmrKSwu56GzrUdgY4LBBUNTqVSCcox4KyYh6rsN+by/OItrj+9Cj+Tg7vPHmHDRsWVjhnVO5P3F2UE3RrwVkxBUXFbB3R+soENiPL8b093tOMaYOjh3SHs25u9j8ZbdbkepEysmIeip2RvIzN/H387sS2y0dZliTDA5pV9b4qIjeW9RtttR6sSKSYjJyC3iqdkbmDgwhVE9ktyOY4ypoyaNoji5bxs+WbqN4rLg6fzRikkIUVXu+mA5sdER3H2qdZliTLA6d0h7CkvKg2oURismIeTdRVks2LiLu05JI6lpI7fjGGOO0oguLWmXEMf7i4OnqcuKSYjYWVTCQ9NXM7RTC85Pt3tKjAlmERHC2YPbMWd9HjsKit2OUysBWUxE5EURyRWRFV7TEkXkcxFZ7/xs4UwXEXlCRDKcQbsGu5fcPQ9OX82+knIeOsvuKTEmFJwzuD2VClN+DI57TgKymAAv4xk7xdsdwJeq2h340nkNnm7xuzuPScBTfsoYML7LyGfK4myuHdWV7nZPiTEhoVOrxgzt1IL3F2UFxT0nAVlMVPUbYFeVyRPxjK2C8/NMr+mvqsc8IEFEwqZr3OKyCv784Qo6toznxjHd3I5jjGlA5wxuz4a8fSzZusftKEcUkMWkBsmqut15vgNIdp63A7yHKMtypoWFJ2dvYKPdU2JMSDqlf1tioyN4Lwi6VwmmYvIT9Rzz1fm4T0QmichCEVmYlxec3Tx7y8gt5KnZGZw5MIXjuts9JcaEmmax0ZzUpw0fL91GSXlg33MSTMUk52DzlfMz15meDXhfvtTemfYLqvqsqqaranpSUnB/+aoqf/5gBfExUdx9mt1TYkyoOnNQO/YWlzM7wMc5CaZiMhW4wnl+BfCR1/TLnau6RgAFXs1hIeuDH7OZv3EXf5rQi1ZN7J4SY0LVcd1a0bJxDFOXbHM7ymEFZDERkTeB74GeIpIlIlcDDwPjRWQ9MM55DTAdyAQygOeA612I7FcF+8t4aPpqBqYm2DglxoS4qMgITuvfli9W51BYXHbkFVwS5XaA6qjqRTXMGlvNsgrc4NtEgeXRz9aya18pL181zO4pMSYMTBzUjle+38zMFTs4L0BvSg7IIxNTs2VZe3ht/mYu/1Un+rZr7nYcY4wfDEpNoGPLeD4K4KYuKyZBpKJSufvDFbRq0ohbT+zhdhxjjJ+ICBMHpPDdhnxy9wZm9ypWTILIGwu2sCyrgLtPTaNZbLTbcYwxfjRxUDsqFaYuDcyjEysmQSKvsIRHZq7hmK4tOWNAittxjDF+1jWpCf3aNQ/Ypi4rJkHi7zNWU1xWwQMT+yJiJ92NCUcTB6awPLuADXlFbkf5BSsmQWBe5k6mLM5m0qgudGvdxO04xhiXnD4gBREC8ujEikmAKy2v5J4PV9C+RRw3ju7udhxjjIuSm8VyTNeWfLQkO+B6ErZiEuBenLuR9blF3H9GH+JirCNHY8LdxIHt2Lxzf8D1JGzFJIBl7znAf75Yz/jeyYxNSz7yCsaYkDehbxtioiL4eGlg9RplxSSAPfDxShTl3tOtI0djjEez2GiO75HE9OXbqawMnKYuKyYB6qs1OXy6MoebxnanfYt4t+MYYwLIaf3bsmNvMYu27HY7yk+smASg4rIK7p26kq5Jjbnm2C5uxzHGBJixack0iopg2rLAaeqyYhKAnpyVwdZdB/jrmX2JibKPyBhzqCaNohjdszXTlm+nIkCauuybKsBk5hXx9NeZnDkwhWO6tnI7jjEmQJ02oC15hSUs2LjL7SiAFZOAoqr85aOVNIqK4K5T09yOY4wJYGN6tSYuOpJpywPjBkYrJgFk2vLtzMnI57aTetK6aazbcYwxASw+Jooxaa2ZsXwH5RWVbsexYhIoikrK+esnq+jbrhmXjujodhxjTBA4rV9bdu4rZX4ANHVZMQkQT3y5npy9Jfx1Yl8ibfREY0wtjO7VmviYSD5Z5n5TlxWTAJCRW8iLczZyQXoqgzq0cDuOMSZIxEZHMi4tmRkrdlDmclOXFROXqSr3Tl1JfEwkt0/o6XYcY0yQOa1/W/bsL+O7DTtdzWHFxGXTl+9gbsZObjupJy2bNHI7jjEmyIzqkUTTRlF84vIIjEFXTERkk4gsF5ElIrLQmZYoIp+LyHrnZ1C0Fe0rKedv01bRu20zLhluJ92NMXUXGx3J+N7JfLrS3aauoCsmjtGqOlBV053XdwBfqmp34EvndcCbPCuD7QXFPDCxj510N8YctQl927C3uJx5me41dQVrMalqIvCK8/wV4EwXs9RKZl4Rz32bydmD25HeKdHtOMaYIDaqRxJx0ZHMXLHDtQzBWEwU+ExEFonIJGdasqoe7PFsBxDQg3+oKvd9vIrYqEjuOLmX23GMMUEuNjqS0b2S+HRljmt9dQVjMTlWVQcDJwM3iMgo75nqGcuy2ndTRCaJyEIRWZiXl+eHqNX7bFUO36zL45bxPexOd2NMg5jQty35RSUsdqlb+qArJqqa7fzMBT4AhgE5ItIWwPmZW8O6z6pquqqmJyUl+SvyIQ6UVvDAx6vomdyUK35lJ92NMQ1jdM8kYiIjmLHcnaauoComItJYRJoefA6cCKwApgJXOItdAXzkTsIje2p2Btl7DvDAxD5ERQbV22+MCWBNY6M5rnsrPl25A08DjX8F27dZMjBHRJYCC4BpqjoTeBgYLyLrgXHO64Czeec+nv4mk4kDUxjepaXbcYwxIeakvm3I3nOAFdl7/b7vKL/vsR5UNRMYUM30ncBY/yeqmwc+XkV0hHDXKda9vDGm4Y1PSyYyQpixYjv92jf3676D7cgkaH25Oocv1+Ry87juJDezk+7GmIbXonEMI7okMnOF/5u6rJj4QXFZBfd/vIpurZtw1cjObscxxoSwCX3akJm/j4zcIr/u14qJHzz7TSZbdu3n/jP6EG0n3Y0xPnRinzYAzPDzDYz2zeZjW3ftZ/KsDE7t15aR3WxMd2OMbyU3i2VIxxZ+vxveiomP/fWTVUSI8Gcb090Y4ycT+rRh1fa9bNm532/7tGLiQ7PX5vLZqhxuHNONlIQ4t+MYY8LEhL6epq6ZK7cfYcmGY8XER0rKPSfdO7dqzDXH2Ul3Y4z/pCbG0yelGZ+uzPHbPq2Y+Mjz325kY/4+7jujD42iIt2OY4wJM+N7J7N4y27yi0r8sj8rJj6wbc8B/vtVBif2Tub4Hu70AWaMCW/j0pJRha9WV9tVYYOzYuIDD05bTaUq95zW2+0oxpgw1SelGSnNY/l8tX+auqyYNLA56/OZtnw7N4zuRmpivNtxjDFhSkQY1zuZb9fnUVxW4fP9WTFpQKXlldw7dQUdEuOZNKqL23GMMWFuXFoyxWWVzFmf7/N9WTFpQC9/t5ENefu49/TexEbbSXdjjLtGdGlJk0ZRfOGHpi4rJg0kZ28x//liPWN7tWZsWkCPGmyMCRMxUREc3zOJL1bnUunj4XytmDSQB6etpqxS+cvpdtLdGBM4TuydTH5RCUuy9vh0P1ZMGsC8zJ1MXbqN60Z1oWPLxm7HMcaYn5zQozWREcIXq3zb1GXFpJ7KKiq596OVtEuI47cndHM7jjHGHKJ5fDTDOyfyuRWTwPbq95tZm1PIX07vTVyMnXQ3xgSecWnJrM8tYlP+Pp/tw4pJPeQWFvP45+sY1SOJE3vbSXdjTGAa73w/+fKqLism9fDwjDUUl1dw3+m9ERG34xhjTLVSE+Pp1aapT5u6rJgcpYWbdjFlcTa/Oa4LXZKauB3HGGMOa3zvZBZu3s3ufaU+2X7IFBMRmSAia0UkQ0Tu8OW+yisqueejlbRtHsuNY+ykuzEm8I1LS6aiUpm11jcdP4ZEMRGRSGAycDLQG7hIRHx2w8cbC7awevte7j61N/ExUb7ajTHGNJh+7ZrTumkjvvRRL8IhUUyAYUCGqmaqainwFjDRFzvaWVTCo5+uZWS3lpzSr40vdmGMMQ0uIkIY06s136zLo7S8suG33+BbdEc7YKvX6yxnWoP7x8w17C+t4P4z+thJd2NMUBmblkx0VARbdjX8JcJh1UYjIpOASQAdOnSo8/qqSseWjbn+hK50a920oeMZY4xPje6ZxA9/HkdkRMP/IRwqxSQbSPV63d6ZdghVfRZ4FiA9Pb3OvZ6JCDeMthPuxpjgFBXpu8aoUGnm+gHoLiKdRSQGuBCY6nImY4wJGyFxZKKq5SJyI/ApEAm8qKorXY5ljDFhIySKCYCqTgemu53DGGPCUag0cxljjHGRFRNjjDH1ZsXEGGNMvVkxMcYYU29WTIwxxtSbqNb53r2QICJ5wOajXL0VkN+AcYKN/f7h/fuDvQfh/Pt3VNWkqhPDtpjUh4gsVNV0t3O4xX7/8P79wd6DcP/9q2PNXMYYY+rNiokxxph6s2JydJ51O4DL7Pc34f4ehPvv/wt2zsQYY0y92ZGJMcaYerNiUkciMkFE1opIhojc4XYefxKRVBGZJSKrRGSliNzsdiY3iEikiPwoIp+4ncXfRCRBRN4TkTUislpEfuV2Jn8Skd87//ZXiMibIhLrdqZAYcWkDkQkEpgMnAz0Bi4Skd7upvKrcuAPqtobGAHcEGa//0E3A6vdDuGS/wAzVbUXMIAweh9EpB1wE5Cuqn3xDHdxobupAocVk7oZBmSoaqaqlgJvARNdzuQ3qrpdVRc7zwvxfJG0czeVf4lIe+BU4Hm3s/ibiDQHRgEvAKhqqarucTeV30UBcSISBcQD21zOEzCsmNRNO2Cr1+sswuzL9CAR6QQMAua7m8TvHgduByrdDuKCzkAe8JLTzPe8iDR2O5S/qGo28CiwBdgOFKjqZ+6mChxWTEydiUgT4H3gFlXd63YefxGR04BcVV3kdhaXRAGDgadUdRCwDwib84Yi0gJPS0RnIAVoLCKXupsqcFgxqZtsINXrdXtnWtgQkWg8heR1VZ3idh4/GwmcISKb8DRxjhGR19yN5FdZQJaqHjwafQ9PcQkX44CNqpqnqmXAFOAYlzMFDCsmdfMD0F1EOotIDJ6Tb1NdzuQ3IiJ42stXq+pjbufxN1W9U1Xbq2onPJ/9V6oaNn+ZquoOYKuI9HQmjQVWuRjJ37YAI0Qk3vm/MJYwugDhSEJmDHh/UNVyEbkR+BTPlRwvqupKl2P500jgMmC5iCxxpt2lqtNdzGT863fA684fU5nAVS7n8RtVnS8i7wGL8VzZ+CN2J/xP7A54Y4wx9WbNXMYYY+rNiokxxph6s2JijDGm3qyYGGOMqTcrJsYYY+rNiokxhyEi7UXkIxFZLyKZIvJfEWnkg/2cICLHeL2+TkQud55fKSIpR7HN90SkyxGWuVFEfl33xMYcyoqJMTVwbkybAnyoqt2B7kAc8IgPdncCXndTq+rTqvqq8/JKPN131JqI9AEiVTXzCIu+iOfeEWPqxYqJMTUbAxSr6ksAqloB/B64XESaOEcM/z24sIh8IiInOM+fEpGFztgX93sts0lE7heRxSKyXER6OZ1mXgf8XkSWiMhxInKfiNwmIucC6XhuFFwiIqeKyIde2xsvIh9Uk/0S4COv5a4SkXUiskBEnjuYW1X3A5tEZFgDvWcmTFkxMaZmfYBDOnV0OrbcBHQ7wrp/VtV0oD9wvIj095qXr6qDgaeA21R1E/A08G9VHaiq33rt7z1gIXCJqg4EpgO9RCTJWeQqPEcXVY08mF1E2gL3O9OOxTMWj7eFwHFH+H2MOSwrJsb4xvkishhPlxt9OPQL/GAHmYuATnXZqHq6rPgfcKmIJAC/AmZUs2hbPN3FAwwHZjsdFJYCb1dZNpc6NqMZU5X1zWVMzVYB53pPEJFmQBtgLdCXQ/8gi3WW6QzcBgxV1d0i8vLBeY4S52cFR/d/8CXgY6AYeFdVy6tZ5kCVfR5OrLO8MUfNjkyMqdmXQLzXVVWRwL+A/6rqATzNXQNFJEJEUvGMxAnQDM9YHwUikoxnmOcjKQSa1maeqm7DM8Lf3XgKS3VW83NT3Hw8TW0tnSEEzquybA9gRS0yGlMjKybG1MBpUjoLOFdE1gM7gUpVfdBZZC6wEc8RzBN4epNFVZfiad5aA7zhLHckHwNnHTwBX2Xey8DTzrw4Z9rrwFZVrakL9Gl4rhBDVbcD9wHfO1mqrjMS+LwWGY2pkfUabEwtOfeBvAmcpaqLXc7yX+BHVX2hhvlxwCxgpHMVmve8K4F0Vb1RRAYBt6rqZb7ObEKbnTMxppZU9Tugo9s5RGQRnma0P9S0jKoeEJF7gXZ4BnWqSSvgnoZNaMKRHZkYY4ypNztnYowxpt6smBhjjKk3KybGGGPqzYqJMcaYerNiYowxpt6smBhjjKm3/wfDhyWFMOcdnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(0, 3 * np.sqrt(10), 0.01)\n",
    "\n",
    "plt.plot(x, pi(x))\n",
    "plt.title(\"Profits ($\\pi$) v.s Quantity ($q$) \\n $\\pi(q) = 90q - q^{3}$\" \"\")\n",
    "plt.xlabel(\"Quantity (q)\")\n",
    "plt.ylabel(\"Profits ($\\pi$) \\n\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To be certain that we have found the global maximum, we really should not rely on the plot and should check if the function is concave for values of $q$ greater than 0.  \n",
    "A real-valued function is concave over an interval if its second derivative is negative over that set.  \n",
    "The second derivative of profits is $\\pi''(q) = - 6q$, which is negative for all positive values of q.\n",
    "So, we have in fact found the global maximum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_descent(\n",
    "    df, gamma, x0, max_iter=1000, tol=10 ** -6, exact_iter=None\n",
    "):\n",
    "    \"\"\"\n",
    "    Performs Gradient Descent for a real-valued function of one variable.\n",
    "    Returns the minimizer, the number of iterations and the optimization path\n",
    "    as dictionary with keys 'argmin', 'n_iter' and 'path'.\n",
    "    \n",
    "    Arguments:\n",
    "    -----------\n",
    "    df (function): The derivative of the objective function, a real-valued function of one variable\n",
    "    gamma (float): The learning rate\n",
    "    x0 (float): The initial guess\n",
    "    max_iter (int): The maximum number of iterations\n",
    "    tol = 10**-6: The tolerance. The simulation stops when |x_{n} - x_{n-1}| <= tol and returns xn\n",
    "    exact_iter (int): The exact number of iterations for the algorithm to run for\n",
    "    \"\"\"\n",
    "\n",
    "    n_iter = 0\n",
    "    x_last = np.inf\n",
    "    x_current = x0\n",
    "    x_path = [x0]\n",
    "\n",
    "    while (\n",
    "        abs(x_current - x_last) > tol\n",
    "        if not exact_iter\n",
    "        else n_iter < exact_iter\n",
    "    ):\n",
    "\n",
    "        x_current, x_last = x_current - gamma * df(x_current), x_current\n",
    "\n",
    "        x_path.append(x_current)\n",
    "\n",
    "        n_iter += 1\n",
    "\n",
    "        if n_iter >= max_iter:\n",
    "            warnings.warn(\n",
    "                \"Algorithm did not converge within the tolerance of {tol} in {max_iter} iterations.\"\n",
    "            )\n",
    "            break\n",
    "\n",
    "    return {\"argmin\": x_current, \"n_iter\": n_iter, \"path\": x_path}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def d_pi(q):\n",
    "    \"\"\"First devivative of profit.\"\"\"\n",
    "    return 90 - 3 * q ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Profit is convex, so to maximize it, we must minimize its negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'argmin': 5.477244902393321,\n",
       " 'n_iter': 10,\n",
       " 'path': [10,\n",
       "  5.8,\n",
       "  5.5816,\n",
       "  5.5123444864,\n",
       "  5.4891879821953315,\n",
       "  5.481316899962672,\n",
       "  5.478626802493688,\n",
       "  5.477705704033555,\n",
       "  5.47739011723345,\n",
       "  5.47728196745145,\n",
       "  5.477244902393321]}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient_descent(df=lambda q: -d_pi(q), gamma=0.02, x0=10, exact_iter=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'argmin': 5.477225574495832,\n",
       " 'n_iter': 7,\n",
       " 'path': [10,\n",
       "  3.7,\n",
       "  5.1678999999999995,\n",
       "  5.4642528631,\n",
       "  5.477027521429311,\n",
       "  5.477222779084842,\n",
       "  5.477225535629513,\n",
       "  5.477225574495832]}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient_descent(\n",
    "    df=lambda q: -d_pi(q), gamma=0.03, x0=10,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'argmin': 5.477225714304883,\n",
       " 'n_iter': 15,\n",
       " 'path': [10,\n",
       "  1.5999999999999996,\n",
       "  4.892799999999999,\n",
       "  5.6200609792,\n",
       "  5.429850730008825,\n",
       "  5.491857255987541,\n",
       "  5.472597721570381,\n",
       "  5.478678622913902,\n",
       "  5.47676828853305,\n",
       "  5.477369382179302,\n",
       "  5.47718034031907,\n",
       "  5.477239802673743,\n",
       "  5.477221099954523,\n",
       "  5.4772269826200795,\n",
       "  5.477225132323104,\n",
       "  5.477225714304883]}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient_descent(\n",
    "    df=lambda q: -d_pi(q), gamma=0.04, x0=10,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the above learning rates produce a solution within tolerance, but you can see that picking a learning rate that is either too small or too large will slow or, in the worst case, prevent convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = [\n",
    "    {\n",
    "        \"gamma\": gamma,\n",
    "        \"n_iter\": gradient_descent(df=lambda q: -d_pi(q), gamma=gamma, x0=10)[\n",
    "            \"n_iter\"\n",
    "        ],\n",
    "    }\n",
    "    for gamma in np.arange(0.001, 0.06, 0.001)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEaCAYAAADKaorCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deXxcdb3/8dc7S5O2SZMuaZuWbpQWKFspZRFQQDbZcUHlqnABRS8oeNXfFfWqeL14Fa/7giIodQMRFasXFKgsAgJ2YymlpS3d0zZtkzRpm7RJPr8/znfSaZqZTJPMksnn+XjMIzNnvnPO58y085nvcr5fmRnOOedcJhRkOwDnnHMDhycd55xzGeNJxznnXMZ40nHOOZcxnnScc85ljCcd55xzGeNJx+UNSW+WtCzbcTjnEvOk4/qEpNWSzslmDGb2dzM7PB37lvSEpGZJTZK2Svq9pOoUX3umpPXpiKsnJE2X9NtwHg2SXpL0CUmF2Y7N5T9POq7fyIEvxY+aWRlwGFAG/G+W4zlokqYCzwPrgGPMrAK4ApgNlGczthhF/LspT/kH69JKUoGkWyStlLRN0v2SRsQ9/1tJm8Iv7qckHRX33D2S7pD0kKSdwFmhRvWp8Ou8QdJvJJWG8vvVKJKVDc//h6QaSRslfVCSSTqsu3Mys3rgQWBm3L6ukbRUUqOkVZI+HLYPBR4GxoVaUpOkcd29L128jx+StELSdklzJY2Le84kfUTS65LqJf1AkhLs6kvAs2b2CTOrCeezzMz+JZwXki6VtCTs6wlJR6bynobzvziubJGkWkmzwuNTJD0b9vuipDPjyj4h6TZJzwC7gEMlnSdpWTjODyU9KemDca+5NhyzTtJfJU1K9T0J72fs83o1LsZxkn4X4n5D0k2JPhPXQ2bmN7/1+gasBs7pYvvNwHPAIUAJ8GPg3rjnryX6hV0CfBtYHPfcPUADcBrRD6TScJwXgHHACGAp8JFQ/kxgfaeYEpV9G7AJOAoYAvwSMOCwBOf3BPDBcH8k8Bjwx7jnLwKmAgLOIPrinNVVXKm8L53KvhXYCswKZb8HPBX3vAF/BiqBiUAt8LYE+9oEXJPkc5wO7ATOBYqB/wBWAINSeE+/APyq03uyNNwfD2wDLgyf5bnhcVXc+7s2fB5FQBWwA3hHeHwzsDfuM7gsxHVkeP4/iZJpt+8JUc1uA3Bi+LwOAyaFuBaE8xgEHAqsAs7P9v+vfLplPQC/5ceNxElnKXB23OPq8OVR1EXZyvBlUREe3wP8vIvjvD/u8e3Aj8L9/b7cuyn7U+B/4p47jO6Tzi6iJGjAYmBikvfjQeDmruLqwftyN3B73OOyUHZyeGzA6XHP3w/ckiCuvSRISOH5zwP3xz0uCF/QZ6bwnh4GNAJDwuNfAV8I9z8N/KLTsf4KXB33/v5X3HNXAf+IeyyiJsFY0nkYuK5TnLuASd29J+G4N3dx7icDaztt+wzws2z//8qnmzevuXSbBPwhNHHUE33ZtgFjJBVK+mpoYtpB9IUGMCru9eu62OemuPu7iL6EE0lUdlynfXd1nM5usqgP5FhgOFEtBQBJF0h6LjR/1RP9oh+VYD+Q5H3pouw4YE3sgZk1EdUSxseVSfU92UaU4BLpfKx2ovem22OZ2YpwHpdIGgJcCvw6lJsEXBE733DOp3eKJf4z2O/zsSgDxA/GmAR8J25f24kSUyrvyQRgZRfnPomoGTQ+xs/S9Wfieqgo2wG4vLcOuNbMnun8hKQPEDWTnEOUcCqAOqIvj5h0TYNeQ1zSIPoiSomZvSzpv4EfhL6AQcDviH6d/9HM9kp6kH3n0dU5JHxfurCR6AsR6OgnGklUAzlYjwHvBH6W5FjHxB1LRO9Nqse6F7iSqObxakhEEJ3vL8zsQ0leG/8+7ff5hDjiP691wG1m9qsU44q3jqgptKvtb5jZtB7s06XIazquLxVLKo27FQE/Am6LdfJKqpJ0WShfDrQQ/foeAnwlg7HeD1wj6cjwq/zzB/n6OUS/gC8lSjolRP0GrZIuAM6LK7sZGCmpIm5bsvels3tDrDMllRC9T8+b2eqDjBngi8Cpkr4uaWw49mGSfimpkuh9uUjS2ZKKgU8SfUbPprj/+4jO/d/YV8uBqM/sEknnhxpuqaKBH4d0uRf4P+AYSZeHf0c3AmPjnv8R8BmFgSeSKiRdkWKMdwGfknSCIoeFz+EFoFHSpyUNDnEeLenEFPfrUuBJx/Wlh4Ddcbdbge8Ac4FHJDUSdZ6fHMr/nKgpZwPwanguI8zsYeC7wONEHdKxY7ek+Po9ROf2eTNrBG4i+sKuA/6F6JxjZV8jShyrQrPNOJK/L52P9RhRUvwdUQ1gKvDegznfuH2tBN4ETAaWSGoI+50PNJrZMuD9RIMVtgKXAJeE801l/zXAP4BTgd/EbV9HVKv9LFFyXgf8PxJ8B5nZVqIO/9uJfpTMCDG2hOf/AHwNuC80zb4CXJBijL8FbiNKio1E/W8jzKwNuJhoVOIb4fzvIqqBuz6i0Fnm3ICmaFjwK0CJmbVmOx63P0XX7awH3mdmj2c7HtdzXtNxA5akt0sqkTSc6Ffznzzh5I7QFFcZmhQ/S9RHlrHasEsPTzpuIPswsIVoJFMbUT+Eyx1vIvpsYs18l5vZ7uyG5HrLm9ecc85ljNd0nHPOZYxfp5PAqFGjbPLkydkOwznn+pUFCxZsNbOqRM970klg8uTJzJ8/P9thOOdcvyJpTbLnvXnNOedcxnjScc45lzGedJxzzmVMv0w6Yd6mF8JCUEskfSlsvycsvLQ43GaG7ZL0XUWLYL0UW7DJOedcZvXXgQQtwFvNrClMSvi0pIfDc//PzB7oVP4CYFq4nQzcQYJ5rpxzzqVPv6zpWKQpPCwOt2RXuV5GtBiYmdlzQKWkZGuKOOecS4N+mXQAwrTji4mmMXnUzJ4PT90WmtC+FeZsgmhhp/gFotaz/2JPsX1eL2m+pPm1tbVpjd855waifpt0zKzNzGYSLex0kqSjiZaWPYJo7fMRREvkHsw+7zSz2WY2u6oq4bVNSW3fuYdvPrqcJRsbevR655zLZ/026cSYWT3RmihvM7Oa0ITWQrQy4kmh2Ab2XxnyEHq26mK3du1p5bvzXmfJhh3p2L1zzvVr/TLphFUWK8P9wcC5wGuxfpqwtO3lROujQLRY1lVhFNspQENYbKrPlZcUA9DY4jPkO+dcZ/119Fo1MEdSIVHivN/M/izpb5KqiNbdWAx8JJR/CLiQaIXIXcA16QpsaEkhAE3NnnScc66zfpl0zOwl4Pgutr81QXkjWmM97YoKCxhcXEhTy95MHM455/qVftm8luvKSoto8uY155w7gCedNCgrKaLRm9ecc+4AnnTSoKzEazrOOdcVTzppUFZS5AMJnHOuC5500sD7dJxzrmuedNKg3JvXnHOuS5500sBrOs451zVPOmkQ69OJLg9yzjkX40knDcpKi2htN1pa27MdinPO5RRPOmlQXhJN9ODX6jjn3P486aRBWWmUdLxfxznn9udJJw3KwkzTfq2Oc87tz5NOGpTFmtd80k/nnNuPJ500KI81r3lNxznn9uNJJw2GlnifjnPOdcWTThqUedJxzrkuedJJg3Ifveacc13ypJMGJUUFFBXI+3Scc64TTzppIMnnX3POuS540kkTX1PHOecO1C+TjqRSSS9IelHSEklfCtunSHpe0gpJv5E0KGwvCY9XhOcnpzvGspIiGr2m45xz++mXSQdoAd5qZscBM4G3SToF+BrwLTM7DKgDrgvlrwPqwvZvhXJpVV7qNR3nnOusXyYdizSFh8XhZsBbgQfC9jnA5eH+ZeEx4fmzJSmdMZb5Qm7OOXeAfpl0ACQVSloMbAEeBVYC9WYW+6ZfD4wP98cD6wDC8w3AyHTGV1Za7EnHOec66bdJx8zazGwmcAhwEnBEb/cp6XpJ8yXNr62t7dW+ykqKfGkD55zrpN8mnRgzqwceB94EVEoqCk8dAmwI9zcAEwDC8xXAti72daeZzTaz2VVVVb2Kq7y0iCaf8NM55/bTL5OOpCpJleH+YOBcYClR8nlXKHY18Mdwf254THj+b5bmtaTLSopo3tvO3jZfPdQ552KKui+Sk6qBOZIKiRLn/Wb2Z0mvAvdJ+m9gEXB3KH838AtJK4DtwHvTHWBs0s+dLa1UDhmU7sM551y/0C+Tjpm9BBzfxfZVRP07nbc3A1dkILQO5XGTfnrScc65SL9sXusPfMlq55w7kCedNOlY3sBHsDnnXAdPOmkSq+n4VDjOObePJ500KfeajnPOHcCTTpp4n45zzh2oz5OOpDGS7pb0cHg8Q9J13b0u33ifjnPOHSgdNZ17gL8C48Lj5cDH03CcnDZ0kPfpOOdcZ+lIOqPM7H6gHTom2GxLw3FyWkGBfCE355zrJB1JZ6ekkURLDRDWuWlIw3FyXrS8gc+/5pxzMemYkeATRHOdTZX0DFDFvvnQBpSyUl9Txznn4vV50jGzhZLOAA4HBCwzswH5c9+XN3DOuf2lY/TajUCZmS0xs1eAMkk39PVx+oOykiJ2ek3HOec6pKNP50NhjRsAzKwO+FAajpPzfMlq55zbXzqSTqEkxR6E5QcG5DTLZaU+es055+KlYyDBX4DfSPpxePzhsG3AKSsp8ut0nHMuTjqSzqeJEs2/hcePAnel4Tg5rzyMXjMz4ip/zjk3YKVj9Fo7cEe4DWhlJUWYwa49bR0riTrn3EDWbZ+OpJslDVPkbkkLJZ2XpPxpkh6VtFzSKklvSFrVt2H3Dz7pp3PO7S+VgQTXmtkO4DxgOPAB4KtJyt8NfBM4HTgRmB3+DjixST/9Wh3nnIuk0uYT64y4EPiFmS1R8g6KBjN7uPeh9X/lXtNxzrn9pFLTWSDpEaKk81dJ5YTJPBN4XNLXJb1J0qzYrU+iDSRNkPS4pFclLZF0c9h+q6QNkhaH24Vxr/mMpBWSlkk6vy/jSaSspBjw5Q2ccy4mlZrOdcBMYJWZ7QqTeV6TpPzJ4e/suG0GvLVnIXapFfhkmHKnnCgxPhqe+5aZ/W98YUkzgPcCRxEtufCYpOlmltbZrzvW1PFJP51zDkgh6ZhZu6TNwAxJqZQ/q08iS36MGqAm3G+UtBQYn+QllwH3mVkL8IakFcBJwD/SGWesec37dJxzLtJtEpH0NeA9wKvsWxfHgKcSlB8DfAUYZ2YXhFrGm8zs7r4J+YDjTQaOB54HTgM+KukqYD5RbaiOKCE9F/ey9SRPUn0iNkza519zzrlIKn06lwOHm9mFZnZJuF2apPw9ZGjlUEllwO+Aj4cRdncAU4maA2uAbxzk/q6XNF/S/Nra2l7HN7SkEPCBBM45F5NK0lkFFB/EPjOycqikYqKE8ysz+3041mYzawsXqP6EqAkNYAMwIe7lh4Rt+zGzO81stpnNrqqq6nWMJUWFDCoq8KlwnHMuSGUgwS5gsaR5QEtso5ndlKB82lcODUO27waWmtk347ZXh/4egLcDr4T7c4FfS/omUQ1sGvBCX8aUSLkvWe2ccx1SSTpzwy1VmVg59DSii1RflrQ4bPsscKWkmUQJbzXRHHCEa4vuJ+qXagVuTPfItRhfPdQ55/ZJZTTaHEmDgOlhU9KVQDOxcqiZPc2+i1bjPZTkNbcBt/VlHKko85qOc851SGX02pnAHKKag4AJkq42s0Sj197RadN0SQ3Ay2a2pXfh9j++vIFzzu2TSvPaN4DzzGwZgKTpwL3ACQnKXwe8CXg8PD4TWABMkfRfZvaLXkXcz5SXFrGxvjnbYTjnXE5IJekUxxIOgJktDyPHku3zSDPbDB3X7fycaKaCp4ABlXR8yWrnnNsnlaQzX9JdwC/D4/cRXXiZyIRYwgm2hG3bJQ24+WB8IIFzzu2TStL5N+BGIDZE+u/AD5OUf0LSn4HfhsfvDNuGAvU9DbS/Kisp9oEEzrl+41uPLmdQUQE3nnVYWvafyui1FqL1cb7ZXdlQ/gZJ7yRaTweiprXfmZkBaZ+XLdeUlxaxp62dltY2SooKsx2Oc84l9dclmzhk+JC07T9h0pF0v5m9W9LLhAs945nZsV28phBYYmZHEM0WMOCVdcy/5knHOZf7ahqaOWnKiLTtP1lN5+bw9+JUd2ZmbWG9molmtrZ3oeWH2KSfTc2tjBg6KMvROOdcYrv2tNKwey9jK0rTdoyESSduOpkbzOzT8c+Fmac/feCrgGhJ6yWSXgB2xu0v2SSheatjyWpfU8c5l+NqGqLLO8ZVDE7bMVIZSHAuByaYC7rYFvP5XkWUZzqWrPbBBM65HFcTrinMSk1H0r8BNwCHSnop7qly4JlErzOzJyVNAqaZ2WOShgADtjNj3+qhnnScc7mtpmE3kL2azq+Bh4H/AW6J295oZtsTvUjSh4DrgRFEa9uMB34EnN3raPuhslJPOs65/iHWvDamoiRtx0i4no6ZNZjZajO70szWALuJRrGVSZqYZJ83Es0CvSPs53VgdB/G3K+Ul/iS1c65/qGmYTejygaldaRtt4u4SbpE0uvAG8CTRBN/PpzkJS1mtifu9UV0MeR6oPCajnOuv6hpaKY6jU1rkNrKof8NnAIsN7MpRM1kzyUp/6SkzwKDJZ1LNDPBn3odaT81uLiQAvlAAudc7qupb07rIAJILensNbNtQIGkAjN7HJidpPwtQC3wMtEiag8B/9nrSPspST7pp3OuX6hp2M24NCedVIZM10sqI5oh+leSthB3/U0XLgd+bmY/6YsA80F5abH36TjnctrOllZ2NLcyNgea1y4DdgH/DvwFWAlckqT8JcBySb+QdHHo0xnQopqOXxzqnMtdHReGVmaxeS3MpfZnM2s3s1Yzm2Nm3w3NbV0ys2uAw4j6cq4EVoalEQasstIidra0ZTsM55xLKHaNzthhWUw6ZtYGtEuqOJidmtleohFu9xGtGnp5jyPMA75ktXMu1+2r6aS3eS2Vpq8m4GVJj7L/XGo3dVVY0gXAe4iWqX4CuAt4d28D7c/KSopYX7cr22E451xCsSlwRg9L34WhkFrS+X24peoq4DfAh8NaPH1O0gSidXrGEF0DdKeZfUfSiHDsyUTXE73bzOokCfgOcCFR/9S/mtnCdMTWFR+95pzLdZt27GZUWUnal2BJZRG3OZIGAxPNbFkK5a/sk8iSawU+aWYLJZUDC0JN7F+BeWb2VUm3EA3f/jTRBKXTwu1k4I7wNyPKSov8Oh3nXE7bWN9MdZqHS0OKMxIAi4lGriFppqS5Scq/Q9Lrkhok7ZDUKGlH34UcLbsQq6mYWSOwlGiOt8uAOaHYHPb1JV1GNIzbzOw5oFJSdV/GlExZSRE797TR1j5gJ2ZwzuW4TQ05knSAW4GTgHoAM1sMHJqk/O3ApWZWYWbDzKzczIb1OtIEJE0GjgeeB8bErQO0iaj5DaKEtC7uZevDts77ul7SfEnza2tr+yzG2PIGO/d4bcc5l5s2NuzOmaSz18waOm1rT1J+s5kt7UVMKQsXrf4O+LiZ7VebMjPjIOd8M7M7zWy2mc2uqqrqszjLSnxNHedc7mpqaaWxuZXqNI9cg9QGEiyR9C9AoaRpwE3As0nKz5f0G+BBoGMggZkdzGCEbkkqJko4v4rb92ZJ1WZWE5rPtoTtG4AJcS8/JGzLCJ/00zmXyzaFa3RypabzMeAoogTya6ABuDlJ+WFEI8TOI5qd4BLg4t6Fub8wGu1uYKmZfTPuqbnA1eH+1cAf47ZfpcgpQENcM1zalfnyBs65HLYxDJdO9wzTkFpN5yIz+xzwudgGSVcQzThwgDAjQbqdBnyA6PqhxWHbZ4GvAvdLug5Yw77rgx4iGi69gighZiLGDuVe03HO5bBNDbGkk/6aTipJ5zMcmGC62gaApEOA7xElBoC/Azeb2fqeBtmZmT0NKMHTB6xQGvp3buyr4x+sspJiwPt0nHO5aWNoXhuT5ilwIEnSCTMLXAiMl/TduKeGEV0nk8jPiJrhrgiP3x+2ndu7UPuvWJ/OTq/pOOdy0KaGZkaVlTCoKJUel95JdoSNwHygmWj+tNhtLnB+ktdVmdnPwgShrWZ2D9B3Q8H6oY4+HU86zrkctLGhOe2zS8ckrOmY2YvAi5J+ZWYH8225TdL7gXvD4yuBhLNSDwRDB0XTSnjzmnMuF21q2M3kkUMzcqyENR1J94e7iyS91PmWZJ/XEnXgbwJqgHeR4Y77XFNUWMDg4kJfU8c5l5Nq6pvTPrt0TLKBBLFh0Qc13NnM1gCX9jiiPFVWWsSO3V7Tcc7llsbmvTS2tDI2AyPXIElNJ3Ydi5mt6eqW6HWS5kiqjHs8XNJP+zbs/md85WDW+fIGzrkck8nh0pDaxaEH61gzq489MLM6ornRBrSpVWWsrG3KdhjOObefmobMXRgK6Uk6BZKGxx6ENW5SuR4or00dPZTNO1pobPZ+Hedc7qjJ4BQ4kHwgwbzw92sHuc9vAP+Q9GVJXyaap+32noeYH6ZWlQGwqnZnNyWdcy5zNtY3I2XmwlBIXgOplnQqcKmk++g0A0CilTfN7OeS5gNvDZveYWav9km0/Vgs6aysbeK4CZXdlHbOuczI5IWhkDzpfAH4PNGMzN/s9JyxL6kcICSZAZ9o4k0aOYSiAnm/jnMup2RqHZ2YZBeHPgA8IOnzZvbljEWUp4oLC5g4Yog3rznncsqmhmYOrcrMhaGQQge/mX1Z0qXAW8KmJ8zsz+kNKz8d6iPYnHM5pqahmdMOG5Wx43XbiCfpf4guFI01md0s6SsJyhZKerxvQ8wfU0cPZfXWXbS2JVt41TnnMqOxeS9NLa250bwW5yJgppm1Q3TxJ7CIaP2a/ZhZm6R2SRVdLHE94E2tKmNPWzvr63YzeVTmqrPOOdeV2DU6mZqNAFK/fqYS2B7uV3RTtolocbVHgY4ODDO76eDDyy/xI9g86Tjnsi2WdDI17xqklnT+h2jSz8eJhk2/BbglSfnfh5vrZGrorFtZ28TZR47JcjTOuYGupj6zF4ZCagMJ7pX0BHBi2PRpM9uUpPwcSYOA6WHTMjPzy/CByiGDGDl0ECu3+Ag251z21TRk9sJQSLF5LUz+OTeVspLOBOYAq4lqRhMkXW1mT/UwxrwytaqMVVt9BJtzLvtqGnZTVVZCcWFmLgyF9MyJ9g3gPDNbBiBpOtGCbiek4Vj9ztTRQ/nrks3ZDsM556hpaM5o0xqkZ8LP4ljCATCz5UBxXx9E0k8lbZH0Sty2WyVtkLQ43C6Me+4zklZIWiYp2XLbaTW1qoztO/ewfeeebIXgnHNALOlkbhABdJN0wnU3rx3kPudLukvSmeH2E2B+z0NM6B7gbV1s/5aZzQy3hwAkzQDeCxwVXvNDSYVpiKlb+yb+9CY251x2bWpozuhwaegm6ZhZG7BM0sSD2Oe/EV1EelO4vRq29anQR7S924KRy4D7zKzFzN4AVgAn9XVMqYgfNu2cc9lS07CbppZWJo0cktHjptKnMxxYIukF9r/uZr8lqSXNM7Ozgf8ys09z4CShmfJRSVcR1a4+GRaRGw88F1dmfdi2H0nXA9cDTJx4MHk2deOHD2ZQUQErfQ4251wWLVwTrbU5a+Lwbkr2rVSSzudT3FePlkLoY3cAXyaaBfvLRIMark31xWZ2J3AnwOzZsy0dARYWiCkjh3rzmnMuqxasqaO0uIAZ44Zl9LipXKfzpKRJwDQze0zSEKCr/pAeL4XQV8ysY1hY6EuKTUy6AZgQV/SQsC0rpo4eytKaxmwd3jnnWLC2jmMPqczocGlIbcLPDwEPAD8Om8YDD3YuZ2YPmNkFwO1mdlanW9oTToi1Ou7h24HYyLa5wHsllUiaAkwDXshETF2ZWlXG2u27aGlty1YIzrkBrHlvG69ubMh40xqk1rx2I1Gn+/MAZva6pNGJCmdq7R1J9wJnAqMkrQe+CJwpaSZRzWo18OEQ0xJJ9xMNamgFbgyDJLJialUZbe3G2m27mDamPFthOOcGqJc3NLC3zThhUm4mnRYz2yNFXTSSioi+1LPKzK7sYvPdScrfBtyWvohSFz+CzZOOcy7TFqypA2DWxMqMHzuVxrwnJX0WGCzpXOC3wJ/SG1Z+m9Ix8aePYHPOZd6CNXVMGTWUkWUlGT92KjWdW4DrgJeJmqseAu7qXEjSiGQ7MbNUr6nJe2UlRYwdVurX6jjnMs7MWLimjjMOr8rK8VMZvdYeFm57nqhZbZmZddW8tiA8L2AiUBfuVwJrgSl9FXQ+mDp6qNd0nHMZt3b7Lrbt3JOV/hxIbfTaRcBK4LvA94EVki7oXM7MppjZocBjwCVmNsrMRgIXA4/0bdj939SqMlZtaaLr/O2cc+kR68/J2aRDdIHlWWZ2ppmdAZwFfCtJ+VNic54BmNnDwKm9CzP/TK0qo7GlldrGlmyH4pwbQBasqaO8pIhpo7MziCmVpNNoZiviHq8Ckl3ZuFHSf0qaHG6fAzb2Kso8FBvBtsL7dZxzGbRgTR0zJ1ZSWKDuC6dBwqQj6R2S3kE0a/RDkv5V0tVEI9f+mWSfVwJVwB+Ilq2uCttcnEN9BJtzLsMam/eyfHNjVi4KjUk2kOCSuPubgTPC/Vog4QIMYZTazZKGmpl/oyYwdlgpQwYVsnKL13Scc5nx4roG2i17/TmQJOmY2TU92WGY9PMuoAyYKOk44MNmdkPPQsxPBQXi0KqhrNrqedk5lxkL1tQhwcwsXBQa0+2Q6TBX2ceAyfHlOy9tEOdbwPlE851hZi9KekuvI81D00aX8/fXa2lvNwqy1L7qnBs4FqytY/rocoaV9vlizilL5eLQB4mml/kT0J7KTs1sXWzanMBntuzC6YeN4g+LNvDyhgaOm5C9Xx7OufzX3m4sWlvHxceOy2ocqSSdZjP77kHsc11oYjNJxcDNwNIeRZfnzjpiNAWCeUs3e9JxzqXVitomGptbs9qfA6kNmf6OpC9KepOkWbFbkvIfIZqZejzRmjUzw2PXyYihg5g1cTiPLd2S7VCcc3ku2xeFxqRS0zkG+ADRImyx5rUuF2WTVAh8x8ze12cR5rlzZozhqw+/xsb63YyrTDgo0DnnemXBmjpGDB3E5JFDshpHKgqsra4AABqYSURBVDWdK4BDzeyM7hZlC2vUTJI0qE+jzGPnHBktTTTvNa/tOOfSZ+HaOmZNrKRTf3vGpVLTeYVo0s5UvxVXAc9Imgt0jAc2s85LWDuimQkmjRzCvKWb+cApk7IdjnMuD23fuYdVtTt51wmHZDuUlJJOJfCapH8CHROFJRkyvTLcCgBfoawbkjj7iDH88vk17NrTypBBqXwkzjmXumdWbAVg9qSkK9BkRCrfcF88mB2a2Zd6GMuAdc6Ro/npM2/w99e3cv5RY7MdjnMuz/xu4XrGVZRmfRABpLaezpOp7EjSt83s45L+RBfLWSepGQ14J04ZQXlpEfOWbvak45zrU5samnlqeS03nHlY1ib5jJfKjASN7Esig4BiYKeZDetU9Bfh7//2XXgDQ3FhAWdMr+Jvr/nsBM65vvWHRRtoN3hnDvTnQGo1nY5+GUXDHi4DTumi3ILwN6WaUW9J+inRAnFbzOzosG0E8BuiKXtWA+82s7oQ93eAC4FdwL+a2cJMxJmqc44cw59fquHF9fUcn8UZYJ1z+cPMeGDBOk6cPJwpo4ZmOxwgtSHTHSzyINHcal2S9IakVZ1vvY70QPcAb+u07RZgnplNA+aFxwAXANPC7XrgjjTE0ytnHl5FYYGY5xeKOuf6yOJ19azMkVFrMak0r70j7mEBMBtoTvKS2XH3S4mu8+nzIRNm9pSkyZ02XwacGe7PAZ4APh22/9yitaGfk1QpqdrMavo6rp6qHDKIEyYN57Glm/nU+YdnOxznXB54YMF6SosLuPCY6myH0iGVms4lcbfziVYNvSxRYTPbFnfbYGbfBi7qk2i7NyYukWwCxoT744F1ceXWh237kXS9pPmS5tfW1qY30i6cc+RoXtvUyPq6XRk/tnMuvzTvbWPuixu54OhqyrM4q3RnqfTpHNS6Op3mZYvVjDJ+8YmZmaQDRtF185o7gTsBZs+efVCv7QtnHzmGrzz0Gn97bQtXvWlypg/vnMsjj7y6mcbmVq7IoaY1SJIMJH0hyevMzL6c4LlvxN1vJXToH3xoPbI51mwmqZp9syhsACbElTskbMspU6vKmDJqKI8t9aTjnOudBxasZ3zlYE45dGS2Q9lPsua1nV3cAK4j6ifpUtz8bGeZ2blm9iEzW9ZnESc3F7g63L8a+GPc9qsUOQVoyKX+nHhnHzGa51ZuY0fz3myH4pzrpzY1NPP067W8c9b4nLsEI9ly1R01FknlROviXAPcx/61mf1I+kSyA/bVHGyS7iUaNDBK0nqimRO+Ctwv6TpgDftqWA8RDZdeQTRkukdLcWfC5ceP566n3+CXz63hhjMPy3Y4zrl+6HcL1+fUtTnxkva1hOtePgG8j2g02Cwzq+tmn7OBEwnLVRMNQHgBeL13oe7PzK5M8NTZXZQ1+smaPkePr+At06v46dNvcO1pUygtLsx2SM65fsTM+N2C9Zw0ZQSTRubGtTnxEjavSfo68E+i0WrHmNmtKSQciPpLZpnZJ83sk8AJwEQz+5LPy5aaG86cytamPdw/f133hZ1zLs7CtfWs2ppb1+bES9an80lgHPCfwEZJO8KtUdKOJK8bA+yJe7yHfUOXXQpOnjKCEyYN58dPrmJvW3v3L3DOueCOJ1YwdFBhTl2bEy9h0jGzAjMbbGblZjYs7lbexbxr8X4OvCDpVkm3As8TNc25FEnihjOnsqF+N3MXb8x2OM65fuKxVzfz2NIt3HT2NMpKcnOZlIOaBicVZnYbUUd9XbhdY2Zf6evj5Lu3HjGaI8aWc8eTK2lvz/glQ865fqZ5bxtf+vMSpo0u49rTp2Q7nIT6POkEQ4AdZvYdYL2k3H0HcpQkbjjrMFZsaeKRVzdnOxznXI6744mVrNu+my9ddhTFhen6au+9Po9M0heJruP5TNhUDPyyr48zEFx49FgmjRzCD59YQTQAzznnDrRm207ueHIllx43jlOnjsp2OEmlIx2+HbiUcDGpmW3El63ukaLCAj5yxlReWt/AMyu2ZTsc51wOMjNunbuE4gLxuYuOzHY43UpH0tkTrosxAEm5N1C8H3nHrPGMGVbCDx5fke1QnHM56NFXN/P4slr+/dzpjBlWmu1wupWOpHO/pB8DlZI+BDwG/CQNxxkQSooK+dCbD+Ufq7bx3Cqv7Tjn9tm9p40v/elVpo8p4+pTJ2c7nJT0adIJK3T+BngA+B1wOPAFM/teXx5noLnypIlMGDGYT97/IvW79nT/AufcgPC9v73Ohvrd/NdlR+f04IF4fRplaFZ7yMweNbP/Z2afMrNH+/IYA9HQkiK+f+UstjQ288n7X/Qh1M45Hnt1M3c8uZJ3nXBIzs0knUw6UuNCSSemYb8D2nETKvnchUcy77Ut3PV0Olb/ds71F8s3N3LzfYs4elwFX77s6GyHc1DSkXROBv4haaWklyS9LOmlNBxnwLn61MlccPRYvvaXZSxYsz3b4TjnsqBu5x4+OGc+gwcVcedVJzB4UP+aFDgdSed8YCrwVqIZpi8Of10vSeJr7zqW8ZWD+eivF1G30/t3nBtI9ra1c+OvF7KpoZkff+AEqisGZzukg5aOaXDWdHXr6+MMVMNKi/nh+2axrWkPn7h/sffvODeA3PZ/S3l25Ta+8o5jOGHS8GyH0yP9Y7iD28/R4yv4/MVH8viyWr49r0+XKXLO5ah7X1jLPc+u5oOnT8nZZQtSkZvTkLpuvf+USby4voHvznud4gLxsbOnZTsk51yaPPxyDZ9/8BXeMr2KWy44Itvh9IonnX5KEl9757G0txvfeHQ5gCce5/LQwy/X8NF7FzFzQiU/+JfjKeon1+Mk4kmnHyssEF+/4jgAvvHocgy4yROPc3kjPuHMufaknF0j52D0/zMY4DoSj+Cbjy7HDG4+xxOPc/3dQy/X8LE8SziQp0lH0mqgEWgDWs1stqQRRFP0TAZWA+82s7psxdiXCgvE1991HEJ867HltLa38+/nTKegQNkOzTnXA7GEc/yESu7Jo4QD+T167Swzm2lms8PjW4B5ZjYNmBce543CAnH7u47l3bMP4Xt/W8H1v1hAw+692Q7LOXcQ9ra18915r+dtwoH8TjqdXQbMCffnAJdnMZa0KCyIBhd88ZIZPLFsC5d+/2le3bgj22E551Lw6sYdXP6DZ/jmo8u58JjqvEw4kL9Jx4BHJC2QdH3YNsbMasL9TcCYzi+SdL2k+ZLm19bWZirWPiWJa06bwn3Xn0Lz3jbe/sNneGDB+myH5ZxLYE9rO99+bDmXfv9pNu9o4UfvP4HvXXl8XiYcAOXjMsiSxpvZBkmjgUeBjwFzzawyrkydmSW8pHf27Nk2f/78DESbPrWNLXzs3oU8t2o7V540kS9cPKPfzdPkXD5bsrGBT/32JZbW7OCymeO49ZKjGD50ULbD6hVJC+K6NQ6Ql6nUzDaEv1sk/QE4CdgsqdrMaiRVA1uyGmQGVJWX8MvrTuZ/H1nOj55cybMrt/L1dx3HSVNGZDs05wa0Pa3t/ODxFfzg8RVUDhnEjz9wAucfNTbbYWVE3jWvSRoqqTx2HzgPeAWYC1wdil0N/DE7EWZWUWEBt1xwBPd+6BTM4D13/oNb5y5h157WbIfm3ID0yoYGLv3+03xn3utcctw4HvvEWwZMwoH8rOmMAf4QLWJKEfBrM/uLpH8SLaV9HbAGeHcWY8y4N00dyV8+/mZu/8sy7nl2NY8v28Lt7zyWk/vR4k/O9Wd7Wtv5/t9e54dPrGT40EH85KrZnDvjgK7lvJeXfTp9IR/6dBJ5btU2/uOBl1i7fRfnHDma604/lFMOHUFI1M65PtTU0spv56/jZ8+sZu32Xbxj1ni+cPEMKof0776bRLrr0/Gkk0A+Jx2AXXta+fGTq/jFc2vYvnMPM6qH8cE3T+HiY8cxqCjvWl2dy7gN9buZ8+xq7n1hLY3NrcyaWMnHzp7GWYePznZoaeVJp4fyPenENO9t4w+LNnD302+wYksTo8tLuPrUyfzLSRP7/Sga57Jh8bp67vr7Kh5+ZRMAbzt6LNedPoVZE/vn+jcHy5NODw2UpBNjZjy5vJa7n36Dv7++ldLiAt456xCuPX0KU6vKsh2eczmtta2dR17dzN1Pv8GCNXWUlxTx3pMmcPWpkzlk+JBsh5dRnnR6aKAlnXjLNjXy06ff4A+LN7CntZ2zDq/iqlMnc/phoyju59OqO9eXtjQ2M3fxRu55djXr63YzccQQrjltMlfMnpC3F3d2x5NODw3kpBOztamFXz63hl8+t4atTXuoHFLMeTPGcOEx1ZzmCcgNUFsam/nrK5v4v5dreP6N7ZjBSZNHcO3pUzh3xhgKB/hEu550esiTzj4trW08tXwr//fSRh5buoWmllYqBkcJ6KJjPQG5/LazpZUX19WzaF09Ty2v5YXVUaKZWjWUi46p5qJjx3H42PJsh5kzPOn0kCedrjXvbePvr2/loZdreOzVzTTGJaALj63mtKmjfPSb69d27Wnl8ddqeWblVhatrWfZph20h6/Jw8eUc/7RY7n42Gqmj/FE0xVPOj3kSad7La1t/H15lIAeDQloWGkRsyePYNbESmZNHM5xEyoZOkDbtl3/sbOllb+9toWHXq7h8WVbaN7bTnlJETMnVnL8xOEcP7GS4ydU5u21NX1pQM695jKjpKiQc2aM4ZwZY2hpbePp17fyyJLNzF+znb+9Fk1tVyA4fOwwZoX/vLMmVjJl1FC/ENVlVfPeNl5a38DCtXXMX13H0ytqad7bzqiyEq44YQIXHVvNiZNHDPj+mXTwmk4CXtPpnfpde1i8rp6Fa+tZtLaOxWvraWyJ5nsbPqS4IwEdH2pDA3Wkj0ufvW3tbGlsoaZ+NzUNzWxqaGbt9l0sXlfP0podtIY2s0kjh3DG9CouPMYTTV/w5rUe8qTTt9rajRVbmli0to6Fa+tYuLaeFVuagKg2NH1MeUczxqyJwzl01FBfbtulrL3dWLW1iYVr6sO/rzpWbGnq6IuJKSsp4ujxw5g1cXjHv7dRZSXZCTpPedLpIU866deway+L19ezcE0di9ZFNaLG5qg2VDG4uCMBHT+xkkkjhjKmooSSIl8PaKBqbWtn+eYm1m7fxaaGqPYS3XazbFMjO8K/nWGlRRw/cTjHjK9g/PDBjK0oZVxF9HdYaZE37aaZ9+m4nFUxpJgzpldxxvQqIPq1urK2iUVr9/1afXJ5LfG/i0YOHUR1ZSljhw3m8LFlHb9YR/iUPXnFzNjaFGuirWPhmjpeWt/A7r1tHWUGFRVQXVHK2GGlXHRsdUeT7aGjyryWnMO8ppOA13Ryw47mvbyyvoEN9ft+2W5q2M3G+mZW1DbRFtpPJo8cwqyJwzmyehjVlaVUVwymuqKU0eUlFPk1RDnFzKjftbejlhLrb9nYsJtNcbWX5r3tABQViBnjhnXUeqdWlVFdUcqIoYO81pKDvKbj+rVhpcWcetioLp/bvaeNl9ZHF+0tXFPHU69v5feLNuxXpkBQXTGYmRMqo2GvE4dz9Phh3kyXJrGE0jmB1DQ0U1PfzKYd+yeUmMICMaa8hLEVpcwYN4yzjxhNdeVgjhlfwTHjK3yZ9TziNZ0EvKbT/5gZO3a3UrMj7kuuYTertu5k0dp6NtTvBmBQYQEzxg1j0sgh+7X3V1eUMiTFL7eiggLGDCsdUF+GZsb2nXs6aiY1O5qpqT8wubS0HphQxg4rZWxFdKseVkp1ZVQTjb3/o8oGeY00T3hNxw0YkqgYUkzFkGKOGDvsgOe37GjuGML94vqor2BzQwt72tq72FtqKocUdzTljRlW0uV0QJWDixlbMTg0+5VSPWwwwwbnVod2fELpaMKMNXvV7w41lGb2dEooRQViTEgoR4+v4LyjxjJ2WCnjKksZM6yUcZWDGVVW4sOQXQev6STgNZ2Bob3d2L5rDzX10S/1zr/SE2lpbWdzaCqKvpib2dLY3NHHFGPAjt17Dxi625XiwugLvLpiX5/U2Lj71ZWljBpactCd5GbGtp17umjuCv0p3SSU6NghhtBxX105mHEVpYz0hOI68ZqOc0kUFIhRZSWMKivhmEMq0nKM1thFirEv/Prmjgtl47W0trE5JIbF6+r5yyvNB9TCigvF6PLUmwGbW9vYvKPlgIQSn+COPaSS848q7UgqsSQ3quzgE5xz3fGk41yaFRUWMK5yMOMqBwOprx7ZuYayKW6kV3NrW/c7AIoLCzr6UcZWDGZcZVRT8YTismVAJR1JbwO+AxQCd5nZV7McknMJSftqYUePT08tzLlMGzDDRSQVAj8ALgBmAFdKmpHdqJxzbmAZMEkHOAlYYWarzGwPcB9wWZZjcs65AWUgJZ3xwLq4x+vDtg6Srpc0X9L82trajAbnnHMDwUBKOt0yszvNbLaZza6qqsp2OM45l3cGUtLZAEyIe3xI2Oaccy5DBlLS+ScwTdIUSYOA9wJzsxyTc84NKANmyLSZtUr6KPBXoiHTPzWzJVkOyznnBpQBk3QAzOwh4KFsx+GccwOVz72WgKRaYE0KRUcBW9McTqbl2znl2/lA/p1Tvp0P5N85pXo+k8ws4UgsTzq9JGl+ssnt+qN8O6d8Ox/Iv3PKt/OB/DunvjqfgTSQwDnnXJZ50nHOOZcxnnR6785sB5AG+XZO+XY+kH/nlG/nA/l3Tn1yPt6n45xzLmO8puOccy5jPOk455zLGE86SUh6m6RlklZIuqWL50sk/SY8/7ykyXHPfSZsXybp/EzGnUhPz0fSSEmPS2qS9P1Mx51ML87pXEkLJL0c/r4107F3pRfnc5KkxeH2oqS3Zzr2RHrz/yg8PzH82/tUpmJOphef0WRJu+M+px9lOvZEevldd6ykf0haEv4/lSY9mJn5rYsb0VQ5K4FDgUHAi8CMTmVuAH4U7r8X+E24PyOULwGmhP0U9uPzGQqcDnwE+H62P5s+OqfjgXHh/tHAhn5+PkOAonC/GtgSe9xfzynu+QeA3wKf6s/nA0wGXsn2OfTxORUBLwHHhccju/uu85pOYqks+nYZMCfcfwA4W5LC9vvMrMXM3gBWhP1lU4/Px8x2mtnTQHPmwk1Jb85pkZltDNuXAIMllWQk6sR6cz67zKw1bC8FcmWEUG/+HyHpcuANos8oF/TqfHJUb87pPOAlM3sRwMy2mVlbsoN50kms20Xf4suE//ANRJk+lddmWm/OJ1f11Tm9E1hoZi1pijNVvTofSSdLWgK8DHwkLgllU4/PSVIZ8GngSxmIM1W9/Tc3RdIiSU9KenO6g01Rb85pOmCS/ippoaT/6O5gA2rCT+c6k3QU8DWiX2z9mpk9Dxwl6UhgjqSHzSzXaqcH41bgW2bWlNsVhZTVABPNbJukE4AHJR1lZjuyHVgvFBE1vZ8I7ALmSVpgZvMSvcBrOomlsuhbRxlJRUAFsC3F12Zab84nV/XqnCQdAvwBuMrMVqY92u71yWdkZkuBJqK+qmzrzTmdDNwuaTXwceCzipYnyaYen09obt8GYGYLiPpRpqc94u715jNaDzxlZlvNbBfRLP6zkh3Mk05iqSz6Nhe4Otx/F/A3i3rT5gLvDSM+pgDTgBcyFHcivTmfXNXjc5JUCfwfcIuZPZOxiJPrzflMCV8GSJoEHAGszkzYSfX4nMzszWY22cwmA98GvmJm2R492ZvPqEpSIYCkQ4m+F1ZlKO5kevPd8FfgGElDwr+/M4BXkx4t2yMncvkGXAgsJ/pF8rmw7b+AS8P9UqJRNSuIksqhca/9XHjdMuCCbJ9LH5zPamA70S/o9XQa3dLfzgn4T2AnsDjuNrofn88HiDrbFwMLgcuzfS598e8ubh+3kgOj13r5Gb2z02d0SbbPpS8+I+D94bxeAW7v7lg+DY5zzrmM8eY155xzGeNJxznnXMZ40nHOOZcxnnScc85ljCcd55xzGeNJxznnXMZ40nHOOZcxnnSciyOpKcPHe7aP9tMW1mh5RdKfwowL3b2mUtINPTjW4DBhZWGSMoMkPRWbJcG5GE86zqWJIkn/j5nZqX10uN1mNtPMjiaaOeLGFF5TSbROysG6Fvi9JZnC3qIp8ucB7+nB/l0e86TjXDckvV/SC6Em8eO4+bMeVLTq6BJJ14dtk8MKjD8nmhbkzZKWSvpJKPeIpMFx+26Ke12X5SR9PuzzaUn3qvsVNP9B3NT0XcUZfBWYGs7r64nOswvvA/4Y9j1N0mpJh4XHxeH1E4AHQ1nnOnjScS6JsEzAe4DTzGwm0Ma+L9JrzewEYDZwk6TYminTgB+a2VHAmvD4B+FxPdEcXF05oJykE0P544ALwrGSxVsInM3+EzYmivMWYGU4r58mOc/4/Q8imndrNYCZvQ7cCcSWZP8oMNfM1hEl3ROTxesGHm9vdS65s4ETgH+GNV0GEy0FDdEX+NvD/QlESWMTsMbMnovbxxtmtjjcX0C0bHFXuio3CvijReviNEv6U4LXDpa0mKiGsxR4NO65ruLsvGRFsvOMN4ooIcZ7BThH0gjgOqIlCTCzNkl7JJWbWWOCuN0A40nHueQEzDGzz+y3UToTOAd4k5ntkvQE0Uy8EM1eHS9+RdI2oi/0rqRariu7zWympCFE083fCHy3mzjjdXmeXR2ni9cvD8e7FfhfM4s//xJyb5lzl0XevOZccvOAd0kaDSBpRFivpgKoC1/kRwCnpOn4zwCXSCpVtHzzxckKW7SQ1k3AJ+MW20oUZyNQHu4nOs/O+68DCiXFJ56VRAt3nQT8PLYxNONtNbO9B3XGLq950nFuf0MkrY/dgLcRrb3ziKSXiJqtqoG/AEWSlhJ1yD+XcI+9YGb/JOqfeQl4GHiZaH36ZK9ZFMpfmSxOi1axfEbSK8A1dH2eXXmEaIni2H72AjuIFsRrjyt3FtFCec518PV0nMtxksrMrCk0nT0FXG9mC7MYzyzg383sA3Hb1gKTLO4LRdLviRLR8iyE6XKU9+k4l/vulDSDqC9lTjYTDoCZLZT0uKTCMFhgMtHgifiEMwh40BOO68xrOs455zLG+3Scc85ljCcd55xzGeNJxznnXMZ40nHOOZcxnnScc85ljCcd55xzGeNJxznnXMb8f0g+1Vy82IHXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Learning Rate on Convergence\")\n",
    "plt.xlabel(\"Learning Rate ($\\gamma$)\")\n",
    "plt.ylabel(\"Number of iterations\\n\" \"required for convergence\\n\")\n",
    "plt.plot(\n",
    "    [result[\"gamma\"] for result in results],\n",
    "    [result[\"n_iter\"] for result in results],\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'gamma': 0.030000000000000002, 'n_iter': 7}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(results, key=lambda result: result[\"n_iter\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton_raphson_minimization(\n",
    "    df, d2f, x0, max_iter=1000, tol=10 ** -6, exact_iter=None\n",
    "):\n",
    "    \"\"\"\n",
    "    Performs Newton-Raphson minimization for a real-valued function of one variable.\n",
    "    Returns the minimizer, the number of iterations and the optimization path\n",
    "    as dictionary with keys 'argmin', 'n_iter' and 'path'.\n",
    "    \n",
    "    \n",
    "    Arguments:\n",
    "    -----------\n",
    "    df (function): The derivative of the objective function, a real-valued function of one variable\n",
    "    d2f (function): The second derivative of the objective function\n",
    "    gamma (float): The learning rate\n",
    "    x0 (float): The initial guess\n",
    "    max_iter (int): The maximum number of iterations\n",
    "    tol = 10**-6: The tolerance. The simulation stops when |x_{n} - x_{n-1}| <= tol and returns xn\n",
    "    exact_iter (int): The exact number of iterations for the algorithm to run for\n",
    "    \"\"\"\n",
    "\n",
    "    n_iter = 0\n",
    "    x_last = np.inf\n",
    "    x_current = x0\n",
    "    x_path = [x0]\n",
    "\n",
    "    while (\n",
    "        abs(x_current - x_last) > tol\n",
    "        if not exact_iter\n",
    "        else n_iter < exact_iter\n",
    "    ):\n",
    "\n",
    "        x_current, x_last = (\n",
    "            x_current - df(x_current) / d2f(x_current),\n",
    "            x_current,\n",
    "        )\n",
    "\n",
    "        x_path.append(x_current)\n",
    "\n",
    "        n_iter += 1\n",
    "\n",
    "        if n_iter >= max_iter:\n",
    "            warnings.warn(\n",
    "                \"Algorithm did not converge within the tolerance of {tol} in {max_iter} iterations.\"\n",
    "            )\n",
    "            break\n",
    "\n",
    "    return {\"argmin\": x_current, \"n_iter\": n_iter, \"path\": x_path}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def d2_pi(q):\n",
    "    \"\"\"Second devivative of profit.\"\"\"\n",
    "    return -6 * q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Profit is convex, so to maximize it, we must minimize the negitive of its derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'argmin': 5.477225575051661,\n",
       " 'n_iter': 5,\n",
       " 'path': [10,\n",
       "  6.5,\n",
       "  5.5576923076923075,\n",
       "  5.477808091562417,\n",
       "  5.477225606024404,\n",
       "  5.477225575051661]}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newton_raphson_minimization(\n",
    "    df=lambda q: -d_pi(q), d2f=lambda q: -d2_pi(q), x0=10\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice that Newton-Raphson almost always outperforms Gradient Descent regardless of the learning rate or initial guess. For example, try setting the initial guess to 100. For some learning rates, you may get an overflow error when computing $x_{n+1}$, this is because the learning rate is too large and makes the expression explode. Alternatively, if it runs, you may find that Gradient Descent does not converge within a reasonable number of iterations.\n",
    "\n",
    "One issue you may have with Newton-Raphson is a `ZeroDivisionError` error. This will occur because at 0, $\\pi''(0)=0$ and the iterative equation require division by $\\pi''(0)$ in this problem. You can solve this issue by perturbing the iterative equation slightly:\n",
    "$$x_{n+1} = x_{n} - \\frac{f(x_{n})}{f'(x_{n})+\\epsilon}$$\n",
    "where $\\epsilon$ is a very small number. When coding this, epsilon should be able to be defined by the user."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_iter = 15\n",
    "x = np.arange(0, exact_iter + 1)\n",
    "data = [\n",
    "    {\n",
    "        \"label\": \"Newton-Raphson\",\n",
    "        \"path\": newton_raphson_minimization(\n",
    "            df=lambda q: -d_pi(q),\n",
    "            d2f=lambda q: -d2_pi(q),\n",
    "            x0=10,\n",
    "            exact_iter=exact_iter,\n",
    "        )[\"path\"],\n",
    "    }\n",
    "]\n",
    "for gamma in np.arange(0.01, 0.06, 0.01):\n",
    "    data.append(\n",
    "        {\n",
    "            \"label\": f\"Gradient Descent ($\\gamma={gamma}$)\",\n",
    "            \"path\": gradient_descent(\n",
    "                df=lambda q: -d_pi(q),\n",
    "                gamma=gamma,\n",
    "                x0=10,\n",
    "                exact_iter=exact_iter,\n",
    "            )[\"path\"],\n",
    "        }\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'label': 'Newton-Raphson',\n",
       "  'path': [10,\n",
       "   6.5,\n",
       "   5.5576923076923075,\n",
       "   5.477808091562417,\n",
       "   5.477225606024404,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661]},\n",
       " {'label': 'Gradient Descent ($\\\\gamma=0.01$)',\n",
       "  'path': [10,\n",
       "   7.9,\n",
       "   6.9277,\n",
       "   6.3879091812999995,\n",
       "   6.063747670043894,\n",
       "   5.860676595865011,\n",
       "   5.7302506910254145,\n",
       "   5.645177501565497,\n",
       "   5.589136630840059,\n",
       "   5.551983182494114,\n",
       "   5.52724766473319,\n",
       "   5.510733662314236,\n",
       "   5.499688097405338,\n",
       "   5.492291022343079,\n",
       "   5.487333202119768,\n",
       "   5.4840084319871885]},\n",
       " {'label': 'Gradient Descent ($\\\\gamma=0.02$)',\n",
       "  'path': [10,\n",
       "   5.8,\n",
       "   5.5816,\n",
       "   5.5123444864,\n",
       "   5.4891879821953315,\n",
       "   5.481316899962672,\n",
       "   5.478626802493688,\n",
       "   5.477705704033555,\n",
       "   5.47739011723345,\n",
       "   5.47728196745145,\n",
       "   5.477244902393321,\n",
       "   5.477232199145703,\n",
       "   5.477227845344195,\n",
       "   5.477226353155366,\n",
       "   5.477225841733389,\n",
       "   5.477225666452267]},\n",
       " {'label': 'Gradient Descent ($\\\\gamma=0.03$)',\n",
       "  'path': [10,\n",
       "   3.7,\n",
       "   5.1678999999999995,\n",
       "   5.4642528631,\n",
       "   5.477027521429311,\n",
       "   5.477222779084842,\n",
       "   5.477225535629513,\n",
       "   5.477225574495832,\n",
       "   5.477225575043824,\n",
       "   5.477225575051551,\n",
       "   5.477225575051659,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661,\n",
       "   5.477225575051661]},\n",
       " {'label': 'Gradient Descent ($\\\\gamma=0.04$)',\n",
       "  'path': [10,\n",
       "   1.5999999999999996,\n",
       "   4.892799999999999,\n",
       "   5.6200609792,\n",
       "   5.429850730008825,\n",
       "   5.491857255987541,\n",
       "   5.472597721570381,\n",
       "   5.478678622913902,\n",
       "   5.47676828853305,\n",
       "   5.477369382179302,\n",
       "   5.47718034031907,\n",
       "   5.477239802673743,\n",
       "   5.477221099954523,\n",
       "   5.4772269826200795,\n",
       "   5.477225132323104,\n",
       "   5.477225714304883]},\n",
       " {'label': 'Gradient Descent ($\\\\gamma=0.05$)',\n",
       "  'path': [10,\n",
       "   -0.5,\n",
       "   3.9625000000000004,\n",
       "   6.1072890625,\n",
       "   5.012442108560181,\n",
       "   5.743755724810096,\n",
       "   5.295146250866783,\n",
       "   5.589360178156501,\n",
       "   5.40321809798225,\n",
       "   5.524003225827812,\n",
       "   5.446811479984401,\n",
       "   5.496648185210922,\n",
       "   5.464676994413544,\n",
       "   5.485272796422647,\n",
       "   5.472040148746498,\n",
       "   5.480556640322459]}]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOydd3hVRdrAf3N7kpteCR0SgtTQgqIgK3Z07bp2ZO24nwVdy9rL6u7qWlYs2F1dRbGigAo2itQQWgADkRIgkEbKzb25bb4/zsnNTSWQ3ETI/J7nPOecOXNm3nPLO3PeeecdIaVEoVAoFF0HQ2cLoFAoFIqORSl+hUKh6GIoxa9QKBRdDKX4FQqFoouhFL9CoVB0MZTiVygUii6GUvwKRYgRQiQLIX4WQlQKIZ45jPunCCEWh0i2XkKIKiGEsYU8UgiRFor6FZ2DUvxdHCHEZUKIVfqff68QYp4Q4oTOluso43qgGIiSUk5vLpMQ4mFdyY7tKMGklDullHYppU+X4UchxLUdVb+ic1CKvwsjhLgDeA74O5AM9AJeAs7pTLmCEUKYOluGdqA3kCtbmC0phBDAVUCpvg85R8lnqzgcpJRq64IbEA1UARe1kMeK1jDs0bfnAKt+bSJQAEwH9gN7gWv0a2OBQsAYVNZ5wDr92ADcA2wDSoCPgDj9Wh9AAn8GdgI/A0bgGbRe82/ALXoeU9CzvKHLsBt4vLZuYAqwGHgaKNPvPyNIrjjgLf35yoDPg66dBeQAB4ClwLAWPqtxwEqgXN+P09PfBjyAW/+8T27m/gmAE7hc/0wsQdemAIuDzk8Ftuh1vQT8BFwb9NneD+zQv5d3gegWPtvaNBPwBOADXLqsL+r3SeBGIE//LGYAIki2JcCz+rV8/bOYAuzSZbg6SPYzgVygUv+u7uzs/0JX3DpdALV10hcPpwPeWuXZTJ5HgWVAEpCoK7/H9GsT9fsfBcz6H7oaiNWvbwNOCSrrY+Ae/fhWvdweaI3Lq8AH+rVaRfQuEAGE6UonV88fCyygvuL/TC8jQpd1BXCDfm2KrnivQ2tAbkJT8rWK62tgll6uGThRTx+hK62x+n1XA9vRG74Gn1McWqNxpa5AL9XP4/XrbwOPH+T7eAOtATSjKf4Lgq5NQVf8QAJQAZyv13Wr/ny1in8qsBXoB9iBT4H/tvDZ1qbVfpY/1pYVVL8EvgJi0N4Ki4DTg2TzAtfon9PjaI3KDP27PRVNydv1/HuB8fpxLDCys/8LXXHrdAHU1klfvNazLDxInm3AmUHnpwHb9eOJaD1UU9D1/cCx+vHjwJv6cSTgAHrr55uASUH3ddOVlylIEfULuv49uiLXz0+mrpeaDNQAYUHXLwV+0I+nAFuDroXr96bo9frRG6sGz/4yeiMXlLYFvWFokH4lsKJB2i/AFP34bVpQ/LpMFcC5+vmrwBdB16dQp/ivAn4JuibQeta1in8hcHPQ9YyDfLa1aQdT/CcEnX9EXSM+BcgLujZUz58clFYCZOrHO4Eb0MY7Ov1/0FU3ZePvupQACQex86aimQxq2aGnBcqQUnqDzqvRepkA/wPOF0JY0Xqn2VLK2rJ6A58JIQ4IIQ6gNQQ+NCVey64Gcuxq5lpvtF7y3qDyXkXr+ddSWHsgpazWD+1AT6BUSlnWxLP3BqbXlqmX27PB8wfLt6NB2g6gexN5m+I8tF7zXP38feAMIURiM3UFnl9q2rSgBVl2UNdA1hL8+bWWwqDj4O8ZYF/QsVOXq2Fabf4L0N4OdwghfhJCHHcYsijaiFL8XZdf0HrK57aQZw+aAqyll552UKSUuWhK5wzgMrSGoJZdaHb2mKDNJqXcHVxE0PFeNDNPLT0blFUDJASVFSWlHNwKMXcBcUKImGauPdFAxnAp5QdN5G34OYH2We1uIm9TXI2mGHcKIQrRzGJmtM+tIfU+C31QOPizaeo781JfObcUkjek4XqllCullOegNcyfo709KDoYpfi7KFLKcuBBYIYQ4lwhRLgQwiyEOEMI8U892wfA/UKIRCFEgp7/vUOo5n9oNugJaMqslleAJ4QQvQH08lvyJPoIuFUI0V1X0ncHPcde4FvgGSFElBDCIIToL4Q48WDC6ffOA14SQsTqzz9Bv/wacKMQYqzQiBBCTBZCRDZR1FxggO4aaxJCXAIMQrOLt4gQojswCW0gOVPfhgP/oGnvnq+Bofp3ZgKmoZmtavkAuF0I0VcIYUfz2JrV4M2sJfahjQ+0O0IIixDiciFEtJTSg2be8oeiLkXLKMXfhZFSPgPcgeYFUoTWy70FrScGmp1+FbAOWA9k62mt5QPgROB7KWVxUPrzwJfAt0KISrSB3pZ8119DU+7rgDVoitaLZh4CTUFa0AaAy4DZaPb71nAlmg18M9oYxW0AUspVaAPCL+plbkWzZzdCSlmCprino5nQ/gqc1eCZW6o/R0r5rZSysHYDXgCGCSGGNKirGLgI+Kde1yC076hGz/Im8F80j53f0Dx0/tIKOWp5HrhQCFEmhHjhEO5rLVcC24UQFWiD9peHoA7FQaj1bFAojhiEEGcAr0gpG5pXuhxCCAOajf9yKeUPnS2P4shA9fgVv3uEEGFCiDN1M0p34CE0F84uiRDiNCFEjD5wfh+aZ8+yThZLcQShFL/iSEAAj6CZXNageQE92KkSdS7HobnaFgNno7mBOjtXJMWRhDL1KBQKRRdD9fgVCoWii3FEBGlKSEiQffr06WwxFAqF4ohi9erVxVLKRhMBjwjF36dPH1atWtXZYigUCsURhRCi4YxyQJl6FAqFosuhFL9CoVB0MZTiVygUii6GUvwKhULRxVCKX6FQKLoYIVP8Qog3hRD7hRAbgtLihBDfCSHy9H1sqOpXKBQKRdOEssf/NtryfsHcAyyUUqajrRR0TwjrVygUCkUThEzxSyl/BkobJJ8DvKMfv0PLi4C0mddunsYr194QyioUCoXiiKOjbfzJ+uIXoC3lltxcRiHE9UKIVUKIVUVFRYdVWXVZNc6qKnyVrQmLrlAoFF2DThvc1dcKbTZCnJRyppRytJRydGJiU0uPHhyj2YKfGnLmdOVAjgqFQlGfjlb8+4QQ3QD0/f5QVhadkAjSRX7OYqgKaVUKhUJxxNDRiv9LtIWl0fdfhLKybn37A7C/IgUWPxfKqhQKheKIIZTunB8AvwAZQogCIcSfgaeAU4QQecDJ+nnISMzQFH9NdTz+la9DxZ5QVqdQKBRHBCGLzimlvLSZS5NCVWdDEvumAhBWE8FWo2DAomdg8jMdVb1CoVD8LjmqZ+7GpiQBYPJbeC9yOKx+Bw7s7GSpFAqFonM5qhW/zR4JGPEJI995wpDCAD/9s7PFUigUik7lqFb8QgjMlgg8RkiuLuDXXhdCzv+gZFtni6ZQKBSdxlGt+AGsETF4hYdeJTX82z8ajBb46R+dLZZCoVB0Gke94rfHJSKlg94lUSwqzcM7+lpY9xEUbels0RQKhaJTOOoVf2RCPNJfRWpVN6R1K790uxwsEfDjk50tmkKhUHQKR73ij0lJAjxEOWMxRvzGV1udcOxNsPEzKFzf2eIpFApFh3PUK/74VM2ls0bYiHS5WJCfjf/YaWCLhh9Ur1+hUHQ9jnrFH5WYAEC1NYKexVDJZnKKgeP+Alu+ht2rO1dAhUKh6GCOesUfERsPgNNiYWRVAqaIbXyXuw+OvRHC4uCHv3eyhAqFQtGxHPWK3x4XB4DbbGRYVSym8B18l1sA1kg44TbYugB2LutkKRUKhaLjOOoVv8UWhtFswyO8dNvvRwoPv1VuYnuxA8ZcBxFJ8P3jnS2mQqFQdBhHveIHsNljkNKBLPRgkAJjrbnHEg7jp8P2RZD/U2eLqVAoFB1Cl1D8EbFxSL8Dp9fKGEs6kTHbNcUPMGoKRHWHH54A2eyCYAqFQnHU0CUUf1RCAlJW4bTFM97dG695B6t2FlLqcIPZpvX6dy2HrQs7W1SFQqEIOV1C8cekJIK/CqctjsEVUfjxIsJ2sHCT3usfcSXE9IIfHle9foVCcdTTJRR/ZFw84Kc6KpHkQhcmYSI6djsLahW/yQIn3g171sCWuZ0qq0KhUISaLqH47bGaS6fTHol3228MTRyKPWY7P/9ajMvj0zIN+xPE9df8+v3+TpRWoVAoQkuXUPyBSVxC4M7bSlbyGA74fsPpc7Bka7GWyWiCiffCvg2Q+3knSqtQKBShpUsofs3UA26fC4/Ly1jRH4kfe/TOOu8egCHnQ+JA+PEp8Ps6SVqFQqEILV1C8YfHxGoHfgdOWzxpZWYsBgs9UnezYNN+/H59QNdghD/cB8VbYP3szhNYoVAoQkiXUPwmsxlrRCTSX4UrLB7ftu2MSBqB37KV4qoa1uw6UJd54NmQMlSL1+/zdJ7QCoVCESK6hOIHsMfGI2UV7qS+1OTlMSZlDHtd2zCZquu8ewAMBvjD36DsN1j7QecJrFAoFCGiyyj+yIR4kA7ciX2oydvK2G5jARjYd399Oz/AgNOh+yj46Z/grekEaRUKhSJ0dBnFb4+NBxy47Mm4t21jUOxAwkxhxMTvZOv+Kn4rdtRlFkLr9Zfvgux3O01mhUKhCAVdR/HHxeH3Oqg22JFuN7JgL6OSR1HqywXgu9zC+jf0Pwl6jYNFz4DH2QkSKxQKRWjoOoo/Nh6QVDk1001NXh5ZKVnsqtrOgFTZ2NwjBJz0N6jcC6ve7HiBFQqFIkR0HcVfuyBLTRUec7im+LtlAZDeex+rd5RRUtXAnt/nBOh7Iix+FtyOhkUqFArFEUnXUfz67F3pr8LbezA1eVsZGDuQSEskxvBt+CV8v3l/4xtPuh8cRbBiZgdLrFAoFKGhUxS/EOJ2IcRGIcQGIcQHQghbqOu0x9Upfk/PDGq25mE0GBmdPJq8ijV0i7Y1NvcA9MyC9FNhyfPgqgi1mAqFQhFyOlzxCyG6A/8HjJZSDgGMwJ9CXW9YVBTCYEBKB+6E3ri378DvdjO221gKqgoYl2FgUV5Q0LZg/nAfOMtg2cuhFlOhUChCTmeZekxAmBDCBIQDe0JdocFgJCImFoOoxhWeBF4v7t+2k5Wi2fmTkgpwenwszitufHPqCBh4FvzyIlSXhlpUhUKhCCkdrvillLuBp4GdwF6gXEr5bcN8QojrhRCrhBCrioqK2qVue2wcBqPm0gmaZ09aTBpxtjjK/LlEWk1Nm3tA6/XXVGrKX6FQKI5gOsPUEwucA/QFUoEIIcQVDfNJKWdKKUdLKUcnJia2S932OG32rqNagMlETV4eQgjGpIxh5b4VTBiQwMLN+/D5m1iFK3mwFr1z2Suwa0W7yKNQKBSdQWeYek4GfpNSFkkpPcCnwLiOqDgiNh6fp5LK0hrMvXtTs3UrAFkpWeyv3s/INC/FVW5ygoO2BXPKoxCZDO/8EbbM7wiRFQqFot3pDMW/EzhWCBEuhBDAJGBTR1Rsj43D667G43ZD2hBq8vIAAnF7TBH5mAyieXNPdA+Y+i0kZsCHl8Ga9zpCbIVCoWhXOsPGvxyYDWQD63UZOsRJvp5LZ/cBeHbtwl9dTa/IXiSFJ7G+ZDVj+8U1Dt9Qr5BEmPIV9J0AX0yDn59WC7QrFIojik7x6pFSPiSlHCilHCKlvFJK2SEhMGvX3kU6cMf3BCmp2ZaPEIKxKWNZWbiSkwcmsa3IQX5RVfMFWSPhso9g6EXw/WMw769qxS6FQnHE0GVm7kKd4tcWZEkACJh7srplUeoqJb2nFpqhWXNPLSYLnDcTjrtFm9U7e6oK4axQKI4IupTij9BNPSazE4fXirBY6g3wAuyoXsegblEHV/ygLdpy2hNw6uPaAu3vXQCu8pDJr1AoFO1Bl1L8tgg7JrMFs8VFZWkNlv79Az3+VHsqPew9WL53OacMSmb1ziaCtjXHuL9ovf+dv8Bbk6GyhTEChUKh6GS6lOIXQhARF4cwVFNR7MSanhZQ/KB596wqXMWkYxKQEhY2FbStOYZfApfOgtJ8eOMUKN4agidQKBSKttOlFD9odn7pr6KqtAZzWjrewkJ8FVrwtayULCo9lRhte0ltLmhbS6SfDFPmaCGc3zwVdq8OwRMoFApF2+iCij8eb00Ffr/E1z0NIGDnH5MyBoAVhSs4eVAyi/KKcLoP0Vun+yj483dgscPbZ8PWBe0qv0KhULSVrqf44+JwV5cjpcQd0wOAmjxN8SeGJ9Ivuh/LCzU7v8vjZ/HWJoK2HYz4/pryj+sH/7sE1s5qz0dQKBSKNtEqxS+E6C2EOFk/DhNCRIZWrNARERuP11MDuKkmHEN4eD07f1ZKFtn7shnZO0oP2naYA7WRyXDN19DrOPjseljyQvs8gEKhULSRgyp+IcR1aDNtX9WTegCfh1KoUBLw5ZdVVJTWYGkwwJvVLQun10negU1MHJjEwk37mw7a1hps0XDFJzD4PPjuAfjmb+D3t8djKBQKxWHTmh7/NOB4oAJASpkHJIVSqFBSG7bBFu6mstiFNT29nuIfk6zZ+WvdOkscbnJ2lR1+hSYrXPAmZN2ghXT+7Hrwutv0DAqFQtEWWqP4a6SUAU2lL55yxAanqe3xW8NqqChxYktPx1dairekBIAYWwwD4waysnAlJw5IxGQQfHuo3j0NMRjgjH/ApAdh/cfwwSVabH+FQqHoBFqj+H8SQtyHtmLWKcDHwJzQihU6ahddN5mcVBS7sKTpnj15dX73Y1LGsGb/GmwWP8f2iz90t86mEALGT4dzZkD+T/D2WVDVPgvMKBQKxaHQGsV/D1CEFknzBmCulPJvIZUqhJhtNqzhEQjhwFFeg6lvf4D6E7lSxuL2u1m7fy2nDEomv8jBtpaCth0KI66AP/0PirZoE71K89unXIVCoWglrVH8f5FSvialvEhKeaGU8jUhxK0hlyyERMTG4fdVgQSnMRJjdHQ9xT8qeRRGYWR54XJOHpQMtCJo26GQcTpc/SW4DsAbp0Lul2rQV6FQdBitUfxXN5E2pZ3l6FDssXF4a7TZupUljQd47RY7g+MHs7JwJd1jwhicGsWC9lT8AD2zYOo3EBYLH10Jr06ATXNUA6BQKEJOs4pfCHGpEGIO0FcI8WXQ9gNQ2nEitj/22DhcDm15xYpiF9YB6dRs3YoMWlBlTMoY1hetp9pTHQjaVtzaoG2tJTEDbvpFC/DmqYZZV6gGQKFQhJyWevxLgWeAzfq+dpsOnBZ60UKHPS6e6vIDCANUljixpKXhr6zEu6+uV5/VLQuv9JK9P5uTj0lGSvh+0yEEbWstRpMW4G3aivoNwMwJsOkrtbqXQqFod5pV/FLKHVLKH6WUx0kpfwrasqWU3o4Usr2JiI3H7/Nij/ZTUezClp4O1B/gHZE0ApPBxIrCFQxOjaJ7TFjb3Tpbol4D8Cq4q2HW5fDqeNUAKBSKdqU1M3ePFUKsFEJUCSHcQgifEKKiI4QLFfY4zZffZndTUeysc+n8tU7xh5nCGJYwjBV7VyCE4ORjkli89TCCth0qRhMM/1NQA+CoawA2f60aAIVC0WZaM7j7InApkAeEAdcCM0IpVKip9eW32GqoLHVhio3FlJhYr8cPWnz+TaWbqHBXcMqgFFweP4vyOsj3PtAArIRzX9EagA8vUw2AQqFoM60K0ial3AoYpZQ+KeVbwOmhFSu01Pb4jWYnzkoPbpdX8+zZWn/xlKyULPzSz+rC1YztF0ekzcSCTSE09zSF0QSZlzbRAExQDYBCoTgsWqP4q4UQFiBHCPFPIcTtrbzvd0tETCwAQmoLq2sunWmaZ0+QN82wxGFYjVZWFK7AbDTwhwwtaJujphOGOBo1AFVBDcBc1QAoFIpW0xoFfiVgBG4BHEBP4IJQChVqjCYzYVHR+HxavJwK3Zdfulx4CgoC+SxGCyOSRrC8cDkAl4zpSVm1m0tm/sL+SlenyF6/AXhZi/nz4aWqAVAoFK3moIpf9+5xSikrpJSPSCnv0E0/RzT2uHg8Tm2MWlt/t7FnD2h2/ryyPEpdpRyflsDrV49m234H581YSt6+Tgy0ZjRB5mVwy6r6DcDMEyHnA6g+oqdaKBSKENIar56zhBBrhBClQogKIUTlke7VA9okLmdlGSaLgcpiF5b+tcHa6iv+rJQsAFYWrgTgpIHJzLrhWGq8fi54eSnL8ks6VvCGBDcA57wErgr4/Eb4V39483RY/Bzs36zeBBQKRYDWmHqeQwvbEC+ljJJSRkopo0IsV8ixx8bhKCslMj6MihInRnsE5u7d60XpBBgUP4gIcwQr9q4IpA3rEcNnN48jKcrGVW+s4Iuc3R0tfmOMJhhxOfwlG677HsbfqY0DLHgIXhoLL2TCvHsg/0fweTpbWoVC0YmYWpFnF7BByqOryxgRG4+j/ADdh5ipKNHs9da0tEY9fpPBxKjkUawoXFEvvWdcOJ/cOI7r/7uKWz/MYfcBJzed2B8hRIc9Q5MYDNqC791HwUl/g/IC+HU+bJkPq96E5S+DNQrSJsGAMyD9FAiP61yZFQpFh9Iaxf9XYK4Q4icgEKxGSvnvkEnVAUTGxYOUhNm9FG5zIqXEOiCdqqVLkR4PwmwO5M1KyeLngp/Z59hHckRyID063My7f87ir7PX8c/5W9hV6uSxcwZjMv6OnJ6ie8CYa7XN7YBtP8Cv8+DXb2HjZyAM0PNYLWLogDMgIV1bO0ChUBy1tEbxPwFUATbA0h6VCiFigNeBIWireU2VUv7SHmW3lgh9JS6z1YXb5aOmWvPlx+PBvWMHVn02L9TZ+VcUruDs/mfXK8dqMvLsxZn0iA1jxg/bKCx38uJlI4mwtuajbX+cVZWYLBbMFmvji5YIOOYsbfP7Yc8arRHYMh++e1Db4vppDUDG6dpC8UZz43IUCsURTWu0U6qUckg71/s8MF9KeaE+RyC8ncs/KLVLMBqNTsBMRbGTqCDPnmDFnxGXQZQlqknFD2AwCO46bSDdY8J54IsNXDLzF968egxJUbYOeZZavB4P791zG0aTiYsfejLwjE1iMECPUdp20v1wYJdmEvp1Pqx8DZbN0BaLTztZ27oNh4QBqiFQKI4CWqP45wohTpVSftseFQohooEJ6DH99fV8O3z18dpF16WsAmKpKHaRMKQfGAzaAO8ZdXkNwkBWSla9Ad6muGxsL7rF2Jj2fjbnvbSUt64Zw4DkyBA+RX02fP8tFUX7MJpMfPzofVz80JOByWoHJaYnZF2nbTVVkP+D9iaQ9w1s+ETLY7RC0jGQMhRShkG3YZA8GKwd94wKhaLttMYYfRMwXwjhbCd3zr5oSzm+pbuJvi6EiGhDeYdFeFQ0wmDA59F9+UucGKxWLL16NRrgBS1M8x7HHgoqCxpdC+YPGUl8dMNxuH2au+fSbcUhkb8hXreb5Z9/RGrGIC782+NUlBTx8WN/o7r8wKEXZrXDMWfDuTNg+ha4eRmc/xqMvR7CYrRQEfPugjdPgyd7wAsj4KOr4OentbGDir3KfVSh+B1z0B6/lLK9u3MmYCTako7LhRDPo63r+0BwJiHE9cD1AL169WpnEUAYDETExuGqLMcabqKyWPfsabAaVy3Bdv4ekT1aLHtI92g+u3kc17y1kqvfXMG/LhzOuSO6t/szBLP+h2+pKi3h9Jtvp8egIZx/90N8+tQjfPz4/Vz0wBOER0UfXsEGo9bLTzoGhl2spUkJlXth7zooXA+F67Tj3C/q7gtP0N4Iat8OUoZCfJpWnkKh6FSaVfxCiIFSys1CiJFNXZdSZh9mnQVAgZRyuX4+G03xNyx/JjATYPTo0SHpPtpj46gqKyEy3lbn0pmeTuXChfhrajBY6wZI+0X3I94Wz/K9yzk//fyDlt0jNpzZN43jhv+u4rZZmrvnzRND4+7pdbtZ8dlHdB84iF5DhgPQc/Awzv3rA3z+j0eZ/fj9XPTg3wmzt1MbLgREpWpbRlC8PlcF7NtQ1xgUrodlL4NPt+SZwjTTUPJgiO0DMb0guqdmZrKnaOMOCoUi5LTU478Drcf9TBPXJHDS4VQopSwUQuwSQmRIKbcAk4Dcwymrrdhj4ynbu5vkAWGU7dUCtlkHpIPfjzs/H9sxxwTyCiHI6pbF8r3LWVu0lozYDGymlgdvo8PMvDM1i7tnr+Nf32yhoKyax84Z0u7unusWfkNVWSln3DK9XsPSe2gm59x1P5//6zFN+d//BDa7vV3rroctCnqP07ZafB4o2qI3BnqDsGkOOBuElDCYIbq73hAENQi1+6geYGoXpzKFosvTrOKXUl6vH54hpawXkUwI0VZ3lb8A7+sePfnANW0s77Cwx8VRkLueqHgbOzaUaL78QZ49wYof4KSeJzHvt3lcMfcKjMJIWkwagxMGMzhe29Jj07EY6ysnq8nIs5dk0iM2nBd/2MqeAy5mXD4Sezu5e3rdblZ88TE9jhlCz8HDGl3vM3wkf5x+H18+/QSzn3iAC+9/DFtECJV/Q4xmSBmibVyK1+1mV+56uvXqjs1bpnkTle/U97u0/bbvobIQrX9Ri4DIlMYNQnQv/Y0hGWwx6q1BoWgFrdE+S9Fs8gdLazVSyhxg9OHe317YY+NxOaqIiDHi8/iprnAT3qsXmM1N2vlP73s6mUmZbCzZyMbijeSW5PL9zu/5NO9TQJvlmx6TXq8xSItNw2wwc+dpGfSIDeNvn2/g4ld+4a1rxpDcDu6e6xbOx1FWyuS/3NmsGanfiDGcfce9fPnMk3z694e44G+PYQ3vcA9aCrflMf+lZykp2InJYmXg8SeSedpkksec2jiz1w0Vu+sag8B+J+zOhtwvwd8g9IQwQFgchMfrWxxEJASdx2tjD+FBeSwRasKaosvRko0/BegOhAkhRgC1/44oOsHvPhTUunSazNoLTUWxi4j+0Vj79q23DGMwKREppESkMKnXJACklOxx7CG3JJeNxRvZWLKRb7Z/w+xfZwNgMVjIiMtgUPwgBscP5tELU/j756WcN2MJb0/NapO7p8ddw4rPP6bHoLrevpQSJAhDfWXWf9RYzrr9br569ik+ffIhLrjvESxhHfM1+rweln36Ecs/m0VETCyn33w7u7fksmnxj2z44Vu6pWeQedpZDBh7PCaL/sZkskBcX21rCr8fqpyQYvAAACAASURBVPbpDcJOcBRBdUnQVgol22DXCu1cNrNkpslW10gENw5hMWCxa66q1kj92B50rO/NYarhUBxxiOZC8AghrkbztR8NrKRO8VcCb0spP+0IAUEb3F21alW7l7t93Ro+eeIBzrjlYX54v4KTrxlExtgUdk+/E+eaNaR9v/CwypVSUlBZoL0Z6FtuSS4OjzaOYDFY8Tq74Xf14OpRxzOuTz+irFFEWiIDm9lw8IlS2XO/4Id3XuPih56k56ChFGwu5Yf3NlNd6SGuWwTx3SOIT7Vr+x52wuwWfl2+hK+e+wepAwZy/r2PYLGFHdYztpainduZP+NZ9m/fxqAJJ/GHKdcHTE0uRxW5Py0k59u5lO3dTVhkFENPOpVhJ59BdFLyQUo+BKQEV3mDhqEEHMV1jYSeVl1eys79bmpcbuKt1cRbqwkztrDwjjBqDYJFbyCs9qBGIqruuLaRMNm0zWzTBrvr7WuvhdXt1YQ5RRsQQqyWUjayrjSr+INuvEBK+UnIJGsFoVL8xbt28M6d0zh92p38+D8/Y//Yj9Fn9qH4lVcpeu45BqxahdHePlMM/NLPzoqdgYZgzb71bCzORYqm566ZhY0wo50Is51ISxTR1ijiwqKJC4smxhqNXYRR+sJcwpITGH7DNez7QVKwzEFUoo1eg+Mp2+ugZLcDV1WdOSQ8ykJ89wikP4+ty/9Lct+BXHDvw4RFtf80Cr/Px8o5n7L0o/ex2e2cct0tpI05tsm8Ukp2rl9LzrdfsW3VCiSSfiPHkHnqZPoMG4EIod3e6/Gw99dNbF+3hh3r1rDvt22N5iDYo6OIT0kkITGGhPhIEuLCiI80YPY7tQioNVXaegjuyqBjPd1dqZ1LfzMSHARhrNdQ+Iw2HP5wqnxWTGYLNpuJMJsFk8WCMFnAaNEaC2PQscHcdHqjY7NWn8EIBlPQ1tT5QfIIoxpv+R3QnOJvjY2/hxAiCq2n/xqabf+e9prJ25nULrrurCgjLCqRihInoHv2AO5tWwkbPrxd6jIIA32i+9Anug+T+00GoNTh4snvFrOzvIhyVwWVngocnkqc/ircwkm10UmJwYUwViGMRQiDC2F0IowuBv0WSVZVHAv77qbk39uJqoljfbefyUmZQ8Q+D0YMmFINhHtjiK3pTlRNKvaaFIrykwlzJ2IOP519+fN5+aY7IeZEvJYDeM2l+MylSHMpwlCFEQNGacAoBQiBQATGEYQ+909PBQwIaUb4rcgaP969hUhXNYaIeHzRGXzz8XK++SgbpBnh1392wg9CotmmJIgErKkT8VbtIH9NDvmrVyBM4ZiiemCOTEGYTIG8grp7hZAghDYULOvMLiLwkioC48RSgt9Tjc9Rhre6FJ/zgK6UBUZbNJbYvpjC4hFGK353FX53FS6Pg13b97Pj19/qKXBhDsNoicBgsWOwRGOwpGKwhCNEY4Un8CGkBPwY/D5NbukH6Ufg1675ffh9NUivG7+3Bp/Pjd/rwe/34vN58Pm8+P0NTVZ+wAW4MAoDBoPAIARGUbsHowBD0LnBUHddQNDYUP0GT0qJBPxSap9b0LlfgkTfNzivzS+1Tx4hBELovxV9j9B+S7W/KS2foU4eIRCG4OuGwFcrJdr3rTfQfi0BqcuvJUskQn8GTRj9qv47CX5a/ZeiVR6UIqi9oO1EUH6B9rcI+r0Ffe91fQcZlCYD99dJUydr/Q5HnZx/uPE60kcdR3vSGsU/VUr5vBDiNCAebSnG/wJHvOK3RkRgslipKislKr4nFcV14ZlB8+xpL8XfFHERNv517smN0qWUVLi8HKh2U1btoazarR07PByodlNaXklk/r9wh8cyoeAiXAYXpVVzOWXhKm4pqsEkJVCrIPbpW920C4mgOiyRnYmJ/BazC+v+r4kOOwNXWJo2QAoYfTVEOPYS4diNzXUAr8mK1xiGzxSG12TDa7ThNYXh1XuhPpMNKSW+mjV4nYtBGDFHnInBnIGo1P9Mfh8mnxOjr0b7UwoDUhi1zWAMnJsMmRijfPg9eXhr1uIp/RVP6TaMloEYrZkYTIdmBpJ+J37vTvye7fg8O0BWafIYYjBahmIw9cFg7oEQVpDgr667VwAmE2ACGeZH+suRvmKkrxi/rxifswSvYwd1f3ADwhCLMCZgMCYgjAkIYzzCEAXSgfRXIf2V+hZ8XAnSQUPlC2aEIRJhiEOYIjEaIhEGO8JgB+lDShdIF9LvBOnCJ114pQv8dWl1v4WmMIKwIYRNU2zSg5QekB7gUNdtsIAwI4QZhAUtMIAP/D7Ah5RewKuPt/iaeFZFU+RlZ3eK4q9t0s4E3pVSbhSdHnS+fRBCBBZkiUoIY99v5QCYe/RA2GxNevZ0lFzRYWaiw8z0jm98feHLr5NT48Bin0zPouX02/wRJl8NlkGDMJx5NjK1JxgN2qu20aSZSoyaUtXSjYQZDMQbDcRtzWX1T3OJSV7GuNMuxuUwU1nup/KAjcoyOyVlfXG7JEaTwGQRmCwGTBaB2WIgTN+bLAK/r5w9uV9S49xOTLcBDDjuXMKiI3ULgsBgAYPBDwL8gT+8DPSCtL67fuz34/dLpH8EUmZSuW8fezZks//XDbgrN2JPTCVp4Ajieg8EYcLvl/rnppXq9/twFO+mYu92KvZsp7qsEACjxUpsr95Ed+9DVGofrJHRmq5DondCg76EgIhBGND6PvFARuCaz+vFVV5KdWkxzrJinAeKcZbtxV21pcXv2WA0YY6IxBIRiSWit7YPtwelRWK0Wts86c/v9eCtcQU2X40z6NyJr8aF1+VEoq1HbTBb9L0Zg8mC0WzGoJ8bTRZ9r+UzmMwYzWaE0XTIckq/H7/Pi/T56vZ+H9Lrxe/zIf36vva636+/CRi0Teh7g7EuLfhcGBqk69f09Nq3Ds0hwo/0+/VjiZT6sT/oWPqR/obXtb2Ufj1daj3/2rcB/a0Gg6HurUHUvjmLQJ5AvqD70HIwasLxbfr+m6I1in+1EOJbtBg79wohItHeL48KIvTZu7172di6ej9+nx+D0dDkoiydhd/ppHrlSioWLyV7vY8tci1GQzfG7P6eXpndiJjyOBHHj8MUd+gLqvQ9dSJxGf35buZ/2Lx6IWffcR+moLUIpJRIv8TQzKQzKSXrF37Dj/99AyHgtBtvZfDEk9t/hvIF51FT7WDjT9+z9tuvyV/0NXtzFjHkpFMZefLpeD0eduh2+l0b1+OpcWEwGumWPpDMSZPoM2wkyf3TMHRgyIia6mpKCnZQvGsHlSUlRMTEEhmfENhs9sjOX7hH0SVpjeL/M5AJ5Espq4UQ8XTShKtQYI+LZ19+HpHxNqRfUlVWQ1RCGNb0dByLF3eKTFJKarZswbF4MVVLluBctZpyawqbBl1NuXU3OKs588prSD/z9HZRHMMmnYb0+1nw+gy+eu4fnH373RhNmvIXQiCMTddRWVrMt6+8wPa12fQaMozTbryNqMSkNsvTHNbwCEaecTYjTj+LXRvXkfPN16ya8ykrv5gdyBPbLZXBEyfRe9hIeg4a2inzFerkDSd1wDGkDjjm4JkVig6kJT/+K6SU70kp/UKIMCnlAQApZYkQ4lJgXYdJGULssXFsKyslMl6bTFVR4goo/vLPPsNbVoYptpWhjduAt7gYx9KlOJYsoWrJUnzFWlRPc3oGeyb/lc3lqYTZwVj8Fd2HjWDA5DMOUuKhMfyUM/D7fXz/5it8/fy/mHzrXzGamv55SCnZtPhHvn/rFXxeLydNvZHMU84MqfdNMEIIeg0ZTq8hw6ksKWbT4h+x2e30Hjqifd1AFYqjlIPF6nlPP/4P9WfqTgVeDJVQHYk9Ng5vTQ1hds1gW1niBGKxpmsDvO6tWzGNGROSuqXHQ/GrM6lcuJCaTZsAMMbGEjFuHBHHH4+r/yh++qqQ4l1VZIxNITxyA798XMlxF14WEnlGnHYW0ufjh3deY+5/nmby/92FwVjfNOI4UMaC12ewdeUyUjMGcfrNtxGbkhoSeVpDZHwCWedc2OZyPB4PBQUFuFyug2dWKH5n2Gw2evTogdncunkfLSn+poa6mjs/YomIqx09dSAEdZ49esweV14e4SFS/EUvzqDk1VcJHz2axNtvJ+L447ENOgYpIfubnaycuQ1ruIkzbhxKz4GRvHbL4/QeNoLuGaEzHYw88xz8fj8//fcNhMHAmbdMDyj/X5ctZsHrL+F2OZlwxVRGTT6nQ23moaSgoIDIyEj69Omj7O6KIwopJSUlJRQUFNC3bzMz3RvQkuKXzRw3dX7EEqn78lcfKCMi1hrw5TclJ2OIjAzZAK9j6VJKZs4k+sILSH388UB66R4HC9/JZf+OStJGJzHhTwMIs1tY8cVsnJUVjLsoNL39YEafdR5+n49F/3sbg9HIxKuu5Ye3Z7J5yU8k90vjjGl3EN+j/ddI6ExcLpdS+oojEiEE8fHxFBUVtfqelhT/QCHEOrTefX/9GP283+GL+fsiQveEqSorISo+NrAgixACa3o67ryt7V6nt7iY3XffjaVfP1Luuw8Av1+Ss2AnK778DbPVyGnXDSFtlDZQ6nY5WTXnU/oMH9lhA4VZ51yI9PtZ/OG7/LpsMdLvZ9zFl5N1zkXN2v6PdJTSVxypHOpvt6V/cJdwRahdkLyqrJSohG7syq2LE29NT6dy/nzdN7d9lIL0+9lzz734yyvo9fobGMLDObCvmoXv5FKYX0Hf4QlMvHwg4VF14Z1zvvla7+1f3i4ytJax512MMBjIz17JH6ZcT3Lf/h1av0KhCA3NumFIKXe0tHWkkKHEbLVhjYigqrSEyPgwHOVuvB5tpqM1LQ1feTneQ3iFOhilb72FY/Fiku+7F1vGANb/WMCsx1dQVljNydcM4owbh9ZT+m5nNSvnfErfzFF0S89oNzlaS9Y5F/KnR/6hlH4HIIRg+vTpgfOnn36ahx9+uE1l5uTkMHfu3DZKVsfbb79NYmIimZmZDBw4kGefffawy3r44Yd5+umn2002RetRUZTQYvZos3c1l87KkvoDvO1l53euXcv+Z58j8tRTibnkEop2VvLzh7+SOiCGSx8cS8bYlEZvFmvmf4WrsoLjOsC2r+hcrFYrn376KcW6K2970N6KH+CSSy4hJyeHJUuW8MQTT7Br1652LV8Reo5OY+0hUjt7NypeC1FcWeIiNiUiEKytJi8P+/Ftmzbtq6hg9x3TMScl0e2xRxFCsHX1fgwGwSlTB2OLaOyG5XZWs+qrz+g7YjTd0jq+t99VeWTORnL3VLRrmYNSo3jo7MEt5jGZTFx//fU8++yzPPHEE/WuFRUVceONN7Jz504AnnvuOY4//niGDh3KokWLiI6OJiEhgWeffZarrrqKq666issvv5wHH3wQp9PJ4sWLuffeeznllFOYOnUq+fn5hIeHM3PmTIYNG8bDDz/Mzp07yc/PZ+fOndx222383//9X4vyxsfHk5aWxt69e+nZsyePPvooc+bMwel0Mm7cOF599VWEEEycOJHhw4fz008/4fV6efPNN8nKygIgNzeXiRMn1qvT4XBw8cUXU1BQgM/n44EHHuCSSy5h4cKF3HnnnXi9XsaMGcPLL7+M1WqlT58+XH311cyZMwePx8PHH3/MwIED2/BtHf202OMXQhiFEO93lDCdhbboel2Pv3bhdVNcHMb4eGq2tm2AV0rJ3gcfwlNYSOozT2OMjkZKybbs/XQfGNuk0ge9t19VybgQ+e0rfn9MmzaN999/n/Ly8nrpt956K7fffjsrV67kk08+4dprrwXg+OOPZ8mSJWzcuJF+/fqxaNEiAH755RdOOOEEHn300UAP/ZJLLuGhhx5ixIgRrFu3jr///e9cddVVgTo2b97MN998w4oVK3jkkUfweFoO0rZz505cLhfDhmmLAN1yyy2sXLmSDRs24HQ6+eqrrwJ5q6urycnJ4aWXXmLq1Kkt1jl//nxSU1NZu3YtGzZs4PTTT8flcjFlyhRmzZrF+vXr8Xq9vPzyy4FyEhISyM7O5qabblLmo1bQYo9fSukTQvQWQliklE0Hjj8KsMdppp7wSDMGk6Ci2Bm4Zk1Pb7Op58DHH1M5fz6Jd9xB+IgRAJTsdlBe5GTEqU27RdZUV7Nqzqf0GzmGlLQBbapfcWgcrGceSqKiorjqqqt44YUXCAurWyRnwYIF5ObmBs4rKiqoqqpi/Pjx/Pzzz/Tu3ZubbrqJmTNnsnv3bmJjY4mIaLzOwuLFi/nkE215jZNOOomSkhIqKrS3m8mTJ2O1WrFarSQlJbFv3z569OjRqIxZs2bx888/s3nzZl588UVsNq3D9MMPP/DPf/6T6upqSktLGTx4MGeffTYAl156KQATJkygoqKCAwcONFvn0KFDmT59OnfffTdnnXUW48ePZ+3atfTt25cBA7T/wtVXX82MGTO47bbbADj//PMBGDVqFJ9+2mFrRB2xtMbGnw8sEUI8IIS4o3YLtWAdiT02Dr/Ph8tRSWScLTCJC7QBXnfeVi0y4GFQk5fHvif+TsS444i/9s+B9G3Z+xEC+mUmNnnfmvlzcDmqQjZLV/H75bbbbuONN97A4XAE0vx+P8uWLSMnJ4ecnBx2796N3W5nwoQJLFq0iEWLFjFx4kQSExOZPXs248ePP+R6rVZr4NhoNOL1epkxYwaZmZlkZmayZ88eQLPxr1u3jqVLl3LPPfdQWFiIy+Xi5ptvZvbs2axfv57rrruu3izohmNXtedN1TlgwACys7MZOnQo999/P48++mirZa8tQ9EyrVH824Cv9LyRQdtRQ+2CLFV6eObKkvo9fn91NZ49ew+5XL/Tye477sBgt5P6j3/Ui2WzLXs/qQNiCIu0NLqvptrB6q8+o9+oLFL6px/GEymOZOLi4rj44ot54403Ammnnnoq//nPfwLnOTk5APTs2ZPi4mLy8vLo168fJ5xwAk8//TQTJkwAIDIyksrKysB948eP5/33Nevtjz/+SEJCAlFRUc3KMm3atEBjk5paPzTH6NGjufLKK3n++ecDSj4hIYGqqipmz55dL++sWbMA7Y0jOjqa6OjoZuvcs2cP4eHhXHHFFdx1111kZ2eTkZHB9u3b2aqbXf/73/9y4oknNluGomUOOrgrpXwEQAgRLqWsPlj+I5HaRdc1l85IinbU/VECnj1b87D06H5I5e578ilq8rbS8/XXMSXW9exL9zgoK6xm6MTGr9EA2fO+xOWoUrb9Lsz06dN58cW6cFgvvPAC06ZNY9iwYXi9XiZMmMArr7wCwNixY/H5NBfk8ePHc++993LCCScA8Ic//IGnnnqKzMxM7r33Xh5++GGmTp3KsGHDCA8P55133mmTnHfffTcjR47kvvvu47rrrmPIkCGkpKQwpkGYE5vNxogRI/B4PLz55pstlrl+/XruuusuDAYDZrOZl19+GZvNxltvvcVFF10UGNy98cYb2yR7l0bqiwc0twHHAbnATv18OPDSwe5rz23UqFEylJQX7ZdPXzxZrl0wT66a95t88YaFssbpkVJK6a2okLkZA2XRzJmHVua8eTI3Y6Dc9/TTja6t+CpfvnjjQll1wNXomstRJf9zzcXys38+engPozgscnNzO1uEo5YTTzxRrly5srPFOOpp6jcMrJJN6NTWmHqeA04DSvSGYi0wIRSNUGcREaOFXa4q1Uw9UOfLb4yMxNSt2yEN8LoLCtj7wIPYhg8jsQmXuG3ZRXTrH01EtLXRtey5X1LjcCjbvkKhCBmtmsAlpWw4Q6OlRTyPOIwmE+HRMdokLt2Xv55nT1oaNa2M2SM9Hnbrsy+7P/MMokGY1AP7qinZXUX/kY0XLHE5qlj99ef0H32smimrOGr48ccfGT16dGeLoQiiNYp/lxBiHCCFEGYhxJ3AphDL1eHYY+O1SVy1vvzBnj3p6bi3bUP6Dt7eFb3wAq616+j22KNYmnCF27ZmPwD9RzT25sme+yU11Y4OicCpUCi6Lq1R/DcC04DuwG60ZRhvDqVQnYE9Lo6q0lJsdjMmqzEQnhk0xS/dbtz6rMnmqFq8hJLXXifmkkuIOv30JvNsyy4iuW8U9lhbvXSXo4rsuV+QNuY4kvocNcFPFQrF75DWKP4MKeXlUspkKWWSlPIKjsLInbVhG4QQRMXbGvX4oeWYPd6iIvbcfTfW9HSS772nyTzlRU6KdlY2aeZZ/fUX1FQ7OO7CS9v4JAqFQtEyrVH8/2ll2hGNPTae6opyfF4vUfG2+r78/fuBEM0qfun3s+fuu/E7HHT/9zMYbLYm8+Wv0aJ8NjTzuKq03n561jjV21coFCGnpcXWjwPGAYkNZupGAW1eb08IYQRWAbullGe1tby2Yo+LAylxHCgjMiGM3XkHAnH4DWFhmHv2bHaAt+T1N3As/YWUxx4NvB00xbY1+0nqHRnwHKpl9dzPcTurVW9foVB0CC31+C2AHa1xCJ6xWwG0fXVruJXf0SBx7exdzbPHhsflo8ZRN/W7uZg91WvWUPT880SdeQYxFzb/sVSWutj3WwX9GvT2nVWVZM/9kvSx40js3br1MhUKhaIttLQQy09Sm7V7rL7/l5TyESnlv6WUbYpaJoToAUwGXm9LOe1JRGzQEox6j7z+AG8a7h078LvrYtX5ysvZM/1OzN26kfLIIy2u0lVn5qlv38/+ura3rzx5FLBv3z4uu+wy+vXrx6hRozjuuOP47LPPDru8houdjBs37rDKOXDgAC+99FKz141GI5mZmQwePJjhw4fzzDPP4D/M+FbtwcHkdTqdnHjiiYEZz+3F/PnzycjIIC0tjaeeeuqQ802dOpWkpCSGDBlSL7/b7WbChAntFoeoNTb+VCFELrAZQAgxXAjR/CfaOp4D/go0+8sQQlwvhFglhFh1KIsIHy6RQWEbmnPpxOvF/dt2QA+1/MCDePbvp/u/n8EY2XL4om3Z+4nvYScmOTyQVlNdTfa8Lxlw7Akk9urTvg+kOOKQUnLuuecyYcIE8vPzWb16NR9++CEFBQWN8h2uUl26dOlh3XcwRRoWFkZOTg4bN27ku+++Y968eTzyyCOHVVd7cDB533zzTc4//3yMxjZbrQP4fD6mTZvGvHnzyM3N5YMPPqgXUbU1+aZMmcL8+fMb3WOxWJg0aVIg5lFb6fCZu0KIs4D9UsrVLeWTUs6UUo6WUo5OTGw6gmV7EhYZhcFopKqslMimJnE18Ow5MGsWld9+S9LttxOmxyNvDseBGvbml5M2sv5z7N68EbfTyfBTzmzPR1G0lXn3wFuT23eb17SnVzDff/89FoulXgya3r1785e//IXt27eTkZHBVVddxZAhQ9i1axfnnnsuo0aNYvDgwcycOTNwzxNPPMGAAQM44YQT2LJlS7067HZ74Pi9994jKyuLzMxMbrjhBnw+H9u3b+eYY47huuuuY/DgwZx66qk4nU7uuecetm3bRmZmJnfddVeLz5GUlMTMmTN58cUXa8O+NFmXw+Fg8uTJDB8+nCFDhgSU2rvvvsuwYcMYPnw4V155ZUjkff/99znnnHMAyMvLo0+fPoEAcB6Ph8zMzENeWWzFihWkpaXRr18/LBYLf/rTn/jiiy8OKd+ECROIi4trsvxzzz03EGCvrbRqBS4p5a4GZoy2vB8dD/xRCHEmYAOihBDv6W6inYYwGIiIicNRVoo1zIQ13BRYkAXA2qcPmEzU5OXh2rKFfX9/kojx44m7ZspBy87PKQJJIzfOgk0bMBhNdBugVtdSwMaNGxk5cmSz1/Py8njnnXc49thjAa3XGhcXh9PpZMyYMVxwwQVs376dDz/8kJycHLxeLyNHjmTUqFGNytq0aROzZs1iyZIlmM1mbr75Zt5//30mTJhAXl4eH3zwAa+99hoXX3wxn3zyCU899RQbNmwIRAU9GP369cPn87F//35KS0ubrCsiIoLU1FS+/vprAMrLy9m4cSOPP/44S5cuJSEhgdLS0naX1+12k5+fT58+fQBIT0/n+uuv55tvviEtLY0XX3yRP/7xj/Ts2TNwz/jx4+tFOa3l6aef5uSTTwZg9+7d9e7p0aMHy5cvb3RPa/M1ZMiQIaxcufKg+VpDaxR/vZm7tHFQVkp5L3AvgBBiInBnZyv9Wuxx2kpcQKPwzMJiwdKnN851a6lcsABDdBSpTz1ZL9Ryc2zL3k9stwhiU+ovjFGwaQMpaQMwWxrH7FF0Imc0b5vtSKZNm8bixYuxWCx8/PHH9O7dO6D0QYvYWWv/37VrF3l5eSxbtozzzjuP8HDNpPjHP/6xybIXLlzI6tWrA1E0nU4nSUlJTJgwgb59+5KZmQloC5ts3749EO3zcGiurssuu6zRgivvvvsuF110EQkJCQCB3m97yltcXExMTEy9tCFDhrBgwQJKS0t54403Gini2pXNOhOj0YjFYqGyspLIg5iWD0ZrFP+NwPPUzdz9Fm0m71GHPTae0j2aPTUq3kbJHke969b0dCrnzQch6PXmG5ji4w9aZnWFmz15Bxh1Zp966R6Xi335Wxl99vntJr/iyGbw4MGB1bEAZsyYQXFxcSDOTfCKWj/++CMLFizgl19+ITw8nIkTJ9Zb+ORgSCm5+uqrefLJJ+ulb9++vdHiKE6ns+HtByU/Px+j0UhSUlKzdQFkZ2czd+5c7r//fiZNmkRsbGzI5Q0LC2v0WQ0YMIAZM2bw8MMPc+eddzZavaw1Pf7u3bvXMw8VFBTQvXvjUO6tzdcUNTU1gRXP2sJBu6tSyuKGM3ellCVtrlkr+8ffgw9/LbWzdwEiE8KoLHEh/TJwvdbOH3/99UQcd1yryszPKUJKSGtg5tm7dQt+n48eAztvmT/F74uTTjoJl8tVby3Z6uqml8AoLy8nNjaW8PBwNm/ezLJlywDNRvz555/jdDqprKxkzpw5Td4/adIkZs+ezf79Wuyo0tJSduzY0axsDRd0aYnaheFvueUW9m+J2AAAHQNJREFUhBDN1tXUgisnnXQSH3/8MSUlJYG87S1vbGwsPp+vnvLv378/2dnZrFixot46xLUsWrQosCBN8Pb/7d17eFT1ve/x9zc3AhkCuRiRhJugkAACAtIePVrsBkHdVXQXa9VtD56j3cdatd3PqdbT2u7dx+Ktais9rUd3vbbWpzcvVesNusXWoxJBAglyEWQCSkiAkJiQZOZ7/lhrhkkyk0xgMmsm8309zzwkM2vWfIHkN2t+67e+n9CgDzB//ny2bt3KRx99REdHB08//XTUT1zxbtdTY2MjpaWl5OZGz+geiH6P+EVkEnADMDFye1Xtv9I04ysq5khrK51H2iksySfQFeSz5g4KRjtHFKMvvhjJzqFkxX+Le5873t/H6BNHUDy29zSPSBZjpw657hfmGIkIf/rTn7j55pu56667OOGEEygoKODOO+/ste2SJUv4xS9+QWVlJVOnTg1PAZ1++ulcdtllzJo1i7Kysl6BKCFVVVX86Ec/YvHixQSDQXJzc1m1ahVjxoyJun1JSQlnnnkmM2bMYOnSpdx9993dHm9ra2P27Nl0dnaSk5PDVVddxbe+9a0+X+vQoUO9AlemT5/ObbfdxjnnnEN2djZz5szh0UcfTXi9ixcvZu3ateGBOzc3l8LCQlauXElWHNO30eTk5PDggw9y3nnnEQgEWLFiBdOnHz2wO//883n44YcZO3ZszO0uv/xy1qxZw/79+6moqOCHP/wh11zjRLauXr2aCy644Jhq6yVak/7IG7AB+CawEDgndOvveYm8DXYQS0jNmtf0nuUXaNPeet25cb8+eN3rumfrgWPeX9vhDl31L2/o3/+4rddjv/3hrfr4d755POWaBLIglsyybt06vfLKK7vdN27cOA0Ggx5V1L9ly5bpli1bYj4+kCCWeOb421X1p4l5m0lt4at3m5oYWeJcRdvc2M5JU45tfzs2NKBB7bWaJ9DVyd4P6zht0dLjqtcYc2xOP/10Fi5cSCAQIDs7m507dzJhwoQ+L8L0UkdHBxdffDGnnnpqQvYXz8D/gIjcjnNS90joTlWtTkgFKcRXfPTq3TFTnCmYyLX8A7W9uoHC0nxKx/m63f/J9m10dXZQUWnz+8Z4ZcWKFeGvJ06cmBIrd2LJy8uLeu7hWMUz8M8ErgLO5eiVtup+P6QcbdvQRE5eNiMK88IRjAN15LNO/HVNzDp3XK+jCH9tDQDldmLXGOOBeAb+LwMnq2pHv1umuWEjCsgZNoyWJmdFQWFpfrd+PQOx84P9BAO9p3nAuWK3uHwcIwpHHVe9xhhzLOI5fV0DjO53qyFARPAVHb2Ia2TJ8G79egZiW3UDvqJhlE3sfqFFMBigvm6zTfMYYzwTzxH/aKBORN6l+xz/kFvOCc4J3tbw1bv5bFu3j2AgSFZ2/Eu8Otq62L25iRlnl/ea5mnYtZOOts+oqJwR49nGGDO44hn4bx/0KlKIr7iET7Z9CDhtGzSotBw40is8pS87a/YT6Aoy+fTezeXqbX7fGOOxfgd+Vf1rMgpJFQXuVI+qMrIk1J65bUAD//bqBkaMymPMyb3n8P21mxhVdiKFpYPfcdQYY6Lpd/5CRD4nIu+KSIuIdIhIQESak1GcF3xFxXR1HOHIZ60UhtozD2BlT+eRAB/XNDJ5ThmS1X2aR1Xx122yaR5jjKfimbh+ELgc2AoMB/47sGowi/KSz13S2XqgCV/xMEQY0JLOXTWNdHVGn+Zp2uOnrfmQTfMYYzwV1xlLVd0GZKtqQFV/BSwZ3LK843OTuA43NZKdnYWvKH9AF3Ftf38fw0fmctKU3guh6ms3AdiKHhOTRS8mRjpGL+7evZuFCxdSVVXF9OnTeeCBB8LbexG9+JmI5AHrReQuEbk5zuelpcjQdXDX8se5pLOrI8DOjY2cPPsEsrJ6X/rtr62hYHQRo8eMTVzBZshQi15MmHSMXszJyeHee+9l8+bNvP3226xatSr8/ERHL8azqucqnIH+G8DNwDjg0oS8egoqCLVtaDranvnjTfF1of54cxNdRwJRL9pSVXbX1lBeOSNl+4EYx53v3EldU11C9zmteBrfOeM7fW7TX/Tieeedx4IFC1i3bh0vvvgiN954I7t376a9vZ0bb7yRa6+9FnCiFx977DHKysoYN25ctwQun89HS0sL4EQZ/vSnP6Wjo4MFCxbw85//nN27d7N06VLOOuss/va3v1FeXs6zzz7bLcpw0aJFvbpdRgpFL86fP58f/OAHiEjU12pvb2f58uX4/X4CgQDf+973uOyyy3j88ce55557EBFOO+00nnjiiYTX+9RTT/HrX/8acJLNFi1axGuvvcaUKVPo7Oxk/vz5PP/8892SsvoTGakIhCMVq6qq4tru1ltv5aSTTgKcttKVlZXU19eHn3/xxRdz6623csUVV8RdUyx9Dvwikg3coapXAO2Ad2/hSZKbN4z8Al+4L39hST6fHeqgqyNATl7fRwfb399HfkEuY0/tPc3T3LCPlsb9Ns1jYrLoRYteDNm5cyfvv/8+CxYsCN+XtOhFVQ2IyAQRycuElg0hBUXFtLg/cIXuks7DTe29ohMjBTqD7Nywn8lzy8iOcrFXqD+PrehJff0dmSeLRS9mZvRiS0sLl156Kffffz+FhYXh+5MdvbgDeEtEngPCWYSq+pPjeuUU5is+evXuyNKjSzr7Gvh31zXR0R5g8pze0zzg9OfJL/BRWjE+8QWbIcGiFy16sbOzk0svvZQrrriCSy7pHcuatOhFYDvwgrvtyIjbkBXZrye0lv9wPyt7tlfvI294DhXTov/g+ms3MXZaVVzh7CYzWfRiZkcvqirXXHMNlZWV4fSySEmNXlTVIT+v35OvuITWg01oMEjBqDyyc7L6XNkTCAT5aMN+Js0qJTun98DeevAAB/bWM/OL5w1m2SbNWfRiZkcvrl27lieeeIKZM2eGp67uuOMOzj//fCD50YsnAHcDLwJvhG79PS+Rt2RFL4ZUv/y83rP8Am050KSqqk9+/+/60i8/iLn9rk1OTOOODQ1RH6/725t6z/ILdM/WukGp1xw/i17MLJkevRjPW9tTQB0wCWdVz04gMaeWU5QvIpAFnBO8fR3xb69uIDc/m3GVsaZ5asgdlk/ZxMmJL9YYM2CR0YtAxkUvxjPwl6jqI0Cnqv5VVVcwBNO3IoWu3g2f4C2JHcgSDATZsb6BiTNLycmNvtyzvraGk06dRnZOPOfSjTHJsGLFivAFXJkWvRjPwN/p/rlXRC4QkTlAccIqSEGhq3ePJnEN50hrFx1tvS+X3rP1IO0tnVF78wC0t7TQsHuXrd83xqSMeA5BfyQio4BvAz8DCnGu4B2yRowaDSLhi7jC7Zkb2ymt6B6cvr26gZy8LMZPL4m6r/otm0HV1u8bY1JGzIFfRPKBrwNTgHLgEVVdmKzCvJSdk8OIwlFH5/hDa/n3t3Ub+INBZcf6BibMKCU3xlW9/toasnNyGDMlMXNzxhhzvPqa6nkMmAdsBJYC9yalohThKy7pFroOvdszf7L9EJ81d8Sc5gGnI+eYKaeSmzcs5jbGGJNMfU31VKnqTAAReQR4JzklpQZfUTGH3YE/vyCX3GHZvdozb6/eR3ZuFhNmRJ/m6Wxv59OPtjH/S0O2p50xJg31dcQfOqmLqiamCXQaiQxdFxGnPXPEEb8Gle3vNzC+qpi8/Ojvn3u21hEMBCx4xRiTUvoa+GeJSLN7OwycFvr6eKIXRWSciKwWkc0isklEbjzWfQ0mX3EJnx06SMANPhhZMrzbEf+nO5tpPXgkagvmEH/tJkSyGHtq5aDXa4wx8Yo51aOqiUso6K4L+LaqVovISGCdiLyqqr0TCzxUEIpgPNhEYWkZhSX5+LccQFUREbZX7yMrR5h4WmnMfdTX1lA26WSGuZ0SjTEmFSS9Y5iq7lXVavfrw0AtzqqhlOILB7IcXdnTdSRAe2snqsr26gbGVxYzbHj0986uzk72bt1i6/fNgFj0YmKkY/Rie3s7Z5xxBrNmzWL69Oncfvvt4e29iF4cNCIyEZgD9EorEJFrReQ9EXmvoaEh2aX1imAMr+Xf307Dx4c53NTOyTFaMAN8un0rXZ0dlNv6fRMntejFhEnH6MVhw4bxxhtvsGHDBtavX8/LL78c7rrqRfTioBARH/B74CZV7XXOQFUfAh4CmDdvnia5vG6h69B9Lf/+3YfJyhImzYo9zRMKXimfWhVzG5OaPrnjDo7UJjZ6cVjlNMZ897t9bmPRi5kdvVhVVYXP51wn1NnZSWdnZ7feQUmLXhwsIpKLM+g/pap/8KKG/gz3jSQrO4fWiAhGcNbyb6tuoGJaEfkFsfti19dtoqRiPCMKRyWlXpP+LHrRohcDgQBz585l27ZtXH/99d5ELw4Gcd7CHgFqNYVTvCQri4KiovDVu3nDcxhWkMOO9Q00N7Qx97wJMZ8bDAao37KZyrO+kKRqTSL1d2SeLBa9mHnRi9nZ2axfv56DBw+ybNkyampqmDFjRvixZEYvJtqZwFXARhEJvR1/V1Vf9KCWPvmKS8IDPzhpXJ9+1IxkCZNmx57madj5ER1tbTa/bwbEohctejFk9OjRLFy4kJdffjk88ENyoxcTSlXXqqqo6mmqOtu9pdygD24EozvHD0dbN4w9ZTTDfXkxn+ev3QRAhV24ZQbAohczO3qxoaGBgwcPAs4nmldffZVp06aFn5fU6MVM5isq4eONG8Lfh/J3p/TRmwecE7ujThzDyJLYnwqM6cmiFzM7evGDDz7g6quvJhAIEAwGWb58ORdeeGH4+UmNXkyFW7KjF0Pe/uMzes/yC7SjrU1VVT989xP95Y1rtOVge8znBINBXXXN5frSqvuSVaZJAItezCyZHr1oR/x9GOku6Ww50EjRSeVMmVvGxNNit2AGaKr303a4mfJKW8ZpTKqKjF7Mzs626EVzVEGP7F0R6XPQh6Pr9y14xZjUZtGLJqpwBGPEyp7++GtrKCgqZvSJJw1WWcYYc1xs4O9DOHQ9YmVPX1QVf20NFdOmp+xHRmOMsYG/D3nDh5M7LD+cvduf5oZPaWlqtGkeY0xKs4G/DyKCr7g43KGzP6H1++XWkdMYk8Js4O9HQVFx3HP8/tpN5Bf4KK0YP8hVGWPMsbOBvx+REYz9qa+robxyOnKMF4AYY0wy2AjVD19xCS1NjTjXQsTWcqCJA3v3WJsGY0zKs4G/H76iYro6OzjS2trndvV1bn8eO7FrjElxNvD3wxdx9W5f/LU15A7Lp2zS5GSUZYYoi15MjHSMXgwJBALMmTOnW5+eIRW9mA7CV+/2s5bfX7uJsVMryUpglJvJLGrRiwmTjtGLIQ888ACVlZXdnjNkohfTRTxX77a1HGb/7l1M/dyxh1WY1PHmMx+yf3dLQvdZOs7Hf13ed58Vi1606EW/38+f//xnbrvtNn7yk+45VWkfvZhOCtxgiL5W9uzZshlUbX7fHBeLXrToxZtuuom77ror6muldfRiusnNG0a+b2Sfc/z+2k1k5+QwZkpiOucZb/V3ZJ4sFr2YWdGLL7zwAmVlZcydO5c1a9b0ejzdoxfTTs8krp78tTWMmTKVnLzYqVzG9MeiFzM7evGtt97iueee48UXX6S9vZ3m5mauvPJKnnzyyfC2aRu9mI76unq3o72NT3dss2kec9wsejGzoxd//OMf4/f7w9N15557brdB36IXk8xXVEJj/e6oj+39cAsaDFIxzYJXzPGx6MXMjl7sj0UvJtmbv3lc7/3KP2og0NXrsbW/fULvvewf9chnrR5UZhLFohczS6ZHL9pUTxx8RcVoMEhbc3Ovx/y1NZRNmkze8BEeVGaMORaR0YuARS+a3gqKo1/E1dXZyd6tW6iwNszGpB2LXjR9GhnjIq5Ptn9IoLPTTuwaY9KKDfxxiHXEXx8KXrETu8aYNGIDfxwKRhWBSK8jfn/dJkoqxjN8ZKFHlRljzMDZwB+HrOxsCkaNpjXi6t1gIMCeLZttmscYk3Zs4I9TKJAlpGHXR3S0tdmJXWNM2vFk4BeRJSKyRUS2icgtXtQwUD2v3vXX1gAWrG6MST9JH/hFJBtYBSwFqoDLRSTlz476ogz8o088iZHFpR5WZYwxA+fFEf8ZwDZV3aGqHcDTwEUe1DEgvuIS2poPEejqRINB/HWbKbd8XWNMGvJi4C8HIhvf+N37UlookKX1wAGa9vhpP9xs8/sm4Sx6MTHSNXpx4sSJzJw5k9mzZ4e7skIGRS+KyLUi8p6IvNfQ0OB1OfhCEYwHGsPz+7aixySSWvRiwqRz9OLq1atZv3497733Xvi+oRC9WA9E5plVuPd1o6oPAQ8BzJs3T5NTWmzh7N0DTfhrN+ErKmbUidE7A5r0tvrRh9i3a0dC91k24WQWfu3aPrex6EWLXuxLukcvvgucIiKTcAb8rwBf9aCOAfEVu20bmpwj/vLKGSnb0MmkJ4tetOhFEWHx4sWICNddd134zRzSPHpRVbtE5BvAX4Bs4D9UdVOy6xio4SMLycrOob5uMy1NjTbNM4T1d2SeLBa9mFnRiwBr166lvLycffv2sWjRIqZNm8bZZ58NDIHoRVV9EXjRi9c+ViKCr7iYHdXOO64Fr5hEs+jFzI5eDD0GznTZsmXLeOedd8IDP1j0oid8RSV0dRwh3zeSkorxXpdjhhiLXszs6MXW1tZwza2trbzyyivMmHF0ZsGiFz0SWtlTPm06cozxbMbEYtGLmR29uGPHDpYtWwZAV1cXX/3qV1myZEn4+YmMXhQnnSu1zZs3TyOXNnnljUd/yfsvPc85V13DvAuXeV2OSaDa2loqKyu9LsMkSXV1Nffdd194xRDA+PHj2bVrV8ou2rjkkktYuXJlzBSuaD/DIrJOVef13NYOWwcgdBGXndg1Jr1levSiTfUMwNTPn0Vnextlk072uhRjzHFasWJF+OtMi160gX8ARpWN4czLrvK6DGOMOS421WOMMRnGBn5jXOmw0MGYaAb6s2sDvzFAfn4+jY2NNvibtKOqNDY2DujCLpvjNwanX4rf7ycVOsEaM1D5+flUVFTEvb0N/MbgXMAzadIkr8swJilsqscYYzKMDfzGGJNhbOA3xpgMkxa9ekSkAYjdiq9vpcD+BJYzGFK9xlSvD1K/xlSvD6zGREi1+iao6gk970yLgf94iMh70ZoUpZJUrzHV64PUrzHV6wOrMRFSvb4Qm+oxxpgMYwO/McZkmEwY+B/yuoA4pHqNqV4fpH6NqV4fWI2JkOr1ARkwx2+MMaa7TDjiN8YYE8EGfmOMyTBDeuAXkSUiskVEtonILV7XE0lExonIahHZLCKbRORGr2uKRUSyReR9EXnB61p6EpHRIvI7EakTkVoR+bzXNfUkIje7/8c1IvIbEYm/jeLg1fQfIrJPRGoi7isWkVdFZKv7Z1GK1Xe3+//8gYj8UURGe1VfrBojHvu2iKiIlHpRW3+G7MAvItnAKmApUAVcLiJV3lbVTRfwbVWtAj4HXJ9i9UW6Eaj1uogYHgBeVtVpwCxSrE4RKQe+CcxT1RlANvAVb6sC4FFgSY/7bgFeV9VTgNfd773yKL3rexWYoaqnAR8Ctya7qB4epXeNiMg4YDHwcbILiteQHfiBM4BtqrpDVTuAp4GLPK4pTFX3qmq1+/VhnAGr3NuqehORCuAC4GGva+lJREYBZwOPAKhqh6oe9LaqqHKA4SKSA4wA9nhcD6r6n0BTj7svAh5zv34MuDipRUWIVp+qvqKqXe63bwPx9yEeBDH+DQHuA/4XkLIrZ4bywF8O7I743k8KDqwAIjIRmAP8P28riep+nB/ioNeFRDEJaAB+5U5FPSwiBV4XFUlV64F7cI7+9gKHVPUVb6uK6URV3et+/QlwopfF9GMF8JLXRfQkIhcB9aq6weta+jKUB/60ICI+4PfATara7HU9kUTkQmCfqq7zupYYcoDTgf+jqnOAVrydnujFnSe/COdNaixQICJXeltV/9RZ552SR6wichvOVOlTXtcSSURGAN8Fvu91Lf0ZygN/PTAu4vsK976UISK5OIP+U6r6B6/rieJM4EsishNnquxcEXnS25K68QN+VQ19UvodzhtBKvkH4CNVbVDVTuAPwH/xuKZYPhWRkwDcP/d5XE8vIvI14ELgCk29i5Am47zBb3B/ZyqAahEZ42lVUQzlgf9d4BQRmSQieTgn1J7zuKYwERGcuelaVf2J1/VEo6q3qmqFqk7E+fd7Q1VT5mhVVT8BdovIVPeuLwKbPSwpmo+Bz4nICPf//Iuk2AnoCM8BV7tfXw0862EtvYjIEpxpxy+p6mde19OTqm5U1TJVnej+zviB092f05QyZAd+9yTQN4C/4PyiPaOqm7ytqpszgatwjqLXu7fzvS4qDd0APCUiHwCzgTs8rqcb99PI74BqYCPO75znl/WLyG+AvwNTRcQvItcAK4FFIrIV55PKyhSr70FgJPCq+/vyC6/q66PGtGAtG4wxJsMM2SN+Y4wx0dnAb4wxGcYGfmOMyTA28BtjTIaxgd8YYzKMDfzGE27nwnsjvv9XEflBgvb9qIj8UyL21c/rfNntCLq6x/0TQx0bRWR2Ipfput1I/2fE92NF5HeJ2r/JDDbwG68cAS5Jtba1biO1eF0D/A9VXdjHNrOBAQ38/dQwGggP/Kq6R1UH/U3ODC028BuvdOFcyHRzzwd6HrGLSIv75xdE5K8i8qyI7BCRlSJyhYi8IyIbRWRyxG7+QUTeE5EP3Z5DoVyBu0XkXben+3UR+31TRJ4jypW/InK5u/8aEbnTve/7wFnAIyJyd7S/oHvF+L8Bl7kXHF0mIgVuH/d33MZyF7nbfk1EnhORN4DXRcQnIq+LSLX72qHOsiuBye7+7u7x6SJfRH7lbv++iCyM2PcfRORlcXrt3xXx7/Go+/faKCK9/i/M0DSQoxtjEm0V8EFoIIrTLKASpx3uDuBhVT1DnCCbG4Cb3O0m4rTmngysFpEpwD/jdMecLyLDgLdEJNQp83ScXu8fRb6YiIwF7gTmAgeAV0TkYlX9NxE5F/hXVX0vWqGq2uG+QcxT1W+4+7sDp/XFCnGCRN4RkdciajhNVZvco/5lqtrsfip6231jusWtc7a7v4kRL3m987I6U0SmubWe6j42G6cD7BFgi4j8DCgDyt2cAMTjYBOTPHbEbzzjdiN9HCeoJF7vulkGR4DtQGjg3ogz2Ic8o6pBVd2K8wYxDScc459FZD1OC+wS4BR3+3d6Dvqu+cAat8laqCPk2QOot6fFwC1uDWuAfGC8+9irqhrq7y7AHW4ritdwWor31yb5LOBJAFWtA3YBoYH/dVU9pKrtOJ9qJuD8u5wsIj9z++CkVHdYM3jsiN947X6cPja/irivC/egRESygLyIx45EfB2M+D5I95/nnr1IFGcwvUFV/xL5gIh8AaelczIIcKmqbulRw4IeNVwBnADMVdVOcbo9Hk9kY+S/WwDIUdUDIjILOA/4OrAcp8+9GeLsiN94yj3CfQbnRGnITpypFYAvAbnHsOsvi0iWO+9/MrAFp2Hfv4jTDhsROVX6D255BzhHRErFifO8HPjrAOo4jNNYLOQvwA0iIm4Nc2I8bxROFkKnO1c/Icb+Ir2J84aBO8UzHufvHZU7hZSlqr8H/jep19LaDBIb+E0quBeIXN3zf3EG2w3A5zm2o/GPcQbtl4Cvu1McD+NMc1S7J0R/ST+fet1EqluA1cAGYJ2qDqRd8WqgKnRyF/h3nDeyD0Rkk/t9NE8B80RkI865iTq3nkaccxM1UU4q/xzIcp/zW+Br7pRYLOXAGnfa6Um8z7A1SWLdOY0xJsPYEb8xxmQYG/iNMSbD2MBvjDEZxgZ+Y4zJMDbwG2NMhrGB3xhjMowN/MYYk2H+P1sJW3nXJRPRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Convergence of Algorithms\")\n",
    "plt.xlabel(\"Number of Iterations\")\n",
    "plt.ylabel(\"Parameter Estimate\")\n",
    "for e in data:\n",
    "    plt.plot(x, e[\"path\"], label=e[\"label\"])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum profit is 328.63353450309955 where the quantity is 5.477225490171372.\n",
      "Scipy's mimimizer with its implementation of the BFGS solver took 3 iterations.\n"
     ]
    }
   ],
   "source": [
    "results = minimize(fun=lambda x: -pi(x), x0=5)\n",
    "print(\n",
    "    f\"The maximum profit is {-results.fun} where the quantity is {results.x[0]}.\"\n",
    ")\n",
    "print(\n",
    "    f\"Scipy's mimimizer with its implementation of the BFGS solver took {results.nit} iterations.\"\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gradient-descent-newton-raphson",
   "language": "python",
   "name": "gradient-descent-newton-raphson"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}