PS1/problem_set_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem set 1: Solving the consumer problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[<img src=\"https://mybinder.org/badge_logo.svg\">](https://mybinder.org/v2/gh/NumEconCopenhagen/exercises-2020/master?urlpath=lab/tree/PS1/problem_set_1.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this first problem set, we will take a look at solving the canonical utility maximization problem for the consumer.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem set structure:** Each problem set consists of tasks and problems. _Tasks_ train you in using specific techniques, while _problems_ train you in solving actual economic problems. Each problem set also contains solutions in hidden cells. *You should really try to solve the tasks and problems on your own before looking at the answers!* You goal should, however, not be to write everything from scratch. Finding similar code from the lectures and adjusting it is completely ok. I rarely begin  completely from scratch, I figure out when I last did something similar and copy in the code to begin with. A quick peak at the solution, and then trying to write the solution yourself is also a very beneficial approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Multiple solutions:** Within the field of numerical analysis there is often many more than one way of solving a specific problem. So the solution provided is just one example. If you get the same result, but use another approach, that might be just as good (or even better)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Extra problems:** Solutions to the extra problems are not provided, but we encourage you to take a look at them if you have the time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Download guide:**\n",
    "\n",
    "1. Follow the [installation guide](https://numeconcopenhagen.netlify.com/guides/python-setup/) in detail\n",
    "2. Open VScode\n",
    "3. Pres <kbd>Ctrl</kbd>+<kbd>Shift</kbd>+<kbd>P</kbd> \n",
    "4. Write `git: clone` + <kbd>Enter</kbd>\n",
    "5. Write `https://github.com/NumEconCopenhagen/exercises-2020` + <kbd>Enter</kbd>\n",
    "6. You can always update to the newest version of the code with `git: sync` + <kbd>Enter</kbd>\n",
    "7. Create a copy of the cloned folder, where you work with the code (otherwise you can not sync with updates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement a Python version of this function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \n",
    "u(x_1,x_2) = (\\alpha x_1^{-\\beta} + (1-\\alpha) x_2^{-\\beta})^{-1/\\beta} \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def u(x1,x2,alpha=0.5,beta=1):\n",
    "    return (alpha*x1**(-beta) + (1-alpha)*x2**(-beta))**(-1/beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## print"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1_vec = [1.05,1.3,2.3,2.5,3.1]\n",
    "x2_vec = [1.05,1.3,2.3,2.5,3.1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct a Python function `print_table(x1_vec,x2_vec)` to print values of `u(x1,x2)` in the table form shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update this code\n",
    "\n",
    "def print_table(x1_vec,x2_vec):\n",
    "    \n",
    "    # a. empty text\n",
    "    text = ''\n",
    "    \n",
    "    # b. top header\n",
    "    text += f'{\"\":3s}'\n",
    "    for j, x2 in enumerate(x2_vec):\n",
    "       text += f'{j:6d}' \n",
    "    text += '\\n' # line shift\n",
    "    \n",
    "    # c. body\n",
    "    # missing lines\n",
    "    \n",
    "    # d. print\n",
    "    print(text) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        0     1     2     3     4\n",
      "  0  1.050 1.162 1.442 1.479 1.569\n",
      "  1  1.162 1.300 1.661 1.711 1.832\n",
      "  2  1.442 1.661 2.300 2.396 2.641\n",
      "  3  1.479 1.711 2.396 2.500 2.768\n",
      "  4  1.569 1.832 2.641 2.768 3.100\n"
     ]
    }
   ],
   "source": [
    "def print_table(x1_vec,x2_vec):\n",
    "    \n",
    "    # a. empty text\n",
    "    text = ''\n",
    "    \n",
    "    # b. top header\n",
    "    text += f'{\"\":3s}'\n",
    "    for j, x2 in enumerate(x2_vec):\n",
    "       text += f'{j:6d}' \n",
    "    text += '\\n' # line shift\n",
    "    \n",
    "    # c. body\n",
    "    for i,x1 in enumerate(x1_vec):\n",
    "        if i > 0:\n",
    "            text += '\\n' # line shift\n",
    "        text += f'{i:3d} ' # left header\n",
    "        for j, x2 in enumerate(x2_vec):\n",
    "            text += f'{u(x1,x2):6.3f}'\n",
    "    \n",
    "    # d. print\n",
    "    print(text)\n",
    "\n",
    "print_table(x1_vec,x2_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reproduce the figure below of $u(x_1,x_2)$ using the `meshgrid` function from _numpy_ and the `plot_surface` function from _matplotlib_. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate utility function\n",
    "import numpy as np\n",
    "x1_grid,x2_grid = np.meshgrid(x1_vec,x2_vec,indexing='ij')\n",
    "u_grid = u(x1_grid,x2_grid)\n",
    "\n",
    "# import plot modules\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm # for colormaps\n",
    "\n",
    "# write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeXxTdbr/P2fJnrZpm6QrXSibgKxuoAKOK4qOioK7XjfGGXdFB2VcRkeGUa/3OooodwYd/XmdEbyIgMOoKC4giOLGTvc1zdKk2ZOz/P5gTkzS7Dltk3rer1df0OTke75Jk0++5/k+z+cheJ7nISEhISExJJDDPQEJCQmJnxOS6EpISEgMIZLoSkhISAwhkuhKSEhIDCGS6EpISEgMIZLoSkhISAwhkuhKSEhIDCGS6EpISEgMIZLoSkhISAwhkuhKSEhIDCGS6EpISEgMIZLoSkhISAwhkuhKSEhIDCGS6EpISEgMIZLoSkhISAwhkuhKSEhIDCGS6EpISEgMIZLoSkhISAwhkuhKSEhIDCH0cE9AIn/heR4cx8Hj8YDnedA0DZIkQVEUSJIESZIgCGK4pykhkVNIoiuRNjzPg2VZMAwT+j/HceA4DjzPRwitIMLCjyTGEj93JNGVSJlosSUIAiR5LEJFEAQoihpwPM/zCAaDCAQCkhhLSEASXYkU4HkeDMOAZdkBYpsIgiBiiijP8wAAhmEQDAYBABaLBXq9XhJjiRGPJLoScRHElmEYAEhZbJMhCGi4kPp8vtDYghiHhyqElbQQNxbEWRJjiXxDEl2JAXAcB5Zl0dvbC61WC4VCkVDcCIIIrV6zIZYYAz+tjKNDGzzPJ1wZS4IskYtIoisRguO4UBgBAPx+P9RqdUriJYboxiNVMY5+DEmSAzIqJDGWGG4k0f2ZE77ZxXEcgJ9isSRJJhXTYDAIi8UChmGgUqkgl8vRuG4dGK8X026+GSQ9eG+xZGLMcRz8fv+Ax8RaHUtiLDFUSKL7M0XIsWUYZoDYCiQKGwQCAVgsFni9XhQWFkKtVoPjODR++CG23XsvuGAQ37zyCk54+GHUnHEGFAoF5HI5aJoedHGLJ8bAT89bWM2HI3zJKJVKSYwlBg2CH8zrQomcI1psEwlUT08PNBoNCgoKQrf5/X6YzWYEAgHo9XoUFBTA7/eDZVm4Wlrw5plnwm+3R4xTc845mPHb30JRUQGGYUAQBORyeehHoVCgs7MTdXV1g/rcEyFsGnZ3d6OysjLi9YgVM5YyKiQyRRLdnwmxcmyB2GIrYDKZoFKpUFhYCK/XC7PZDJZlYTAYoNFoQo/1+XxwW614+5xzYG9sjDkWJZdj+m234ZQHHgCt0SAQCIR+/H4/PB5PaDUc/iOTyYZM3FiWRXd3N6qrq0O3CeEX4UfKNZbIFkl0RzixxDZVUejt7QUAeDweEAQBg8EAtVo94DiPy4V3LrkEHZ99lnRMtdGI0x59FJOuuSZiHq2trRg1ahSCwSD8fn9IkIPBIAiCgEwmixDlwRDjWKIbD0mMJTJFEt0RSqyChlQ/9DzPw+VyoaurCzRNo6qqCkqlMu7xm5YswY+vvZbW/MqmT8e8lStRdcopAI6Jbm1tbcxjOY4LVbUJgiwUVYSHKLIVY4Zh0NPTk5LoxkP4OAliHI4kxhKAJLojjlgFDemIrdPphMViCW16KRQKFBcXx33M7v/+b2x78MGM5zt+4ULMeeIJ2Fg2rugmmm94iCKWGAuCnIoYiyG6ieYq/BtLjMPT26TCj5GNJLojBGFHPlOxdTgcsFqtUKlU0Ov1kMvlsFqtIEkyruge3bIF6y+7DPy/sx8yhVarMfaGG3DWo49CplJlNRYQKcaCIKcixgzDwGQyoaqqKus5pDNX4d/wUIXb7YZCoQhlUkiFHyMHSXTzHCETwePxwGq1oqqqKuUPJMdxsNvtsNls0Gq10Ov1oMPyam02GwCgpKRkwGN7f/wRb5xxBgJOpzhPBEDBqFE4/fHHMeGyy0QbMxwhHzk8ZhwIBAAAMpkMNE3D6/WivLwccrl8WIVNqAZUKpUDVsZS4Ud+I4luHhKroCEYDMJkMqGmpibp4zmOg81mg91uR2FhIUpLSwc4hAFAX18fOI5DaWlpxO1ukwmvnX46+tvaxHlCUVTNmoV5K1eibNq0QRk/GuG19Hg86Ovrg1KpjBDj8JjxUImxyWRCYWEhVFEr/0RhCqnwIz+QiiPyiEQFDalUj7EsC6vViv7+fuh0OtTX18cUW4FYxRGMz4f1ixYNmuACQOfOnXhz3jxMuuYanPboo1AbDIN2LgChvGGCIODxeFBRUQHgJzEWQhQulyuuGMtkMlHMgASiMyLC5xr+b/Rj4hV+hAuxJMbDiyS6eUCsgoZY1WNcnNgqwzCwWCxwuVwoKSnB6NGjU7ZmjBbdLUuWoGvXruyeUArwHIcf//Y3HN6wAScvXYoZt90GSi4f3HPGWDkKq1utVhtxnCDGgUAALpcr5IomiHF4zDgTMY4nuolIZKUpiHG4WZDg3CYVfgwtkujmMPFybGN9IGKtdAVfBLfbjdLSUpSVlaX9YQof8/M//AH7//73zJ5MhgT6+/HZ736HH159FXOfegoN8+cP6fljES7G4USLscfjQSAQCLUyii7+SCTGmYhuovmmKsYCUnrb4CGJbg6SqENDPMJXpX6/HxaLBT6fD3q9HuXl5Rl9WMLHPLBuHT5/8sn0n4xI2Bsb8e7ixag980zM++MfUTp+/KCcJxtRSSTGDMOENvBSEWMxRTfRfOOJcfgXiNADT6vVSmIsAtJGWg6RTUEDABw+fBgqlQrBYBAGgwFarTarD4PT6YTb7QbX3o43zzkHjNeb8VhiQtI0ptx0E2Y/9BCUCXKI0yUQCMBqtYZiuoON8PcOzzMWxJhhGGi12pBzW7KV8WBi/7eXRlFRUWjeUuFH5kiimwNkU9AAHCvTNZvNcLvdqK2thUajEWVeLpcLPYcO4f2FC+Hu6RFlTDFRlpRg9sMPY8qNN4JMsCGYKn6/HzabbchENx48z6O9vR3FxcURK2RhZRxdhTfYYmyz2UDTNAoLC+POV/g3UeGHJMbHkER3GMm2oMHtdsNisYAkSRgMBnR2dmLMmDGizc/e24t/nHcebPv3izbmYKCfNAnz/vhH1Mydm9U4fr8ffX19KC8vF2lmmdPe3o6qqqoIQQ1fGYevjqPFWBBkscTYarVCoVBEbCamQrzCj/COH+FVeD+XjAoppjsMRHdoyMQXwWw2QyaToby8PKEvQqbwHId/3nJLzgsuAFj27cO6Cy/EmAsvxJwnn4Suvn64p5Q1sWK6gvGPTCaLuJoR9gAEEe7v70cgEADHcaAoakDMOFGaYCzCLUDTQer4ERtJdIeI8M0Jm80GgiCg0+nSEtv+/n5YLBYolUpUV1cP2LARC4/Hg38tXYqWrVsHZfzB4uh776Fp61Ycf8stOOn++6EpLk5rtZdLF33pzIUgCNA0DZqmY4qxsCqOFuPooo94YiysTMXi597xQxLdQSZeQUOqqwee52G322G1WqHRaFBTUwOZTDYoc3W73TCbzWhavx4H//KXQTnHYMMFAvjuxRdx+O23Memuu1Bz4YUDYqDDXeKbKtnOMVyMwy05UxHj8Ncr05VuJvMN/zec8PQ2oWpQqNYLF+LwlXKukrszy3MSFTSQJDngsioajuPQ19eHvr4+FBQUoK6uLqU3UrqpRkJs2Gw2g6IoMEePYs9jj6X8+FzF29uLPQ8/jI533sFpf/gDCqdMCRUyhFeVCcKiUCiGJE0rF4gnxgAiYsaCGAthC6VSGSHG6YYpsp2z8LdhWTYkstG5xk6nEzzPp+1YN5RIoisyqRQ0kCQZt3qMZVnYbDY4HA4UFRUlLdUNJ538znCxpWkaFRUV8HZ24m/XXw/u345cI4Ger7/GuvnzMeHyy3H673+PkspKAAPNb1wuF3w+HziOA8dxonn05huxxLizsxOlpaXgOA6BQABOpxN+vx8cx4EkyYgrCIVCMehiLJwXGLgfInYoZDCQRFck0iloiCW6DMPAarXC6XSiuLg45VLdcGKV7caap7ARJ5fLUVlZCYVCAZ/djnWXXgqv1ZrWOfMCnsfBf/wDjZs348R77sEJd94J+t+rtvC4uNfrhcPhQHFxMQKBAHw+H/r7+yNsIcNXxkPRZDMX4DgOMpkMFEUNWBmHb+C5XC7YbDawLAuSJAeEdMS65A8X3XhzzWUk0c2SWAUNqVSPhbuDWa3WtH0REo0ba6UhXHqZzeYBG3Ecw2DDVVfBeuhQRufNF4JuN3Y8+SR+fP11XPHfT0P7i/MGHCOs3BQKRURDTqF7hd/vh8/ng8PhGNBkcyg7Hg8liVaPghDHEuNoo6B4YpyuYXu897hwXy7HcwFJdDMmVkFDqmJJkiRYlkVXVxe8Xm/Gvgixxo1e6YZnPahUKowaNWpA1sMH99yDlm3bsjp3PnHauApU3LsIwXkXwL/sWfBlP4Uc4hEuxuEIl9xCuazdbo8Q4/CVcb52g8gk1k1RFFQq1QBrSkGMw+Pr4WIcHTOOdd5EG3uJBDlXkEQ3TbIpaACOJeD39vbC7XajqqoKFRUVohqbhCekOxwOWCyWhFkPe158EXvXrBHl/PnA3DkzcULLVwAA2SebQe/5DP47HkVw0c0ZjUeSJJRK5YBcaUGMhU7HfX19ccUl11dmYpKKGLvd7gFhivDXLNGXgCS6I4hsChoARLQwLykpAcuyccsqM4UgiFBKjdVqhVarTZj10PSvf+GjBx4QdQ65zAknHo+55n0gwvxmCVc/lCvuA/3+2/A/8CcQBnFa9cQT40QrvfBL7lzKGR4KUhFj4csrGAyira1tQNEHTdOhzIZcRioDTkCsDg2Z+iIACLUwZ1kWra2tGD16tGhz5TgOLS0tYBgGRUVFKC0tTbiCMu/fjzfmzYO/v1+0OeQyEyaNxeVkD0h3/PZCvEwO+xW/An3nI4BscL17ownfkPL7/XA6naHNq1SLGAaLRJ2ahxqe59HW1oZRo0ZFpLN5PB4sX74cLMti0qRJmDBhAk488URMmjQp4Xgsy2L58uVobm4GRVFYsWJFRPeVbdu24cUXXwRN01i4cCEWLVoEn8+HpUuXhnLnV65cGbOlVTwk0Y2BkPvndDrBcRxUKlVGvghC7qvBYIj4Bud5Hk1NTWhoaMh6ruH5vARBwGAwJF1BszYLXp13BsyHj2R9/nygdvQoXKvzgLKnlpnBjp4A/++eBzv9lEGeWWw4jkNnZydGjRoV04Usurx3MI1vBJHLFdENf21isW/fPsjlcrS2tkKpVGJ+Ev/lDz/8EB999BFWrFiBXbt24dVXX8VLL70E4Ngm9/nnn49169ZBpVLhyiuvxOrVq7Fp0ya4XC7ccccd2Lx5M/bu3Yvly5en/Byk8EIYQtoXy7LgOA4+nw8Mw6Ts2iVkCAgtzIV0rGhSSe1KRnifMyGf12w2J/9i4Dgo778WS3RO7D//NGz6aDcC/kBWc8lljOUGXFPGgTKlngpHNR2E6sZzEbzsRvjvehzQihsGSkZ4zDJRRZkgxHa7fYDxTfjKOBsxzrWCESEUEw+FQoEpU6ZgxowZKY131llnYd68eQCArq4u6PX60H2NjY2oqakJWVrOnDkTe/bswddff42bbz62BzBnzhysWrUqrecgiS7iFzRQFBWqXkr2+PAW5oPpiyAUT9jt9gH5vKn0SaOefwzklx8DAI639GDi1HLsK23Apo92IxgYOUURAFBYXIgbxxeCbm9M+7EEz0P+9l9Ab38f/t8+A+YXCwZhhrFJJnSJvBbCV8bxzNLTaSOUa8UGiXJ0hSvUdOdL0zQefPBBfPDBB3j++edDt7tcroi0QY1GA5fLFXG7RqOBM82O2D9r0U1W0JCocgwY2MJ8MH0RwptKFhcXo6GhYcCbK9kKmvxkC6g1T0fcRll6MMXSg0nTK/BjyWhs+mg3mBEgvgqVEjdNLYei7XBW45C9XVDdexWCZ14E/2+fAW8YfNvHTFeXiVzIwn153W73gJ5u8Xwphsp3IVWSFUZkmpa3cuVK3H///Vi0aBE2b94MtVoNrVYLt9sdOsbtdqOgoCDidrfbnfaG+M9SdFPNsaUoKqboRrcwT9UXIRPSqVRL1JwSHc2gf3sTiDiiTJm7MdXcjeOmlmFf8Whs+fgrMMHE/hC5Ci2j8et5x6Hw0F7RxpR9tBGe+n54774TGuYs0caNhdiX9OFiHH2eeN2Oo0MTuRJmSJQSlknmwoYNG2AymbBkyZLQ3o0wRkNDA1pbW2G326FWq7Fnzx7cdNNN6Orqwvbt2zFlyhR8+umnmDlzZlrn/FltpKXboUHIqRWC9tEtzIuLizPeST569CgaGhrinj+6g69Op0t62WSz2QBg4E6q3wfZ1WeA3J+6CLHGKnyvq8PmD78Eywxs6Z2zEMBvFpwM/T5xOxa7fjkVzc+Q4Ek3CgOLUeK/B8QgrVmG20w93JfC4/HA7XaH3ufRq+Kh9qVwOp0IBoMxswV8Ph/MZjMmTpyY8ngejwfLli2DxWIBwzC45ZZb4PV64fF4sHjx4lD2As/zWLhwIa6++mp4vV48+OCDIU/rZ599FgaDIeVz/ixEN9OChmAwiM7OTlRXV0esNovT9GmNRVNTE+rq6gaME93BNx3P3b6+PnAch9LS0ojb6d/9CtT6VzOap7+0DPtLx+SN+N580amo+vELUcfs/8U4tL6kBU+7Qrcp2ONh9K4AzYsvjD6fD3a7PSc6WHi9XjidThiNRvA8P6BrRTAYjCiFHmxfCqFfm06nG3CfUBE4fpCalorFiA4vZFvQwLIsfD4fWlpaUFpamnBlmi5CvFgQ3WAwCLPZHCoLzqSDb6zwArn+1YwFFwAUVhOmW004/sRqfF9Ygy0f7cpZ8b1qwWmo+vFzUcf0zGpA20s68LQ94nY/9QPaVVdBa3sIWv5UUdvj5MqlPBAZ0yUIIq4vhSDGsXwpok1vsnluibwV8qEwAhiBoitGQUN4C3OCIEQVWwFBdAOBACwWC7xeL/R6fVZlwdHZC8SB70A/ebco86VNHZhh6sCUk0bhu4JReD/HxPei82Zj7H5xBdc3rQbNa8vA0ZbYB1D9cOkfgq9vMbiOhQBPDMgSyMQwPddEN9mXSaql0IIvRTo+C+nMJx9KgIERJLrxOjSk8+YVYkJCC/PKyko0NjYOygeA4zh0d3eDZdmsxVYgInuh3w767itA+H0izPYn6J52zOxpx9STR+FbTTXe/2gXODZ+hsdQcOYZJ2Ha0S9FHdM/rhxNb9aBlSXpgkxwYEr+F8rCJhi8T4APFkRkCURvTKXi0ZtLoptNylgiMQ63g4znQBbLmzffHcaAESC60QUNsUzDk+H1etHb2wue50OluoP1pvf7/TCbzfD5fDAYDCgpKRHV8IbjOPAcB2Lp9SDbm0UZNxZ0dztOQDsmn1CJvZpqfLh9z7CI7ymnTMWpXd+CSJDaly6B2lI0vXMcGEVnyo/x0bvQpbkWRu9T0MqmRtwXvjEltMeJjoWGm9/kkugORsoYSZIJfRaETAqr1TrAKN3v90eENqIfH6sYKdfIW9EVxNbn88Fms8FoNKZdqiv4Igjls9GeoGIirKIZhoHBYAitAsRODQoEAmgyvYjC+SwqvlGAcvuTPzALlOYuzDJ3YeYpNfhGWY4PPhk68Z0ybQLOcR4BERSvoi5oLETzxhkIqlrTfixLmtCtXoIS/x0oCl4duj1cXMOJZwspiG54ddlwXTbzPD9k505kehO+cWc2m0PzEl6fvr4+8Dw/qJ9hscg70Y1V0ODxeNIS23RamAuX7JmKo+AuxnEcDAZDKGnd7XYnLLxIB+E5mUwmsHQniNq/wVbnh+Oy6dBvM0L/0DZQNlfygbJA3t2GU9CGmbNq8bWiHP/6eDd4bvASYxrG1eGXMIHwekQbk9GpcXTzyQhq0xfcEAQDm/I5+KjvYPA9AhLauIfGu/y22+2hbrjhK76h8lsIJxcu2cON0p1OJyorK0Oe1MLK+N1338VXX30FlmVRUVGB6dOn49e//nXCcYPBIB566CF0dnYiEAjgtttuw5lnngkAMJvNuPfee0PHHjhwAPfddx+uuOIKzJkzB3V1dQCAadOm4b777kvr+eRdyhjDMKFLM0EQUzGPiW5hrtfrU7oUaW5uRk1NTdrf9rHcxcLp7e2FQqEI1XVnQnTrncKiQpg0vwKn+jHiOJIrROnOahiWfQq605bx+dIhWFWHPfIyfDAI4ltRVYYbKwHaahJtTFajwJEPzkCgrEm0MWmuBmXelZBzY9N6nMPhAMdxKC4ujrg92vzG7/fHLPEVs9uxYH6fqv/IYNPa2oqampqYz6+rqwtGoxGBQABmsxmTJ09OONb69etx8OBBPPzww+jr68Mll1yCTz75ZMBxe/fuxXPPPYe1a9eio6MDK1aswOrVqzN+Dnm30iVJMuLbPVnpa7YtzIVv1FRFV3AXIwgCRqNxwKVS+LiZrnSjxbaqqgoKhQIW/nVwyh8HHM+R/TCfuh+W7TUo+XYuDA99Cfnh7ozOnSqyzhbMQgtmnlqPPbQBH27/ShTxLS7V4boaOWhTuwizPAanoNGy9VwEyg6KNiYAMGQbutQ3oNT3WxQwF6b8uHhXVvHMb6JLfDPZvItHrpUBA/Fb07MsC5qmUVxcjLKysqTjnHfeeTj33HNDv8drc/XEE0/gmWeeAUVR2LdvH0wmE6699loolUosW7YsbYvWvBPdVMm0hXk0qYhjtJVjspCFMG66FxnxxBYAAuiATfFC4scTHlin74NtSymKDs9C2SPfQLGnJa05pIu8sxmz0YwTTqvHV5QBH2UhvuoCDW6dWgZli3i93DiKROv7F8BdcUC0McPhCT8sqsfhC+xFqf9BkEhtoydVoUtW4iv0dIvevEu1jVAmBjLDRbopY8Lq3eVy4c4778Tddw9Mr9y2bRvGjh0bElaDwYBbb70V8+fPx549e7B06VKsX78+rXmOGNEVVgdC5wTBFyGdFuaxiOe/IJwzPD5cUVGRVGwFSJIMVcglI5HYCphkj4EnvKmNRwRgH/8j7G+pUdh2CYxPHID6Y3FXedHIO5pxKppx4mmjsZsqxUeffAWkob1yhRy3nToGysPfiTqvtvcugrN2v6hjxsIl34gAdRBG758g46sTHiuGs1f45l28QoZYubPRYpxLLmPJFimZ5Ol2d3fjN7/5Da666ipceOHAq5GNGzfiuuuuC/0+efLk0DlOOOEEmEymtPd88k50Yz05kiQRDAZht9tDvgjZdNWNHjtadFMRwUzGjSb6PPEsI+3kOnjI3WmdHwBAMOiv3Yf+NQTo1nNR9VwnijYNDE+IibyjCaehCSfNacAuogTbticXX4Ik8Kuzp0K7/ytR53L0f+fDM37wBVcgQB1Gl+Ya6H2PQsOcEfe4wUwZS9ZGKDpdS6gsU6lUQ7Z5F49kXwDpVqRZLBbceOONeOSRRzBr1qyYx+zbty/Cm/eFF16ATqfDLbfcgoMHD6KysjL9ytF820gT6r8FGIZBU1MTCIIIeRWI+aYI3/CKbmNuMBgy9s11Op1wu90x6+ujxdZoNMY9TxA9aJFfCo4QJztBbRsH4wtmFL4mrmFMPPw1Y7CL1+HjT/fEFd8lF81G+Y87RD1v5+u/hHX24IQUUqEwcDVK/HfENM2xWq1QKBTQauNnPgwV7e3t0Ol0EXHj8M27eJaQgwHDMDCZTKiqGtjHjud5HDlyBDNmzEh5Hk8++STef//9iJjs5ZdfDq/Xi8WLF8Nms+E//uM/8O6774budzgcWLp0KTweDyiKwiOPPJJ2B5i8E10AoT++UD5LEAQqKiriblplg8ViAUmSoCgqtJOr1+uzNil3u93o7+9HRUVF6LZwsVUoFCmJervsV/CQ4goSAKicDTD81Q3d85+KPnYsPFX12MXr8OkXeyPE9/qLTkOdyH4K3asvhPls8eLCmaJgpsHoWwGaj3SoyqWMgba2NlRXV0csZKI374R/AXE27+Lh9/ths9kiPjMCLMuiubkZ06dPF+Vcg0nehRd4nkdXVxd8Pl+ofLanp0e0nNfocwmXW4WFhaKalIeHF6LFNtXOEw5yw6AILgB4CxrRdhfQeeM8lP2dRMmfPgE5iEUP6s5mnAFg9tyx2MkVYvunX+Oy808VXXBNz8zPCcEFAD/9LTrVV8Po+wNU7Imh23OpIi3WXAZ78y4eI8F3AchD0RXCCOGXM0Jal1iEp5nJZDLodLqUUlDSQZizEK5IR2wBgEEveumnkx+YJWxBF7puBszXnQ7DRjWKH/8IlGfweqop2o5gHoATl58FpdkMiLi/Z3nsbJguFi8PVww40oYe1W+gC9wKXeAmEMi+f57YpJNJkcrmXV9f3wCvhXDjm3gkEt18cRgD8lB0AUCpVEa8MRNlGKRDePsdIc3M6/VGtOwQAyHFzOVygaKojHqq9dBPgCPS682UDUF5N7ouA3ovPgGlH5ZCv3wbqD5xXxcB9rSJUJd+BqLYD3b6dOCfdlDfZ+cjYb3/DHRd0wEQuSVoAACCg12xGn7qBxi8v8+pla4YpLJ553Q64ff7E1beJTO7yXXRFb408lJ0owsiKIrKaqUbntMbnWYmlqADA1PMFApFzE2BZPSTm+Gmtosyp3RhaAtM51lgOnM8dJ9VoPLRHaC7+kQbnxtbBeIkE4h/l8FS7r3g51BgzpwF8u9HQHbFsVlMQN9tp6HzV70AkTtWlLHw0l+gXXkFQNwHhf8EyGSyYS/BHUxieS0IZf6CGId3OhYeI/R1C7/azWXRje6/OCL+ohRFIRhMv5mi0Ousr68POp0uZk5vNpVjArFitjRNo6WlJe2xGFhhov+Y1XxEQeaC/RdH0D+vDsV758Gw7AvIG3uzGpIt0oJZQELhj2yXTvAs6MBO8JdrwfhOA/XGHhDu1CwrHdechI77HACRH802eZkZqF8Op+UmuHrmg2GYCOOb4U7bGmzCOx1HV95ZLJZQ3Di8n9uePXvQ3t6OhoYGKJVKVFVVJX19EvkuAMDatWuxbt26UFugxx9/HJWVlVi6dGmounXlypUx2waFI6xut23bhvUtah8AACAASURBVIMHD+LEE0/MT9GNvvRKN6Yb3cY8Vmfd8LHFKNeNjtkK/r/pYqKfBEc4MprPYMCRblhn7oNtaxkKf5yNsse+gfLbtrTH4WkK7C31UHh/iHsMwbhA05+DXVIGrmMq6Ld3JYwWOC+ZhrZHA+CJwXVaEx2SAWN8GZriFlT6HgbHyEMZAuErv/BOvkJbdTHDErkU5hC8VlQqVUQqnfA6cByHo0ePYseOHTCbzVizZk3CLr0bN26ETqfD008/HfJdCBfdffv2YeXKlRH+DWvXrsW4ceNwxx13YPPmzVi1ahWWL1+ecN6CrqhUKvT19WH79u35KbrRpBoCSKWNeTSZbNJF5/PGitlm8mbuJ7fCRX2U9uOGAp7ww3H8fjjWaVHQeDHKntgH9edHUn68/1czofSmVuBB+U2gDCb47q4Dv0cJ1WcDd9vc50xE658o8OTQxb3Fxi3bigB5CEbfn6ChRw9oqy5kCgQCgVDDxliZApmGKHLNdyFWCIEgCFRWVuKMM84IFSqlQjLfhX379uGVV16B2WzGvHnzsGTJEnz99de4+eabAQBz5szBqlWrUp77rFmzMHXqVBw9ejQ/RTf6jZAsphvexrykpCStarV0VrrRYjtq1Kis83kFGPShl35KlLEGFYKBc8x+OF+loO28CMaVzdBuib96BQDmhlOh5NNvJqkMtABTgMDJ00BssUK2/5gJjvPkOjQ/rwbI/kyeQU4RpFrQpb4eet9D0DLzQ7cn8+gVWuSEZwqEhydS6emWSyXAQOLsBaFgI1WS+S5ccMEFuOqqq6DVanH77bfj448/hsvlCmVlaDQaOJ3pfaGr1WpMmTIlP0U3mniiG93GPJPS4FSMaQZTbAV66RVgCfE2rAYdgoOr+iBcfwbUjyyA8b97UPi/ewYcxs6fDkr7ZVo+DNHIvd+C/wUN5pzZCP7Qj7YXCwFZ+htuuQpPeGFW/Q6+wLco9d8HAvHfW4kyBYTOFQ6HI2ZlWXSIItdWusKXR7z7xPJd4Hke119/fUhg586di/3790Or1YYymdxud8LwhYDwRbFq1Sp0d3ejoaFhZIhu9Go0uo15WVnZoLx5hkJsAcASfB9OxT9FH3eo8BgOo+VJQLn0PBjWOFH80rFVLTu9AeS4IyBEaHBJ8Az8qia8++cbUUU4UIwPsh4z13DK18NPHYDR+0fI+Mq0HhtuBC4ghCjC07aEDWlhJSz4L2RSzCA2YhZHJPJdcLlcWLBgAbZs2QK1Wo1du3Zh4cKFUKlU2L59O6ZMmYJPP/0UM2fOTHoeYb7jx49HMBhER0dHfpYBcxw3IFvh6NGjqK2tjWhjXlRUJMob5ejRoxgzZkzod7E8GBobGzF69Oi4c/R6vTBZGuGv+w1A59EqNwlydw1K3pKhNPA9KG/qfcgSEZSrsf6mu2BSHct8qGbLUE98BJocOSteAZIvhMH7ONTs6YMyvhCicLlccLlcoGk6IkQRvnk3lOGHRAbmHR0dqKqqSmn1CST3XdiwYQNef/11yOVyzJo1C3feeSe8Xi8efPDBUMrns88+C4PBkOAssRkRohsIBNDY2Ai5XA69Xo/CwkJRv5WPHj0aMrUQy/AGAJqamlBbWzvgG9rn86G3t/fYt3fdS/Ao3s9q/rkGzxM4yF2OIBPErN0H0PDVxqzG40Bg042/RbMuMtVMyakxkfOiiP44q/FzEp5AUeB6FAduA4HByU91u93weDwhYQnvVSaEKga7c0U4ra2tqK2tjXlfW1sb6uvrc7ZHGsdxob5ueRleEP6gQmfdQCAAkiRRX18/KN+8BEGEWv2oVCrRwgjR8eLw5pVGoxF8wTfolI0swQUAM3sReuhegAI2nV6O8pnLcPInO1B3MLOCj+2L7xsguADgIz34hgRqmMWooz4ARQxNq6IhgeDhULwKP/UDjL6nQPGlop8ieiMtXogiUeeK6CyKwQpR5HoZcPjrmJeiy7Is2tvbEQwGYTAYoNVq0dzcLHrNuhBGEOJdYsdshZiZ3+9Hb29vSGw1Gg1YONFCPy7auXIFHzsZB6hIG8oetRXvnj8elbNPxpztn6CsMXVv4D3zl+D7qsRZCm20CWZuNiZyDhRSn2U071zFR3+NTvXVMPiegoqdkfwBaZDKRloi8xthRez1euFwOBAMBiM6+KbitxA+XiJyXXTD55eXokuSJEpKSqBWq0NvCiFXV4wXPjxmq1KpoFarsw4lxDtPd3c3WJYNia3wfMz0M2CI7Cq8cg2WU+N7rh4cFTvVpktnx1sXTceYvjMx+4N3UdyZ2GD80OzLsGNCauf2ki58DQp1zGLUUv8EmUMFJtnCkhb0qG5DceA2FAWuBwFxVpPZpIwRBBES1nDCS3z7+/sRCATAsuyAQo/oEEWiQg2h0CiXRTd8bnkrutFeo0LaWDbWi9FiK6xsOzo6RLWOFLqVer1e6PV6lJaWRryh3MQOOKj/E+18uUIrewHcMnPigwgeR0usaFo0FxN6F+CUrX9HgWVgS/TOiXPxwcll4NMs722hTejl5mASZ4WWGhxbzGGBYNGneAF+8nvofY+DQkHyxyRhMPqjxfNbSBaiEEQrlvgK8xzu7Ip4mEwmvPfee6HCirzcSAMQYYIBHMu5KywszMj4OVpso03Ksxk7nGAwGBJbg8EQ8ukNL2vk4Eaz/FIwxOB26x1qHMwv8A2d/hcXxckwqUOOU/71BlT9xwTbVj0Jb192Fnxkaj3h4jGaKUM1tRmUSF03cgWaq4LRuxIKLsXLgDjYbDbIZLIIm8ahJDpEIWRShK+iCYKA3++H2+3GlClThmWeybBarXjhhReg1Wpx9tlnjxzRNZlMUKlUKaeMAMnFNnxstVqd8ZtPyBsWdoILCgpAEETMcXvoJ+Cg3s7oPLlKkKvGbmIyAkRqJjWxoBk5Jh3yYcrX2/B/V1wCl1yc8l4tV4iJXDc0tLj914YbgpejxH8/CoOXZjyGxWKBUqnMibZBXq8X/f39KCsriwhRCOW6LpcLNTU1GD9+PK666qqE5cDJzG42bdqE1157DRRFYdy4cXjsscdAkiQuvvji0Ge1uroaK1asSHn+nZ2deOedd9Df35+/ohsMBiMu+S0WC2iahk6nS/rYVMVWQMjLS2XscBiGgdlshtvthsFgGJDKJvRAKyoqAgB4iN1ol92Sm56vGcLzFPZxl8JMJQkrpDIWR+Gf5odwvHI3qou2ijA7YWACDawBVeR7oLJcPeca2uD5KPU9BBKpdakOp7e3F1qtNifSsKLT18LxeDyw2WwoKSnBoUOHMHHiROj1+rhjrV+/HgcPHsTDDz8cMrv55JNPABzLIFqwYAHee+89qFQq3Hvvvbjgggtw2mmnhfJ3s2HDhg35GdONRSqeuvFitqmMnU5MN7z8WK/Xo7y8PG4XY2FcDh70yB4bUYILAD3sL2GmsxdcAPjWcRu2sgZsdV+AswOzcI5uNWSynuwHJng00r3oDp6FyVw7NPS32Y+ZI7hkW+Dm9kFt/h3U1Ji0UrdyyXshWTUaTdMoKytLqcNLIrMbuVyOt956KxRzZhgGCoUCBw8ehNfrxY033giGYXDvvfdi2rRpKc15z549WLduHYqKilBZWZm/ohvL9Caep2602Kbb64wkSTAMk/Q4lmVhsVjgdDpTKj8mCCIkumbqeQSJjpTnlA94mBk4RNlFGcviOQOve8eGfv8gWIKvzL/FDQVfoq7gLVHO4ZE58RVfjDHMZaii3gORb5aQceAVrfBU3gnWcjccvSeDYZgBrXJiVZflkveCmF0jEpndkCQZWiW//vrr8Hg8OPXUU3H48GHcdNNNuPzyy9HS0oJbbrkF//znPxOa7AivnUajCTX3VKlU+Su60cSyYOR5PqKoIdPGksnsHcMtI0tKStDQ0JDSm5UkSQSDQXiIb2Cn/jfteeUyHF+EH8lq8ET2Tl/BYDledvwSXFQqlB0k/ss5G7O8U/BL3V+glDdmfS6e4HGEtqCXOx/HcUehohI7pOULPOmBz/gUCnWLUeK/BxxLRKRuRVeXKRQKsCybF6LLsmza9pXxzG6Ecz399NNobm7Gn//8ZxAEgfr6etTW1ob+r9PpYDabY3YmFhBeu+bmZpx11lmYMOHYxuaIEd3wEIBYYisQz94xWmzTdTEjSRIcfOihHx1xYYUmdj7cdPaX/jxH4p2+O2Hm47+uOxktvrPciWvVP2BC4VoQZPYGOg6yD7t4A0Z7L0K18n2QedJ5Ihn98r/DT+2H0bsCKqo8buqW8NPV1RVzVTzUObGJ0kE5jkvr853I7AYAHnnkEcjlcqxatSr0eV63bh0OHz6Mxx57DCaTCS6XK6nvwqFDh9Dd3Y21a9di0aJFIdHN2400lmUjLvkDgQC6u7uh0+lCYmswGERpme71emGz2UI7ohzHwWq1wm63o6SkBMXFxRnFvlwuF8z0fyJQtC7rOeYSNuY8fEdnnqkQznd9v8Zab+qpT8dTPiwuehNapXhx2cKADhP4fdAoUjdlz3VIrggG35NQswNFR6C9vT30ng/3W/D7/aHVZbgQD5bnApB4U89sNkOtVqO8vDylsRKZ3UyePBkLFy7ECSecEHou1113HebOnYtly5ahq6sLBEHg/vvvx4wZiSsAOzs7sXv3bjz33HOYOHEiJk+efKwAaiSIrtAyXRBdscRWQCjTraqqCvVUKy4uRklJSVYbDX3+3egtuBUgxCu8GG58TD2+osaCIbJv097nOR1P2i8Dm2aFFQ0eV6oaMb3oZZCkOHFZkqcwni1EGbURRI43uEwZnoQucCN0gVtBYOD7OJGrl9BAMlyIowsawlfF2YpxT08PiouLB1S4AcdSOnU6XcKMheGCYRi88cYbkMvlsNlsMJlM+Su6gv1ceBjB4/FEWDCKhd/vR1vbsb5fOp0OpaWlWe/qcgigmb4MDNUiwgxzA46X4dvgAjjk2dtQBhkjnjUvQw+f+WVsAxnENYX/h2L151nPR6CU1WM8voOCGjmrXiVzEoy+P4DiiyNuT+TqFY/wggbhR/DjjV4Vp/MZ6urqiruY6urqgtFoRHFxcYxHDh9vv/02FixYAKfTCaPRGLo9b2O6fr8fzc3NETHbo0ePinoOjuNgt9thtVrBcRzGjBkjWizLSr00ogQXADqZi+CQD3T7SheeI/Cu7a6sBBcAGjkZfm+/HJd4T8Vs3WrQVPZ+C1bKgl18LSYwE2Gg3gMxAq5SfPRudKqvgdH7FJTc1KzGSuS5IIhwdHPNcCGO11xTTAPzoaKqqgoqlQp/+ctf0NXVheLiYowdOzZ/RVcmk2W9QRYPIVxhtVpRWFiI+vp6tLS0iPaH9RH7YaNeFWWsXMHFnoKjtDjWiT/Yb8bnjDilpzwIvOOvws7ex3B9wQco127JekyWYLCPdsDALsQ4Yg/kZLMIMx1eWNKEbvUSlPjvQFHwatHHT9S5wu/3w+fzweFwgGGYCOEWxDhRq55cFd0dO3bA7XajvLwcZWVl8Hg8aG9vz1/RpSgqpuBm0zY6XGwLCgpQV1eXcSfVeAQZL9plDwEjJS4IgOX1+JHQAyJ4GFicJ+Nv3skQySgrRDdP4Y/95+Ec70k4W6SiCjNlRh8/Bscxk1FKbQKR7xkoBAOb8jn4qG9h8D06+KcLa64ZXgov2J2GO5EFAgF0dXVFCLGwcZerto6VlZX44YdjKYcKhQI6nQ4VFRX5G9MVYkfhNDU1oa6uLu14K8/zcDgcsFgs0Gq10Ov1A8Q2umVPugiFEw7lGvDGf2Q8Ti5ymLsCnWT2IhYIluK/LMvRlWVYIRnF4HB9wU7UFfxdtDGNrBHjiJ2Qke2ijTmc0GwN2Ka7UFc+d7inAgBoaWlBVVVVxKad3W7HihUroNfrcfLJJ2PixImYMmVKQq+IZL4L27Ztw4svvgiaprFw4UIsWrQIPp8PS5cuhdVqhUajwcqVK1FSUpLSvNvb29HZ2Yl9+/ahu7v72KIwX0UXOBbXDae1tRWVlZUphxyEfF6z2QyNRgODwRB3ZZup6ArpZQ6HA4VGG/r0vwKI5NVt+YKVuxDfk9mbz/Acgbd7HsEOiN8BIR6zaBcu0v0PVPImUcaT8QpMZCmU0NmHMHICTg69fxkKmAuTHzvIxNvU83g8+Pzzz0EQBA4fPowpU6bg/PPPjztOIt+FYDCI888/H+vWrYNKpcKVV16J1atXY9OmTXC5XLjjjjuwefNm7N27F8uXL084X5fLhQ8++ACXXHIJgsFghCblbXgBOHZ5Ev6dkaqnbnRZcG1tbdLHCOdKNXTBcRz6+vrQ19cHnU6H+tE1aFc8NKIE188dh32EOAYxP9pvGFLBBYSiirtwjfp7HFf4atZFFUHCj+9ooIJdjAbiM8jILpFmOkyQAVhUj8MX2ItS/4MgMTBdayhIVI6sVCpRX1+fNGdWIJHvQmNjI2pqakIGVDNnzsSePXvw9ddfh7xw58yZg1WrVqU055KSEjidTixfvhzFxcUwGo2YNGlSjOS8PCZZua4gtk1NTXC5XBg1alTKK+N4VWmxzmGz2dDU1ASWZVFfXw+9Xg+77FX4yYPpPJ2chuOV2I/jwIrwJWL3noLXfInNQwYLDwi84pmKteYV6HNNFGXMbsqE3cQU9DHniDLecOOSb0S3+j8QJIYndJIocyHdeK5Go4FWq43pu+ByuSJiyxqNJtQRWbhdo9HA6Ux+ZVdQUIDZs2eDIAicd955qK6uRn9/Pz777LP8XulGE88NjOf5Y9VfZjMUCkVGvc4EQY/3Bw6PCxcUFKC+vj50rJ9ohJV6Jf0nlMN0chfDTmUfx2XYEvyPfTEYsXfO0uR7Von9jltxmecQTtT/FVSWRRUBwodvaaCSWYwG8hPQpEmkmQ4PAeowujTXQu97FBrmjCE9t9jpYvF8F7RaLdxud+h3t9uNgoKCiNvdbndKnt1C77j29na43e7QSrm3tze/V7qxnMaiV7oulwvNzc1wOByoqqpCVVVVRr3O4q10hbhwU1MTvF4v6urqUFZW9lN7EbDooX8HXoQKrVzBxc7DUREEFwA22+5GxyBvnKUKQ5B4izkOf+5+AibHSaKM2UWbsBsz4GDOTH5wjsMRLvSqlsKqeA48hi5MJqbDmOC7sHTpUlx22WUR9zU0NKC1tTWUR7xnzx5Mnz4dM2bMwPbtxzpVf/rpp5g5c2bS8/T29uKjjz7CmjVrsGPHT22hjEZjfm+kRRuZ9/X1geM4lJaWwu12o7e3FzRNw2g0xiwfTIeuri7odLpQnmH06tloNMYMU9iotTDTz2V17lyC4SuxG9PgJzxZj3XIcSNecg9PWCEZBHhcIu/CbN1LoOnsndIAYBRTjjryA9Bk9gUkww3hOQ4q88NQ0ZWD7rvg8XhCjQCicblccDqdGDduXEpjJfJdWLx4cSh7ged5LFy4EFdffTW8Xi8efPDBUDODZ599NqnZjcvlwjfffIOXX34Zvb29mD9/PkiSxLRp0/JbdBmGiVjZOhwOuFwuBAIBUBQFo9EIpTJ9x/xY9PT0QKvVhi41ent7IZPJYDQa466cA0QzWmSLwI8UX1aexEH+cvSIcKns8J2IJ23XIDjMYYVkVBIsriv8F8o174synpJTYxLnQSH9iSjjDRs8gZ0tyzDKdzbGUMc+c0LerdhuZE6nE4FAAKWlAzdaBVvKcBHNFTiOw65du2CxWEBRFA4cOACHwzFyRNfj8YTamdfU1IgmtgK9vb0gSRIulwskSSYVdB4c2mQ3wEeOnC4EZvYS/Ehl76vAMDo8b34EbXz+bCmcI7PiLN1qyGXixGZrmTLUUltBEeKYvA8pvAJ/P7gCTx24ADqax7rpHpyk40J+KNFuZEKpryDI6a6KHQ4HOI6L6a1gt9vBcVzaHhGDjRCH3rBhA5qbm7FgwQJUV1cfe+75LLosy8LpdKK3txcAUFRUBLfbnbApXSb4fD60tx/bua2uro7wII1HH/UGeuk/iTqP4cTHTcUuQg9OhEq6dzofxqdE8rYquUYxONxQsAO1BeIUtyhZNSZydhTJ8qcVPMEV46mvX8Rb7T/FNbUUj/+d5sXckoHvjWiPXqG6LJ1VcV9fHyiKirmBZbVaQdM0qqurxXuSIrJnzx589tlnoUIutVqd36JrtVphsVhgNBqhVqsRCATQ09ODmpoaUcYXLB0ZhoFKpYJMJot5iRNNAO1okS8En0X321yCYdX4ljgHTjL7VdkRx3V40X2CCLMaPmbTLlwoYlFFhVuHMcoPQFPu5AcPIzxTg9u+WIOd1roB9ylJHq9N8WK+IbUv5WSrYuFHJpPBZrNBoVDErDRL10t3OPB4PGhtbUVzczM+/PDD/Bbd6JguwzBob29HfX19VuMGAgGYzWb4/X4YjUZotVo4HA4EAoGkAXQePNplN8FL7slqDrnEQffF6NZkL7hO33Q8absB/hyP46aCBjyu1XyP8QVrQZDZO41puEIcx/agQLZbhNmJT8A3HZd/8iJaPPEXHTKCxyuTfVhYnllmQ6zOFcFgEAzDQKlUQq1WD1gVm0wmFBUVJf1cDhex0t3yWnQ5jotoRsnzPJqamtDQ0JDReELLdKHVc0FBQSj25HQ6Q45Biegj30Kv7KmMzp+L9AXPxLey7EMKLFuIP5sfQwuXP3HcVJhK+XB50RvQKr/PfjCeQI1Ph3rFv0CSubPqtTnPxYUfPw0XkzwDiASP5yf6cV2VeO2Nuru7oVarwfN8KETBsizeeutYQ9IZM2Zg5syZqKmpScl35bvvvsMzzzyD119/PXSb2WzGvffeG/r9wIEDuO+++3DFFVdgzpw5qKurAwBMmzYN9913X9JzCGK7cuVKdHZ2Yty4cRgzZgzGjh07skQXyMwjIbxlusFgQGFh4YBAv9vthsPhQGVlZdxxguhCs/xS8CKkU+UCvmA1viIngqGyzzF+3/wotgaHtsx3qJCBxxUidqrQsIWYzHdCTX8twuyy43DvjVj0+VLwaaT0E+CxYrwfv64RR3i7urqg1+sjsoR4nkdXVxd27tyJvr4+tLS0wGg0JvVEWLNmDTZu3AiVSoV//CN2bH7v3r147rnnsHbtWnR0dGDFihVYvXp1RnPfsWMHjhw5ApvNBp/Pd8y6Mp9FN5bTWDqiG95YsrS0FDqdLu6uqs/ng8ViSRiwb5ctgYfcmfoTyGF4nsQP7KWw0pasx2p0XIM/u8UpNshlxpBBXF34DorVX2Q9FsETGMPqUUm9B3I49gZ4Eh80PYT7v7s24yGWN/jxwOjsv7A7OjpQUVERc6Otra0NdXV1obbqydi6dSvGjx+PBx54IKboCvm5zzzzDEaPHo0tW7ZgzZo10Gq1UCqVWLZsWcbpaZ2dnbBarfldkRaPZN8jHMfBYrGgubkZFEVh9OjRKC4uTpjGksx7wU6uHzGCCwAmThzBdfmmYY37RBFmlPsc5WT4vX0RPrM+CJZNXiqaiGOt4M3Yy8+Hl82um0O6ELwaL+3976wEFwCebFTgkSPpV39GI2YZ8LnnnpvQI3vbtm0YO3ZsSFgNBgNuvfVWvP7661iyZAmWLl2a3uTxkx5VVVVhypQpI8t7ATgmjvHcwKKdv9JpmZ7ITCeIHpjpZ7Oady7h5U7EQTL7LhAsW4C/2q+FbwRsnKUKDwLr/VXY2fs4ri/cijLNP7Mar5/swy6+FGOZy1BBbQQ5yOXkXLAUD+9+GVtMx4sy3n+1KOBkCPznBD8yLVhL5O4ntoH5xo0bcd1114V+nzx5cmj8E044ASaTKe1GCdHH5vVKN9YTj+W/wPM8+vr6Bjh/pWN2nmil20P/HpwIXRNyAZbX4QeiArwIvb8+6rsbTZz47ZTygS6ewgrH+Xjf/CgCwexyknmCw2Hagu/4BfCxk0Sa4UDYwBhct/1t0QRX4C8dcizZpwQrciCT53nRW/Xs27cvwibyhRdewGuvvQYAOHjwICorK7Mud877lW4iT91wk3KtVptV+x1hBR2NKfgPeBTidZsdblr4+XCT3VmP0+S4AlsCuZnGM5RsDZZit3mZKEUVdtKGXXw5xjHjUUa9C1LElk8e7yxc8vGf0eMTpzddNG91y+BmgbXH+yAXaaknfB6z6cz93nvvwePxYPHixbDZbNBoNBGieuutt2Lp0qXYvn07KIrCihUrsp53Xm+kAQh1FRXo6upCYWEhOI4LJU4n6giRDuGbdF6vFz2WAwjU3w5QI2OV28edj2/J7DMvXP7j8QfrzfD+jMIKqXAq7cSFuv+BUp59I8sSthTjsQ9KKnuP5h77xbjokyfhH4Krkl+UMnhzqhfqFBenPM+jvb09ZsETwzBoa2vD1KlDG/POlrwOLwCRIQae58GyLLq7u+F0OjFq1ChUVFSI2lzS7/ejra0NJpMJVO1fR4zgBrkx+JHIPsWHZbV4re8GSXBj8AVTgN9b7sYB+43guew+ejbKil1kFUzMxeCzsMb8tut2nLtt5ZAILgBss9K49BsV+lOsnxDTwDxXyHvRFfB4PGhpaYHf70dxcXHGvrnxCAQCCAaD6OzsRGlpKUpG74dX9qlo4w8nPC/HfkwDI4LofmK/G0d+pnHcVPCAwMueaXil6/dwuLPrVMERLPbTdvzAXYoAm2b/Pl6GDYf/iOu/vCOrOWTCDjuN83bSONJjCzmIxbvgzsfW68nI+/CC0+lEd3c3CIKA0WiE1+sNeeqKQXiVGsuyaGhoAE/Z0SK/BGw+OkTFoJNbjMMi2DW2OhfjOeepIszo54EMPC6XHcLMkv8BlWUBCsXLcByrgp7aBCLJJijBFeI/v3serzbPyuqc2TJeFcTrY3pQyHkQDAZBEESE74JCoQh1/Y1VCep2u9Hf35+yl26ukPcbaYFAAEajMeT8JdRrZ0t44YRer0d5eTlaW1vB8zxM9FMjRnDd3Kk4TPRmPY4nMBGvOGeLMKOfD0EQeDM4AbvNK3B14XoUqzN3G2OJIH6k66QFAQAAIABJREFUgzCwCzGO2A052RrnwCrcufNlfNI7NuNziYWLoxDUlKCyQAfg2GdOMMERfHKFTCTB9EbwXSAIIm/DC3m/0o02vREayWXqOsRxHGw2G/r6+lBSUoLi4uLQ5U1bWxsOykqxj/0Ep1evQoFieBr1iQXLG/AVToKXyK7On+VUeKX3CRzixAvn/NwgwONSRSdOKVwFmSy7fQKal+M4VoZSagsI4qePdzAwGVd+shpHXMOfVXJiEYs3p3pRpkgsP06nEx6PByqVKsKN7IsvvkBTUxPGjx+P008/HfX19Snt3cTyXQCAtWvXYt26dSgpKQEAPP7446isrMTSpUthtVqh0WiwcuXK0P3ZkPeiy7IsGOanqLzX64XNZkvbU5fnedjtdlitVhQVFaG0tHRALKm1vQMXHq3DYQ8FGcFjrsGK86r+hVOrVkMtz6/GgzxP4DC/GF1k9r3O/tV9P7bw4thp/twp4wK4qfgDGDVbsx7LyBgwjtwJGdmBftcZuHDbc7Azyb2gB5srK4J4fqIPihR2lOIZmNvtdnzxxRfo6OhAd3c3zGYz/vrXvyYU3kS+C/fffz9uuOEGTJ48OXTb2rVr4XK5cMcdd2Dz5s3Yu3dvUm+HVBhxopuup250Lq/BYIh7yfLs90483K4fcLuC5HGm0YTzqv+JkyteglKW+6EHC3sRfqCy7/3V7lqIZ/vnijAjiXDOlVlwVvFqyOjsQj8UK0Ng/0l44MjVYDG8l+IkeDw+1o+76lIP//X19YEkSRQVFQ24z2w2Q6VSoaKiIqWxEvkuzJ8/H2PHjoXZbMa8efOwZMkS3H777bj55psxbdo0OJ1OXHHFFdi8eXPKc49H3sd0B5TYJSjXjcblcsFkMkGlUqG2tjZmY0kBNwP8V8/AdiEA4OcIbOkpx5aeG6ChrsfZZZ04p2ojTixfA7ks9xzH/Nwk7BchH9cbGI+X++eIMCOJaLYG9fjK/BCu136B2oK3Mx7nq2234287J2DWaD92qkhkXIubJYU0j/+Z7MV5KZqcC3AcF3f1mui+WJx77rno6OiIed8FF1yAq666ClqtFrfffjs+/vhjuFwuFBQcKxbRaDRwOp1pzT0eeS+60VAUldCYBjiWXmYymUJtPlLpFPxfzRTMweQrBTdLYENXNTZ0/RqF1BKcW9GO86r/D9PK1oKmhr9BJcersI8YD5bIzluB41R4ve8WuKR83EHDxpN4znk6TvNOw4I0iyp4jsTnm5Zh/bfHrvh2Nilw4igWPxSRQ+6FUa/i8PdpXkzQpl9aLmb79XjwPI/rr78+JLBz587F/v37Q01ogWOZErHaBWXCiMnTFYguCw7H5/OhtbUVZrMZFRUVGDVqVEqC2+sHnmtM/4/bz1J4u6MON315D87Z8i2e3rMJ35luAMcN32VeB/dLOLIUXAD4wn4X9rPSxtlQ8DlTgCcsd+OQI7WiCp6RYeu6x0OCK/BVO4W6Hg56EXw1UmVOMYOPT3ZnJLjA0OTpulwuLFiwAG63GzzPY9euXZg8eTJmzJiB7du3AwA+/fRTzJw5M8lIqZH3K91UzCcCgQB6e3sRDAZhNBpT9t4UePIwDReb3erAGiTxRttYvNG2DGWKB3BBxSGcPeotHFf6DxDk0ITV+9kz0Ehlv3HW5boE633xzdwlxMcNAi+5p2GqbwUW6d6ARvFDzOO4gBrvvPUoPm/Rxbz/oIVChY/D6NEsmgb5y/+m6gCeHu8HncXSTkxbx2jCfRfuueceXHfddZDL5Zg1axbmzp2Lk046CQ8++CCuvPJKyGQyPPusOE6Ceb+RFs/IvKGhASzLore3F16vN9TrLF2HoCMuAtM/lYHhB+eSrFrJ4MySr7Fg7NsYV/reoJwDABiuCrswBQHSm9U4vsBY/MFyO5xSWGHYkIHHVaqjmFr0Mkjyp/c+4y3Bm/9vOb7pSr6oKJDzaBjP4dssSojjQRM8nqx3Ykld9ivRjo4OlJeXx4zdtrS0YMyYMSl1584l8l50gWMFEeE0NjZCo9EkbL+TKld8TWNDz9CEA+rVAZxj+Apnj3oNY43bRRuX50nsZy9Db5Y74RynwF/NT+JHNnlIRmLwGUMGcU3ROuhUOxHor8Rf3/gtDlpS/9tQBI+Tj2OxgxTvgreE5rB6jBUzFMeKG4TNLqVSGSpuoGk65c9jW1sbRo0aFfP4pqYmHHfccaKW+w8FI0J0hdptjuNgtVphNpthMBig1+uz8r78so/AvB3D8wedoPXhgqpvcOaoNagqzLxSCQB6mItxgM4+je0L2wN42xe/XZHE0EOAxzX2Dnz5rxoc6c/svXpqA4MdSgp8lpkNEzQs/j7Ni3r1T5LC8zyCwWBEh1+GYUCSZIQQy+XymJ/V1tZW1NbWxjzfkSNHMG3atKysHYeDESG6fr8fNpsNNpsNOp0OXq8XBoMBSqUyq3HP2CHDzr7h/4MeX+jBgqpd+MWol2HU7k3rsa7A8dgjM2RtSt7lugh/6j8rqzEkxEdv4aG6Wg6NloftXAK9GQrQjCoWh4pJuDMMG52nZ/CX470oSHHRzLIsfD5fSIgDgQAIgoBcLo8Q4/b29piiy3EcGhsbMX369KxNxYeaESG6R48ehVwuh16vB0VR6Orqgk6ng1qtznjMjT0kFn2de25ZJ+hcOL/qC5wx6iWUqA8kPJbltNjDnQEPnV1+ocdfh6csd8FF5F+d+0imyM6j9HoZmg4d+7sYDRxKr+RxQJ7Z32lMCQtXNYGeNG0n767z47ExAZBZah/HcRErYuFHo9EMCE8wDIPW1lZMmzYtu5MOAyNCdKOt4UwmE9RqdSjvLu3xGA7Tt9No9OWe6AqQ4HFKST/Or96OOdWrUKQcmMPZxF6B1iyzFVhWjr/2PI59RHoZHxKDi8bNY9RNNA5+G7m0VMh5TL+awZelmb13jWoOhXUMjpLJQxUKksfzx/lwZWWK5rhpwvM82traUFFREbEqPnLkCN58803U1tZi7ty5mDBhAurr65OueOP5LmzatAmvvfYaKIrCuHHj8Nhjj4EkSVx88cUhDamurhalawQwQkQ3GAxGFERYLBbQNA2dLnbaTDyEXmqrm3j83pQ/KVEUweO0EivmV2/D6dWroFV0w86eg71Z2gUCwOfme7AuWC/CLCXEQuHnMf42Gb7fEX9FO/vCIL4cR4PL4NJbRXGYMJ7B3gTCWybn8OY0L04sGrycX6EhQXV15D4Cx3Fob2/H3r174Xa7cfDgQdx+++2hri6xiOe74PP5sGDBArz33ntQqVS49957ccEFF+C0007D4sWLsWHDBtGfV97n6QIDc3VjNadMBM/zcDqd6O3tBaUuxMv2/NosYnkC2616bLcugvz7y3GxvhOjiw7BeNxGkDJfxuN2Os6TBDfHoBgex99PY08CwQWAHe/JMPUEFi1zCTiI9MIFXpbEt/tlOHXC/2/vzMOjKLM1/lb13tk76e6EhCTskc2wiEZxFBER0RmQTWAEhMs4Xh3HuS6M6+OMC+N1xrmKiqKAIyiLQREERAF1BNkim4AR2bek9yS9d1fVd//IVFmddGftJenU73nyKN2d6i9J1dunznfOexjsClOqNSSNxYelXuSrYxuvRarRpWkaOTk5GDFiBPr169eiYxUWFmLRokV47LHHQh5XKpVYvXq1UHbGMAxUKhUqKyvh9Xoxd+5cMAyD//mf/4laKiMpRLchNE232FPX7XbDZDJBpVKhqKgIL51Rw+TvXIl5MQFC4eypPKw9X4DMz0dh9BWXcOWV25FeuLtVx/EHi/GOZ1yMVinRJgjBVU/LsefLll22hytk6G7mkH0nh9Ot7FAgoLCrUo6yHgz2p8jA/GeD7U5jEIsH+KCJQ3o/mt1okXwXeAEHgBUrVsDj8eC6667DiRMnMG/ePEyZMgVnz57F/Pnz8fnnn0dl9FdSiG64SLc5/wWfzweTyQSKopCfnw+VSgWTH/i/0517s2gwxWLf+fqfocZHY93B7lh3cA5662Zg1OAf0XvwRqgym/YB5jg51jj+GzUxagiRaBvXvSDHrk9ad8leOE8j9V2C4b9lUZHa+nN79xk5BucyOJtD44GeAfy5V/sHBLSU5uajRWv2IcdxePnll3HmzBksWrQIFEWhR48eKCoqEv4/MzNTsA9oL0khug1pKr0gbgk2Go0hFQ4vRKHdN5HIQOC6GH79J+1KnPz6StDfDEZZkRNXD9qD/AGfQaZsbMLzfe2DOMC0r9xOIrqMfJXGzpVtu1xdLgrfv03j+qksvi1ovdPYCYsMbw9z4fbCIILB+nQcRVGgaVr4byyIZQuwmGeeeQZKpRJvvvmm8H7l5eU4ceIEnn32WZhMJqHRKhokpeiGs3cUzzoL1xL8s4vCsguJr8ltD9ewHHY5mj4ROUJh19l07Dp7C9K33ozR/S/jyiu3I7OovgHD7BmLD73hi9ElEsNVi4PYuTi1XccghMK3a2QYcQOLo8NpeFoovDkagrUT/bi6mwwcRwlNSPx/AYRca9EU4uYcxpqyYm0O3ndh4MCBKC8vx/DhwzF79mwAwKxZszB58mQ8/vjjmD59OiiKwosvvhi1yDopqhc4jgvJ4TIMgwsXLqBHjx7gOA5Wq1WYdZaRkRG2tMTkBz66LMNGE41ddipmXguxIoci8P8EOANtW3fPrCDGlvyIrfm9cTqzc6dYkonryinseiq6XZF9+nHw3AZckjUtjL2zOHwyyY+emZElQizAhBDhi4cQIsw0A9AqMW7KwNxsNiMtLQ0Gg6HFx+soJKXoEkJw6tQp6HQ62Gw2ZGVlQafTtfgPbg8AW8w0PjPR+MJCw90JUg7XuFjsOds+sRxRweLIHhrDJrG4OIvFOSngTShXb6Gw/2EFOC76519mJgvj1AB+Sg9vFlOWz2LNBD+y2+Alw0fAHMeFiLEYiqKEr0jXpc1mg1KpDFtvX11dDZ1OF7Wp3/EkKURX7DTGj9+5ePEicnJyhC61tuJjga+sNDaaaGwy0TC1MZKMJf1pFj8eoUHa4fw1xMPi4Ju//J5ommDozSy897A4NiQaq5RoDcP+TeHwfyvAMLE73+RygqtnctilD70+JvVjsGRcAOooJx95AebTEQ03uxsKscVigVarDWvFeunSJeTm5ra6Fr8jkFSi63K5YDaboVar4Xa70adPdMZME0LgcrlQbTLjiE+DPawRWx1qnHAnPgdMEYI+JoIT1ravRQUC41qC8+fDH6P/UBYp81hUjOJA2tvrKdEsgysonPgvBXy++PyuR4zx4cBgFRiKwp+uCuK5XwXjNtVHLMTiiJgQArPZjMzMTGg0mkYR8YULF1BYWIjU1PbluhNBUogux3E4efIkKIqC0WiESqUSPHXba4bh9XpRXV0NhUIBg8GAmpoaqNVqpKen44SLwkYTjY3VNPbVUOAS4DF7Lcfiu+PtSyuMPM9i59rmj9G9mEPhHAYHfsPBq5HENxaUHAMuzlbC5Yrv77f/YD8mP2bHlJ42KBQKqNVq4as1VozRwO/3w2QyQSaTIScnJ+wkmIsXL6JPnz6dzksXSBLRBYC6urqQ0TunT59GcXFxm3dRA4EATCYTWJZFbm6u4FhmtVohk8kajYQ2+YHPTPV54K+sNHwxyMM1JI1jQJ+kUBtou+gWshxMb1Lwt6IhJCuLYMBvGfw4nYVNJ4lvtOh5CrD/VokaR3x/pzIZwduvBjB9MitYMfp8PuGLYRjBE5f/UigUURdiQghqamrgcDhgNBpD0gp8KiIQCMBms8HlcmHgwIHtqmBIFEkjug1Nb86dO4du3bq1+o8iLi0zGo2Nbl/sdjsIIU0m8N0M8IWlPg/8uZmGPRibi+haL4vv2jC7TcyQb1kc3Nu2Y6jVBEMnsrg8m8XZ4nYto8tTcBEIzFTAbIpvykomI3j39QCmTmy6bZ73xOWFOBAIQCaThQhxJE/clhAIBFBdXQ2lUgmDwRA2WHK5XKiqqkJmZiYKCgqiVqcbb5JWdC9cuNAqT13eAL22trbJ0rLa2loEAoEWFUpzHAd/kMW3NmCzVYEtFgXO+aJzUfWmGJw+KgPXjtK2EXUs9i1p/4lLUQTDbmbhu4fF0aHtPlyXw2AmUP5WiYsRcuqxQi4nWPp6AJMntG4sOg/vict/+f1+0DQNlUolCLFKpWrybpM3maqtrW3UrCR+H7PZDI/Hg6KioqhN5U0USSO6DZ3GWuqpy9/S2Gw2ZGZmNlta5nQ6hSi4KViWFcrYZDKZcMwjdRQ2mWXYaJbjcF3bL7I+l3z42dH2rjEty0K7gsBqje4W9RVDWKTNZVFxEwdOJqUemiOrhiBrlhKnT8RfcN9bHMDEO9omuJHgOK6REAMIEWK1Wg2aphEIBFBVVQW1Wg29Xh/2unM6naiuru700a2YpBXdlnjq8s5iWq0Wer2+RR0nbrcbtbW16NYtvPUjXzPMtzDyX+G44OUFWIaddrrFDRnDGT8qKts3p+y6nxns+jR2DYkFRfymGwufNvFVHh2RNBdB/jwFKg/HV0gUCoL33w7g17dFV3AjwZuTi4WYv0YyMjKQlpYGtVodIqgsy8JkMsHr9aK4uLjN3tgdkaQRXYZhQtoRm/LU9Xq9wu6o0Whs1WA7r9cLm80W1uOTX0NzYhuOmiCw1SLDZ2YZvrTI4IzQkJEKAu1pArOn7ULWO8jhzOsU2Dg0fWRkcBj0WxaVM1hYs6XIl0ftI+hzrwI/tDGf3laUSoIVSwK4/db4CG5D/H4/qqqqoNFokJ6ejkAgIIgxy7LYtm0bXC4X8vPzMWTIEAwaNCgpolsxSSu64Ta8xGY3ubm5bSo34RP+hYWFAELFlqKokFRCWwlwwNc2Gp9W1UfCFuaXiHSkn8XOn9txEhKC/l8SHD8S3+hTpSIYNoHF5TksznZxi155kGDQH2Q4+HV8d95VKoIP3vFj3C2xMx6PBCEENpsNTqcz4rXHMAyOHj2KI0eOwOFw4OTJk8jIyMArr7wS9/XGkqQV3draWvj9fhgMBrAsC4vFArfbHdbsprXvww/LY1lWeM9oiK34PaxWK3w+H/R6A34MpmKjWYZ9Zhq799EItqMc7Vobi++WJy5yoCiCoaM5BOYw+GF4wpaRMCiOYOgjwPeb4+viplIRrFrmx9jR8Rdcn8+H6upqpKSkRJzQzedus7KykJ+fn3TRrZikEV2WZcEwv8xqcrlccDqdUCgUqKmpQXZ2NjIzM9tdW8hxHM6cOYP8/HwA0RVbfieXX296enqj9f5opfDxT3J8/JMMlfbWva+OEHBLgZqajnGbX3Ili/S5LCpu7jqbbtf+RYbvVsXX3E+tJli93I8xo+IruIQQWK1WuN3ukFp3MQzDwGQywefzJV3uNhJJKbqEEFgsFthsNuTk5CA7OztqVnOBQACXLl0Cx3FQq9XQaDTQaDRQqVRtFnS+zdhqtSItLa3F5jzHLPUC/MkJGX5qgQCXHWWx+/OOF0EUFHIomsPg4AQOHm3yiu/IV2jsXBLflIJGQ7DmPT9G3xBfweU7OdPS0pCdnR322qirq4PJZIJOp0N+fn7MfHk7GkkjunzVgMvlEsbvMAyD4uLiqBy74SYZUL8p4PV6hR1ZiqIEIW5psbjP54PZbIZcLoder29zh81Ryy8R8M+OxifvAD+L46/TIB3YsjI9g8PAGUGcmElgzem462wLI9+hsfMf8RVcrYbgoxV+3DgyfoLLW6l6PB7k5eWFdIny8NGt3+9HcXFxp/RPaA9JI7o+nw/nzp0DTdMwGo2QyWSCp25bae0mmbhY3Ov1Cl07vAhrNBqhj51hGGFTz2AwRLWH/IiZwicn6gX4pIOGDAQ9PiU4+XPniCSUSg4Dx3lg/S8K5/t0fp/961ZT2PVsdD1xm0Or4fDhUidG3xi99Fdz8NFteno6dDqdFN1GIGlEl2EYOJ1OoRmCEILTp0+jV69erT4Wx3FR2yRjGEYQYZ/Ph2AwKLgqZWZmIisrK6b944fNFLbvluH1pxUwmTtX9EhRBIOuDyB4TxA/lnW+HnsAuOYzCnsfVcT1DiMlheBfi60oHVQHv98PQkjY5oRowXEcLBYLfD4f8vLywpZgdvXoVkzSiK7YU5fn5MmT6N27d6uOI+4ko2k6aruovM8vn7dVq9WCGLMsC6VSGZKaiPbubTAIbNoqw7KVcuz4pmOnGcJRciWLtHsYfD+GdJpNt+FfUTj4gCIu9dA8aakEn3zoR9mIX1IKhJCQ5gSfzweO4xoJcVvOOY/HA5PJhIyMDGRlZTWKbgkhcDqdXT66FSOJ7n9oTSdZa/F6vTCbzVCpVMjJyWnU+cavXRwR8xcFL8LRjE7OnqOw/AM53l8lh9nSOQSMJ787h+LZDA5N5OBO6bhrv3Ifhcr5ila5t7WX9DSC9av8uHp48znc9goxx3Ewm80IBALIzc2NGN1WV1cjEAigR48eYc3IuyJJI7oAhD5vnpZ46ra3k6wpgsEgzGYzWJYVfH5bCn9R8CLs8/kAQLgY2lsxUb8+4LPPZVi6Qo6vv+1c0W9GJsGg6Qx+msHCou9Y6+5/BDh3jxJud/zWlZFO8OlqP64a2vZNM/7DX3zOcRwHpVIZ8uHv8/lgMpmQlZUVtgyTv6szm83Izs5Gt27dunx0KyapRLeh01hTnrqx6CQTH5uvT9Tr9VHLX/E97O2tmAjH6bMUliwn+GCtEnZH59m8UioJBo3zwjyX4EK/xK+79wnAMkuJ2jjWQmdmEGxY7cewIdGvUhDfhXk8HjidThBCoNFooNVqQ4zOgV+i22AwiOLiYim6DUNSi244T91obpI1hBCC2tpa2O32iFFAtGlNxUQk+Iic4zjodEZs3a7B0hVyfLOz80S/FEVQeiMHZi6LH65KzCldeAHwzFDCGseUTVYmwcY1Pgy5MrY/Mz8KS6fTIT09vZHR+RdffIEDBw6gqKgIpaWluP7661tkf9paWJbFU089hTNnzkAmk2HhwoVCSz4A7NixA2+88QbkcjkmTZqEqVOnRn0N7SWpRLeh01hDT91YbZIB9e5jFosFGo2m3cMw20u4igm5XN5IiDmOg91uh9PpDBuRnzpDYflKOVaslsNq6xziCwB9B3HImsugYgwHVh6fdeeaCOiZSly+GL/b6GwdwYY1PpQOit0lzLt98RNUwlXaBINBXLp0CRcvXoTL5cKpU6dw+vRp/POf/2yxn3VL2bZtG7Zv346FCxdi7969eO+997B48WJhHbfddhvKy8uh0Wgwffp0vPXWWzER//aQ+PuxGCKTyYQx0LHaJONNdACgW7durXIsixVyuRypqamCiBJCwDAMvF4vvF4v7HY7gsEgWJaFRqOBXq8PWyfcqwfB808H8cyCIDZsrq98+Peujh/9nviBBv6kRLcCDj3mxH7TTecgUM9T4mwcBTdHR/BZuQ+D+sdOcJ1OJywWC3JycpCWlhYxd2symaDX61FSUhLz3O3NN9+MG2+8EUC9Z3ZOTo7w3KlTp1BYWIiMjAwAwLBhw1BRUYFx48bFdE2tJalFl6Zp+Hw+YZKoXC6P2knBsiysViu8Xi8MBkOzZumJhKIoKBQKKBQKqFQq+Hw+aDQaZGZmgmEYoQU5UsWEUglMnsBi8gQWJ0/XR78r13T86PfyRRqXn1ciYxHByLsYnJjJwmyI7ppT6zhkzKFx5mQcBTebYFO5DwOviI3g8jW1AFBYWBjWZzoYDKK6uhoMw6Bv375xPf/lcjkWLFiAL7/8Eq+99prwuMvlCvFuSElJgcvlitu6WkpSpRf4jTF+k8zpdMJut4OmacEjoSU5zqYQm9LodLqIY306GuIPCaPRGDaybU3FRCAAfLpZhmXvy/Hv7zqen0M4lEqCob9mYZ7D4nTryrfDovES9LpXgaP74vfzG/QEmz7yoX9JbC5bvpY8Jycn7Fgcft/CbDZDr9cjLy8vYZUJFosFU6dOxaZNm6DValFZWYl//OMfeOeddwAAL774IoYOHYpbb701IeuLRFJFunwagd8kS09PR2ZmJliWFW6ta2trEQwGoVAoQnKcLcnBulwuWCwWpKamtmvScDwRb+7pdDoYDIaIHxJ8JYQ4DyeumLDb7SEVE7eM0uDX49Q4d0GF5SsV+HCtHFZ7x/0ACgQo7CmXA+VyDLmBBTuXxZGr2yZeigBByUMKHIyj4BoNBJvLfSjpG33B5asOaJpGUVFR2OshkdEtz/r162EymXDvvfdCo9EIlUcA0KtXL5w7dw41NTXQarWoqKjAvHnz4r7G5kiqSPf06dOC+1dKSgoUCkVYYWyY4+Rd6/lbaz6iExvb8JMmDAZDpxn7zE/IiPbmXqSKCVqmwY5/Z+KDtanYtadzfJ73HfifTbdbWr7pRrMEwx+VY9/m+P2MebkcNpf70bd3dC9XPi9rs9mg1+vDWiuKo1uDwYC8vLyE3d15PB48/vjjsFqtYBgG8+fPh9frhcfjwbRp04TqBUIIJk2ahJkzZyZknU2RVKLrdDpRV1cHj8cDj8cDjuMEAeXFNJJgRrq15qsh+NutzpBK4M10GIZpdVNGe95TXDFx4iSFjzfm4NPNGXDUdPz0Q143Dr3uYXFoIgtXahN/Y0JQ9owcuz+Kp+Cy2FzuQ9/e0T33+MhVLpfDYDBEjG6rqqrAcRyKi4s79N5FZyGpRLchgUAALpcLLpcLHo8HXq8XMpksJEcZrtWR4zghb8sLrc/nQyAQgFwuD/nejhT1ivPNkXac47kWhmFQU+vFJxtl+GCtFhUHO36hfHo6weDpDH6eycIUZtNt5P/KsHNZ/AS3Wx6DFW9fglFfvyEk3uhsbrx5JMQpJ36SSrjX1NTUwGKxJDy6TTaSWnQbwo+HdrlccLvd8Hg88Pv9gtmMWq3Gt99+i5SUFJSWloY1E+fTEnxUxzBMozbJRNTout1umM1mpKamRs20PdpUngDe/ReN1euUcNR0vPWJUSgIht3BwjoBhtAGAAAT9klEQVSHxcm+9Y+NfIvGzv+L34ds93wOWz72o0dR/SUqnqrr9XqFtnexV4JarW62EaaqqgpKpRJ6vb7Z6LZHjx5RtR2V6GKiGw6WZeF2u3Ho0CEsWrQIer0e06ZNE5oq+LREpPpbQgiCwWBIfpi30ovGVInmENcJG43GDhV5R8LnA9ZtoLBkOVBxsONf0KW/YmG8xoat/2uI23sWdeewZZ0fRYVNX56RWsMbRsQAUFNTA4fDAaPRGLY9VxzdGo1G5ObmStFtDOjyosuzevVqlJaWoqSkRJhAIU5LiHf2+ZM5XP0iUH/yivOb/IUgrpZQKBTtOqE5joPNZoPL5YLBYOg0Pe7iW9vs7GxcqsrA8g8UWPWRHI4OMrutIYNLTuHYgQ9QOrQ3LljGwGyLbYdTj6J6we1e0LZLk7+jE5+DDMNAoVAgKysLWq22kUcHP+WaEILi4mIpuo0hkui2AD6aEKclfD4fFApFozlpkW7r+R1/PiIWt+Y2J+JieH9Sq9UqmKB3lmiEd6dSqVSNbm19PuDjjTIsWyHH7jiWYTWHIacWnHMJrFYPAECtlmP41WXYf+x6+P3Rv6voWVwvuAX57b8s+Rx/bW0t9Ho9aJoOCQQA4KOPPkL37t2Rm5uL0tJS5Ofnx+R8CgaDeOKJJ3Dp0iUEAgHcd999GD16tPD88uXLUV5eDp1OBwD4y1/+gp49e0Z9HR0BSXTbCMdx8Hg8IUIcDAZDcmt8RBtJiPm0hNjMnPcz5YVY/L28aPH5uJaIdEdA3JgRaSqsmOOVFJatlGPVR3LU1CbuA0UmY9G321L8eLyq0XP5+RnIK74ZFUcGRu39evesF9xuee2/JP1+P6qrq4VywXDnoNfrxY4dO/DTTz/BZDLh/PnzuPXWW3Hfffe1+/0bsm7dOlRWVuLJJ5+Ew+HAxIkT8fXXXwvPP/LII5gzZw4GDoze77OjIoluFOFbasVpCUJIiAjzHXHhEPuZ8mIM1O9Y894RRqOx05TtiKPytnTveb3Aug0yLF8px5798Y9+h/f/BBV7jzT5mtIhRajzj8PpC8Z2vVff3hw2r/Mhr32HASEEdrsddXV1yM3Njdh56HA4YLVakZubC6PRKPxdCCExiXTdbjcIIUhNTYXD4cDkyZOxfft24flx48ahT58+sFgsuPHGG3HvvfdGfQ0dBUl0Ywg/sr1hWoIvWxOnJcLtIvMXh91uFy4ev9/f6Pvb09YcK/gcoUKhiEpUfuzHes+HVeXxiX4H96vAkYpNLXqtXE7j6muH4YefR6HO1fpcaEkfDpvW+ZDbzn06v9+PqqoqpKSkRKxgCQQCqKqqAkVRKC4ujroLWHO4XC7cd999mDp1Ku644w7h8ddffx0zZsxAamoqHnjgAUyfPh2jRo2K69rihSS6cYbjOHi93hAhDgQCIWNStFotjh49CrVaDb1ej+zs7BBRDmfd2Ja25lj9fPwGXyyicq8XKP+0Pve77/vY/IwFuVWoqXoPLleg+ReLyM7Rom//Udh9cCiAlpXEXdGPw6ZyH4zt2JsjhAi/80jpm6ai23hRVVWF+++/HzNmzMDkyZND1iY2q/nggw9QU1OD+++/P67rixeS6HYAGIaB2+2Gy+XChQsX8O6778LlcuGhhx4SzGma66bjTaXFbc0NjWpiXbvLG13Ha4Pvh+MUlq+QY/U6OWrrovNeGrUPBu3bOHeups3HKLkiD5T6Vvx4srDJ1w24gsOmj3zQ5zT5sibx+Xyorq4W6rPD/c4THd0CgNVqxd13341nnnkGZWVlIc85nU7cfvvt2Lx5M7RaLf74xz9i0qRJuOGGG+K+znggiW4H4+mnn8bo0aNx7bXXNsoP0zTdKD8cKaIN19YcrdE+DQkGgzCZTKAoKiHeFB7PL9Hv/gPti36vGVSOPd8da/eaKAq4umwgTl68BVZHYz+DQQM4fLbWh5zsth2fECKMhMrLywvb6i2ObvPy8po0O4o1zz//PLZs2RJSkTBlyhR4vV5MmzYN69evx4oVK6BUKlFWVoYHH3wwIeuMB5LodhL42ks+Ig7XTddcRMsfg9+oa9jW3NQmXzjEmzYdpVb4yLFfot86Z+sE5rqhe7Drm61RXU9amhJXDrseew5fA4ap/91eOYjDxjU+ZOvadkyv14vq6mqkp6dDp9OFFVK+eiGR0a1EeCTR7cTw3XT8F9+WLC4744U0khA3dFtr2Nas0WjCfi/fdpyWlha2XTrRuN2/RL8VB5uPfkt6ncfp4+8jEGBjsp7iHjpkGsaCkvfCxrV+ZGW2/hj8wFO+9C5SdGu322Gz2RIe3UqERxLdJIPvphMLMYAQb4jmytaaamtWKBSoqakBIQRGo7FDjCdqjsNH/xP9fiyHM0z0m5XhhoZbgsuX62K6jh490rBh453INaY2qsFuDo/HA5PJhIyMjIj5cr56QSaTobi4OC7uchKtRxLdJIfvphNHxOJuOr5aoqm0BCFEMDF3u92QyWSNBl22t605lvC7+yazC7v2FmLFag0OHK6PfimKQ2nvlTh44ExM16DTqfH551PQrZuq2akcYjiOg8Vigc/nQ15eXtgPOXF0261bN+j1+g77t5CQRLdLwnfT8SIsLltrKKQ0TcPj8cBsNkOr1QrdTdFqa441fISYlpYWsrt/6AcKy1bIcebnndjxxb9jugaFgsaGDVPwq1+FVjM0NKvx+XzCZqlGoxE+LHQ6HTIzM6XoNklIWtFlWRZPPfUUzpw5A5lMhoULF6Kw8JeTnneYl8vlmDRpEqZOnZrA1SYevptOnJaoq6vD2rVrIZPJ8PDDDyM9Pb3JsrWGtpfNtTXHEpZlYTabEQwGkZubGzENYrV68Mor+7BkyUF4vUxM1vLmm2Mxe/bgFr2WHy1ltVqFiRwN7yr4DzObzQa73S5Ft52MpBXdbdu2Yfv27Vi4cCH27t2L9957D4sXLwZQn/e87bbbUF5eDo1Gg+nTp+Ott96CXh9b96jORGVlJR599FHMmDEDZWVlgpiKu+Ga8w+O1NYs/v5Y2F7ywxWzs7NbPO2jqsqFv/99D5YtOxLVzbQ//vEqvPjijS1+vdvthslkCmmbFjfD1NbW4sknn0RKSgr69euHkSNHYsSIEcLY8WjTnFGNFLy0nqQVXaA+epPL5fjkk09w4MABPPfccwDqBeXll1/G0qVLAdRPDR0yZAjGjRuXyOV2KPx+P1iWDekoa8oEXhyJKZXKZsvWxEbcMpmsUSTXFiHmx8/ws+zakt64cKEOf/vbbqxceRQMw7X6+8XcfntvrFo1ATTd/M/CR+YMwyA3NzfsHQWfbrBarZDL5bh8+TKOHj2KgoIC3HXXXe1aaySaMqqRgpe2kfikWwyRy+VYsGABvvzyS7z22mvC4+KWQwBISUmBy+VKxBI7LOFygzRNQ6vVhggxX7bG1w7bbLaQbjh+o44XkXDHaBjJiduaI41UEiMeUxTJoLuldO+ejjfeGIuHHx6BF174DmvX/giOa31cMniwAcuWjW+R4PKdfE1F5j6fD1VVVVAoFBgwYABUKhUGDx4c8/Hit956K8aOHSv8W/x3OHXqFAoLC4Uoe9iwYaioqJCCl2ZIatEFgJdeegmPPPIIpk6dik2bNkGr1SI1NRVut1t4jdvtDjsFVaJ5ZDIZ0tPTkZ6eLjwmNoGvq6sTRntH6qaTy+VITU0VZnWJ25r5sff8lOeGTSB8G6xWq0VxcXHUcsY9e2Zh6dLxePTRa/D88zuxfv0JtPSeMDc3BeXldyIlpelyOpZlYTKZwHEcCgsLw0bmfHRrt9uRn5+PnJycuOZu+Q8wl8uFBx98EA899JDwnBS8tI2kFd3169fDZDLh3nvvhUajAUVRwkXeq1cvnDt3DjU1NdBqtaioqMC8efMSvOLkgZ9QkJWVBaCxCbzVam3SBJ6iKCiVSiiVSkHMxW3NDocDfr8fDMOAEIKsrKyYTWouKcnGypW/weHDJjz33C5s2XKqyddrNHKsXTsR+flNf4g7nU5YLBZhynQ4xNFt//79E1YTLTaqETuDScFL20janK7H48Hjjz8Oq9UKhmEwf/58eL1eeDweTJs2TdgAIIRg0qRJmDlzZqKX3KXgd+mdTic8Hk9YE3itVhu2m47fbEpPT4dGoxFSE+1ta24J+/dfxl//uhM7dpxr9BxFAe+//2vceWe/iN/PMAxMJhOA+pl2kaJbq9UKh8ORkOhWTFNGNcFgEOPHj8fatWuh1Wpx1113YfHixTAaW24KHCv/3o5M0opuR0EaU9JyGprAezz8iJx6EfZ6vdi0aRPuuOOOiJtNbW1rbi07d17AX/+6E7t2XRQee+aZkViwoCzs68WG7nq9PmJEyEe3SqUSRUVFCe/4a86oJlrBi9/vh1wuT5glaTyRRDfGSGNK2g5vAu90OrFhwwasWbMGd955J8rKyoRItikTeOCXsjWx7WU0pzV/+eUZPPfcTvTpo8PSpePDvoZhGCGvbTQaIxrW89FtQUFBRJvGZKCqqgqrVq3CH/7wBygUCqxZswaVlZUwGo34/e9/n+jlxZykzel2FJra/QWAY8eOYcmSJV1iTElr4TffGIaBw+HAunXrkJKSEmICX1tbG9YEnu+m48eRq1QqYZdd3Almt9tDpjXzEXFL25rHjOmBMWN6gGUbl5cRQlBXVwebzQaDwSBsFDbE5/Ph8uXLUKlUCc3dxovDhw/DarVCoVBg27ZtsFgsePjhh/Hhhx+ipqam1WOdOhuS6MaYpnZ/AWD8+PEhY0q++uqrpB1T0lZSU1PxyCOPCP9OSUkJKQsTm8B7PB5YrdZG1Q5iE3iapoXHeMRtzXzZWmvammWy0JQFXzMsl8tRVFTUpaNbQgguXbqEd999F9dccw1omkbv3r0BAKWlpbj55psB1OfqMzMzEQgEkvqDR0ovxAFpTEn84WfTNTSBb+i21lRaIlJbs/gYDfPDhBDU1tbCbrc3WTPs9XpRVVUFtVqNoqKiuBu/J4JVq1bh8OHDOHDgAG655ZaQD9Lz589jw4YNyMrKgsFgwJgxYxK40tgiiW6MkcaUdAya6qYTR8NNddM119Ysl8ths9mgUqlgMBjCHof3xK2trUVBQUFEE/JkQlyh4HK5MGXKFLAsi9mzZ2PAgAEoLS3Fvn378Nxzz2H27NkhgUkyIolujJHGlHRcGprAezyeRiY9vBBHgh806nA4BNvLhkNC+bbmrhjdiuHb8l944QWMGzcO33zzDQ4cOIDFixfj7Nmz0Gq1wnWSzKVkkuhKSIgQd9PxaQl+tlw4E3h+6CM/uZmm6ZC2Zq/XizfeeAMOhwO9evXCtddei7KysogNEdHk8OHD+Pvf/44VK1aEPJ7oMsVp06Zh6dKlSE1NbZS/5eUoWQUXkDbSugySW1TLaK6bzmazCb63W7duxfnz5/HEE08gIyNDSCeI25q9Xi/mz5+PmpoaOBwOHDhwAJs2bcKrr74aU2F55513sGHDhpDNQp5jx47hpZdeSkiZ4t69e5GVlSV4byiVSnAcJ/zukllseaRIt4sguUVFD5vNhgceeAD9+/fH5MmTwTAMgsFgSFpCpVKhrq4OdXV16N69e1xG0ovZunUr+vXrh8ceewxr164NeW7cuHHo06dPwsoUvV5v2A+DroIU6XYRJLeo6KHVavH888+jV69ewmNiE3iXywWTyQStVov+/fsnJHc7duxYXLx4MexziS5T1Gg0YFm2S3SfhUMS3S6C5BYVPTQaTYjgAvUphczMTGRmtmHMbxwhhGD27NnC3/uGG27A8ePH414b3lUFFwA61txsiZhSVVWFWbNm4Te/+Y3kFtVFcblcuP322+F2u0EIwd69e6UW9DgjRbpdBKvVirlz54atF5asLpOfjRs3Cg57f/rTnzBr1iyhTFGqC48v0kZaFyFeblESEhJNI4muRMLoqHWkEhKxREovdFI6e8dOR60jlZCINdJGWieEF9zKykoEAoFEL6dNFBYWYtGiRWGf4+0up0+fjrfffjvOK5OQiC2S6HZC+Aj39ddfx5YtWwAAu3fvxrZt28AwTCKX1mLGjh0b0Spx/PjxePbZZ/Gvf/0L33//Pb766qs4r05CInZIotuJ+d3vfoePP/4Y33zzDd588004nU5ByDiusal2Z4CvI9XpdFAqlUIdqYREsiCJbieEF1SWZXHq1Cl8/fXXePLJJzFx4kQ4nU4AiNoo8njT1etIDx8+jLvvvrvR4zt27MCkSZMwbdq0Rm29Ep0LaSOtk0EIAU3T+O6777B9+3Z4vV6MGDECJSUl8Hg82L9/P1asWIGZM2cKjvydAamONPLmYjAYxMKFC0O8MUaNGiV5Y3RSJNHtZFAUhVdffRWHDh3CggULEAwG4ff7AQAqlQo33XQTDh06JES6HbnKoaCgQIjaxB1yEyZMwIQJExK1rITBby4+9thjIY9L3hjJRee8B+3CcByHoqIi/O1vf0NJSQkmTpyI1atXw+v1Cv3sJpMJffr0AdA1rPKiRaJv7SNtLkreGMmFFOl2MmiaDokCS0tL8fTTT0Oj0SAQCKC2thYMw6B79+4JXGXnoyPf2kveGMmFFOl2QsSVCRRFYcCAAQDqh/s9/fTTOH78OPbv35+o5XVKItUNi2/tlUqlcGsfT8TeGIFAABUVFRgyZEhc1yARPaRItxMSqTKhd+/eeOutt2Cz2TptyViiiOQ/m8hbe/Hm4p///GfMmzdP8MYwGo1xWYNE9JFEN4ngN82ys7MTvZSkId639pE2F2+66SbcdNNNMXtfifghpReSCGnTLPpIt/YS0UaKdCUkwiDd2kvECsnaUUJCQiKO/D8ucwiVBY0slAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# a. plot\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1,projection='3d')\n",
    "ax.plot_surface(x1_grid,x2_grid,u_grid,cmap=cm.jet)\n",
    "\n",
    "# b. add labels\n",
    "ax.set_xlabel('$x_1$')\n",
    "ax.set_ylabel('$x_2$')\n",
    "ax.set_zlabel('$utility,u$')\n",
    "\n",
    "# c. invert xaxis\n",
    "ax.invert_xaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following minimization problem:\n",
    "\n",
    "$$\n",
    "\\min_x  f(x) = \\min_x \\sin(x) + 0.05 \\cdot x^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve this problem and illustrate your results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update this code\n",
    "\n",
    "# a. define function\n",
    "def f(x):\n",
    "    return 0 # wrong line\n",
    "\n",
    "# b. solution using a loop\n",
    "import numpy as np\n",
    "N = 100\n",
    "x_vec = np.linspace(-10,10,N)\n",
    "f_vec = np.empty(N)\n",
    "\n",
    "f_best = np.inf # initial maximum\n",
    "x_best = np.nan # not-a-number\n",
    "\n",
    "for i,x in enumerate(x_vec):\n",
    "    f_now = f_vec[i] = f(x)\n",
    "    # missing lines\n",
    "\n",
    "# c. solution using scipy optmize\n",
    "from scipy import optimize\n",
    "x_guess = [0]      \n",
    "# missing line, hint: objective_function = lambda x: ?\n",
    "# missing line, hint: res = optimize.minimize(?)\n",
    "# x_best_scipy = res.x[0]\n",
    "# f_best_scipy = res.fun\n",
    "\n",
    "# d. print\n",
    "# missing lines\n",
    "\n",
    "# e. figure\n",
    "# missing lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best with loop is           -0.88366802 at x = -1.51515152\n",
      "best with scipy.optimize is -0.88786283 at x = -1.42756250\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAECCAYAAAD9z2x7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd1gU19vG8e/SBEVF7JXYa4yK2HsLGgsxFkAwb0wEjYrBSCxRY9fYY4xRjCW22DXG2EuQWIhiFyRW7IpYERBh5/2DsD9RLMDuzsI+n+vywp2ZnbmdXR5nz545R6MoioIQQohsz0LtAEIIIYxDCr4QQpgJKfhCCGEmpOALIYSZkIIvhBBmQgq+EEKYCSu1A7xOaGio2hGEECJLcnZ2TnO5yRZ8eH3otwkPD6dy5cp6TpN5ppoLTDeb5EofyZU+2THXmy6WpUlHCCHMhBR8IYQwE1LwhRDCTEjBF0IIMyEFXwghzIQUfCGEMBNS8IUQwkxku4L/9OlTgoODCQkJUTuKEEKYlGxX8JcsWYKvry+zZs1SO4oQQpiUbFfwW7ZsCcCePXvQarUqpxGmbPLkybRr144RI0bg5eVFUlLSa7dNSEigR48eJCYmGjGhEPqV7Qp+xYoVKVy4MFFRUZw5c0btOMJEXbt2jWPHjrF161YqV65M69atsbS0fO32NjY21K9fn61btxoxpRD6le0KvkajoV69egDs3r1b5TTCFF26dAkvLy9u3ryJm5sba9eu1X0yBPD29ubAgQMAzJw5k/HjxwPQqlUr/vjjD1UyC6EP2a7gA7qCv2fPHpWTiLfRaDSv/RMYGKjbLjAw8I3bpkeZMmVwc3Nj4MCBrFmzhqioKEqUKKFb7+fnx7x589i8eTPh4eEMGzYMgPLly3P69Gn9/MOFUEG2LvhBQUEkJCSonEaYon///ZeKFSvy4MEDcufOnWqdi4sLiqKwZMkSZsyYoWvqsbS0xNrampiYGDUiC5Fp2bLgFy5cmEqVKlG8eHGuX7+udhzxBoqivPaPj4+PbjsfH583bpteFy5coHz58tja2r5yURAREUFUVBQ2NjbY29unWpeQkECOHDky9o8VQmXZsuAD/PPPP0RERFCmTBm1owgTExMTg5WVFXZ2duTNm5ekpCSePXsGwN27dxk8eDBz587Fzs6O4OBg3fMePHiAo6Mj1tbWakUXIlOybcF/+WO6ECnOnz9P+fLldY8bNmxIaGgocXFxDBgwgKFDh1K2bFm+/PJL5syZo9suJCSEpk2bqhFZCL0w6Rmv9OHatWsULFgQW1tbtaMIE1GzZk1q1qype+zl5cXixYtp0KABq1ev1i13cXFJ9XjLli0MGjTIqFmF0Kdse4UP0L17d0qVKsW+ffvUjiJMWJUqVahbt+5bb7xq1aqVNBGKLC1bF/xy5coB0h9fvF2XLl3eeuOVm5ubERMJoX/ZuuC3aNECQK7whRACldvw3dzcdF+ulihRgkmTJul1/w0aNMDGxoYTJ05w//59HB0d9bp/IYTISlQr+Cnd4JYtW2awY9jZ2VG/fn2CgoIICgri448/NtixhBDC1KnWpHPu3Dni4uLo1asXPXv25MSJEwY5Tkqzzt69ew2yfyGEyCo0SkZuU9SDiIgITp48SdeuXbly5Qq9e/dm+/btWFklf+gIDQ0lZ86cGdp3fHy8rhtmaGgo3t7eVKxYkY0bN+otf2ZzmRpTzSa50kdypU92zBUbG4uzs3Oa61Rr0ildujROTk5oNBpKly6Ng4MDUVFRFC1aVLdN5cqVM7Tv8PBw3XPLli1LwYIFadiw4Su3yRvbi7lMjSlmK1IE7tx5dXnhwnD7dsb3u2HDBi5dusTgwYMzvA9TPF8gudIrO+YKDQ197TrVmnTWrVvH5MmTAbhz5w4xMTEULFhQ78exsbHhww8/VL3Yi/RLq9i/abkQ4s1Uu8Lv0qULw4YNw8PDA41Gw8SJE3XNOUIYw6JFi/jzzz+xsrKidu3aBAQE8PjxYwICAoiJiSEpKYmBAwdSv3592rVrR+3atTl//jx58+bF19dX7fhCpJtqFdbGxobp06cb5Vj37t3D19eXGzducPjwYaMcU5i2yMhIQkJCWLVqFVZWVgwYMIB9+/bxzz//0KBBAz799FPu3LmDh4cHu3fvJj4+ng4dOuDi4sKUKVPYsWNHquEZhMgKsvWNVykcHBzYvXs3ISEhXLt2Te04wgSEh4fzwQcfYG1tjUaj0V29X7x4ERcXFyB5mG17e3vu37+PlZWVbnmtWrW4ceOGmvGFyBCzKPhWVlY0adIEkLtuRbLKlStz6tQpEhMTURSFI0eOULp0acqWLcvRo0eB5O+WHj9+jIODA4mJiZw7dw5I/lKsVKlSasYXIkPMptG8RYsWbNmyhb1799KzZ0+144h3ULjw63vpZJaTkxO1atXCw8MDrVaLs7MzrVq1wsXFheHDh7Njxw7i4+MZO3as7rulBQsWcPPmTYoVK8ann36a+RBCGJnZFPzmzZsDyVf4iqKkex5UYXy3bxum21znzp11f//ss89SrXNwcGDu3LlpPm/ixIm62a7Cw8P1mkkIYzCLJh2A6tWr4+joyNWrV7l8+bLacYQQwujMpuBbWFjoZiuSdnyRXnv37pW5bEWWZzZNOgC9evWiQYMGMk2dEMIsmVXBb9++Pe3bt1c7hhBCqMJsmnREFnTxIoUmT07ulmNpmfzT3x8uXlQ7mRBZktkV/DNnzvDdd9+xefNmtaOIN9m2DerVQ8mRAw4ehGfPkn/a2UG9esnrhRDpYnYF/++//2bs2LGsWLFC7SjidS5ehJ49YfNmovz9oWxZsLJK/jlxImzenLzeQFf64eHhzJkzxyD7ftnNmzd1czVMmDCBmzdvvvU5UVFRjB492sDJRHZkVm34IP3xs4Q5c6B3b6hfH9Lq716/PnzxBfz0E8yYoffDV65c2WhD5h4+fJhLly7RokULvv3223d6TsGCBaXgm4hr166xZ88e9uzZw/Hjx/nwww91Y4TduHGDRo0akT9/fkqWLImTkxNOTk4UK1aMIkWKULNmTfLkyQPA5cuX2b9/P0eOHKFgwYJ069bNIHnNruBXqFCBokWLcuvWLc6ePUu1atXUjiRetnJlcvPNm3zxBTRsmO6Cf/nyZYYNG4aVlRWWlpZMmTKF+fPnc+rUKZ4/f86AAQPInTs3q1atYubMmbRs2ZIPPviAq1evUr58eSZMmICnpyefffYZlStXJigoiL/++ovvvvtOd4ywsDDGjRuHpaUlOXLkYNy4cWi1WgYOHEjBggW5c+cOTZo0wc/Pj8DAQOLj46lZsyZLlixh9OjRbN26lcjISB48eMCjR4/w9PRk586dXL58me+//54CBQowaNAgli5dSu/evQFITEzk5MmTzJkzhydPnjBz5kwsLS0pWbIkY8eOxdraOt0vg3i9Q4cO8e23377SxbtChQq6v0dHR3PlyhWuXLmS5hj1QUFBuiFfJk+eTGBgoG4fUvD1RKPR0Lx5c1auXMm+ffuk4Juie/fAyenN25QqlbxdOh08eJCqVasydOhQjh49yrp163jw4AHr1q0jKiqK5cuX06BBA932d+7cYeDAgTg5OTFw4EB2795N165d2bt3Lx9++CHr169/ZajkESNGMGHCBCpXrszu3buZPHky33zzDTdu3GDhwoXkzp0bT09P2rRpg4+PD5cuXaJly5YsWbJEtw9bW1sWLlxIYGAgQUFBzJs3j/Xr1/Pnn3/qhnWwtbVl2bJlKIpCQEAAbm5uFClSBH9/f1auXEn+/PmZNWsWGzduNFgBMUffffcdY8eOBcDe3p7mzZvTsmVLGjVqRMmSJXXbVapUiYsXLxIVFcXVq1eJjIwkMjKS27dvc/v2bUqUKKHb9r333uOTTz7BxcVFN0ifIZhdwYfkcXVWrlzJ3r17GTBggNpxxMsKFIDIyOQ2+9e5ejV5u3Tq0qULCxYs4IsvviB37txUr16dGjVqAMlNJf7+/oSEhOi2L1q0KE7//edTs2ZNLl++TM+ePfnpp5+Ijo7m9u3bVK1aNdUx7t69q2sScnFx0X3Er1SpEg4ODkDynd9vuuO7SpUqAOTOnZty5coBkDdvXp49e/bKtuPGjaN06dJ0796dQ4cOcffuXb766isgeaq8hg0bpvs8iddzdXVl+vTp+Pv78/XXX+te05fZ2NhQpkwZypQpQ926dd+4z2HDhqV6bKihO8zuS1v438TmQUFBJCUlqZxGvMLTExYufPM2v/ySvF067dmzB2dnZ3799VdcXV1ZtWoVp0+fBuDJkyd8/vnnqba/c+cOUVFRABw7doxy5cphZ2fH+++/z4QJE+jUqdMrxyhUqJBuZM0jR47w3nvvAXDx4kXi4uJISkri1KlTlCtXDgsLC7Ra7Sv7eNfvlmbNmoWiKPTr1w+APHnyUKRIEebOncuyZcvo06fPW4uNeDNFUQgODtY9rl+/Pjdu3GDcuHGvLfamyiyv8EuXLk2jRo0oW7YsT548yXIvWrbXv39y18sOHSCt1+bQoeSCn4HJbKpVq0ZAQAA//vgjFhYWzJ49m40bN+Lh4UFSUpKucKawsbFh3Lhx3Lp1iw8++EB3sdC6dWuGDx+u+/L0woULLF++nNGjRzN+/HjGjRuHoihYWloyceJEAKytrRk4cCD37t3D1dWVSpUqodVq+fnnn1/5lPAuTp06RWBgIHXq1MHb2xtIvrnw22+/xcfHB0VRyJUrF1OmTEn3vkUyRVHo06cPgYGBLF++nB49egDJn7ayJMVEHT16NMPPDQsL02MS/THVXIpigtm2blWUAgWUqC++UJQLFxQlISH559ChilKgQPJ6I2jQoEGay3///XclICDgnfdz7do1pWvXrvqK9Vom9zr+J6vmGjdunAIodnZ2yoYNG4yUKnPn60210yybdEQW0LYtHD6M5vnz5N44dnbJP589S76yb9tWtWjLly/n559/lu9/srnVq1czcuRINBoNq1ev5uOPP1Y7UqaZZZNOiqioKP766y/c3Nyk25opKluWu0OGkP+F3ivGduDAgVeWeXl54ezsnKpHxtuUKFGCNWvW6DOaMKDDhw/rekNNmzaNDh06qJxIP8z6Cr9Zs2Z069aNI0eOqB1FCGEiHjx4QOfOnXn27Bk+Pj74+/urHUlvzLrgp9x1u2fPHpWTCCFMxaVLl9BoNDRo0IA5c+Zkq7vxzbrgp/S4SBnLRAghnJ2dOXv2LL/99lu2a+pVteBHR0fTtGlTLqo03G3Tpk3RaDQcPHiQuLg4VTIIIUyDoii6vzs4OFCqVCkV0xiGagX/+fPnjBo1CltbW7UikD9/fmrUqEFCQgIH3zZ2ixAiWxs8eDBDhw7N1hd/qhX877//Hnd3dwoVKqRWBECadYQQyQPezZo1i2nTpvHvv/+qHcdgVOmWuWHDBhwdHWncuLFuhLi0ZHQ8ifj4+Hd+brly5dBoNISHhxts/IqM5DI2U80mudJHcqVPSq6BAwei1Wrx8PDAxsZG9ayGOl8a5cWGKyPp0aMHGo1GV2jfe+89fv75ZwoWLKjbJjQ0FGdn5wztPzw8/J3HM3/27BlPnz7F0dExQ8dKj/TkMjZTzSa50kdypU94eDgxMTHUqVMHOzs7Ll68SNGiRdWOlanz9abaqcoV/ouzTXl7ezN69OhUxd6YcuTIQY4cOVQ5thBCfcOHDwfAz8/PJIq9IZl1t8wXKYrC3bt31Y4hhDCiw4cPs3v3bvLmzcs333yjdhyDU31ohWXLlqkdgUuXLtG4cWPs7e2JiIhQO44QwkiCgoIACAgIMEqzrtpUL/imoFSpUsTExHDz5k2uXr2aLfvfCiFeNWTIELy9vc1mzgBp0gGsrKx0wyzs3r1b5TRCCGNq2bIl9vb2ascwCin4/2nVqhUAu3btUjmJEMLQrl69ytmzZ9WOYXRS8P/TunVrIHkgtbSmnBNCZB+TJ0+mWrVqLF26VO0oRiUF/z8VKlSgRIkSREVFcerUKbXjCCEMJDo6miX/zbHQoEEDdcMYmRT8/2g0mlRX+UKI7CkwMJC4uDhcXV0pV66c2nGMSnrpvOCrr77i888/p06dOmpHEUIYQEJCAj/++CMAgwYNUjmN8UnBf0H16tXVjiCEMKDff/+dW7duUbVqVVq1asW5c+fUjmRU0qQjhDAbCxcuBMDHxydbzWT1rqTgv+TgwYO0adPGLG6zFsKcaLVaihcvjqOjIz169FA7jiqk4Kdh165dbNq0Se0YQgg9srCwYOHChdy6dYv8+fOrHUcVUvBfUqdOHRwcHDh//rxqUy8KIQzHxsZG7QiqkYL/EisrK9q0aQPA9u3bVU4jhNCHo0ePsnr1ap49e6Z2FFVJwU+Dq6srANu2bVM5iRBCH6ZOnYq7uzszZsxQO4qqpOCnIaXg79u3j/j4eJXTCCEy4969e2zatAmNRoOXl5facVQlBT8NRYsW5YMPPiA2Npbg4GC14wghMmHNmjUkJCTQpk0bSpYsqXYcVcmNV68REBDA06dP+eCDD9SOIoTIhJUrVwLJ06maOyn4r2Gu/XSFyE4iIyM5cOAAdnZ2dOrUSe04qpMmHSFEtrVq1SoAOnbsaDaTnLyJFPw3OHPmDF999RWLFy9WO4oQIgOKFStGjRo18PT0VDuKSZCC/wZhYWH88MMPLFq0SO0oQogM8Pb25vjx43To0EHtKCZBCv4buLq6Ym1tzcGDB4mKilI7jhAig8xxoLS0SMF/gzx58tCiRQu0Wi1btmxRO44Q4h0pisK0adOIiIhQO4pJUa3gJyUlMWzYMNzd3enRowdXr15VK8obpXyz//vvv6ucRAjxrkJCQggICNBdsIlkqhX8ffv2Acnfovv5+TFp0iS1orxRx44dAdi5cyexsbEqpxFCvIu1a9cC0K1bNywspCEjhWpnolWrVowbNw6AmzdvUqBAAbWivFHx4sWpXbs2cXFx7N69W+04Qoi3UBSF9evXA9ClSxeV05gWjaIoipoBhgwZwq5du5g9ezaNGjXSLQ8NDSVnzpwZ2md8fDy2trb6isjq1au5ePEiXbt2pXz58hnej75z6ZOpZpNc6SO5krtTd+vWjYIFC7Jv3743XuFnx/MVGxuLs7Nz2isVE3D37l2lWbNmytOnT3XLjh49muH9hYWF6SOW3plqLkUx3WySK30kl6IMHTpUAZR+/fq9ddvseL7eVDtVa9LZtGkT8+fPB8DOzg6NRoOlpaVacYQQ2YCiKKxbtw6Q5py0qDaWTps2bRg2bBg9evQgMTGR4cOHkyNHDrXivNWjR49Ys2YNMTEx+Pv7qx1HCJGG2NhYXFxcUBSFxo0bqx3H5KhW8HPmzMkPP/yg1uHT7caNG/j4+ODg4EC/fv3Mepo0IUxVrly5WLlyJVqtVnrnpEHOyDuqUqUK1atX5+HDh+zYsUPtOEKIN5BinzY5K+ng7u4O/G8EPiGE6YiMjGTjxo1yv8wbSMFPh+7duwPJd93Km0oI07JixQo6d+6Mn5+f2lFMlhT8dChTpgx169bl6dOn/Pnnn2rHEUK8YNOmTcD/7o4Xr5KCn07SrCOE6blx4wZHjhwhZ86ctG7dWu04JksKfjp17dqV5s2by/jaQpiQzZs3A8ndve3s7FROY7pkTtt0Kl68OHv37lU7hhDiBSnNOW5ubionMW1yhS+EyNIePnzI3r17sbCw4KOPPlI7jkmTgp9Bp0+fpl+/foSEhKgdRQizdv36dapUqULjxo1NdtRdUyFNOhm0YsUK5s6dy9OnT6lbt67acYQwW9WqVePkyZPSVfodyBV+BvXu3RtIHjr5wYMHKqcRQmR0OHVzIgU/g8qWLUurVq2Ij49nxYoVascRwixdunSJCxcuqB0jy5CCnwk+Pj4AzJ8/H0XdeWSEMEvTp0+nfPnyzJo1S+0oWYIU/Ezo1KkThQoV4syZM+zZs0ftOEKYFUVRdP3vX5wtT7zeWwv++PHjgeQpt0RqNjY2DBw4EICJEyeqnEYI83LixAmuX79OsWLFqFWrltpxsoS3FvzDhw8D4OnpafAwWVG/fv3o3bs3P/30k9pRhDArKVf3HTp0kOGQ39Fbu2U2btyY7t27ExUVxbp166hUqRLly5c36dmpjClv3rwEBgaqHUMIs5NS8GWwtHf31oI/ZMgQrl27hre3N9evX2fv3r1cuHABa2tr+bLkJYqiEBsbS65cudSOkoqiKNy6dYtixYrpli1atIjExEQ6depE4cKFVUwnRPpdv36dY8eOkTNnTlq0aKF2nCzjnW68KlmyJIsXL6Z06dK6ZU+fPuX8+fMGC5bVnDp1it69e+Pk5MSaNWvUjgPA1atXmTp1Khs3buThw4c8ePAAa2trAObMmcPx48fp27cvrVq1okWLFlSoUEEmkhdZwunTp8mZMydt2rTB1tZW7ThZxjvfaftisYfkuSNr1Kih90BZlaOjIydOnODIkSOcOHFC1XNz+fJlJk2axJIlS3j+/DkABQsW5NKlS1SsWBGAzz77jOLFi7N9+3Z27tzJzp072b59O0uXLqVkyZKqZRfiXbRt25bo6Giio6PVjpKlyDcdelKiRAn69u2Loih8+eWXaLVao2dITExkypQpVKlShQULFpCUlISHhwdHjx7l9u3bumIPMGDAAP744w/u3LnDvHnzyJ8/P3/99RfVq1fnzJkzRs8uRHrZ2tpSvHhxtWNkKVLw9WjMmDEUKVKEQ4cOsXjxYqMfX6vVsmzZMuLj4/Hw8CAsLIyVK1fi7Oz82l4Mjo6O+Pr6smnTJtq3b0+VKlWoVKmSkZML8e7u3r3Ls2fP1I6RJUnB16O8efMyY8YMIPnLbmN83FQURffmt7GxYcmSJWzbto2VK1emuqJ/m/z587N582a2bt2KlZWVbt9CmJpvvvmGAgUKsGHDBrWjZDmqFPznz58TEBCAp6cnXbp0yVZ3qbq7u9OiRQuio6MZNmyYQY/14MEDPvnkE3x9fXXLnJ2dcXV1zdD+NBoNefPmBZJvtPvkk09M5gtoIQCSkpLYsmULMTExVK5cWe04WY4qBX/z5s04ODiwcuVKFixYwLhx49SIYRAajYaffvqJAgUKUL16dYMdJyQkhJo1a7Jx40Y2btzI1atX9br/tWvXsnHjRry8vAgKCtLrvoXIqEOHDhEdHU25cuWk6TEDVCn4rq6uuiEJgGzXFbBSpUpcvXqV/v37633fiqIwc+ZMGjVqRGRkJC4uLhw/fpxSpUrp9TheXl74+/vz/PlzvL29efjwoV73L0RGvHizlUajUTlN1qNRVGyojYmJoW/fvnTr1u2VScFDQ0MzPL51fHy8SfXNDQsLo0iRIuTMmTNTue7fv8+oUaN0c+r27NmTQYMGYWNjk+mMaZ2zxMREevTowenTp2nfvj1TpkzJ9HH0kcsUSK700Veujz76iMuXL/Prr7/i4uJiMrn0LTO5YmNjcXZ2TnulopKbN28qH3/8sbJ27do01x89ejTD+w4LC8vwc/Vtx44dip2dndKwYUMlNDQ0U/sKCAhQACVv3rzKxo0b9ZQw2evO2b///qvkzJlTAZRVq1bp9ZjvwpReyxdJrvTRR65z584pgJIvXz7l+fPnekiVPc/Xm2qnKk069+7do1evXgQEBNClSxc1IhjN+++/j6OjIwcOHMDd3Z2IiIh0PT/lximAUaNG4e3tzYkTJ3Bzc9N31DSVL19e1/OoT58+XL9+3SjHFeJl+/btA5Kv8lN6kon0UaXgz5s3j8ePHzN37ly8vb3x9vbOtsMvFy1alF27dlGxYkXOnz9P7dq1WbVq1Vu7PEZHR/PNN99QoUIFnj59CoC9vT1Lly7lvffeM0Ly//Hx8aFjx464u7ub5MdfYR769OnDuXPn+Pbbb9WOkmWp8t/kiBEjGDFihBqHVkXlypU5cuQI3bt3Z9u2bXh4eDB+/HhWr15N1apVgeR2t0ePHrFr1y7+/PNPtm7dSkxMDABbt26la9euquXXaDRs2rRJviQTqkvPvSXiVfK5yEhy587NtGnTaN++PaNGjeLy5cs4OTnp1r///vtcunQp1XNat27NhAkT9PLlVGa9WOyfP3+OpaWljEEujCYhIUEvnRPMnfzGGpFGo+HLL7/k1q1bBAcHY29vr1tna2tL7ty5adOmDT/88AMXLlxg586dJlHsX7Rz506qVavGypUr1Y4izIinpyc1a9bkn3/+UTtKliZX+CqwtrZ+ZUq2s2fPqpQmfW7evMm///7LsGHD6Ny5c4a7zgrxruLi4ti2bRuxsbEULVpU7ThZmlzhi3Tp2bMnNWvW5Pr16zL5jTCK3bt3ExsbS+3atWXo7kySgi/SxcLCgqlTpwIwY8YM3RfLQhjKxo0bAYzWFTk7k4Iv0q1FixbUr1+f6Oho5s+fr3YckY0lJibqhlOQgp95UvBFumk0Gl1f6GnTpmXbeyiE+g4ePKgbLK1KlSpqx8nypOCLDGnXrh01atRAo9Gk++5hId7V77//DiRf3ct9IJknvXREhmg0GtatW0fx4sXl7lthMGPGjKFhw4Zyda8nUvBFhpUtW1btCCKbs7e3p3PnzmrHyDakSUdk2oMHDwgMDJQpEYUwcVLwRaYoikKdOnXw9fXVjdMvhD40bdoULy8v7ty5o3aUbEMKvsgUjUbDp59+CsDs2bNVTiOyi4iICPbv38+WLVvIly+f2nGyDSn4ItN8fHywsbHhjz/+4OLFi2rHEdnAunXrAOjUqZMMmqZHUvBFphUqVAgPDw8UReGnn35SO47IBtauXQug6rDg2ZEUfKEXfn5+ACxcuFCGWxCZcv78eU6ePEmePHlo3bq12nGyFSn4Qi9q1apFo0aNePz4Mb/++qvacUQWlnJ136lTJ3LkyKFymuxF+uELvRk4cCC2traUL19e7SgiC0sp+Nl9vms1SMEXetOlSxf5JRWZtnz5ctavX0+bNm3UjpLtSMEXQpiUqlWr6uZ6FvolbfhC7w4fPsznn3/OgwcP1I6S5dy7d49NmzYxePBgmjRpgqurK48ePVI7lsgmpOALvRs1ahSLFi1i+fLlakfJMh49eoS/vz9FihTh448/Zvr06QQHB7Nv3z5y56M5icYAABmhSURBVM6t227atGns2bMnWw5jERYWRs2aNWWOBQNSteCfPHkSb29vNSMIA/Dx8QGQ8XXeUVhYGJUqVWLWrFkoikKzZs0YOXIkW7duZevWrVhYJP+aXrlyhaFDh9KqVStatWrF7du3VU6uX8uWLePEiRMcO3ZM7SjZlmpt+AsWLGDz5s3Y2dmpFUEYSMeOHSlUqBBnzpzh8OHD1K9fX+1IJq18+fLky5eP0qVL89NPP1GzZs00t8uXLx9jx45l5syZ7N27l1q1arFmzRoaNWpk5MT6p9VqWbFiBQBeXl4qp8m+VLvCL1WqFD/++KNahxcGZGNjw2effQYgH89fIzExkYSEBACsra3Zs2cPf//992uLPUDevHkZPnw4p0+fpmnTpty6dYtmzZpli9+j4OBgrl27hpOTEw0bNlQ7TralUVT8zH39+nUGDRrEmjVrXlkXGhpKzpw5M7Tf+Ph4k5yUw1Rzgf6zRUZG0rZtW+zs7AgKCsLe3t4kculLZnIlJCQQEBCgmxDeyir9H7QTExOZNWsWixYtAmD06NF069Yty56vkSNHsn79enx8fPjqq69MJpdaMpMrNjYWZ2fntFcqKrp27ZrStWvXNNcdPXo0w/sNCwvL8HMNyVRzKYphsjVt2lQBlAULFmR4H6Z6zjKaS6vVKr169VIAJW/evMq5c+cylWPhwoVKhQoVlKtXr2Yql6G9KVdcXJySJ08eBTB6/qx4vt7mTbVT+uELg/Hz88PFxYWmTZuqHcVkLFiwgEWLFmFnZ8eePXuoWLFipvbXq1cvPD09TfIq9V1t2bKFx48f4+zsTOXKldWOk61JwRcG07lzZ5me7gX//PMPAwYMAJJ7ML32Y3c6pRR7RVFYuHAhrVq1wtXVVS/7NoZ27dqxatUq6cBhBKoW/BIlSqTZfi9EdhMVFcUnn3xCQkIC/fv3N0hPlA0bNjB9+nR++eUXQkJCMv3pwVhy5sxJ9+7d1Y5hFuTGK2FQiYmJzJ07l3bt2pGYmKh2HNV89913XL9+nQYNGjB9+nSDHOPjjz+mVatWPHr0iE6dOvHw4UODHEefkpKS1I5gVqTgC4OytLTkhx9+YNu2bezYsUPtOKqZOnUqQ4cOZdWqVQabwcnCwoJJkybx/vvvExERgaenJ1qt1iDH0ofExESqVq3K559/zuPHj9WOYxak4AuD0mg09OrVC0ieHMVc5cqVi0mTJlGyZEmDH+f3338nf/78bNu2jcmTJxv0eJmxbds2IiIiMtVtV6SPFHxhcD179sTCwoI//viDqKgoteMY1YoVK4x+9Vq6dGndOEYjR44kODjYqMd/Vyk35fn6+uqGjxCGJWdZGFzRokVxdXUlMTGR3377Te04RrN37168vLxwdnbW3VVrLK6urgwZMgRnZ2eKFy9u1GO/i8jISLZu3YqNjQ3/93//p3YcsyEFXxhFyi/1kiVLVM1hLImJifTv3x+ATz/91GDt9m8ybtw4/v77b8qUKWP0Y7/NggULUBSFLl26ULBgQbXjmA0p+MIoOnToQL58+Th+/DgnT55UO47BLViwgPDwcMqUKUNAQIAqGaytrXX/0SiKwqlTp1TJ8bL4+Hh++eUXAPr06aNyGvMiN14Jo7C1tWXs2LHkyZOHcuXKqR3HoB49esSoUaMAmDJliuoTcT9//pxPPvmEHTt2cOTIEapXr65qnr1793Lnzh1q1qyZLUb6zEqk4AujSWniyO4mTZrEvXv3aNSokUncaWxtbU2xYsVISEjAw8ODI0eOZHhgQn1o164dZ86c4cmTJ2g0GtVymCNp0hFCj27dusXMmTMBmDFjhskUtBkzZlCpUiXCwsL4+uuv1Y5D1apVqVevntoxzI4UfGFU169fp0+fPrrx8rObIkWKsGnTJkaOHImLi4vacXRy5szJb7/9ho2NDfPmzWPjxo1Gz6DVajl48KDMgqYiKfjCqKysrPjll19Yvnw5d+/eVTuO3mk0Gtq2bcvYsWPVjvKKGjVqMGXKFCB5lM3IyEijHn/Tpk00bNgQDw8Pox5X/I8UfGFURYoU0Y2rk90mOb9x44baEd7Kz8+PDh068PDhQ6PehavVapk4cSIAjRs3NtpxRWpS8IXRvTjUQnb5eH/06FGcnJz48ssv1Y7yRhqNhsWLF/Pdd9/xww8/GO24f/zxB6GhoRQpUkT3+gvjk4IvjO6jjz6iUKFChIWFceTIEbXj6MWIESNISkoiT548akd5q/z58zN69Gij3Qz2+PFjpk2bBsD3338v496rSAq+MDpra2u8vb0BdHOyZmXBwcHs2LGD3Llzq3aTVUY9ePCA7t27ExERYbBjjBkzhujoaBo0aGCQeQDEu5OCL1SR0kvnt99+Iy4uTuU0mZNyk9WgQYPInz+/ymnSZ8yYMaxZs4ZOnTrx6NEjve8/PDyc2bNno9Fo+PHHH2WQNJXJ2ReqqFq1KlOnTiUoKChLf8QPCgrir7/+wsHBAX9/f7XjpNv48eOpVq0aEREReHl56X38fCcnJ4YPH463tze1atXS675F+knBF6oZPHgwNWrUUDtGpowZMwYAf39/8ubNq3Ka9LO3t+f333/H0dGRLVu2MGDAAL1+kZ4zZ07GjBnD0KFD9bZPkXFS8IVJyIrTH8bHx1OsWDEcHR3x8/NTO06GlSlThnXr1pEjRw7mzp1LQEBApot+cHAwN2/e1FNCoS9S8IWqgoODqV+/PiNGjFA7SrrZ2tqyfPlyLl26hIODg9pxMqV58+asX78ea2trpk+fnqk7cS9cuEDHjh1xdnY2+s1d4s2k4AtVWVpacvjwYRYvXmz0SUL0JSs25aTlo48+4rfffqNfv364ubllaB+3bt3Czc2Nhw8fUrduXYNP6SjSRwq+UFX9+vWpVq0ad+/eZdOmTWrHeWdTpkxh+/bt2ebGsRSffPIJc+bM0fWmuXTpEtevX3+n5x4+fBhnZ2fOnj1LpUqVWLp0qfTKMTGqvRparZZRo0bRvXt3vL295aOfmdJoNPj6+gL/m+PU1B06dIglS5bQrVs3g3RlNBWxsbF8/PHH1KpVi8DAwNd2n01MTGT+/Pk0bdqUW7du0aRJE4KCgrLETWjmRrWCv3v3bhISEli9ejVff/21Ucf1EKbFy8sLOzs79u7dy/nz59WO81YpPXP8/PyyfNv9m8THx1O4cGGioqLw9fWlZMmSfPvtt6xfv57du3cTHx8PwJMnT/Dz8yMhIYF+/fqxe/duChUqpHJ6kRbVCn5oaKhuEKUaNWpw5swZtaIIlTk4OODu7g5AYGCgymneLCQkhB07dpAzZ84s2e8+PRwdHdm2bRsrVqygdu3aREdHM3HiRLp06ULr1q11n27y5ctH//79+fXXX5kzZw7W1tYqJxevo1FUaoT89ttvadOmDU2bNgWgWbNm7N69Gyur5Em4QkNDMzwrT3x8PLa2tnrLqi+mmgvUz3by5Ek8PDwoVaoU27Zt000conaul/n6+hIcHMxnn31mksMoGOp8KYrC8ePH2bx5M9HR0cTExPDzzz+/87FM7XVMkR1zxcbG4uzsnOY61aY4tLe35+nTp7rHWq1WV+xTVK5cOUP7Dg8Pz/BzDclUc4H62SpVqkRiYiJubm7Y29ubTK4X/fPPPwQHB5MrVy6++OILk8n1IkOerypVqtCjR48MPdeUXscXZcdcoaGhr12nWpNOrVq12L9/PwAnTpygQoUKakURJkCj0eDl5ZWq2JualKkL+/fvT758+VROI0T6qXaF37p1aw4cOIC7uzuKougmRxAiNjaW+/fvU6JECbWjpDJv3jwqV65M3759uXfvntpxhEg31Qq+hYWFSU4DJ9S1Z88eunXrRuPGjU2uX37evHl1I2NKwRdZkdwVIUxKtWrViImJYfPmzVy8eFHtOABERkZm+SGchQAp+MLEFC5cGE9PTxRF4ccff1Q7Doqi0KNHD8qWLZttZucS5ksKvjA5AwcOBJLnvH38+LGqWXbs2MGBAwdISEigYsWKqmYRIrOk4AuTU6NGDVq0aEFMTAwrVqxQLYdWq9WN4jl06FAZKkBkeVLwhUkaOXIkAL/++qtqV/m//fYboaGhFC1alC+//FKVDELokxR8YZKaNWtGkyZNKFmyJLdu3TL68WNjY3WzNE2YMCHDd30LYUpU65YpxNts3LiR27dvq9J2PmPGDK5fv06NGjXo2bOn0Y8vhCHIFb4wWY6OjroxdYytVq1alC9fnhkzZmBpaalKBiH0TQq+MHnnz5+nT58+Rm3Lb9euHWFhYTRv3txoxxTC0KRJR5g8X19f9u3bR+7cuZk6dapBjxUXF4ednR3AK4P5CZHVyRW+MHlTp05Fo9Ewa9YswsLCDHacZ8+e4eLiQr9+/VKN5CpEdiEFX5g8Z2dnfH19SUxMpH///gabR3bcuHGcPXuWPXv2SLu9yJak4IssYcKECeTPn599+/axZs0ave//xIkTTJ48GY1Gw8KFC01yUgwhMksKvsgSHB0dmTRpEgCDBg3iyZMnett3QkICvXr1Iikpif79+9OwYUO97VsIUyIFX2QZn3/+OXXq1CEqKopDhw7pbb9+fn4cP34cJycnmZdBZGvSDUFkGRYWFixdupSHDx9St25dvexz5cqVzJ8/nxw5crB27VqTnnFLiMySgi+ylJfvuk1MTMxU98mOHTvSuXNn3NzccHFxyWw8IUyaNOmILGv9+vVUq1aNq1evpl5x8SL4+0PhwmBpmfzT3z95+Uvs7e1Zt24d3t7eRkothHqk4IssSVEU5syZQ0REBPXq1ePYsWPJK7Ztg3r1wM4ODh6EZ8+Sf9rZJS/fto2dO3fSoUMHnj17BqDa8A1CGJsUfJElaTQaNmzYQLNmzbh16xaNGzdmT2Ag0R/1pP69zWgmTURTriwaays05cpSZNFElN9/J65bN/q3bcuWLVtYuHCh2v8MIYxKCr7IsvLly8eOHTvo2bMnsbGxnPb1Zb7Sm8PUf2XbO3eg/qBBzIqJoY9Wy9ChQ+nTp48KqYVQj3xpK7I0GxsblixZQrly5fAYNYoGfP7abUNCQnieNy/BikLO//r0C2FO5ApfZHkajYaRI0dSyMKCSJxeu93s2bPZf+UKOWWcHGGmVC34u3bt4uuvv1YzgshGNAUK4ETka9cPGDCAXNHRUKCAEVMJYTpUK/jjx49n+vTpaLVatSKI7MbTk895yxexv/wCnp7GySOEiVGt4NeqVYvRo0erdXiRHfXvj69mAfV4ddiFwoWBQ4eSC36/fsbPJoQJ0CiGGmv2P2vXruXXX39NtWzixIlUr16dkJAQVq1axcyZM195XmhoaIYnjo6PjzfJ0Q5NNReYbrb05sq1fz/Fhg3j4Sef8LBLF54XLYr1rVs4rFuHw/r13Jw0iadNmhg9l7FIrvTJjrliY2NxdnZOe6WiosOHDytfffVVmuuOHj2a4f2GhYVl+LmGZKq5FMV0s2Uo14ULiuLvryiFCyuKpWXyT3//5OVq5jICyZU+2THXm2qndMsU2U/ZsjBjRvIfIYSOdMsUQggzoeoVft26dfU2zK0QQog3kyt8IYQwEwbvpZNRoaGhakcQQogs6XW9dEy24AshhNAvadIRQggzIQVfCCHMRJbvh79r1y62b9/O9OnTAThx4gQTJkzA0tKSRo0a0b9//1Tb379/n8GDBxMfH0+hQoWYNGkSdnZ2BskWGBhIcHAwAI8fP+bevXscOHAg1TZ9+vTh4cOHWFtbkyNHDn755ReDZHmRoig0adKE9957D4AaNWq8MojdnDlz+Ouvv7CysmL48OFUr17d4LmePHlCQEAAMTExPH/+nKFDh1KzZs1U24wfP55jx46RK1cuAObOnUvu3LkNkker1TJ69GgiIiKwsbFh/PjxODn9bzTONWvWsGrVKqysrOjbty/Nmzc3SI6XPX/+nOHDh3Pjxg0SEhLo27cvLVu21K1fvHgx69atw9HREYAxY8ZQpkwZo2Rzc3PTvR4lSpRg0gvDUKt1vjZs2MDGjRsBePbsGeHh4Rw4cIA8efIAxn1PpTh58iTTpk1j2bJlREZGMnToUDQaDeXLl+e7777DwuJ/1+Lx8fEEBAQQHR1Nrly5+P7773Wvbbpl+HYuEzBu3Djlww8/THW3bseOHZXIyEhFq9UqX3zxhXLmzJlXnrN+/XpFURRl/vz5yuLFi42S1cfHR9m/f/8ry9u2batotVqjZEhx5coVxdfX97Xrz5w5o3h7eytarVa5ceOG0rlzZ6Pk+uGHH3Svx8WLFxU3N7dXtnF3d1eio6ONkmfHjh3KkCFDFEVRlOPHjyt9+vTRrbt7967Svn175dmzZ8rjx491fzeGdevWKePHj1cURVHu37+vNG3aNNX6r7/+Wjl9+rRRsrwoPj5e6dSpU5rr1DxfLxo9erSyatWqVMuM+Z5SFEUJDAxU2rdvr3Tt2lVRFEXx9fVVDh8+rCiKoowcOVLZuXNnqu0XLVqkzJ49W1EURdmyZYsybty4DB87SzfpvDwAW0xMDAkJCZQqVQqNRkOjRo04dCj1QFqhoaE0btwYgCZNmnDw4EGD59y5cyd58uTRHTfFvXv3ePz4MX369MHDw4N9+/YZPAvA2bNnuXPnDt7e3vTu3ZtLly6lWh8aGkqjRo3QaDQUK1aMpKQk7t+/b/Bc//d//4e7uzsASUlJ5MiRI9V6rVZLZGQko0aNwt3dnXXr1hk0z4vvlRo1anDmzBndulOnTlGzZk1sbGzInTs3pUqV4ty5cwbNk8LV1ZWBAwfqHltaWqZaf/bsWQIDA/Hw8GD+/PlGyQRw7tw54uLi6NWrFz179uTEiRO6dWqerxSnT5/mwoULdO/eXbfM2O8pgFKlSvHjjz/qHp89e5Y6deoAadekl2vWyzUtPbJEk87rBmBr164dISEhumUxMTHY29vrHufKlYtr166lel5MTIzu41quXLl48uSJQTNWr16d+fPnMyON2/yfP3+u++V49OgRHh4eVK9enfz58+sl0+tyjRo1Ch8fH9q2bcvRo0cJCAhg/fr1uvUxMTE4ODjoHqecpwx/jHzHXCnnKyoqioCAAIYPH55qfWxsLF5eXnz22WckJSXRs2dPqlWrRqVKlfSW60Uvv58sLS1JTEzEysoq1fsIks9RTEyMQXK8LKXpISYmBj8/P7766qtU6z/66CM8PT2xt7enf//+7Nu3zyjNJ7a2tnz++ed07dqVK1eu0Lt3b7Zv3676+Uoxf/58+r00Uqqx31MAH374IdevX9c9VhQFjUYDpF2T9FmzskTB79q1K127dn3rdvb29jx9YTajp0+f6trpXt7G1tY2zfX6znjhwgXy5MmTqu03RYECBXB3d8fKyor8+fNTuXJlLl++rNeCn1auuLg43VVh7dq1uXPnTqo3XVrnUd9tmq87XxEREQwaNIhvvvlGd9WTws7Ojp49e+q+c6lXrx7nzp0z2C/ny+dBq9ViZWWV5jpDnKM3uXXrFv369cPT05MOHTroliuKwqeffqrL0rRpU8LCwoxS8EuXLo2TkxMajYbSpUvj4OBAVFQURYsWVf18PX78mEuXLlGvXr1Uy439nkrLi+31b6pZr1ufrmNl+JkmyN7eHmtra65evYqiKPz999/Url071Ta1atUiKCgIgP37979+GFE9OXjwIE1eMxzvwYMHdVdnT58+5fz580b5cm3OnDm6q+tz585RrFgxXbGH5HP0999/o9VquXnzJlqtVq9X969z4cIFBg4cyPTp02natOkr669cuYKnpydJSUk8f/6cY8eOUbVqVYPlqVWrFvv37weSOwNUqFBBt6569eqEhoby7Nkznjx5wsWLF1OtN6R79+7Rq1cvAgIC6NKlS6p1MTExtG/fnqdPn6IoCiEhIVSrVs0oudatW8fkyZMBuHPnDjExMRQsWBBQ93wBHDlyhAYNGryy3NjvqbRUqVJF11Kxf/9+g9asLHGFnx5jxoxh8ODBJCUl0ahRIz744AMePnzIiBEjmDNnDn379mXIkCGsWbOGfPny6Xr3GMrly5dp2LBhqmVTpkzB1dWVpk2b8vfff9OtWzcsLCwYNGiQUQqrj48PAQEBBAUFYWlpqetJkZKrevXq1K5dm+7du6PVahk1apTBMwFMnz6dhIQEJkyYACT/B/7zzz+zePFiSpUqRcuWLenQoQPdunXD2tqaTp06Ub58eYPlad26NQcOHMDd3R1FUZg4cWKqLN7e3nh6eqIoCv7+/q9852Ao8+bN4/Hjx8ydO5e5c+cCyZ+Y4uLi6N69O/7+/vTs2RMbGxvq16+f5n+ehtClSxeGDRuGh4cHGo2GiRMnsmzZMtXPFyT/HpYoUUL3WK33VFqGDBnCyJEjmTFjBmXKlOHDDz8EoFevXsybNw8PDw+GDBmCh4cH1tbWmapZcqetEEKYiWzVpCOEEOL1pOALIYSZkIIvhBBmQgq+EEKYCSn4QghhJqTgCyGEmZCCL4QQZkIKvhDvyNvbWze89cyZMxk/frzKiYRIn2x3p60QhuLn58fs2bOJjo4mPDycn3/+We1IQqSL3GkrRDp4eXkRGxvL0qVLU42kKURWIE06QryjiIgIoqKisLGxkWIvsiQp+EK8g7t37zJ48GDmzp2LnZ2dbupKIbISKfhCvEVcXBwDBgxg6NChlC1bli+//JI5c+aoHUuIdJM2fCGEMBNyhS+EEGZCCr4QQpgJKfhCCGEmpOALIYSZkIIvhBBmQgq+EEKYCSn4QghhJqTgCyGEmfh/O/NOzJgx2+UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# a. define function\n",
    "def f(x):\n",
    "    return np.sin(x)+0.05*x**2\n",
    "\n",
    "# b. solution using a loop\n",
    "import numpy as np\n",
    "N = 100\n",
    "x_vec = np.linspace(-10,10,N)\n",
    "f_vec = np.empty(N)\n",
    "\n",
    "f_best = np.inf # initial maximum\n",
    "x_best = np.nan # not-a-number\n",
    "\n",
    "for i,x in enumerate(x_vec):\n",
    "    f_now = f_vec[i] = f(x)\n",
    "    if f_now < f_best:\n",
    "        x_best = x\n",
    "        f_best = f_now\n",
    "\n",
    "# c. solution using scipy optmize\n",
    "from scipy import optimize\n",
    "x_guess = [0]      \n",
    "objective_function = lambda x: f(x[0])\n",
    "res = optimize.minimize(objective_function, x_guess, method='Nelder-Mead')\n",
    "x_best_scipy = res.x[0]\n",
    "f_best_scipy = res.fun\n",
    "\n",
    "# d. print\n",
    "print(f'best with loop is           {f_best:.8f} at x = {x_best:.8f}')\n",
    "print(f'best with scipy.optimize is {f_best_scipy:.8f} at x = {x_best_scipy:.8f}')\n",
    "\n",
    "# e. figure\n",
    "import matplotlib.pyplot as plt\n",
    "fig = plt.figure() # dpi = dots-per-inch (resolution)\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "ax.plot(x_vec,f_vec,ls='--',lw=2,color='black',label='$f(x)$')\n",
    "ax.plot(x_best,f_best,ls='',marker='s',color='blue',label='loop')\n",
    "ax.plot(x_best_scipy,f_best_scipy,ls='',marker='o',\n",
    "        markersize=10,markerfacecolor='none',\n",
    "        markeredgecolor='red',label='scipy.optimize')\n",
    "\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$f$')\n",
    "ax.grid(True)\n",
    "ax.legend(loc='upper center');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following $M$-good, $x=(x_1,x_2,\\dots,x_M)$, **utility maximization problem** with exogenous income $I$, and price-vector $p=(p_1,p_2,\\dots,p_M)$,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "V(p_{1},p_{2},\\dots,,p_{M},I) & = \\max_{x_{1},x_{2},\\dots,x_M} x_{1}^{\\alpha_1} x_{2}^{\\alpha_2} \\dots x_{M}^{\\alpha_M} \\\\\n",
    " & \\text{s.t.}\\\\\n",
    "E & = \\sum_{i=1}^{M}p_{i}x_{i} \\leq I,\\,\\,\\,p_{1},p_{2},\\dots,p_M,I>0\\\\\n",
    "x_{1},x_{2},\\dots,x_M & \\geq 0\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem:** Solve the 5-good utility maximization problem for arbitrary preference parameters, $ \\alpha = (\\alpha_1,\\alpha_2,\\dots,\\alpha_5)$, prices and income. First, with a loop, and then with a numerical optimizer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the following functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def utility_function(x,alpha):\n",
    "    # ensure you understand what this function is doing\n",
    "\n",
    "    u = 1\n",
    "    for x_now,alpha_now in zip(x,alpha):\n",
    "        u *= np.max(x_now,0)**alpha_now\n",
    "    return u\n",
    "    \n",
    "def expenditures(x,p):\n",
    "    # ensure you understand what this function is doing\n",
    "\n",
    "    E = 0\n",
    "    for x_now,p_now in zip(x,p):\n",
    "        E += p_now*x_now\n",
    "    return E\n",
    "\n",
    "def print_solution(x,alpha,I,p):\n",
    "    # you can just use this function\n",
    "    \n",
    "    # a. x values\n",
    "    text = 'x = ['\n",
    "    for x_now in x:\n",
    "        text += f'{x_now:.2f} '\n",
    "    text += f']\\n'\n",
    "    \n",
    "    # b. utility\n",
    "    u = utility_function(x,alpha)    \n",
    "    text += f'utility = {u:.3f}\\n'\n",
    "    \n",
    "    # c. expenditure vs. income\n",
    "    E =  expenditures(x,p)\n",
    "    text += f'E = {E:.2f} <= I = {I:.2f}\\n'\n",
    "    \n",
    "    # d. expenditure shares\n",
    "    e = p*x/I\n",
    "    text += 'expenditure shares = ['\n",
    "    for e_now in e:\n",
    "        text += f'{e_now:.2f} '\n",
    "    text += f']'        \n",
    "        \n",
    "    print(text)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can initially use the following parameter choices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = np.ones(5)/5\n",
    "p = np.array([1,2,3,4,5])\n",
    "I = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving with a loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update this code\n",
    "\n",
    "N = 15 # number of points in each dimension\n",
    "fac = np.linspace(0,1,N) # vector betweein 0 and 1\n",
    "x_max = I/p # maximum x so E = I\n",
    "\n",
    "# missing lines\n",
    "for x1 in fac:\n",
    "   for x2 in fac:\n",
    "        for x3 in fac:\n",
    "            for x4 in fac:\n",
    "                for x5 in fac:\n",
    "                    x = np.array([x1,x2,x3,x4,x5])*x_max\n",
    "                    E = expenditures(x,p)\n",
    "                    if E <= I:\n",
    "                        u_now = utility_function(x,alpha)\n",
    "                        # misssing lines\n",
    "\n",
    "# print_solution(x_best,alpha,I,p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Extra:** The above code can be written nicer with the ``product`` function from ``itertools``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving with a numerical optimizer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update this code\n",
    "\n",
    "from scipy import optimize\n",
    "\n",
    "# a. contraint function (negative if violated)\n",
    "# missing line, hint: constraints = ({'type': 'ineq', 'fun': lambda x: ?})\n",
    "# missing line, hint: bounds = [(?,?) for p_now in p]\n",
    "\n",
    "# b. call optimizer\n",
    "initial_guess = (I/p)/6 # some guess, should be feasible\n",
    "# missing line, hint: res = optimize.minimize(?,?,method='SLSQP',bounds=bounds,constraints=constraints)\n",
    "\n",
    "# print(res.message) # check that the solver has terminated correctly\n",
    "\n",
    "# c. print result\n",
    "# print_solution(res.x,alpha,I,p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions using loops"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using **raw loops**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = [2.14 1.07 0.71 0.36 0.43 ]\n",
      "utility = 0.758\n",
      "E = 10.00 <= I = 10.00\n",
      "expenditure shares = [0.21 0.21 0.21 0.14 0.21 ]\n"
     ]
    }
   ],
   "source": [
    "N = 15 # number of points in each dimension\n",
    "fac = np.linspace(0,1,N) # vector betweein 0 and 1\n",
    "x_max = I/p # maximum x so E = I\n",
    "\n",
    "u_best = -np.inf\n",
    "x_best = np.empty(5)\n",
    "for x1 in fac:\n",
    "   for x2 in fac:\n",
    "        for x3 in fac:\n",
    "            for x4 in fac:\n",
    "                for x5 in fac:\n",
    "                    x = np.array([x1,x2,x3,x4,x5])*x_max\n",
    "                    E = expenditures(x,p)\n",
    "                    if E <= I:\n",
    "                        u_now = utility_function(x,alpha)\n",
    "                        if u_now > u_best:\n",
    "                            x_best = x\n",
    "                            u_best = u_now\n",
    "\n",
    "print_solution(x_best,alpha,I,p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using **smart itertools loop:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = [2.14 1.07 0.71 0.36 0.43 ]\n",
      "utility = 0.758\n",
      "E = 10.00 <= I = 10.00\n",
      "expenditure shares = [0.21 0.21 0.21 0.14 0.21 ]\n"
     ]
    }
   ],
   "source": [
    "import itertools as it\n",
    "\n",
    "N = 15 # number of points in each dimension\n",
    "fac = np.linspace(0,1,N) # vector betweein 0 and 1\n",
    "x_max = I/p # maximum x so E = I\n",
    "\n",
    "x_best = np.empty(5)\n",
    "u_best = -np.inf\n",
    "for x in it.product(fac,fac,fac,fac,fac):\n",
    "    x *= x_max\n",
    "    E = expenditures(x,p)\n",
    "    if E <= I:\n",
    "        u_now = utility_function(x,alpha)\n",
    "        if u_now > u_best:\n",
    "            x_best = x\n",
    "            u_best = u_now\n",
    "          \n",
    "print_solution(x_best,alpha,I,p)       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions using solvers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution using a **constrained optimizer:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "x = [2.00 1.00 0.67 0.50 0.40 ]\n",
      "utility = 0.768\n",
      "E = 10.00 <= I = 10.00\n",
      "expenditure shares = [0.20 0.20 0.20 0.20 0.20 ]\n"
     ]
    }
   ],
   "source": [
    "# a. contraint function (negative if violated)\n",
    "constraints = ({'type': 'ineq', 'fun': lambda x:  I-expenditures(x,p)})\n",
    "bounds = [(0,I/p_now) for p_now in p]\n",
    "\n",
    "# b. call optimizer\n",
    "initial_guess = (I/p)/6 # some guess, should be feasible\n",
    "res = optimize.minimize(\n",
    "    lambda x: -utility_function(x,alpha),initial_guess,\n",
    "    method='SLSQP',bounds=bounds,constraints=constraints)\n",
    "\n",
    "print(res.message) # check that the solver has terminated correctly\n",
    "\n",
    "# c. print result\n",
    "print_solution(res.x,alpha,I,p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution using an **unconstrained optimizer:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "x = [2.00 1.00 0.67 0.50 0.40 ]\n",
      "utility = 0.768\n",
      "E = 10.00 <= I = 10.00\n",
      "expenditure shares = [0.20 0.20 0.20 0.20 0.20 ]\n"
     ]
    }
   ],
   "source": [
    "# a. define objective function\n",
    "def unconstrained_objective(x,alpha,I,p):\n",
    "    \n",
    "    penalty = 0\n",
    "    E = expenditures(x,p)\n",
    "    if E >= I:\n",
    "        ratio = I/E\n",
    "        x *= ratio # now p*x = I\n",
    "        penalty = 1000*(E-I)**2\n",
    "    \n",
    "    u = utility_function(x,alpha)\n",
    "    return -u + penalty \n",
    "    # note: \n",
    "    #  \"-u\" because we are minimizing\n",
    "    #  \"+ penalty\" because the minimizer \n",
    "    #   will then avoid the E > I\n",
    "\n",
    "# b. call optimizer\n",
    "initial_guess = (I/p)/6\n",
    "res = optimize.minimize(\n",
    "    unconstrained_objective,initial_guess,\n",
    "    method='Nelder-Mead',args=(alpha,I,p),options={'maxiter':5000},tol=1e-10)\n",
    "\n",
    "print(res.message)\n",
    "\n",
    "# c. print result\n",
    "print_solution(res.x,alpha,I,p)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cost minimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following 2-good **cost minimziation problem** with required utility $u_0$, and price-vector $p=(p_1,p_2)$,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "E(p_{1},p_{2},u_0) & = \\min_{x_{1},x_{2}} p_1 x_1+p_2 x_2\\\\\n",
    " & \\text{s.t.}\\\\\n",
    "x_{1}^{\\alpha}x_{2}^{1-\\alpha} & \\geq u_0 \\\\\n",
    "x_{1},x_{2} & \\geq 0\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem:** Solve the 2-good cost-minimization problem with arbitrary required utility, prices and income. Present your results graphically showing that the optimum is a point, where a budgetline is targent to the indifference curve through $u_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classy solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem:** Implement your solution to the utility maximization problem and/or the cost minimization problem above in a class as seen in Lecture 3. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}