Integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Analysis - IMPA 2020\n",
    "### Professor Dan Marchesin\n",
    "### Hallison Paz, 1st year Phd student"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Integration\n",
    "--------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Polygon\n",
    "from matplotlib.collections import PatchCollection\n",
    "\n",
    "from ipywidgets import interact, interact_manual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "OUTPUT_DIR = os.path.join('images', 'integration')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple analytical functions for test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_polynomial(coefs):\n",
    "    '''Returns a function that represents a polynomial with coefficients coefs'''\n",
    "    def poly(x):\n",
    "        v = 0\n",
    "        for i in range(len(coefs)):\n",
    "            v += coefs[i]*x**i\n",
    "        return v\n",
    "    return poly\n",
    "    \n",
    "def integrate_polynomial(coefs, a, b):\n",
    "    '''Integrates a polynomial with coefficients coefs in the interbal a–b'''\n",
    "    v = 0\n",
    "    for i in range(len(coefs)):\n",
    "        v += coefs[i]*(b**(i+1))/(i+1)\n",
    "        v -= coefs[i]*(a**(i+1))/(i+1)\n",
    "    return v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Checking functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyEAAAIgCAYAAACSxWbiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3SUZeL28euZNAghgdCFQKQTeg29F3HFwtplRaSI9Ca9F+lVpIuIve1a14I0qYFQQgcpkZoAAgkhIW2e94/sb19dCy2Zeyb5fs7hrEsyM99zXNdcPDP3Y9m2LQAAAABwFYfpAAAAAAA5CyMEAAAAgEsxQgAAAAC4FCMEAAAAgEt5mw4AAAAAcOd27dpV2Nvbe4WkKnLPiwtOSQfS0tK61a5d++Kvv8AIAQAAADyQt7f3iqJFi1YqVKjQVYfD4XZH3jqdTuvSpUthMTExKyQ9/OuvueNiAgAAAHBrVQoVKhTvjgNEkhwOh12oUKE4ZVyp+e3XDPQAAAAAuHcOdx0g/+c/fb/bHIwQAAAAAC7FCAEAAADgUowQAAAAAHdl4cKFBcqXLx9WoUKFsEcfffT+230cp2MBAAAAHu6VT6JCjsVc98/M5yxfNG/izMern/mzr0dGRuaaNWtWsW3bth0pVqxYWmxsrNftPjdXQgAAAADcse+++y6wQ4cOV4sVK5YmSUWKFEm/3cdyJQQAAADwcH91xcIdcSUEAAAAwB1r165d/Jdffpk/JibGS5Lu5O1YXAkBAAAAcMfq1Klzc/DgwReaNGlS0eFw2FWqVEn89NNPo2/nsZZtu/X9TQAAAAD8gaioqOjq1atfNt1xK1FRUQWrV68e+uvf4+1YAAAAAFyKEQIAAADApRghAAAAAFyKEQIAAADApRghAAAAAFyKEQIAAADApRghAAAAAFyKEQIAAADApRghAAAAAO7KokWLgqtWrVqpYsWKYc8++2yptLS023qcdxZ3AQAAAMhqn/UO0cVD/pn6nIXDEvXo62f+7Mu7d+/O9cknnwRHRkYe8fPzszt16lRyyZIlBfr06fPLrZ6aEQIAAADgjn377bd5Dxw44F+9evVKknTz5k1H4cKFb+tSCCMEAAAA8HR/ccUiq9i2bT3xxBO/vP766+fu9LF8JgQAAADAHXvggQfiv/rqq/znzp3zlqTY2FivY8eO+d7OYxkhAAAAAO5Y7dq1b44ePfpcq1atypcvXz6sZcuW5c+cOeNzO4/l7VgAAAAA7kr37t2vdu/e/eqdPo4rIQAAAABcihECAAAAwKUYIQAAAABcihECAAAAwKUYIQAAAABcihECAAAAwKUYIQAAAABcihECAAAAwKUYIQAAAADuytGjR31Lly5d+emnny5VtmzZyo0aNSqXkJBg3epx3DEdAAAA8HBjtowJOX71uH9mPmfZ/GUTJzWadOZW33f69Olc77zzzsmGDRv+/OCDD5ZevXp1/l69el35q8dwJQQAAADAXStevHhyw4YNkySpZs2aidHR0X63egxXQgAAAAAPdztXLLKKr6+v/X9/7eXlZSclJd3yQgdXQgAAAAC4FCMEAAAAgEtZtm3f+rsAAAAAuJWoqKjo6tWrXzbdcStRUVEFq1evHvrr3+NKCAAAAACXYoQAAAAAcClGCAAAAACXYoQAAAAAnsnpdDpveXdyk/7T5/zf32eEAAAAAJ7pwKVLl4LcdYg4nU7r0qVLQZIO/O/XuFkhAAAA4IHS0tK6xcTErIiJiaki97y44JR0IC0trdv/foEjegEAAAC4lDsuJgAAAADZGCMEAAAAgEsxQgAAAAC4FCMEAAAAgEsxQgAAAAC4FCMEAAAAgEsxQgAAAAC4FCMEAAAAgEsxQgAAAAC4FCMEAAAAgEsxQgAAAAC4FCMEAAAAgEsxQgAAAAC4FCMEAAAAgEsxQgAAAAC4FCMEAAAAgEt5mw64EwULFrRDQ0NNZwCAx0hISJAkBQQEGC4BAM+ya9euy7ZtFzLdkV151AgJDQ1VZGSk6QwA8BibN2+WJDVu3NhwCQB4FsuyfjbdkJ3xdiwAAAAALsUIAQAAAOBSjBAAAAAALsUIAQAAAOBSjBAAAAAALsUIAQAAAOBSjBAAAAAALsUIAQAAAOBSjBAAAAAALsUIAQAAAOBSjBAAAAAALsUIAQAAAOBSjBAAAAAALsUIAQAAAOBSxkeIZVlelmXtsSzrK9MtAAAAALKe8REiqb+kw6YjAAAAALiG0RFiWVYJSX+TtMJkx62cuX7GdAIAAECOcuXmFSWkJJjOQBYxfSVknqShkpyGO/7UJ8c+0cP/eljHrx43nQIAAJBjTIuYpse/fFyp6ammU5AFjI0Qy7IeknTRtu1dt/i+HpZlRVqWFXnp0iUX1f1/rUq2kr+Pv6btmCbbtl3++gAAADnNzpid+ib6G3Uo00E+Xj6mc5AFTF4JaSTpYcuyoiV9IKmlZVnv/O832ba9zLbtOrZt1ylUqJCrG5U/V371qdlHETERWvPzGpe/PgAAQE6S5kzT1B1TVSxPMb1Y5UXTOcgixkaIbdsjbNsuYdt2qKSnJa2zbbuTqZ6/8kT5J1Q+f3nNipylpLQk0zkAAADZ1sfHPtZPV3/SK3VfUW7v3KZzkEVMfybEI3g7vDW83nBduHFBKw+sNJ0DAACQLV29eVUL9yxUeNFwtS7Z2nQOspBbjBDbtjfYtv2Q6Y6/UrdoXbUPba+V+1fq7PWzpnMAAACynQV7FuhG6g0NrzdclmWZzkEWcosR4ikG1RkkL4eXZu6caToFAAAgWzn4y0F9euxTPVPxGZXNX9Z0DrIYI+QOFM1TVD2q9dC6M+u09dxW0zkAAADZgtN2amrEVOXPlV+9avQynQMXYITcoefDnlfJvCU1dcdUzq0GAADIBF+d/EpRl6I0oNYA5fXNazoHLsAIuUO+Xr4aVm+YouOj9e7hd03nAAAAeLSElATNiZyjagWr6ZGyj5jOgYswQu5C0xJN1bREUy2OWqxLia6/gSIAAEB2sSRqia7cvKIR4SPksPjRNKfg7/RdGlZ3mFKdqZq7a67pFAAAAI908tpJvXv4XT1W7jFVKVjFdA5ciBFyl0oGllTnyp315ckvtffiXtM5AAAAHsW2bU3bMU25vXOrX81+pnPgYoyQe9C9ancV9i+sVyNeVboz3XQOAACAx1h3ep22XdimXjV6qUDuAqZz4GKMkHvg7+OvwbUH6/CVw/rn8X+azgEAAPAIN9NuambkTJXNV1ZPVXzKdA4MYITco/b3t1ftIrW1YPcCxSXHmc4BAABwe28efFPnEs5pRL0R8nH4mM6BAYyQe2RZlkbUG6H4lHgt3LPQdA4AAIBbO59wXm/sf0NtS7VVvWL1TOfAEEZIJqgQXEFPln9SHx37SEevHDWdAwAA4LZmRc6SJUtD6gwxnQKDGCGZpE/NPgr0DdSrEa/Ktm3TOQAAAG5n+4XtWvPzGnWr2k3FAoqZzoFBjJBMEuQXpH61+mn3xd365tQ3pnMAAADcSqozVVMjpqpEQAm9UOUF0zkwjBGSiTqW7ahKwZU0O3K2ElMTTecAAAC4jfcPv6+TcSc1tO5Q+Xn5mc6BYYyQTOTl8NLI8JG6mHRRS/ctNZ0DAADgFi4nXdbiqMVqVLyRmoc0N50DN8AIyWQ1CtfQI2Ue0eqDq3X86nHTOQAAAMbN2DlDyenJGl53uCzLMp0DN8AIyQKD6gySv4+/Jm2fxIfUAQBAjrbt/DZ9c+obda3aVaFBoaZz4CYYIVkgOFewBtUepN0Xd+vzE5+bzgEAADAiOT1ZUyKmKCRviLpV7WY6B26EEZJFHiv3mKoXqq45kXN07eY10zkAAAAut/LASv0c/7NGh4/mw+j4DUZIFnFYDo2pP0bxKfGat3ue6RwAAACXOh1/Wiv2rdADoQ+oYfGGpnPgZhghWahCcAV1qtRJn/70qfZe3Gs6BwAAwCVs29aUiCny8fLRK3VfMZ0DN8QIyWK9avRSEf8imrh9olKdqaZzAAAAstx3P3+nree3qm/NvirsX9h0DtwQIySL+fv4a0S9Efrp6k967/B7pnMAAACyVEJKgmbsmKFKwZX0dIWnTefATTFCXKBlyZZqVqKZXt/7umJuxJjOAQAAyDIL9y7U5aTLGttgrLwcXqZz4KYYIS5gWZZGhI+QbduatmOa6RwAAIAsceiXQ3r/yPt6ssKTqlKwiukcuDFGiIsUDyiul6q/pLWn12rjmY2mcwAAADJVujNdk7ZNUn6//OpXq5/pHLg5RogLdQ7rrDJBZTR1x1QlpSWZzgEAAMg0Hx/7WAd+OaBX6r6iQN9A0zlwc4wQF/Lx8tHo+qN1LuGclu1bZjoHAAAgU1xOuqwFuxcovFi4Hrz/QdM58ACMEBerU7SOHi7zsFYdXKUT106YzgEAALhnsyJn6Wb6TY0OHy3LskznwAMwQgwYXGew/L39NXn7ZNm2bToHAADgrm2/sF1fn/xaXat2VWhQqOkceAhGiAHBuYI1sPZARcZG6suTX5rOAQAAuCsp6Smasn2KQvKGqFvVbqZz4EEYIYZ0LNdR1QtV16ydsxSXHGc6BwAA4I6tPLBS0fHRGhU+Sn5efqZz4EEYIYY4LIfG1B+j+JR4zds9z3QOAADAHTkdf1rL9y1Xu9B2alS8kekceBhGiEEVgivouUrP6ZNjn2jvxb2mcwAAAG6Lbdt6NeJV+Xj5aGjdoaZz4IEYIYb1qtFLRfyLaNL2SUp1pprOAQAAuKXvfv5OW85vUd+afVXYv7DpHHggRohheXzyaET4CB27ekxvHXzLdA4AAMBfikuO09SIqQorEKanKjxlOgceihHiBlqVbKU2pdpo8d7Fio6LNp0DAADwp2ZFZhyqM6HhBHk7vE3nwEMxQtzEiHoj5OftpwnbJshpO03nAAAA/M6289v02fHP9ELlF1QxuKLpHHgwRoibKORfSEPqDFFkbKQ+/elT0zkAAAC/kZSWpInbJqpUYCn1rN7TdA48HCPEjTxW9jGFFw3XnMg5uph40XQOAADAfy3au0hnE85qXINxyuWdy3QOPBwjxI1YlqWxDcYq1ZmqKdunyLZt00kAAAA6ePmgVh9arcfLP666ReuazkE2wAhxMyUDS6p3jd5ad2ad1vy8xnQOAADI4VKdqRq7dawK5CqgQbUHmc5BNsEIcUP/CPuHKgVX0qsRryouOc50DgAAyMHeOviWjl09plH1Rymvb17TOcgmGCFuyNvhrQkNJ+ha8jXNjpxtOgcAAORQp+JOafHexWpTqo1alWxlOgfZCCPETVUqUEmdK3fWv47/S9svbDedAwAAchin7dSEbRPk5+2nkeEjTecgm2GEuLGXq7+sUoGlNGHrBCWlJZnOAQAAOcgnxz7RrthdeqXOKyqYu6DpHGQzjBA3lss7l8Y1GKezCWe1aO8i0zkAACCHiL0Rq7m75iq8aLgeLfuo6RxkQ4wQN1e3aF39vdzftfrQah385aDpHAAAkM3Ztq3JEZOV5kzTuAbjZFmW6SRkQ4wQDzCoziAVyFVA47aMU6oz1XQOAADIxr7/+XttOLNBvWv0VkhgiOkcZFOMEA8Q6BuoUeGjdPTqUb118C3TOQAAIJuKS47TqxGvKqxAmDqFdTKdg2yMEeIhWpVqpTal2mjx3sWKjos2nQMAALKhWZGzFJccpwkNJ8jb4W06B9kYI8SDjKg3Qn7efhq/bbycttN0DgAAyEa2nd+mz45/phcqv6CKwRVN5yCbY4R4kEL+hTSkzhDtit2lT3/61HQOAADIJpLSkjRh2wSVCiylntV7ms5BDsAI8TCPlX1M4UXDNSdyjmJvxJrOAQAA2cDre17XuYRzGtdgnHJ55zKdgxyAEeJhLMvS2AZjleZM08TtE2XbtukkAADgwfZe3Ku3D7+tx8s/rrpF65rOQQ7BCPFAJQNLqn+t/vrx7I/67PhnpnMAAICHSkpL0pgtY1TEv4gG1x5sOgc5iLERYllWLsuydliWFWVZ1kHLsiaYavFEz1Z6VrWL1NaMnTMUcyPGdA4AAPBAr+15TdHx0ZrYaKICfANM5yAHMXklJFlSS9u2q0uqIekBy7LqG+zxKA7LoUmNJindTte4reN4WxYAALgju2J36Z1D7+ipCk+pfjF+BINrGRshdoaE//xXn//84ifpOxCSN0SDaw/W1vNb9clPn5jOAQAAHiIxNVGjN49W8YDiGlR7kOkc5EBGPxNiWZaXZVl7JV2UtMa27QiTPZ7oiQpPKLxYuGbtnKVzCedM5wAAAA8wd9dcnUs4p0mNJsnfx990DnIgoyPEtu1027ZrSCohqZ5lWVX+93ssy+phWVakZVmRly5dcn2km3NYDk1sODHj1KwtY7mJIQAA+EsRFyL0wdEP9Fyl51SnaB3TOcih3OJ0LNu2r0laL+mBP/jaMtu269i2XadQoUKuj/MA9wXcp1fqvKIdMTv0wZEPTOcAAAA3lZCSoLFbxqpUYCn1q9XPdA5yMJOnYxWyLCvff/46t6Q2ko6Y6vF0Hct1VKPijTRv9zydjj9tOgcAALih2btmKyYxRpMbTVZu79ymc5CDmbwSUkzSesuy9knaqYzPhHxlsMejWZal8Q3Gy9vy1pgtY3hbFgAA+I2t57bqk2OfqHNYZ9UoXMN0DnI4k6dj7bNtu6Zt29Vs265i2/ZEUy3ZRdE8RTWs3jDtvrhb7xx6x3QOAABwE/Ep8Rq7daxKB5VW75q9TecA7vGZEGSeh8s8rOYlmmvBngU6FXfKdA4AAHADM3fO1OWky5rSeIr8vPxM5wCMkOzGsiyNbTBWfl5+Gr1ltNKd6aaTAACAQRvPbNRnxz/Ti1VeVJWCvzuIFDCCEZINFfIvpFHho7Tv0j69degt0zkAAMCQuOQ4Tdg2QeXyl1PP6j1N5wD/xQjJptrf316tS7bWwj0LdfzqcdM5AADAgKk7purqzaua0miKfL18TecA/8UIyaYsy9Lo+qMV4BOgUVtGKdWZajoJAAC40Nqf1+rrk1+rR7UeqlSgkukc4DcYIdlYgdwFNLr+aB365ZBW7l9pOgcAALjI1ZtXNXH7RFUKrqRu1bqZzgF+hxGSzbUNbav2oe21ZN8SHb1y1HQOAABwgSkRUxSfEq/JjSfLx+FjOgf4HUZIDjAyfKSCfIM0cvNIpaSnmM4BAABZ6JtT3+i76O/Uq3ovlc9f3nQO8IcYITlAvlz5NLHRRB27ekwLdi8wnQMAALLIhYQLmrRtkqoVqqYuVbqYzgH+FCMkh2haoqmeqvCU3jr0lrZf2G46BwAAZLJ0Z7pGbh6pdDtd05pMk7fD23QS8KcYITnI4DqDdX/Q/Rq1eZTikuNM5wAAgEy06uAqRcZGakT4CIXkDTGdA/wlRkgOkts7t6Y1maYrN69owrYJsm3bdBIAAMgEh345pIV7F6pNqTZ6pMwjpnOAW2KE5DBhBcLUp0Yfrfl5jT4/8bnpHAAAcI+S0pI07MdhCs4VrHENxsmyLNNJwC0xQnKgFyq/oDpF6mhqxFSdiT9jOgcAANyD2ZGzFR0frSmNpyjIL8h0DnBbGCE5kJfDS682flVelpdGbB6hNGea6SQAAHAXNp7ZqA+PfqjOYZ1Vv1h90znAbWOE5FDFAoppTIMxiroUpeX7lpvOAQAAd+hy0mWN3TpWFfJXUL9a/UznAHeEEZKDtb+/vR4q/ZCW7luqqEtRpnMAAMBtsm1bY7eM1Y3UG5rWZJp8vXxNJwF3hBGSw40MH6ki/kU0/MfhupF6w3QOAAC4DR8e/VCbzm3SwNoDVTZ/WdM5wB1jhORweX3zamqTqTp/47ymRkw1nQMAAG7hxLUTmhU5S42KN9KzFZ81nQPcFUYIVKtILXWr2k2fn/hc30d/bzoHAAD8iZT0FA3fNFz+3v6a3Ggyx/HCYzFCIEnqWb2nqhasqgnbJijmRozpHAAA8AcW7lmoI1eOaGKjiSqYu6DpHOCuMUIgSfJx+Ghqk6lKdaZq9JbRctpO00kAAOBXdlzYoVUHV+mJ8k+oeUhz0znAPWGE4L9KBZbSsLrDFHEhQm8fett0DgAA+I+45DiN3DxSpQJLaUidIaZzgHvGCMFvdCzXUS1DWmr+7vk6euWo6RwAAHI827Y1afsk/ZL0i6Y1nSZ/H3/TScA9Y4TgNyzL0viG45XPL5+G/ThMN9Numk4CACBH+/Lkl/ou+jv1rtlblQtUNp0DZApGCH4nf678mtxosk7EndCMnTNM5wAAkGNFx0Vr8vbJql2ktrpU7mI6B8g0jBD8oYbFG+rFKi/q42Mf69tT35rOAQAgx0lOT9aQjUPk5+WnaU2mycvhZToJyDSMEPypPjX7qEahGhq/bbxOx582nQMAQI4yc+dMHb16VFMaT1HRPEVN5wCZihGCP+Xj8NGMpjPkZXlpyMYhSklPMZ0EAECO8G30t/rw6IfqUrmLmpZoajoHyHSMEPylYgHFNLnRZB2+clizImeZzgEAINs7E39G47eOV7VC1dS3Vl/TOUCWYITgllqUbKF/hP1D7x95Xz/8/IPpHAAAsq2U9BQN+XGIHJZDM5vOlI/Dx3QSkCUYIbgtA2sNVJUCVTR2y1idvX7WdA4AANnSnF1zdOiXQ5rUaJLuC7jPdA6QZRghuC0+Xj6a2WymJGnoj0OVmp5quAgAgOxl7em1evfwu+pUqZNalWxlOgfIUowQ3LYSeUtoQqMJ2n95v+btnmc6BwCAbONcwjmN2TJGlQtU1qDag0znAFmOEYI70qZUGz1d4WmtPrRaG85sMJ0DAIDHS01P1dCNQ2XbtmY2mykfLz4HguyPEYI7NqTuEFUKrqRRm0fpQsIF0zkAAHi0BXsWaN/lfRrfcLxC8oaYzgFcghGCO+bn5aeZzWYqzZmW8fkQJ58PAQDgbvx49ketOrhKT1V4Su1C25nOAVyGEYK7UiqwlMY3HK+9l/bq9T2vm84BAMDjxNyI0cjNI1UhfwW9UvcV0zmASzFCcNfa399efy/3d71x4A1tPrfZdA4AAB4jzZmmYT8OU2p6qmY1myU/Lz/TSYBLMUJwT4bXG66y+cpq5KaRir0RazoHAACP8Pre17X74m6NbTBWoUGhpnMAl2OE4J7k8s6l2c1m62b6TQ3bNExpzjTTSQAAuLUt57Zoxf4V6liuo/5W+m+mcwAjGCG4Z6Xzldbo+qO1K3aXFkctNp0DAIDbir0Rq5GbR6psvrIaXm+46RzAGEYIMsXDZR7Wo2Uf1bJ9y7h/CAAAfyA1PVWDNg5SUlqSZjWbpdzeuU0nAcYwQpBpRoWPUqXgShq5aaROx582nQMAgFuZvnO69l3ap0mNJqlMvjKmcwCjGCHINLm8c2lui7lyOBzqv76/ElMTTScBAOAWvjjxhT48+qFeqPwC9wMBxAhBJiseUFwzmszQiWsnNH7beNm2bToJAACjDv9yWBO3TVS9ovXUv1Z/0zmAW2CEINM1LN5QfWv21TenvtE7h98xnQMAgDHXbl7TwA0Dlc8vn2Y0nSFvh7fpJMAtMEKQJbpW7aqWIS01O3K2ImMiTecAAOBy6c50Dd80XBcTL2pO8zkqkLuA6STAbTBCkCUclkNTGk9RSN4QDdk4RBcTL5pOAgDApRZFLdKW81s0InyEqhWqZjoHcCuMEGSZAN8AzWsxT4lpiRq0YZBS01NNJwEA4BLrT6/Xsn3L9FjZx/R4ucdN5wBuhxGCLFUmXxlNajRJUZeiNGPnDNM5AABkuZ/jf9bIzSMVViBMo+qPkmVZppMAt8MIQZZrF9pOL1R+QR8c/UBfnPjCdA4AAFkmMTVRA9YPkLfDW3Obz5Wfl5/pJMAtMULgEv1r9Ve9ovU0cdtEHf7lsOkcAAAynW3bGrd1nE7GndSMpjN0X8B9ppMAt8UIgUt4O7w1o+kM5fPLp4EbBiouOc50EgAAmertQ2/r2+hv1bdmXzW4r4HpHMCtMULgMgVyF9Cc5nN0MfGihm0apnRnuukkAAAyxc6YnZqza45alWylrlW6ms4B3B4jBC5VrVA1jQgfoS3ntmhx1GLTOQAA3LPYG7EasnGIQvKGaHKjyXwQHbgN3LYTLvd4uce1/9J+Ld23VFUKVlHzkOamkwAAuCup6akatHGQbqbd1Jvt3lSAb4DpJMAjGLsSYllWiGVZ6y3LOmRZ1kHLsvqbaoFrWZalUfVHKaxAmEZuGqmTcSdNJwEAcMds29bUHVO179I+TWo0SaXzlTadBHgMk2/HSpM02LbtMEn1JfW2LCvMYA9cyM/LT/Oaz5OPl4/6ru3LB9UBAB7nvSPv6eNjH6trla5qG9rWdA7gUYyNENu2L9i2vfs/f31d0mFJxU31wPWKBRTT/BbzdeHGBQ3eMFipTu6oDgDwDFvPbdWMnTPUIqSF+tXqZzone4o/L51Yb7oCWcQtPphuWVaopJqSIv7gaz0sy4q0LCvy0qVLrk5DFqtRuIbGNxyviJgITd8x3XQOAAC3dDLupIZsHKKy+cpqWpNpclhu8eNU9nFmp/TJi9K8qtI/e0jpaaaLkAWMfzDdsqwASZ9KGmDbdvz/ft227WWSlklSnTp1bBfnwQUeLvOwjl87rjcPvKky+cromYrPmE4CAOAPxSXHqe/avvLx8tFrLV+Tv4+/6aTsIS1FOvSZFLFEOrdL8guSwntK9bpLXsZ/XEUWMPp31bIsH2UMkHdt2/6nyRaY1b9mf526dkrTd0xXqcBSanhfQ9NJAAD8RqozVYM3DNaFGxe0st1K7oieGRIuSbvelHa+ISXESAXKSg/Okqo/I/lx0lh2ZmyEWBmHaL8h6bBt23NMdcA9eDm8NK3pNP3jm39oyMYheu/B9xQaFGo6CwCA/5q+Y7oiYiI0pfEU1Shcw3SOZ7uwL+Oqx/5PpPRkqWxrKfx1qUxLycHb23ICk3+XG0n6h6SWlmXt/c+vBw32wLA8Pnn0WsvX5OPwUd91nJgFAHAf7x95Xx8e/VBdqnTRw2UeNp3jmdLTpEOfSyvbS0ubSAc/k2r9Q+q9U+r0qVSuNQMkBzF2JcS27c2SuKUofqN4QHHNbT5XXb/vqiEbh2hR60XycfiYzgIA5GBbz2/V9B3T1bxEc/WvyW3N7he7QVQAACAASURBVFjSVWn3amnHCinutJSvpNR2ilSzk5Q7n+k6GMLchNupVaSWxjUYp+0XtmvGjhmmcwAAOdipuFMasnGISucrrWlNp8nL4WU6yXNcPCJ9NVCaEyatGSvlLyU99a7Ub6/UsA8DJIfjuAG4pUfLPqoT105o1cFVKpOvjJ6u+LTpJABADhOXHKe+6/rKx5FxElYenzymk9yf0ykdXyNtXyydXC95+UnVnpTCX5KKVjVdBzfCCIHbGlBrgE7FndK0HdNUKrCUGtzXwHQSACCHSHWmavDGwTqXcE5vtH1DxQO4n/JfSr4u7X1PilgqXTkh5S0mtRwj1e4i5Slgug5uiLdjwW15Obw0rck03R90vwZvHKzouGjTSQCAHGL6jumKuBChcQ3GqVaRWqZz3NeVk9I3w6XZlaRvhkr+wdLf35AG7JeaDmGA4E8xQuDWAnwD9FrL1+RteXNiFgDAJT448oE+PPqhXqj8gh4t+6jpHPdj29LJDdJ7T0sLakk7l0sVHpC6rZO6/SBVfVzy4lAZ/DVGCNxeibwlNLfFXJ1NOKshG4cozZlmOgkAkE1tO79N03ZMU7MSzTSg1gDTOe4lJVHatUpa3FBa/Yh0dqfU9BVpwAHp7yukErVNF8KDMELgEWoXqa2x9cdq+4XtmhoxVbZtm04CAGQzJ6+d1OCNg3V/0P2a1oSTsP4r7qy0Zpw0N0z6sr/k8JIeWSQNPCi1HCUFFjNdCA/EB9PhMR4r95ii46O18sBKFQsopm5Vu5lOAgBkExcTL6rnDz3l6/DVwlYLFeAbYDrJLNuWzkRknHJ1+EtJtlTxIan+y1LJBpLFrd5wbxgh8Cj9a/VXzI0Yzd89X0X8i6hDmQ6mkwAAHi4hJUG9fuila8nXtOqBVTn7JKy0ZOnAP6WIJdKFvVKuIKlBb6le94ybDAKZhBECj+KwHJrUaJIuJ13W2C1jVTB3QY7uBQDctVRnqgZtGKTj145rYauFCisQZjrJjOuxUuTKjF83LkoFK0h/myNVf1ry5f4oyHx8JgQex9fLV3NbzFVoUKgGbhioo1eOmk4CAHgg27Y1fut4bbuwTeMajFPj4o1NJ7ne+T3SP1+S5laWNk6T7qshdfqn1DtCqtuVAYIswwiBRwr0DdTi1ouVxyePev3QSxcSLphOAgB4mIV7F+qLE1+oV41eeqzcY6ZzXCc9TTr4L+mNdtKy5tKRr6Q6L0p9dknPfSyVbcVnPpDlGCHwWEXzFNXi1ouVmJaol394mXuIAABu28fHPtayfcvUsVxH9azW03SOayRekTbNkeZXkz5+QUqIkdpNlQYdkh6cIRUsa7oQOQgjBB6tfP7ymt9ivn6+/rMGrB+glPQU00kAADe38cxGTd4+WY2LN9bo+qNlZfc/9Y89JH3RT5oTJq2dIBUoKz3zgdR3t9SgV8aHzwEXY4TA49UrVk+TG01WZGykRm0eJaftNJ0EAHBTBy4f0Cs/vqKKwRU1u9ls+Tiy6Z29nenSkX9Lb3WQFjeQ9n0oVXtSenmb1PkLqUL7jPt9AIZwOhayhb+V/ptiE2M1d9dcFc1TVIPrDDadBABwM2fiz6j32t4KzhWs11u9Ln8ff9NJme9mvLTnHWnHMunqKSmwuNR6vFSrs+QfbLoO+C9GCLKNLpW7KOZGjFYdXKWieYrquUrPmU4CALiJKzevqOcPPeW0nVrSeokK5i5oOilz/XJCilgq7X1XSkmQQupLrcdJFTtIXvy4B/fD/yqRbViWpWF1hyn2Rqym75iuIv5F1LpUa9NZAADDktKS1HddX8UmxmpF2xUKDQo1nZQ5bFs6sS7jxoI/fS85fKQqf5fq95Tuq2m6DvhLjBBkK14OL01vOl3dvu+m4ZuGa3nu5apZmP8jBoCcKt2ZrmE/DtP+S/s1t/lc1Shcw3TSvUu5IUV9kHHl4/JRKU9hqfkIqXYXKW8R03XAbeGD6ch2cnnn0mstX1OxPMXUd11fnYw7aToJAGCAbduaumOq1p9ZrxHhI9SqVCvTSffm2mnp+zHSnErS14Mkn9zSY0ulgQek5sMZIPAojBBkS/lz5dei1ovkZXnp5TUvK+ZGjOkkAICLLdu3TB8e/VBdqnTRMxWfMZ1zd2xbit4ifdhJml9d2va6VLqF9OJ3Uo8NUvWnJW8/05XAHWOEINsKyRuixa0XKz4lXj3W9NCVm1dMJwEAXOTdw+9q4d6F6lC6gwbUGmA6586l3pT2vCstbSKtelCK3iw17CcN2Cc9+ZZUsj53NYdHY4QgWwsrEKaFrRbqQsIF9VzTU/Ep8aaTAABZ7LPjn2najmlqVbKVJjaaKIflQT/uxF+Q1k2W5laWPu8lpadJHeZLAw9JbSZIQSVMFwKZwoP+qQTuTu0itTW3xVz9dO0n9VnbR4mpiaaTAABZZM3PazRu6zg1KNZAM5rOkLfDQ87gObtL+rSbNK+K9OMsqURd6fnPpV7bpNovSL7Z8J4myNEYIcgRGhdvrOlNpivqUpQGrB+glPQU00kAgEy2+dxmDf1xqKoVrKZ5LebJ18vXdNJfS0+V9n8irWgtrWgpHf1WqtdD6rdbevYDqXRz3nKFbMtD/ngAuHdtQ9sqMS1RY7aM0dAfh2pWs1me8ydkAIC/tCt2lwauH6iy+crq9dZufjf0G5elXW9KO9+Qrl+QgktL7WdINZ6V/PKargNcgp/AkKM8WvZR3Ui9oWk7pmnslrGa3HiyZ71XGADwOwd/Oajea3urWEAxLW2zVIG+gaaT/ljMASlisbTvYyk9OeOUqw7zpbJtJAf/LkLOwghBjvNcpeeUkJKghXsXKo9PHo0MHymLy90A4JFOXDuhnmt6Ksg3SMvaLFNwrmDTSb/lTJeOfpNxV/PoTZKPf8YVj/CeUuGKpusAYxghyJF6VOuhhNQErTq4Snl986pfrX6mkwAAd+jM9TPq8X0PeTu8tbztchXNU9R00v+XdE3a87a0Y1nGTQaDQqQ2E6Vaz0u585uuA4xjhCBHsixLg2oPUkJqgpbvX648PnnUtWpX01kAgNsUeyNW3b/vrmRnst5s96ZKBpY0nZTh8k8ZVz32vi+l3pBKNZLaTpEqPCh58WMX8H/4pwE5lmVZGh0+WompiZq3e57y+OTR0xWfNp0FALiFqzevqseaHrqWfE0r2q5QufzlzAY5ndKJdRmf9zj+g+TlK1V9Qgp/SSpW3Wwb4KYYIcjRvBxemtx4shJTEzUlYory+ORRhzIdTGcBAP7E9ZTremnNSzqXcE6LWy9WlYJVzMUkJ0hR70sRS6VffpICikgtRkm1u0gBhcx1AR6AEYIcz8fho1nNZ6n3D701ZssY+fv4q1XJVqazAAD/IyktSX3W9tFPV3/S/JbzVbdoXTMhV6OlHcul3W9LyXHSfbWkjiuksEckbze/NwngJjgPDpDk5+WnBS0XqHKBynpl4yvacm6L6SQAwK8kpydrwPoB2ntpr6Y2naqmJZq6NsC2pVM/Su8/K82vkfG5j3Ktpa4/SD3WS9WeYIAAd4ARAvyHv4+/FrVepDL5yqjfun7adHaT6SQAgKSbaTfVd21fbTu/TeMbjNcDoQ+47sVTk6Tdq6UljaW3OkhntktNBksD9kuPr5RCDF2NATwcIwT4lSC/IK1ou0Jl8pVR//X9tfHMRtNJAJCjJaUlqc+6Ptp+YbsmNpqox8o95poXjjsnrZ0ozQmTvuib8XsPL5QGHpRajZEC73NNB5BNMUKA/xHkF6TlbZerfP7yGrBhgDac2WA6CQBypMTURPVZ20c7LuzQ5MaT9WjZR7P2BW1bOrND+riLNL+atGmOVKqh1PkrqedmqdY/JJ/cWdsA5BCMEOAPBPkFaVnbZaqYv6IGbhiodafXmU4CgBwlMTVRfdb1UWRspKY0nqKHyzycdS+WliLt+0ha3lJ6o410fG3GHc3775Wefle6v4lkWVn3+kAOxAgB/kSgb6CWtV2msOAwDd4wWD/8/IPpJADIERJTE/XyDy9rV+wuvdr41aw7Oj3hkrRhujSvivTP7lLydenBWdKgQ1K7KVL+0Kx5XQCMEOCv5PXNq6VtlqpywcoasnGIvo/+3nQSAGRrN1JvqOcPPRV1KUrTm0zX30r/LfNf5EKU9K+Xpblh0oZXpaJVpec+lXrvkOp1l/wCMv81AfwG9wkBbiHAN0BL2yzVyz+8rKE/DpXTduqB+114MgsA5BAJKQnq+UNPHbh8QNObTle70HaZ9+TpadLRr6XtS6TTWyWfPFKt56V6L0mFymfe6wC4LYwQ4Dbk8cmjxa0Xq9cPvTRs0zA5baceLP2g6SwAyDaup1xXzx966tDlQ5rZbKbalGqTOU+ceEXa83bGzQXjzkj5Skptp0g1O0m582XOawC4Y4wQ4Db9d4is7aURm0fIKaceKv2Q6SwA8HjxKfHquaanDv9yWLOazVKrUq3u/UkvHsm4oeC+D6XURCm0ifTANKlCe8nhde/PD+CeMEKAO+Dv469FrRap77q+GrV5lGzbzroPTAJADhCXHKeX1ryko1ePak7zOWpRssXdP5nTKR1fI21fLJ1cL3n5ZdzJPPxlqWiVzIsGcM8YIcAd8vfx18JWC/87RNLt9Kw/ux4AsqG45Dh1/767jl87rnnN56lZSLO7e6Lk69Ked6UdS6UrJ6W8xaSWY6TaXaQ8BTI3GkCmYIQAdyG3d24tbLlQ/db109gtY+W0nepYrqPpLADwGFdvXtVLa17KGCAt5qlpiaZ3/iRXTkoRy6Q970gp16US9aQWo6SwRyQvn8yPBpBpGCHAXcrlnUsLWi7QgPUDNG7rOMUnx+uFKi+YzgIAtxdzI0Y91vTQuevntKDlAjUu3vj2H2zb0qmNGadcHftWcnhLlR+T6veUitfOumgAmYoRAtyD/xsiIzeP1Oxds3Ul+YoG1hooizvrAsAfOhl3Ui+teUkJKQla2map6hStc3sPTEmU9n8kRSyVLh6S/AtKTV+R6naV8hbN2mgAmY4RAtwjXy9fTW8yXUG+QXrzwJuKS47TmPpj5O3gHy8A+LUDlw/o5R9elpflpTcfeFMVgyve+kFxZzOO1939lpR0NePGgo8skqr8XfLJlfXRALIEPyUBmcDL4aXR9UcrOHewlkQtUVxynKY3nS4/Lz/TaQDgFrad36b+6/srOFewlrdZrpDAkD//ZtuWzkRknHJ1+EtJtlTxIan+y1LJBhJXmwGPxwgBMollWepdo7fy+eXTtB3T1OuHXprfYr4CfANMpwGAUd9Hf6/hm4YrNChUS1svVSH/Qn/8jWnJ0oF/Ztzf48JeKVeQ1KC3VK97xk0GAWQbjBAgkz1X6TkF+QVpzOYxevG7F7W49WIVyM0RkQBypo+OfqTJ2yerRuEaeq3lawryC/r9N12PlSJXZvy6cVEqWEF6aK5U7SnJN4/rowFkOUYIkAUeKv2QAn0DNXjDYHX+trOWtlmq4gHFTWcBgMvYtq3l+5frtT2vqUnxJprdfLZye+f+7Ted35NxytWBTyVnqlSuXcYpV6Vb8JYrIJtzmA4AsqumJZpqWdtlunLzip7/9/M6fvW46SQAcAmn7dSMnTP02p7X9FDphzS/5fz/P0DS0zLecvVGW2lZc+nIVxknXPXdLT33kVSmJQMEyAEYIUAWqlm4plY9sEq2bHX+trP2XtxrOgkAslSqM1WjNo/SO4ffUadKnTSl8RT5OHykxCvSpjnS/GrSJ12khItSu6nSoENS++lSgTKm0wG4ECMEyGLl85fX6varlc8vn3qs6aHN5zabTgKALJGUlqQB6wfoq5NfqW/Nvhpad6gcF49IX/ST5oRJaydIBcpKz3wg9d0lNeiV8eFzADkOIwRwgRJ5S+it9m+pVGAp9V3bV/8++W/TSQCQqeKS4/TSmpe06ewmjQkfrR6+JWStflha3EDa96FU7Unp5W1S5y+kCu0lh5fpZAAGGf1gumVZKyU9JOmibdtVTLYAWa1g7oJa2W6l+q3rp+GbhismMUZdKnfh7uoAPN7Z62fVd11fRcdHa2axNmr33avS1WgpsITUerxUq7PkH2y4EoA7+dMrIZZl/duyrNAsfv1Vkh7I4tcA3EZe37xa0maJ2oW209xdczVu6zilpqeazgKAu7b34l4999VTio2L1pLYK2q3dYUUUFR6YpXUP0pqPJABAuB3/upKyJuSvrcs6y1JM2zbzvSflGzb/tEFQ+eeTPjyoA6djzedgWzG1qMq6O2lfx3/l747ekghqT3lJc7CR+a7cvKgJCn4IG99QSazbQXZHyky91oVTUvVgtgrOu/dWMMLPqpTaeWkzZI27zRdCQ8Xdl+gxnWobDoDWeBPR4ht2x9blvWNpDGSIi3LeluS81dfn+OCPlmW1UNSD0kqWZK7pSJ7sORQ4bRH5OssrAs+b+uU7zSFpPaRn13EdBoA/CU/5001TvpBV72/0DtBUvWbaWqT0FQz8j+iOC+ueAC4Pbf6TEiKpBuS/CTl1a9GiKvYtr1M0jJJqlOnju3q12d9I2s10O7Y5uq/vr8u556lec3nqU7ROqajkI1s3pwuSWrcuIHhEni8a6elHcuUvHu1xuX11tcBefRw/moa126JfP3yqrPpPgAe5U9HiGVZD0iaI+kLSbVs2050WRWQg9QqUkvvPfieeq/rre5rumt8g/F6pOwjprMAQLJt6eetUsRi6cjXuuLwUv/QstrrvKH+Nfupa9VuHK4B4K781ZWQUZKesG37oKtigJwqJDBE7zz4jgZtGKTRW0YrOj5afWv2lcPiFG0ABqTelA58mjE+YvZLufPrRL0u6p2wX5eTr2p2s9lqG9rWdCUAD/anP+HYtt0kqweIZVnvS9omqYJlWWcty+qala8HuLNA30Atbr1Yj5d/XCv2r9CQjUOUlJZkOgtAThJ/QVo3WZpbWfq8l5SeJnWYr61Pv6lOV7cp2U7TqgdWMUAA3DOj9wmxbfsZk68PuBsfh4/G1h+r0MBQzY6crQsJF7Sg5QIV8i9kOg1AdnZ2V8ZVj4P/kpzpUvkHpPo9pfub6aNjH+vVjQNVJl8ZLWy5UMUCipmuBZANGB0hAH7Psix1rtxZJfOW1LBNw/Tsv5/VwpYLVSG4guk0ANlJeqp06HMpYol0dqfkFyjV6yHV6y4Fl1a6M12zds7QO4ffUdMSTTWj6Qzl8eEocQCZgzecA26qRckWWt1+tZy2U89/87zWn15vOglAdnDjsvTjTGleVenTrlLiFan9DGnQIemBqVJwaV1Pua7+6/vrncPvqFOlTlrQYgEDBECm4koI4MYqBlfU+397X/3W9VO/9f3UvWp39a7RW14ObjwH4A7FHMh4y9W+j6X0ZKlMS6nDAqlsa8nx//9M8tjVYxq0YZDOXT+nUeGj9HTFpw1GA8iuGCGAmyvsX1hvtX9Lr0a8quX7l+vA5QOa3nS68ufKbzoNgLtzpktHv8l4y1X0JsnHX6r5nBTeUyr0+7d4fnXyK03cNlF5fPLojXZvqFaRWgaiAeQEjBDAA/h5+WlCwwmqXqi6pmyfoie/elJzms1R1UJVTacBcEdJ16Q9b0s7lmXcZDAoRGozUar1vJT793+AkZqeqpmRM/X+kfdVq3AtzWo2iwMxAGQpRgjgQTqW66gKwRU0eMNgdf62s4bXG64nyj/BzcIAZLh0TNqxVNr7vpR6QyrVSGo7RarwoOT1x//Kj7kRoyEbhyjqUpSeD3teA2oPkI/Dx8XhAHIaRgjgYSoXqKwPH/pQwzYN06TtkxR1KUqj649Wbu/cptMAmOB0SifWStsXZ/ynl69U9Qkp/CWpWPW/fGjEhQgN/XGobqbd1Kxms9QutJ2LogHkdIwQwAMF+QVpUatFWhK1REuilujIlSOa13yeQgJDTKcBcJXkBCnqfSliqfTLT1JAUanFaKn2C1LAX7+VyrZtrTywUgv2LFBoYKjmtpur0vlKu6YbAMQIATyWw3KoV41eqlqwqoZvGq6nvnpKrzZ5Vc1DmptOA5CVrpySdq6Qdq+WkuOl4rWljiuksEckb99bPvx6ynWN3jxa686sU7vQdprQcALH7wJwOUYI4OGalGiijzp8pIHrB6rvur4c4wtkR7adcbrV9iXS0X9LDq+M0RH+shRS97af5tfH7w6tO1SdKnXiM2UAjGCEANlA8YDievvBtznGF8huUpOk/R9nvOUq9oDkX0BqMliq21UKvO+OnorjdwG4E0YIkE387zG+T3z5hKY2maq6RW//T0kBuIm4c1LkG1Lkm1LSFalIFenhhVLVxyWfOzuEIjE1UbMiZ+njYx9z/C4At8EIAbKZjuU6qmJwRQ39cai6ftdVL1Z5Ub1r9JaPF0duAm7NtqWzOzNOuTr8hWQ7M47WDe8phTaW7uJtUwd/OajhPw7Xz/E/q0vlLupbqy/H7wJwC4wQIBsKKxCmjx76SDN2ztAbB97Q1vNbNa3pNJUO4vQbwO2kpUiHPssYH+d3S35BGcOjXncpf+hdPWW6M11vHnxTr+95XcG5g7W87XKFFwvP3G4AuAeMECCb8vfx1/iG49WkRBON3zpeT335lIbUGaInKzzJB1EBd5BwSYpcmfG2q4RYqUA56cFZUvVnJL+Au37a8wnnNXLzSO2K3aV2oe00pv4YBfkFZWI4ANw7RgiQzbUq2UrVClbTmC1jNDlisjad26QJDSeoQO4CptOAnOlCVMYpVwc+kdJTpLJtMq58lGkpORz39NRfn/xaU7ZPkVNOTWk8RR1Kd+APHQC4JUYIkAMU8i+kRa0X6f0j72tO5Bx1/KKjJjWapKYlmppOA3KG9DTp6NcZ4+P0Vsknj1Src8ZdzQuWu+enj0+J15TtU/TvU/9WjUI1NLXJVJXIWyITwgEgazBCgBzCYTn0XKXnVK9oPQ3fNFy91/bWUxWe0uA6g5Xb+85O2wFwmxKvZNxUcOcKKe6MlK+k1HaKVLOTlDtfprxEZEykRm4eqYuJF9W7Rm91q9pN3g7+9Q7AvfH/UkAOUy5/Ob33t/e0YPcCrT60Wjtjdmpak2mqVKCS6TQg+7h4RIpYIkV9IKUlSaFNpPbTpfIPZNxoMBOkpqdqUdQivbH/DZXIW0Kr269WtULVMuW5ASCrMUKAHMjPy0+v1H1FjYs31ujNo/Xsv59V35p91TmsM3daB+6W0ykdX5NxytXJ9ZJ3LqnqExmf9yhaJVNf6lTcKQ3fNFyHfjmkjuU6aljdYfL38c/U1wCArMQIAXKwBvc10KcPf6oJ2yZo7q65Wnd6ncY3GK+y+cuaTgM8R/J1ac+70o6l0pWTUt77pJZjpNpdpDyZewBEqjNVbx96W4v3Lpaft5/mNp+r1qVaZ+prAIArMEKAHC5frnya03yOvjr5lWbsnKEnvnpCXat0VY9qPeTr5Ws6D3BfV05KEcukPe9IKdelEvWklqOlSg9LWXBz0IOXD2rc1nE6evWoWoa01Kj6o1TYv3Cmvw4AuAIjBIAsy1KHMh3UqHgjzdw5U0v3LdV30d9pfMPxql2ktuk8wH3YtnRqY8YpV8e+lRzeUuXHpPo9peJZ889KYmqiFu5dqHcPv6sCuQpoXvN5alWqVZa8FgC4CiMEwH8F5wrW1CZT9VDphzRp+yS98O0Lerz84xpYe6ACfQNN5wHmpCRK+z6UIpZKlw5LeQpJzf5fe/cdXld1oGv8XSpustx7kQtuuAFGLsEmEGODMT30lBuqDblzJ8ykTEhmMsnkmdxMSCa5GRgXDCSBACEEAgTTMiQEkbjjhsHgbgt3427LKuv+sQVxEhsksM4+kt7f8+jROdKxzsdiq3xnr73WV6D4BijsUmdPW1Jawrf/9G3ePvA2Vw24ittOv43CJoV19nySlCmWEEl/Y2z3sTx28WP89+L/5v7X7+eljS9x++jbmVA0wY3P1Ljs2QTz7oZFP4VD70CX4XDpNBjySchvVmdPu/PQTr43/3vMXjubvq378tNJP2VE5xF19nySlGmWEEnH1CK/BV8a+SXO73s+3/rjt/jH3/8jn+j5Cb42+mt0Kai7V36l1MUIG+YkS+y+/hQQYdCFMOZWKPoY1GERjzHy5OonuWPBHRwoP8DnT/k8Nw670euzJDU4lhBJ72tI+yE8eMGDPLDiAe5afBeXPnEpXxjxBa4eeDU5ISfteNKJU1EGyx+DudNg8xJo1gbO+DsYeVOyyWAd27h3I9+a8y3mbp7LaZ1O45sf+yZ92/St8+eVpDRYQiR9oLycPK4beh3n9DqHb//p23xn7nd4es3TLuerhmHfVlhwb/J2YBt0HAQX/hCGXw1NCur86curyvnZaz9j2pJp5Ofk8y9j/oUrBlxhyZfUoFlCJNVYz8KezJg4473lfK946gquHHAlnz/187Rt1jbteFLtlC5KplwtfwyqyqH/eckqV30/UadTrt4VY+QPm/7A9xd8n3V71zGhaAJfHfVVOhd0rvPnlqS0WUIk1cq7y/mO6z6OuxbfxS/f/CVPr32aqcOn8qlBnyK/DvZHkE6YyvLkOo+502HjXGjSEkbeCKOmQPuTMhbjzXfe5I75dzBn8xx6t+rNnePv5KyeZ2Xs+SUpbSHGmHaGGisuLo4LFixIO4ako6zevZrvL/g+JaUlFBUW8Y/F/8j4nuNdRStLlJSUADBu3LiUk6Ts4C5Y+BOYPwv2lkLbPjB6Kpz6aWiWueWndxzawV2L7+Kxtx6jsEkht55yK1cNvIr8HMu7lG1CCAtjjMVp52ioPBMi6SM5qc1JTJswjVdKX+GO+Xdw2+9uY2SXkXy5+Muc3P7ktOOpsdu6IrnQfOkjUHEY+pwFF/wA+p8LObkZi1FWWcb9K+5n1rJZlFWU8alBn+KWU26hddPWGcsgSdnEEiLphBjbfSyju47mV2/+irsW38XVv7maS/pdwt+f9vd0bNEx7XhqTKoq4c3nkvKx9g+Q1xxOuQZGBHbSwwAAIABJREFU3wKdMluMY4w8t/45frTwR5TuL+XsnmfzxdO/SO/WvTOaQ5KyjSVE0gmTl5PH1YOu5vy+53P30rt54PUHeG7dc9w49EY+N+RzNMuru83dJA7vgVcfgHkz4Z110KoHTPgmjPgctGiX8TjLti/je/O/x+LtixnQdgCzzp3F6K6jM55DkrKRJUTSCdeqSSu+WPxFrhpwFT9c9EPuXHwnj771KF8Y8QUm95ns0qM6sXasgnkzYPGDcGQ/9ByTlI9BF0Fu5n/NbTmwhR8t+hFPr3ma9s3a860zvsUlJ11Cbganf0lStrOESKozPVv15D/P/k8WbFnA9+Z/j9tfvp17lt3D1FOmcm6vcy0j+vBihNUvJqtcvfU85OTD0MuTJXa7nZZKpM37N3PP8nt47K3HCARuHnYzNw67kYL8ut9rRJLqG0uIpDpX3KWYhy98mGfWPsOMpTP48ktfZnrr6UwZPoXzep/nK8SquSMHYMnDMHcG7FgJBZ3g7Nvh9OuhMJ39Nd7e/zazls3i8VWPA3BZv8u4adhNdGvZLZU8klQfWEIkZUROyOGCvhcwqfckXlj/AjOWzuCfXv4npi9Nysik3pPIy/FHko5j94bkWo9FP0uu/eh6Klw2A4ZcBnlNU4m0ad8mZi2bxROrniCEwOX9L+fGoTfStWXXVPJIUn3ib3xJGZWbk8ukPpM4t/e5vLD+BaYvmc7tL9/O9CVJGZncZ7JlRIkYYf0fk1Wu3ngaCDD4Yhh9K/QclZFdzY9l496NzFw2k6dWP0VuyOXKgVdyw9Ab6FLQJZU8klQf+ZteUipyQg7n9T6Pib0m8uKGF5m+ZDpfL/k605dM5+ZhN3PhSRe6gVtjVX4Ylj+aXO+xZRk0bwtjvwAjb4LWPVKLtX7vemYuncnTa54mLyePawZdw/VDrqdzQTrTwCSpPrOESEpVTshhQq8JjC8az+82/o4ZS2bwjT9+gxlLZzBl+BQuOukiy0hjsW9LsqP5gvvg4A7oNBgu+jEMuxKatEgt1po9a7h76d3MXjubJjlN+NTJn+L6Ide7/40kfQSWEElZISfkcE7ROYzvOZ6XNr3E9CXT+dc//it3Lb6LqwZcxRUDrqB98/Zpx1Rd2LQwmXL12uPJRoMDz082Fuzz8dSmXMUY+dPbf+KhNx7ipU0v0SyvGZ89+bNcN/Q6OjTvkEomSWpILCGSskoIgbN7ns1ZPc6ipLSE+1fcz52L72TG0hlM6j2Jawddy7COw9KOqY+qshxWPAFzpkHpAmjaCkZNgVE3Q7u+qcXaf2Q/T65+kofeeIh1e9fRrlk7bhp2E58++dOWYEk6gSwhkrJSCIEze5zJmT3OZM2eNTz8xsM8seoJnlrzFMM6DOPaQddyXu/zaJLbJO2oqo0DO2DhfTD/Hti3GdqdBOffAadeC00LU4t19DF2sOIgwzoM4zvjvuMxJkl1JMQY085QY8XFxXHBggVpx5CUkmO9Sn3FgCu4csCVrkx0HCUlJQCMGzcu3SBbliUXmi/9JVSWwUnjk1Wu+k2AnHQ2raysquTl0pd58PUH+dPmP5Gfk+/ZNknvCSEsjDEWp52jobKESKp3qmIVczbP4aHXk/n6715Pcu2gazm98+mElK4jyEaplpCqSlg5G+ZMh/UlkN8CTrkmud6j48DM56m2p2wPj7/1OA+vfJjS/aV0atGJqwdezeX9L3fKlaT3WELqltOxJNU7OSGHM7qdwRndzmDjvo08svIRHnvrMZ5f/zwD2g7g4pMuZlLvSS6dmpZDu+HV+5PNBXdvgNZFMPHbMOKzyXK7KaioqmDe5nk8vfZpnl/3PIcrD3N659P5h9P/gfFF412BTZIyzDMhkhqEQxWHmL1mNo+++SjLdy4nEBjZZSQX9L2ACb0m0KpJq7QjpiKjZ0K2v5lMuVryEJQfhF5jk7MeAydDbuZf84oxsmzHMmavnc2za59l5+GdFOYXcm7vc7l20LUMbJfe2RhJ2c8zIXXLEiKpwVm3Zx3PrH2Gp9c+zfq968nPyefjPT7O5D6TOavnWTTNbZp2xIyp8xJSVQWr/ydZ5Wr1/0Bu02Rfj9FToevwunnOD7B2z1pmr53N7DWz2bBvA01ymnBWz7OY3GcyZ/Y4s1H9/5f04VlC6pYlRFKDFWNkxc4V/GbNb3h23bPsOLSDlvktOafoHC7oewGjuowiNyc37Zh1qs5KSNl+WPwgzJsBO1dByy7JjuanXwctM7+J39YDW3l23bPMXjubFTtXEAiM6jqKC/pcwDm9zmm0Z8IkfXiWkLrlNSGSGqwQAkM6DGFIhyF8qfhLzNsyj9lrZ/Pb9b/lidVP0KF5Byb1nsSEXhM4peMp5OX4I/ED7VoL8+5Orvko2wvdT4dPzoLBl0BeZpey3XZwGy9vepln1j7DvC3ziESGtB/Cl4u/zPl9zndHc0nKYp4JkdToHK44zMulLzN7zWxe2vQS5VXlFOYXMqbbGMZ1H8cZ3c5oMEv+npAzITHCupeTVa5WzoacXBh8aXK9R8+RJyjpByuvKmfxtsWUlJbwSukrrHxnJQBFhUVc0PcCJveZTO/WvTOWR1LD5pmQuuXLfpIanWZ5zZjYayITe01k35F9zN08l5LSEl4ufZkX1r8AQL82/Tiz+5mM7T6W0zqd1jg3rCs/BMt+CXNnwNbl0KI9nPlFGHkjtOqWkQib92+m5O0SSjaVMHfLXA6UHyAv5HFa59O4bcRtjOs+jgFtB7gssyTVM5YQSY1aYZNCJvSawIReE4gxsmr3Kl4pfYWS0hLuf/1+7nvtPprnNWd019GM6zaOcT3G0b1l97Rj1609pTB/Fiz8CRzaBZ2HwsV3wrArIL95nT51WWUZC7cufO//wZo9awDoWtCVyX0mM7b7WEZ3GU3LJi3rNIckqW5ZQiSpWgiB/m37079tf64beh0Hyw8yb8s8SkpLKCkt4fcbfw9zk+k/wzsOZ1iHYQzrMIyB7QbW/zMlMcKm+ckqVyueAGKytO6YW5OlduvoTMPWA1tZvmM5y3Yse+/tUMUh8nPyKe5czOX9L2dc93H0ad3Hsx2S1IBYQiTpOFrkt+Dsnmdzds+ziTGybu86Xil9hXlb5jFn8xx+s+Y3AOTn5DOo3SCGdhj6XjEpalVETshJ+b+gBiqOwIpfJ+Xj7UXQtHVSPEZNgba9TuhT7T+yn9d2vsayHcuS4rF9GdsObQMgLyePgW0Hclm/yxjbfSzFnYtpkd/ihD6/JCl7WEIkqQZCCPRp3Yc+rfvwmcGfIcbI1oNb33v1fvmO5fx61a956I2HgGSa19D2Q98rJv3a9KNry67ZswLX/u2w4F5YcA/s3wrt+8Pk78Mp10LTjz7VaU/ZHjbs3cDru15n6falLN+xnDV71hBJFkPp1aoXI7uO/IuzSe7fIUmNR6q/DUMIk4D/B+QCs2KM300zjyTVVAiBLgVd6FLQhYm9JgJQWVXJmj1rWL5jOUt3JH9437v8XipjJQB5IY9uLbvRs1VPigqL6NWqFz0Lk9vdW3YnPze/7oNvXpKscrX8Uag8Av0mwphboO94yKn5mZsYI7vLdrNh3wY27N3w3vuN+zayYd8G9pTtee+x7Zq1Y2iHoUzqM4lhHYYxtMNQWjdtXRf/dZKkeiK1JXpDCLnAm8BEYBMwH7g2xrjieP/GJXol1TeHKg6xctdK1u1d9xd/rG/Yt4ED5Qfee1xOyKFrQVeKCosoalVEl4IutGnahtZNW9O6SevkffVbs9xmNb4+oqSkBKoqGdf+HZg7HTb8CfIL4NRPJbuad+j/F4+vrKpk35F97C7bzZ4je9hTlrztPbKXnYd2vlcyNu7dyL7yfX+T/91SVdSqiJ6FPRnYbiDdCrp5PYekescleutWmiXkY8A3Y4znVd+/HSDG+H+P928GDRoUZ82alaGEklR3YozsK9/HtgPb2Hpoa/L+4Fa2H9zOloNbOFRx6Lj/Ni/kUZBfQMsmLWmR14KC/AJa5LcgN/zV7u+V5Wxf9joc2E7HTpWQ1xzaFEGrbsTcPI5UHuFA+YH33vaX73/f580hh/bN29O5RWc6tehEpxad6FyQ3O7QvAP5ORk4kyNJGXLmmWdaQupQmtOxugMbj7q/CRj91w8KIUwBpgB07tw5M8kkqY6FEGjVpBWtmrSiX9t+f/P5wxWHOVhxkANHDnCgorokHNmfvK/Yz8EjB9+7vePQDg7uO0hVrEr+cVUllB+GysMcPFQOeXlsKWyT7GheuRfe2QtA09ymFOQXUNikkC4FXSjIL0jKTX7L924ffb9FXgtyajFlS5Kk48mSKySPL8Y4E5gJyXSsj7TrryQ1VFVVsOqFZJWrNb+Dps1g2JWUhNHQrs9H2zFdkqQTLM0SUgr0POp+j+qPSZJq6vBeWPwgzJsBu9ZAYTc45xsw4jooaA8lJWknlCTpb6RZQuYD/UMIfUjKxzXAp1LMI0n1x641MHcmvPoAHNkHPUbB+H+Gky+GTKyyJUnSR5BaCYkxVoQQ/g54jmSJ3ntjjK+llUeSsl6MsOb3ySpXbz4HOXkw5LJkid3up6edTpKkGkv1mpAY42xgdpoZJCnrHTkIS38Bc2fA9tehoCOc9RUovgEKu6SdTpKkWsv6C9MlqdHavRHm3w0LfwqHd0OX4XDpNBjySchvlnY6SZI+NEuIJGWTGGHDHJg7DV7/DRBh0IUw5lYo+hi46Z8kqQGwhEhSNqgog+W/Sq732LwEmrWBM/4ORt4MbXp+8L+XJKkesYRIUpr2bYUF98CCe+HAdug4CC78IQy/GpoUpJ1OkqQ6YQmRpDSULkrOeix/DKrKof95ySpXfT/hlCtJUoNnCZGkTKksh9efSsrHxrnQpCWMvBFGTYH2J6WdTpKkjLGESFJdO7gLFt4H8++BvaXQtg9M+i6c+mlo1irtdJIkZZwlRJLqytbXkrMeSx+BisPQ5yy44AfQ/1zIyU07nSRJqbGESNKJVFUJbz6blI+1f4C85slF5qNvgc6D004nSVJWsIRI0olweA+8+gDMmwnvrINWPWDCN2HE56BFu5TDSZKUXSwhkvRR7FgF82bA4gfhyP5kQ8EJ30o2GMz1R6wkScfib0hJqq2qKljzIsyZDqtegNwmMPRyGD0Vup2WdjpJkrKeJUSSaqpsPyx5KJlyteNNKOgEZ98OxTdAy05pp5Mkqd6whEjSB3lnfVI8Ft0PZXug66lw2UwYchnkNUk7nSRJ9Y4lRJKOJUZYV5KscrVyNhBg8MUw+lboOcpdzSVJ+ggsIZJ0tPLDsOyXMHcGbF0GzdvC2C/AyJugdY+000mS1CBYQiQJYO/byY7mC++Dgzuh02C46Mcw/CrIb552OkmSGhRLiKTGbeN8mDsNVjyRbDQ48PxkY8E+H3fKlSRJdcQSIqnxqTiSlI6506B0ITRtBaOmwqiboF3ftNNJktTgWUIkNR77t8PCn8D8WbB/C7Q7Cc6/A069FpoWpp1OkqRGwxIiqeHbvDRZ5WrZo1BZBieNh4v/C/pNgJyctNNJktToWEIkNUyVFcnSunOnw/pXIL8FnPaZZFfzjgPTTidJUqNmCZHUsBx6Bxb9DObNgj0boHURTPw2jPhsstyuJElKnSVEUsOwfWVy1mPJw1B+EHqNhfP+HQZOhlx/1EmSlE38zSyp/qqqglW/TVa5Wv0i5DaBYVcmS+x2HZ52OkmSdByWEEn1T9k+WPxgsqv5rtXQsgt84p/h9OugZce000mSpA9gCZFUf+xaC/NmwqsPQNle6H46fHIWDL4E8pqknU6SJNWQJURSdosR1v4hud5j5TOQkwuDL4Uxt0KP4rTTSZKkD8ESIik7lR+CpY8k5WPbCmjRHs78Ioy8EVp1SzudJEn6CCwhkrLLnlKYf3eys/mhd6DzMLjkLhh6BeQ3SzudJEk6ASwhktIXI2ycl6xyteJJICZL6465NVlqN4S0E0qSpBPIEiIpPRVl8NrjMGcabF4MTVsnxWPUFGjbK+10kiSpjlhCJGXe/m2w4F6Yfw8c2AYdBsAFP4Dh10DTlmmnkyRJdcwSIilz3l6cXGi+/FdQeQT6TYQxt0Df8ZCTk3Y6SZKUIZYQSXWrsgLeeCrZWHDDnyC/AEZ8DkZPhQ79004nSZJSYAmRVDcO7oJFP4V5s2DvJmjTC879dzjtM9C8TdrpJElSiiwhkk6sba8nU66W/AIqDkHvM2Hy92DApGSjQUmS1OhZQiR9dFVV8NZzySpXa1+CvGYw/CoYNRW6DE07nSRJyjKWEEkf3uG9sPjnyfUe76yFwm5wzjdgxHVQ0D7tdJIkKUtZQiTV3s7VSfFY/HM4sh96jIJz/gVOvhhy89NOJ0mSspwlRFLNxAhrfgdzpsNbz0NOHgy5LFlit/vpaaeTJEn1iCVE0vs7cgCWPJyc+dixEgo6wllfgeIboLBL2ukkSVI9ZAmRdGy7N8C8u2HRz+DwbugyHC6dBkMvh7ymaaeTJEn1mCVE0p/FmGwoOGcavPGb5GMnXwSjb4Gij0EI6eaTJEkNgiVEElSUwfJfJeVjy1Jo1gbO+D8w8mZo0zPtdJIkqYGxhEiN2b4tMP8eWHgfHNgOHQfBhT+E4VdDk4K000mSpAbKEiI1RqULk1WuXnscqipgwHkweir0/YRTriRJUp2zhEiNRWU5rHgiWeVq0zxoUggjb4RRU6D9SWmnkyRJjYglRGroDuxMplvNvwf2vQ1t+8Ck78Kpn4ZmrdJOJ0mSGiFLiNRQbX0tudB82S+h4jD0PTu53qP/uZCTk3Y6SZLUiFlCpIakqhJWPgNzp8O6lyGvOZxyTbLEbqeT004nSZIEWEKkhuHQbnj1AZg3E3avh1Y9YMI3YcTnoEW7tNNJkiT9BUuIVJ/teCu50Hzxg1B+INlQcOK/waALIddvb0mSlJ38K0Wqb6qqYPWLMHcarPot5DaBoZcnU666nZp2OkmSpA9kCZHqi7L9sOSh5MzHzregZWc4+2tQfD207JR2OkmSpBqzhEjZ7p11MO9uWHQ/lO2BbqfBZTNhyGWQ1yTtdJIkSbVmCZGyUYywriRZ5WrlbCDA4Ith9K3Qc5S7mkuSpHotlRISQrgS+CZwMjAqxrggjRxS1ik/BMseTcrH1uXQvB2MvQ1G3gStu6edTpIk6YRI60zIcuCTwIyUnl/KLnvfhvmzYOFP4OBO6DQYLvoxDL8K8punnU6SJOmESqWExBhfBwhOKVFjt3F+ssrViieSjQYHnp+sctXn4065kiRJDVbWXxMSQpgCTAEoKipKOY10AlQcSUrH3GlQuhCatoJRU2HUTdCub9rpJEmS6lydlZAQwm+BLsf41NdjjE/U9OvEGGcCMwGKi4vjCYonZd7+7bDwPph/D+zfAu1OgvPvgFOvhaaFaaeTJEnKmDorITHGCXX1taV6ZfPS5ELzZY9CZRmcNB4u/i/oNwFyctJOJ0mSlHFZPx1LqpcqK5KldedOh/WvQH4LOO0zMHoqdByYdjpJkqRUpbVE72XAfwEdgadDCItjjOelkUU6oQ69A4t+BvNmwZ4N0LoIJn4bRnwWmrdNO50kSVJWSGt1rMeBx9N4bqlObF+ZnPVY8jCUH4Re42DSd2DgZMjJTTudJElSVnE6lvRhVVXBqt8mq1ytfhFym8KwK5MpV12Hp51OkiQpa1lCpNoq2weLH4S5M2DXamjZBT7xz1B8PRR0SDudJElS1rOESDW1ay3MmwmvPgBle6F7MVx+D5x8MeQ1STudJElSvWEJkd5PjLD2D8n1HiufSa7vGHwpjLkVehSnnU6SJKlesoRIx3LkICx7JJlytW0FtGgPZ34RRt4IrbqlnU6SJKles4RIR9uzCebPgoU/SZbb7TwMLrkLhl4B+c3STidJktQgWEKkGGHj3GTK1YongZgsrTvmVug1FkJIO6EkSVKDYglR41VRBq89DnOmwebF0Kw1fOzzMPJmaNsr7XSSJEkNliVEjc++rbDwPph/DxzYBh0GwAU/gOHXQNOWaaeTJElq8CwhajzefhXmTIfXHoPKI9BvIoy5BfqOh5yctNNJkiQ1GpYQNWyVFfDGU0n52DgH8gvg9Otg1BTo0D/tdJIkSY2SJUQN08FdyQpX82fB3lJo0wvO+w6c9pnk2g9JkiSlxhKihmXrimSVq6WPQMUh6H0mTL4DBkxKNhqUJElS6iwhqv+qquCt55JVrta+BHnNYPhVMGoqdBmadjpJkiT9FUuI6q/De+HVB2DeTHhnLRR2g3O+ASOug4L2aaeTJEnScVhCVP/sXA1zZ8Din8OR/dBzNJzzL3DyxZCbn3Y6SZIkfQBLiOqHGGHN75JVrt56HnLyYOgnYfRU6H562ukkSZJUC5YQZbcjB2DJw8mZjx0roaAjnPUVKL4BCruknU6SJEkfgiVE2Wn3Bph3Nyz6GRzeDV2Gw6XTYOjlkNc07XSSJEn6CCwhyh4xwoY/JatcvfGb5GMnXwSjb4WiMRBCuvkkSZJ0QlhClL7yw7D8V8n+HluWQrM2cMb/gZE3Q5ueaaeTJEnSCWYJUXr2bYH598CCe+HgDug4CC78UbLHR5OCtNNJkiSpjlhClHmlC5NVrl57HKoqYMB5MPoW6Hu2U64kSZIaAUuIMqOyHFY8kaxytWkeNCmEkTfCqCnQ/qS000mSJCmDLCGqWwd2wsL7kmlX+96Gtn1g0n/AqZ+CZq3STidJkqQUWEJUN7Ysh7nTYOkvobIsmWp14Q+h/7mQk5N2OkmSJKXIEqITp6oSVj6TrHK17mXIa56c8Rg9FTqdnHY6SZIkZQlLiD66Q7vh1Qdg3kzYvR5a9YAJ34IR/wtatEs7nSRJkrKMJUQf3o63krMeix+C8gNQ9DGY+G8w6ELI9dCSJEnSsfmXomqnqgpWv5hc77Hqt5DbBIZeniyx2+3UtNNJkiSpHrCEqGbK9sOSh5Ildne+BS07w9lfg+LroWWntNNJkiSpHrGE6P29sw7m3Q2L7oeyPdDtNLhsJgy5DPKapJ1OkiRJ9ZAlRH8rRlhXklzvsXI2EGDwJcmUq56j3NVckiRJH4klRH9WfgiWPZqUj63LoXlbGHsbjLwJWndPO50kSZIaCEuIYO/bMH8WLPwJHNwJnQbDRT+G4VdBfvO000mSJKmBsYQ0ZhvnJ6tcrXgi2Whw4GQYcwv0PtMpV5IkSaozlpDGpuIIrPh1MuWqdCE0bQWjpsKom6Fdn7TTSZIkqRGwhDQW+7fDwvtg/j2wfwu07wfn3wGnXgtNC9NOJ0mSpEbEEtLQbV6anPVY9ihUlsFJ58Aldybvc3LSTidJkqRGyBLSEFVWwMqnYc502PBHyG8Bp30GRk+FjgPTTidJkqRGzhLSkBx6Bxb9LNlccM9GaF0EE78NIz6bLLcrSZIkZQFLSEOw7Q2YNwOWPAzlB6HXOJj0f5PVrnJy004nSZIk/QVLSH1VVQWrXoA502DN7yC3KQy7Mply1XV42ukkSZKk47KE1Ddl+2DxgzB3BuxaDYVdYfw/w+nXQ0GHtNNJkiRJH8gSUl/sWpNc6/HqA1C2F7oXw+X3wMkXQ16TtNNJkiRJNWYJyWYxwtqXklWu3nw2ub5j8KUw5lboUZx2OkmSJOlDsYRkoyMHYdkjyZSrbSugRXv4+Jeg+AZo1S3tdJIkSdJHYgnJJns2JVOuFv00WW638zC45C4YegXkN0s7nSRJknRCWELSFiNsnJuscvX6U0BMltYdcyv0GgshpJ1QkiRJOqEsIWmpKIPXHk/Kx+bF0LR1UjxGTYG2vdJOJ0mSJNUZS0im7dsKC+5N3g5sgw4D4IIfwPBroGnLtNNJkiRJdc4Skilvv5qscrX8V1BVDv3PTTYW7DsecnLSTidJkiRljCWkLlVWwBtPJeVj4xzIL4Di62HUVOjQL+10kiRJUiosIXXh4C5Y+BOYPwv2lkKbXnDed+C0z0Cz1mmnkyRJklJlCTmRtq6AudNh6SNQcQj6fBwmfx8GnJdsNChJkiTJEvKRVVXCm8/B3Gmw9g+Q1wyGXwWjb4HOQ9JOJ0mSJGWdVEpICOEO4CLgCLAauD7GuDuNLB/a4T3w6s9h3gx4Zx206g7n/Cucfh20aJd2OkmSJClrpXUm5AXg9hhjRQjhP4DbgX9KKUvt7FwNc2fA4p/Dkf3Qc3RSPk6+CHLz004nSZIkZb1USkiM8fmj7s4BrkgjR43FCKtfTK73eOt5yMmHoZ9Mplx1H5F2OkmSJKleyYZrQm4AfpF2iONa/SI881XYsRIKOsJZX02W2S3sknYySZIkqV6qsxISQvgtcKy/1L8eY3yi+jFfByqAn7/P15kCTAEoKiqqg6QfIK855DeDS6cnZz/ymmY+gyRJktSA1FkJiTFOeL/PhxCuAy4Ezokxxvf5OjOBmQDFxcXHfVydKRoDU16CEDL+1JIkSVJDlNbqWJOArwBnxRgPppGhxiwfkiRJ0gmVk9Lz3gkUAi+EEBaHEKanlEOSJElShqW1Ola/NJ5XkiRJUvrSOhMiSZIkqZGyhEiSJEnKKEuIJEmSpIyyhEiSJEnKKEuIJEmSpIyyhEiSJEnKKEuIJEmSpIyyhEiSJEnKKEuIJEmSpIyyhEiSJEnKKEuIJEmSpIyyhEiSJEnKKEuIJEmSpIyyhEiSJEnKKEuIJEmSpIwKMca0M9RYCGE7sD6Fp+4A7Ejheesrx6v2HLPacbxqx/GqHcerdhyv2nG8aifN8eoVY+yY0nM3ePWqhKQlhLAgxlicdo76wvGqPcesdhyv2nG8asfxqh3Hq3Ycr9pxvBoup2NJkiRJyihLiCRJkqSMsoTUzMy0A9QzjlftOWa143jVjuNVO45X7TheteN41Y7j1UB5TYgkSZKkjPJMiCRJkqSMsoQcQwjhjhDCGyGEpSGEx0MIbY7zuEkhhJUhhFUhhK9mOme2CCFcGUJ4LYRQFUI47gqiSkSAAAAGFUlEQVQWIYR1IYRlIYTFIYQFmcyYbWoxZh5jQAihXQjhhRDCW9Xv2x7ncZXVx9fiEMKTmc6Zpg86VkIITUMIv6j+/NwQQu/Mp8wuNRiz60II2486pm5KI2c2CCHcG0LYFkJYfpzPhxDCj6vHcmkIYUSmM2aTGozX2SGEPUcdW9/IdMZsEkLoGUL4XQhhRfXvxi8c4zEeYw2MJeTYXgCGxhiHA28Ct//1A0IIucBdwPnAYODaEMLgjKbMHsuBTwJ/qMFjPxFjPNXl9j54zDzG/sJXgf+JMfYH/qf6/rEcqj6+To0xXpy5eOmq4bFyI/BOjLEf8EPgPzKbMrvU4vvrF0cdU7MyGjK7/ASY9D6fPx/oX/02BZiWgUzZ7Ce8/3gBvHzUsfVvGciUzSqAL8YYBwNjgP99jO9Hj7EGxhJyDDHG52OMFdV35wA9jvGwUcCqGOOaGOMR4GHgkkxlzCYxxtdjjCvTzlGf1HDMPMb+7BLgp9W3fwpcmmKWbFSTY+XoMXwUOCeEEDKYMdv4/VULMcY/ALve5yGXAD+LiTlAmxBC18ykyz41GC8dJca4Oca4qPr2PuB1oPtfPcxjrIGxhHywG4BnjvHx7sDGo+5v4m+/YfSXIvB8CGFhCGFK2mHqAY+xP+scY9xcfXsL0Pk4j2sWQlgQQpgTQmhMRaUmx8p7j6l+kWUP0D4j6bJTTb+/Lq+e+vFoCKFnZqLVS/68qr2PhRCWhBCeCSEMSTtMtqieKnoaMPevPuUx1sDkpR0gLSGE3wJdjvGpr8cYn6h+zNdJThH+PJPZslFNxqsGxsUYS0MInYAXQghvVL9a1CCdoDFrNN5vvI6+E2OMIYTjLevXq/oY6wu8GEJYFmNcfaKzqtF4CngoxlgWQphKciZpfMqZ1DAsIvl5tT+EMBn4Nck0o0YthNAS+BVwW4xxb9p5VLcabQmJMU54v8+HEK4DLgTOicdex7gUOPpVsR7VH2uQPmi8avg1SqvfbwshPE4yHaLBlpATMGYeY9VCCFtDCF1jjJurT79vO87XePcYWxNC+D3Jq2mNoYTU5Fh59zGbQgh5QGtgZ2biZaUPHLMY49HjMwv4XgZy1VeN6ufVR3X0H9gxxtkhhP8OIXSIMe5IM1eaQgj5JAXk5zHGx47xEI+xBsbpWMcQQpgEfAW4OMZ48DgPmw/0DyH0CSE0Aa4BGtVqPLURQigIIRS+exs4l+TibB2fx9ifPQl8rvr254C/OZMUQmgbQmhafbsDMBZYkbGE6arJsXL0GF4BvHicF1gaiw8cs7+ab34xyTx1HduTwP+qXsFoDLDnqCmU+ishhC7vXpMVQhhF8vdYo31RoHos7gFejzH+53Ee5jHWwDTaMyEf4E6gKcmUIYA5McZbQgjdgFkxxskxxooQwt8BzwG5wL0xxtfSi5yeEMJlwH8BHYGnQwiLY4znHT1eJHP4H68ezzzgwRjjs6mFTllNxsxj7C98F3gkhHAjsB64CiAkyxvfEmO8CTgZmBFCqCL5hf7dGGOjKCHHO1ZCCP8GLIgxPknyC/7+EMIqkgtmr0kvcfpqOGZ/H0K4mGRa7i7gutQCpyyE8BBwNtAhhLAJ+FcgHyDGOB2YDUwGVgEHgevTSZodajBeVwC3hhAqgEPANY38RYGxwGeBZSGExdUf+xpQBB5jDZU7pkuSJEnKKKdjSZIkScooS4gkSZKkjLKESJIkScooS4gkSZKkjLKESJIkScooS4gkNUAhhJ4hhLUhhHbV99tW3++dbjJJkiwhktQgxRg3AtNI9lih+v3MGOO61EJJklTNfUIkqYEKIeQDC4F7gZuBU2OM5emmkiTJHdMlqcGKMZaHEL4MPAucawGRJGULp2NJUsN2PrAZGJp2EEmS3mUJkaQGKoRwKjARGAP8Qwiha8qRJEkCLCGS1CCFEALJhem3xRg3AHcA3083lSRJCUuIJDVMNwMbYowvVN//b+DkEMJZKWaSJAlwdSxJkiRJGeaZEEmSJEkZZQmRJEmSlFGWEEmSJEkZZQmRJEmSlFGWEEmSJEkZZQmRJEmSlFGWEEmSJEkZZQmRJEmSlFH/HxEVczYJzFdLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,8))\n",
    "for i in range(3):\n",
    "    coefs = [0]*(i+1)\n",
    "    coefs[i] = 1\n",
    "    polynomial = make_polynomial(coefs)\n",
    "    x = np.linspace(-2, 2, 50)\n",
    "    plt.plot(x, [polynomial(i) for i in x], label='Polynomial degree {}'.format(i))\n",
    "fig.legend('center right')\n",
    "plt.axvline(0, color='#11111155')\n",
    "plt.axhline(0, color='#11111155')\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integration Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_rect(f, samples):\n",
    "    v = 0\n",
    "    for i in range(len(samples)-1):\n",
    "        mid = (samples[i] + samples[i+1])/2\n",
    "        v += (samples[i+1] - samples[i])*f(mid)\n",
    "    return v\n",
    "\n",
    "def integrate_trapez(f, samples):\n",
    "    v = 0\n",
    "    for i in range(len(samples)-1):\n",
    "        v += (f(samples[i]) + f(samples[i+1]))*(samples[i+1]-samples[i])/2\n",
    "    return v\n",
    "\n",
    "def simpson(f, samples):\n",
    "    v = 0\n",
    "    for i in range(len(samples)-1):\n",
    "        mid = (samples[i]+samples[i+1])/2\n",
    "        v += (f(samples[i]) + 4*f(mid) + f(samples[i+1]))*(samples[i+1]-samples[i])/6\n",
    "    return v\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate_rect(make_polynomial([0, 1]), [0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate_trapez(make_polynomial([0, 1]), [0,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integration error measures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadrature_error(ground_truth, estimated):\n",
    "    return abs(ground_truth - estimated)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_function_integration(f, a=0, b=1, points=2):\n",
    "    '''draw a function graph and a representation of aproximations of it's integral'''\n",
    "    \n",
    "    total_x = np.linspace(-2, 2, 50)\n",
    "    integration_x = np.linspace(a, b, points)\n",
    "\n",
    "    fig, axs = plt.subplots(1, 2,figsize=(18, 8))\n",
    "    rectangles = []\n",
    "    trapezes = []\n",
    "    # compute rectangles position\n",
    "    for i in range(len(integration_x)-1):\n",
    "        x1, x2 = integration_x[i], integration_x[i+1] \n",
    "        mid = (x1 + x2)/2\n",
    "        verts = [(x1, 0), (x1, f(mid)), (x2, f(mid)), (x2, 0)]\n",
    "        rectangles.append(Polygon(verts, facecolor='#eff3c6', edgecolor='#eff3c6'))\n",
    "    \n",
    "    # compute trapezes position\n",
    "    for i in range(len(integration_x)-1):\n",
    "        x1, x2 = integration_x[i], integration_x[i+1] \n",
    "        verts = [(x1, 0), (x1, f(x1)), (x2, f(x2)), (x2, 0)]\n",
    "        trapezes.append(Polygon(verts, facecolor='red', edgecolor='red'))\n",
    "    # drawing commands\n",
    "    axs[0].plot(total_x, [f(i) for i in x], label='Rectangles')\n",
    "    axs[0].add_collection(PatchCollection(rectangles, alpha=0.2))\n",
    "    \n",
    "    axs[1].plot(total_x, [f(i) for i in x], label='Trapezes')\n",
    "    axs[1].add_collection(PatchCollection(trapezes, alpha=0.2))\n",
    "    \n",
    "    fig.legend('center right')\n",
    "    # show coordinates system axis\n",
    "    for ax in axs:\n",
    "        ax.axvline(0, color='#11111155')\n",
    "        ax.axhline(0, color='#11111155')\n",
    "    plt.xlabel('X')\n",
    "    plt.ylabel('Y')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2cfa6ab67eb94938bdebc10a1282f999",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Text(value='0, 1', description='coefficients'), FloatSlider(value=0.0, description='a', …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact_manual(a=(-2, 2, 0.1),\n",
    "                b = (-2, 2, 0.1),\n",
    "                points=(2, 50, 1))\n",
    "def draw_line_integration(coefficients=\"0, 1\", a=0, b=1, points=2):\n",
    "    '''draw a line and a representation of integration aproximations'''\n",
    "\n",
    "    coefs = [float(c) for c in coefficients.split(',')]\n",
    "    line = make_polynomial(coefs)\n",
    "    \n",
    "    print(\"Ground Truth\", integrate_polynomial(coefs, a, b))\n",
    "    samples = np.linspace(a, b, points)\n",
    "    print(\"Rectangle\", integrate_rect(line, samples))\n",
    "    print(\"Trapeze\", integrate_trapez(line, samples))\n",
    "    \n",
    "    draw_function_integration(line, a, b, points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "78e5dc76e084443292cbc06ca28f2091",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Text(value='1, 2, -1', description='coefficients'), FloatSlider(value=0.0, description='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact_manual(a=(-2, 2, 0.1),\n",
    "                b = (-2, 2, 0.1),\n",
    "                points=(2, 50, 1))\n",
    "def draw_parabola_integration(coefficients=\"1, 2, -1\", a=0, b=1, points=2):\n",
    "    '''draw a line and a representation of integration aproximations'''\n",
    "    \n",
    "    coefs = [float(c) for c in coefficients.split(',')]\n",
    "    parabola = make_polynomial(coefs)\n",
    "    \n",
    "    print(\"Ground Truth\", integrate_polynomial(coefs, a, b))\n",
    "    samples = np.linspace(a, b, points)\n",
    "    print(\"Rectangle\", integrate_rect(parabola, samples))\n",
    "    print(\"Trapeze\", integrate_trapez(parabola, samples))\n",
    "    \n",
    "    draw_function_integration(parabola, a, b, points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error comparinson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_polynomial_integration_error(a=0, b=1, maxpoints=10, coefficients=\"0, 1\"):\n",
    "    coefs = [float(c) for c in coefficients.split(',')]\n",
    "    polynomial = make_polynomial(coefs)\n",
    "    \n",
    "    ground_truth = integrate_polynomial(coefs, a, b)\n",
    "\n",
    "    errors_rect = []\n",
    "    errors_trapz = []\n",
    "    \n",
    "    samples = range(2, maxpoints+1)\n",
    "    for points in samples:\n",
    "        x = np.linspace(a, b, points)\n",
    "        errors_rect.append(quadrature_error(ground_truth, integrate_rect(polynomial, x)))\n",
    "        errors_trapz.append(quadrature_error(ground_truth, integrate_trapez(polynomial, x)))\n",
    "\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(16, 16))\n",
    "    total_x = np.linspace(-2, 2, 50)\n",
    "    axs[0, 0].plot(total_x, polynomial(total_x))\n",
    "    #print(errors_rect)\n",
    "    axs[0, 1].plot(samples, errors_rect, '-', linewidth=4, label='Rectangle') \n",
    "    axs[0, 1].plot(samples, errors_trapz, '-o', label='Trapeze', color='red')\n",
    "    axs[1, 0].plot(np.abs(errors_rect), np.abs(errors_trapz), '-o', label='Error Rect x Trapeze', color='green')\n",
    "    axs[1, 0].set_xlabel(\"Error Midpoint rule\")\n",
    "    axs[1, 0].set_ylabel(\"Error Trapeze rule\")\n",
    "    fig.legend()\n",
    "    plt.xlabel('Number of Points')\n",
    "    plt.ylabel('Quadrature Error')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "20cd901b53e94692a4c7b55d712b6c20",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='a', max=2.0, min=-2.0), FloatSlider(value=1.0, descr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact_manual(a=(-2, 2, 0.1),\n",
    "                b=(-2, 2, 0.1))\n",
    "def error_on_line(a=0, b=1, maxpoints=10, coefficients=\"0, 1\"):\n",
    "    draw_polynomial_integration_error(a, b, maxpoints, coefficients)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "28c7d793a7c04389aaf2b64ec357264e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='a', max=2.0, min=-2.0), FloatSlider(value=1.0, descr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact_manual(a=(-2, 2, 0.1),\n",
    "                b=(-2, 2, 0.1),\n",
    "                maxpoints=(5, 50, 1))\n",
    "def error_on_parabola(a=0, \n",
    "                  b=1, \n",
    "                  maxpoints=10,\n",
    "                     coefficients=\"0, 0, 1\"):\n",
    "    draw_polynomial_integration_error(a, b, maxpoints, coefficients)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing Composite Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import sin, cos, tan, exp, pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_simpson(a, b, m, M4):\n",
    "    return M4 * (b-a)**5 /(2880 * m**4)\n",
    "\n",
    "def error_trapezoidal(a, b, m, M2):\n",
    "    return M2 * (b-a)**3 /(12 * m**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Case 2: integrate sin(x) between 0 and pi/2\n",
    "\n",
    "def solve(f, ground_truth, a, b, minpoints=2, maxpoints=100, M2=1, M4=1):\n",
    "\n",
    "    errors_simp = []\n",
    "    expected_error_simp = []\n",
    "    errors_trapz = []\n",
    "    expected_error_trapz = []\n",
    "    \n",
    "    samples = range(minpoints, maxpoints+1)\n",
    "    for points in samples:\n",
    "        x = np.linspace(a, b, points)\n",
    "        \n",
    "        errors_simp.append(quadrature_error(ground_truth, simpson(f, x)))\n",
    "        expected_error_simp.append(error_simpson(a, b, points, M4))\n",
    "        \n",
    "        errors_trapz.append(quadrature_error(ground_truth, integrate_trapez(f, x)))\n",
    "        expected_error_trapz.append(error_trapezoidal(a, b, points, M2))\n",
    "\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(16, 16))\n",
    "    axs[0, 0].plot(samples, errors_trapz, label='Error Trapezoidal', color='red')\n",
    "    axs[0, 0].plot(samples, expected_error_trapz, '--', linewidth=2, label='Expected error T', color='purple')\n",
    "    \n",
    "    axs[0, 1].plot(samples, errors_simp, label='Error Simpson', color='blue')\n",
    "    axs[0, 1].plot(samples, expected_error_simp, '--', linewidth=2, label='Expected error S', color='green')\n",
    "    \n",
    "    axs[1, 0].plot(samples, errors_simp, label='Error Simpson', color='blue')\n",
    "    axs[1, 0].plot(samples, errors_trapz, label='Error Trapezoidal', color='red')\n",
    "#     axs[0, 0].plot(total_x, polynomial(total_x))\n",
    "#     #print(errors_rect)\n",
    "#     axs[0, 1].plot(samples, errors_simp, '-', linewidth=4, label='Simpson') \n",
    "#     axs[0, 1].plot(samples, errors_trapz, '-o', label='Trapeze', color='red')\n",
    "#     axs[1, 0].plot(np.abs(errors_simp), np.abs(errors_trapz), '-o', label='Error Simpson x Trapeze', color='green')\n",
    "#     axs[1, 0].set_xlabel(\"Error Simpson rule\")\n",
    "#     axs[1, 0].set_ylabel(\"Error Trapeze rule\")\n",
    "    fig.legend()\n",
    "    plt.xlabel('Number of Points')\n",
    "    plt.ylabel('Quadrature Error')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f7e20aa5ca2e4ab49cdfad59c4e8fe18",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=3, description='minpoints', max=10, min=2), IntSlider(value=20, descript…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact_manual(minpoints=(2, 10, 1),\n",
    "                maxpoints=(5, 50, 1))\n",
    "def integrate_sin(minpoints=3, maxpoints=20):\n",
    "    solve(sin, 1, 0, pi/2, minpoints=minpoints, maxpoints=maxpoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "957948b8d46d41948752eed234126cba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=3, description='minpoints', max=10, min=2), IntSlider(value=20, descript…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact_manual(minpoints=(2, 10, 1),\n",
    "                maxpoints=(5, 50, 1))\n",
    "def integrate_exp2(minpoints=3, maxpoints=20):\n",
    "    def fc(x):\n",
    "        return exp(-(x**2))\n",
    "    solve(fc, 1.49365, -1, 1, minpoints=3, maxpoints=20, M2=2, M4=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "74f6a2c739e64935a2eb0a68bbcd5393",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=3, description='minpoints', max=10, min=2), IntSlider(value=20, descript…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact_manual(minpoints=(2, 10, 1),\n",
    "                maxpoints=(5, 50, 1))\n",
    "def integrate_tan(minpoints=3, maxpoints=20):\n",
    "    # f: tan(x)\n",
    "    # f(2): 2 sec^2(x) tan(x)\n",
    "    # f(4): 2 sec^2(x) (sec^2(x) + 2 tan^2(x))\n",
    "    solve(tan, 0, -pi/4, pi/4, minpoints=3, maxpoints=20, M2=4, M4=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.7 64-bit ('venv': venv)",
   "language": "python",
   "name": "python37764bitvenvvenv19f441ca6ee64d0c9084fa56e706e572"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}