{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Robert Johansson\n",
    "\n",
    "Source code listings for [Numerical Python - Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib](https://www.apress.com/us/book/9781484242452) (ISBN 978-1-484242-45-2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sympy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, x0 = sympy.symbols(\"x, x_0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = (x+0.5) ** 3 - 3.5 * (x+0.5) **2 + x  + 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{0} + \\left(x - x_{0}\\right) \\left(- 7.0 x_{0} + 3 \\left(x_{0} + 0.5\\right)^{2} - 2.5\\right) + \\left(x_{0} + 0.5\\right)^{3} - 3.5 \\left(x_{0} + 0.5\\right)^{2} + 3$"
      ],
      "text/plain": [
       "x_0 + (x - x_0)*(-7.0*x_0 + 3*(x_0 + 0.5)**2 - 2.5) + (x_0 + 0.5)**3 - 3.5*(x_0 + 0.5)**2 + 3"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f.subs(x,x0) + f.diff(x).subs(x,x0) * (x - x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f_func = sympy.lambdify(x, f, 'numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def arrowify(fig, ax):\n",
    "    xmin, xmax = ax.get_xlim()\n",
    "    ymin, ymax = ax.get_ylim()\n",
    "\n",
    "    # removing the default axis on all sides:\n",
    "    for side in ['bottom','right','top','left']:\n",
    "        ax.spines[side].set_visible(False)\n",
    "\n",
    "    # removing the axis labels and ticks\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.xaxis.set_ticks_position('none')\n",
    "    ax.yaxis.set_ticks_position('none')\n",
    "\n",
    "    # wider figure for demonstration\n",
    "    #fig.set_size_inches(4,2.2)\n",
    "\n",
    "    # get width and height of axes object to compute\n",
    "    # matching arrowhead length and width\n",
    "    dps = fig.dpi_scale_trans.inverted()\n",
    "    bbox = ax.get_window_extent().transformed(dps)\n",
    "    width, height = bbox.width, bbox.height\n",
    "\n",
    "    # manual arrowhead width and length\n",
    "    hw = 1./25.*(ymax-ymin)\n",
    "    hl = 1./25.*(xmax-xmin)\n",
    "    lw = 0.5 # axis line width\n",
    "    ohg = 0.3 # arrow overhang\n",
    "\n",
    "    # compute matching arrowhead length and width\n",
    "    yhw = hw/(ymax-ymin)*(xmax-xmin)* height/width\n",
    "    yhl = hl/(xmax-xmin)*(ymax-ymin)* width/height\n",
    "\n",
    "    # draw x and y axis\n",
    "    ax.arrow(xmin, 0, xmax-xmin, 0., fc='k', ec='k', lw = lw,\n",
    "             head_width=hw, head_length=hl, overhang=ohg,\n",
    "             length_includes_head=True, clip_on=False)\n",
    "\n",
    "    ax.arrow(0, ymin, 0., ymax-ymin, fc='k', ec='k', lw = lw,\n",
    "             head_width=yhw, head_length=yhl, overhang=ohg,\n",
    "             length_includes_head=True, clip_on=False)\n",
    "\n",
    "    \n",
    "    ax.text(xmax*0.95, ymin*0.25, r'$x$', fontsize=22)\n",
    "    ax.text(xmin*0.35, ymax*0.9, r'$f(x)$', fontsize=22)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAADyCAYAAAAWeS+KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XlYFfX7//HnsC+KLIoIKijgviLiUqa471upqbmUmqZlZprmmqWmZX1yX8s0TdwtK1PE3L6uqYA7ZQiCK4rs+5nfHwY/TVSUA3M43I/r8upq5pyZ+yDy4n7Pe96jqKqKEEIIIV6cidYFCCGEEEWdhKkQQgiRTxKmQgghRD5JmAohhBD5JGEqhBBC5JOEqRBCCJFPEqZCCCFEPkmYCiGEEPkkYVqMKYrSXlGUvYqixCiKkqUoiqooykhFURoriqJTFGWOHs5hoyjKDUVRTiqKouijbiGEMDSKrIBUPCmK0g74DcgC/gBuAiowF1gDVAEqqaoaq4dzvQssBAarqromv8cTQghDI2FaTCmKcgRoAvRVVTXgoe39gPXAZ6qqTtPTuSyACEDHg4BO18dxhRDCUMgwbzGkKIoD0Bi4D2z+z+4xPOhQv9PX+f4Nz3WAK9BLX8cVQghDIWFajCiKUlVRFBW4ByiAPZD577XSZEVRGgENgQOqql59wjFm//v6wFz2KYqirP93/2+Kopg/tDt7eHekXj/U4zVUKMjjCyFEbsye8/UyJlyEhYSE8PXXXxMSEkJwcDA+Pj7Url0bADc3N2vg2OzZs5k5c2YLnvB3HR8fj5eXF7dv3269d+9etXXr1jn7Ro0axeLFi2nWrBm7d+/uYG1tnTOcq6oqLi4u3L59u+mdO3fUMmXK6P3zqapKhQoV0Ol0mJjI74lCCL3I28RJVVWf548wAoMGDVIBdcOGDY9sb9q0qQqoQUFBT33/4sWLVUD19fXN2TZ16lQVUBs0aKDGxcXl+r7u3burgLpx48b8f4hcnD17Vi1ZsqQaGhpaIMcXQhRLecpH+fW9GDp16hQADRo0eGR7cHAwANWrV3/q+99++22qVavGn3/+yZYtW5g/fz6fffYZ1atX5/fff8fOzi7X99WoUQOAM2fO5Pcj5CogIICEhAQCAgKe/WIhhNCj5x3mFUVcSkoKFy9exM7ODi8vr5ztSUlJJCcnA+Dk5PTUY5iZmTF37ly6devGO++8w927d/Hw8CAwMJDSpUs/8X2Ojo4A3Lp1Sw+f5HHZvyRk/1cIIQqLdKbFTEhICFlZWdSrV4+H11CIi4sDwNLSEgsLi2cep2vXrtSsWZOYmBjKlCnD3r17cXNze+p7sjvW+/fv5+MT5E6n03H16lUAIiIi0Ol0ej+HEEI8iYRpMXP69GkAfHx8Htlub28PQFpaGmlpac88zoIFCzh//jwAqampTxzafVh8fDwADg4Oz1VzXpw7d47o6GgAoqOjCQ0N1fs5hBDiSSRMi5nsIdD/hqmNjQ22trYA3Lt376nHWLNmDWPGjMHNzY0uXboQHx/PjBkznnnuu3fvAuDs7PwipT/Vhg0bSEhIACAhIYGNGzfq/RxCCPEkEqbFzJM604e3Xbhw4Ynv3759O0OGDMHR0ZHAwEAWL16MlZUVy5cvJyws7Knnzj5ubufOr/9Oasr+nEIIURgkTIuRtLQ0zp8/j42NDdWqVXtsv7+/PwBHjx7N9f179+6lb9++2NjY8Pvvv1O9enUqVKjAu+++S2ZmJhMnTnzq+Y8dO4aiKLRo0SLfn+VhOp2OiIiIR7bJdVMhip+1a9dy584dTc4tYVqMnD17loyMDOrUqYOpqelj+7t37w48CM3/OnbsWM7+n376CV9f35x9H3/8MaVKlWL79u0cPnz4iee+desWTZo0Qd8LNoSEhBAVFfXItujoaEJCQvR6HiGE4Tp79iyDBg2iRo0aZGRkFPr5JUyLkSddL81Wv359GjduzMGDB3NmxsKDb9KOHTuSlpbGxo0bczrYbI6OjkyYMAGAcePG5XrsNWserCY4cqT+VxMMCAggMTHxkW2JiYlyv6kQxcjy5csB6N27N+bm5s94tf4971NjZDlBIxcQEEDfvn2ZOnUqn376qV6OmZ6ejru7OyYmJoSHh+fp1pvn0a5dO/bs2fPY9rZt27J79269nksIYXgSExNxdXUlISGB0NDQnGVS9SRPywlKZyoe0adPH/z8/Fi4cCGxsfl+lCkAK1as4ObNm8yePVvvQfrw9dLs30azzyHXTYUoHrJn87/00kv6DtI8kzAVj1AUhQULFhAXF8fcuXPzfbzk5GRmzZqFr68vAwcO1EOFjwoODs65vzR70f3s/0ZHRxfY0oVCCMOgqipLly4FYMSIEZrVIcO8okgbP3488+bNo0GDBhw8eBBbW1uSkpJ45ZVXOHXqFOPGjePLL7/UukwhRAE5efIkfn5+ODk5ERUVhZWVlb5PIcO8wviFhITg5eXFtm3bsLGxAR4sQLFt2za8vLxkRq8QRi67K33zzTcLIkjzTDpTUWRlZWXh7OzMli1bcmYYK4pC9vf0gQMH6NmzJ7dv3871ViAhRNEWGxuLm5sbKSkphIWF4e3tXRCnkc5UGLfk5GSWL1/+2K062Zo3b87y5ctznoYjhDAuP/zwAykpKbRu3bqggjTPpDMVRuXhzlQIYbxUVaVmzZpcvHiRrVu30rNnz4I6VZ46UwlTYVQkTIUoHvbv34+/vz/lypUjIiKiIBdqkGFeIUTuAgMDURSFyMhIrUsR4oUsXLgQgGHDhmmy4tF/SZgKUQwFBwfj5ORExYoVtS5FiOcWERHBjh07MDc31/Te0odJmApRDAUHBxfIo/CEKAxLlixBp9PRu3dvypUrp3U5gISpEMXSmTNnqFOnDjNnzqRixYo4ODjQv39/4uLitC5NiKdKTk5m5cqVALz33nsaV/P/mWldgBCicGXfkxcXF0fnzp1ZvXo1586dY+LEiZQoUSLn6RtCGKL169cTGxuLn58fjRo10rqcHBKmQhQzZ8+eJSsri+7du7N48WIAWrVqxd9//01AQICEqTBYqqqyYMECAEaPHq1xNY+SYV4hipkzZ85gaWnJlClTHtleq1Yt7t69K0/aEQbrwIEDnDt3DhcXF3r16qV1OY+QMBXCSCQmJubpsXnBwcHUq1fvsYkbt27dwtXVFRMT+bEgDFN2VzpixAi9P84xv+RfjRBF3M2bN+nTpw+lSpXC0dGRrl27EhMT88TXBwcH4+7u/tj2rVu30qZNm4IsVYgXdvXqVX766SfMzc0ZPny41uU8Rq6ZClGEpaen4+/vz6VLl3K27dy5k/79+7N79+7HXq/T6QgNDSUjI+OR7Vu3buXcuXOsXbu2wGsW4kUsXLgQnU5H3759cXFx0bqcx0hnKkQRtn379pwgXbZsGbt27QJgz549XLx48bHXh4WFkZyczN27dxk3bhx//PEHX3zxBQMHDuSTTz6hbt26hVq/EHlx//59VqxYAcDYsWM1riZ30pkKUYQdOXIEADMzMwYPHoylpSV169bFxMSEyMhIqlev/sjrg4ODMTU1JTAwkKFDh7Jo0SLc3d353//+x9tvv63FRxDimVauXEliYiItW7Y02MVGZKF7YVSK20L3vXv3ZvPmzZQvX55r165pXY4Qepeenk7lypWJjo7mt99+o0OHDoVdgix0L4Sxi4+PB6BUqVIaVyJEwdi4cSPR0dHUqFGD9u3ba13OE0mYClGEJSYmAlCyZEmNKxFC/1RVZd68eQCMGzcORclTk6gJCVMhirCEhASgYML0999/p3Xr1pQuXRpTU1MURWHJkiUAHDt2DBMTEyZOnJivcyQnJ1OuXDkaNmyYr+F5RVEM+geteDF79+4lNDQUFxcX+vXrp3U5TyUTkIQowgqqM929ezedOnXC1NQUf39/XFxcUBQFf39/VFVl9OjR2NnZMWHChHydx8bGhsmTJ/Pee++xdu1aBg0apKdPIIxBdlc6evRoLC0tNa7m6WQCkjAqxW0CkrOzM3fu3GHgwIGsWbNGb8dt2rQpR48eZcOGDbz++uuP7Pvxxx/p378/U6dO5dNPP833udLT03F3d8fExITw8PAXWtkmuystTn/3xi40NJS6detia2tLZGQkjo6OWpUiE5CEMHYF0ZnGxsZy7Ngx7O3tc13/9JtvvkFRFN566y29nM/CwoI33niD69evs3nzZr0cUxR9c+fOBeCtt97SMkjzTMJUiCIqKyuLlJQUQD9hevnyZRRFwdHREVVVuX//PmZmZiiKgo2NDTqdjpMnT3Ly5EmaN2+Oh4fHY8eYNGkSiqLkuiyhqqr0798fRVHo2LHjI6swZQ/vZl+Tzc3Zs2fp0aMHjo6O2Nra4uPjw6pVq576mV60HqGtK1euEBAQgJmZGePGjdO6nLxRVfV5/ghh0B58SxcP9+/fV3lw6UWdOXNmvo8XEhKiDho0SK1Xr54KqD4+PuqgQYPUQYMGqZMmTVJVVVUnTZr01PPFx8erzs7OKqAGBgY+sm/UqFEqoDZr1kxNTk5+7L1ly5ZVFUVRb9++/di+/fv3q9bW1iqgVq1aVX399dfV5s2bqyYmJuoHH3yQ83XQZz1CO8OGDVMB9c0339S6FFXNYz5KmAqjUpzC9Nq1azkhMn/+fL0dd9CgQSqgbtiw4bF9TZs2VQE1KCjoie9fvHixCqi+vr4526ZOnaoCaoMGDdS4uLhc39e9e3cVUDdu3PjI9uTkZNXNzU0F1I8//ljV6XQ5+/bv36/a2Ng8MUzzU4/QRlRUlGpubq4qiqJeunRJ63JUVcJUFEfFKUwvXryYEyLfffed3o5bq1YtFVDDwsIe25cdXNevX3/i+zMyMtRq1aqpgLp582b1m2++UQG1evXq6p07d574vuyud+LEiY9sX7t2rQqonp6eamZm5mPvGzt27FPD9EXrEdoYM2aMCqi9e/fWupRsecpHuWYqRBGVfY8p6G8CUkpKChcvXsTOzg4vL69H9iUlJZGcnAyAk5PTE49hZmaWM3nknXfe4YMPPsDDw4PAwEBKly79xPdlTzK5devWI9sPHDgAwOuvv46pqelj7xswYMBTP9OL1iMK3507d1i+fDnw4Hp3USJhKkQRlT2TF6BEiRJ6OWZISAhZWVnUq1fvsUUQ4uLiALC0tHzm7Stdu3alZs2axMTEUKZMGfbu3Yubm9tT32NnZwc8eELIw6KiogCoVKlSru/LbSKUPuoRhW/+/PmkpKTQqVOnIvcEIwlTIYqoguhMT58+DZDrkzns7e0BSEtLIy0t7anHWbBgAefPnwcgNTU1JyifJnudYQcHh+eqOS8rH71IPaJwxcXFsWjRIgAmT56scTXPT8JUiCLq4c5UX2F66tQpIPcwtbGxwdbWFoB79+498Rhr1qxhzJgxuLm50aVLF+Lj45kxY8Yzz3337l3gwUIUD8vuIK9evZrr+8LDw5963BetRxSuxYsXExcXh7+/P02aNNG6nOcmYSpEEfVwZ6qvYd6ndaYPb79w4UKu+7dv386QIUNwdHQkMDCQxYsXY2VlxfLlywkLC3vqubOP+d9zN2/eHICAgACysrIee9/69eufeMz81CMKT3x8fM7SgUWxKwVkNq8wLhSj2bzz5s3LmcWa272Zzys1NVU1NzdXbWxscp01q6qqOm3aNBVQP/vss8f2BQYGqpaWlmrJkiXVkydP5mwfN26cCqg9evR46vmfdJ9pUlKSWq5cORVQp0yZ8sitMYcOHVJtbW1znc2b33pE4fn0009VQH3llVce+fs1EHJrjCh+ilOYTp8+PSdEUlJS8n28kydPqoDauHHjJ77m9OnTKqA2b978ke1Hjx5VbW1tVUtLS3Xfvn2P7Lt7965aqlQpFVAPHTqU63FDQ0NVQG3atGmu+4OCglQrKysVUKtVq6b27dtXbdGixRMXbchvPaLwxMbG5vx97N+/X+tyciNhKoqf4hSmH374oQqoZmZmejnesmXLVEAdOXLkU1/XuHFjVVEUNTw8XFXVB0Ho4OCgmpmZqTt27Mj1PbNnz1YBtVGjRrnuz/4s69ate+J5g4OD1a5du6r29vaqtbW1WrduXXXp0qWqqqqPhKk+6hGFJ3u0o2XLllqX8iR5ykd5aowwKsXpqTHDhw9nxYoV2NvbExsbW2jnDQgIoG/fvgb11BhRNN27dw8PDw8SEhI4dOgQL7/8stYl5UaeGiOEMcu+77MgHgz+NH369MHPz4+FCxfqJcRXrFjBzZs3mT17tgRpMfPVV1+RkJBAmzZtDDVI80zCVIgiIDo6mqZNm+Lu7s61a9cAiIiIAKBcuXKFWouiKCxYsIC4uLiclYVeVHJyMrNmzcLX15eBAwfqqUJRFMTExDB//nwAo7hVScJUiCKgbNmyXL58mcjISFatWkVwcDB//vknAI0bN37m+7du3Urr1q0pW7Ys1tbW1KpVi4CAgBeup1GjRuh0OubMmfPCx4AH967euHGDkydP5mnxBWE8vvjiC5KSkujQoUORvK/0v+SaqTAqxnzNdOrUqcycOfORbSVKlCA4OBhPT8+nvnf69OmUK1cOT09PFEVhx44dLFu2jNDQUGrUqFGQZQvxmKioKLy9vUlNTeXEiRM0bNhQ65KeJk+/5UmYCqNizGGamZnJuHHjWLx4MZmZmTRu3JilS5dSr1695zpOVlYWqqpSsmRJlixZwptvvllAFQuRu6FDh/Ltt9/Sq1cvNm3apHU5zyJhKoofYw7TbCkpKSQlJeX5iSdpaWmsWLGC7777jn/++SdnDVyATZs20atXr4IqVYjHXLhwgdq1a6MoChcvXsTb21vrkp4lT2FqVtBVCCH0y9raGmtr6zy/vlu3bpw+fZoxY8bg4+ODk5MTBw8eZNy4cTLEKwrdpEmT0Ol0vPPOO0UhSPNMOlNhVIpDZ/o8jhw5wksvvcSePXto06ZNzvaRI0eyevVqEhMTc31GqBAFIfv70cbGhitXruDi4qJ1SXkh95kKUdxl30ZTrVq1nG0XLlxg9erV1KpVS4JUFBpVVZkwYQIAY8eOLSpBmmcSpkIYMR8fH0xNTRk7dixBQUEsWrSIrl27Ym1t/dwTl4TIj19++YXDhw/j5OTE+PHjtS5H7yRMhTBi3t7erFy5kqNHj9KtWzd+/fVXtm7dSkpKioSpKDQZGRk5XenUqVON8uHscs1UGBW5ZiqE4Vm4cCGjR4/G09OT8+fPY2lpqXVJz0NujREFIzU1lRs3bnD9+nWuX7/OjRs3iI+PJykpieTkZJKSklBVFVNTU8zMzDA1NcXW1hYHBwccHBxwdHTE2dmZihUr4ubmhrm5ud5qkzAVwrDcvXsXb29vYmNj2bFjB926ddO6pOclt8aI/NHpdFy+fJnjx49z4cKFnD9Xr17VW2CZmJjg6upKpUqVqFGjBjVr1qRmzZrUqlULZ2dnvZxDCKGdGTNmEBsbS8uWLenatavW5RQY6UxFDp1Ox59//smePXs4cuQIx44dy/WpICamJpRyKYV9OXvsXeyxK2uHVSkrLG0ssbC2wNLGEsVEQZelQ5elIzMjk/TkdJLvJ5Mcl0xKbArxt+K5G32XuBtxTwxmd3d3/Pz88PPzo1GjRvj5+T1zeEg6UyEMx4ULF6hTpw6qqhIcHEzt2rW1LulFyDCveLbExER27drFr7/+yq5du7h9+/Yj+0uVK4WnnyduNdxwqepCuSrlcKnsgrmlfoZmM9Mzib0ey60rt4i+HM2NSze4efkmUReiSEtMe+S11tbWvPTSS/j7+9OyZUsaNmz42K0dEqZCGAZVVenQoQO7d+9mxIgRLF26VOuSXpSEqchdZmYmQUFBrFu3jm3btpGcnJyzz7G8I7Xb1MazqSfejbxxdHPExKTwJ33rsnTcuHyDK6euEH46nH+O/8P1i9cfeY2TkxMdOnSgc+fOtGvXDnt7ewlTIQzEr7/+SufOnSlVqhR//fUXZcqU0bqkFyVhKh4VFRXFkiVLWL16NTdv3szZXrlhZep2rEutNrUoX728JuGZF3G34wj7vzAuHrrIpf2XuBN+J2efqakpLVu2JDAwkDt37uR53VohhP6lpaVRp04dwsLC+Oqrrxg7dqzWJeWHhKl4MNRy9OhR5s+fz9atW8nKygLA2dOZRr0b4dfLj7IeZYvcsyRVVeXmXzcJ3RNK6J5Q/j76N7pMHfAgWP39/enXrx+vvfYaJUuW1LhaIYqXmTNnMnXqVKpVq0ZISAgWFhZal5QfEqbFmaqqBAUFMX36dI4cOQI8mDjk09WHFsNa4N3I22A70BeReC+RM7+eYe3otZiYmeQEq42NDT179mTgwIG0bNlSls8TooD9888/1KxZk9TUVPbt24e/v7/WJeWXhGlxdeDAAaZNm8bBgwcBsHWwpdngZjR/qzlObk5Frgt9HsMchvG/K//j9M7THN14lL+P/p2zr0LFCgwbOowhQ4bg6uqqYZVCGCdVVencuTO//fYb/fv3Z926dVqXpA8SpsXN+fPn+eCDDwgMDATAxt6Gtu+1xX+YPzYlbTSurnAMcxjGytiVOf9/O/w2Rzce5VjAMWIiYoAHw8CdO3fmnXfeoU2bNkbVoQuhpW3btvHqq69SqlQpLl++TNmyZbUuSR8kTIuL+/fv88knn7Bo0SKysrKwtrOm9butafV2K2xL2WpdXqH6b5hm0+l0XNx/kQNrDhDyW0jOMHDValUZ/d5oBg4cSIkSJQq7XCGMRmJiItWrVycqKorFixczcuRIrUvSFwlTY6eqKqtXr2bChAnExMSgmCi8MvgVunzchVKlS2ldniaeFKYPu3/zPofXH+bgdweJvf5gUQq7UnYMGzqM999/nwoVKhRGqUIYlfHjxzNv3jx8fX05duyYMc1PkDA1ZhEREQwdOpS9e/cC4N3Um95zeuNey92or4k+S17CNFtmRiZnfjlD0PIgrhy/AoCZmRl9+/Zl3Lhx1KlTpyBLFcJonD59Gj8/P3Q6HSdOnMDX11frkvRJwtQYqarKihUrGDduHImJiZR0LMlrc16j8auN5dofzxemDws/Hc6exXs4teMUqu7Bt3m79u2YMnkKL7/8sr7LFMJoZGRk0LBhQ0JCQnj//ff55ptvtC5J3yRMjc2NGzcYOHBgTjfas2dP/Of4Y1HaAhNFghRePEyz3bl6h8Clgfzfuv8jPTkdgFeav8LUKVNp1apVse76hchN9j2lHpU9OBd6Dltbo5unkad/9PITuIgIDAykXr167N27l9KlS7Nx40a2bNlCKefieW20oJTxKEO/uf2Ye3YuncZ3wqaUDQcPHKRNmzY0btKY33//XZYrFOJf58+f59NPPwXg1f+9aoxBmmcSpgYuMzOTKVOm0K5dO27fvk3Lli05e/YsvXv3li6pAJVwLEH3Sd35PORzuk/tTgmnEpw4foIOHTrQpGkT9uzZI6EqirWsrCzeeustMjIyaDKoCc38m2ldkqYkTA3Y7du3ad26NbNmzUJRFGbMmMGePXtwcXHRurRiw6aUDZ3GdmJOyBxenfEqJZxKcPzYcdq1a8fLzV5m//79WpcohCa++eYbTpw4gb2rPV1ndMVUMZrZuy9EwtRAhYaG0rBhQw4cOICLiwtBQUFMmzbNmKabFymWtpa0H92ez4M/p8f0Htg62nLk/47g7+9P23Zt+fPPP7UuUYhCc+nSJaZMmQJAn6/7YG1nrXFF2pMwNUA//fQTTZs2JTIykkaNGnHmzBlatGihdVkCsCphRccxHZkTPIeuk7piVdKKwD2BNGzYkJ6v9uTSpUtalyhEgUpPT6d///6kpqbS6PVGVG9bXeuSDIKEqQFRVZW5c+fSo0cPkpKS6N+/P/v375dhXQNkVdKKLuO78Hnw57Qd3RZzK3O2b9tOrVq1ePvtt7l+/fqzDyJEETRjxgxOnz6NU0UneszpgZK3ya5GT8LUQGRlZTFy5EgmTpyIqqrMnj2bH374ASsrK61LE09RwrEEvWb0Yvbp2TQb3AwVlZUrV+Lp5cmUKVOIj4/XukQh9ObQoUN8/vnnKCYKbyx9Q4Z3HyJhagDS0tLo168fy5Ytw9LSkq1bt/Lxxx/LbN0ixL6cPQP/N5Dp/zedep3rkZqSyqxZs/D08mTJkiVkZGRoXaIQ+RIXF8eAAQNQVZW2Y9pSqUklrUsyKBKmGktMTKRLly5s2rSJkiVLsnv3bnr27Kl1WeIFuVZ1ZdQPo5jw+wQ8G3kScyeGUaNGUbN2TX766Se5nUYUWaNHjyYiIoIK9SrQ5qM2Mrz7HxKmGrp37x6tWrUiMDAQZ2dnDhw4QPPmzbUuS+iBVyMvJuyawIg1IyhTuQx/Xf6L7t2749/Sn+DgYK3LE+K5rF+/nrVr12JhbcGAZQMwtzDXuiSDI2GqkXv37tG6dWtOnDiBu7s7hw8fpn79+lqXJfRIURQadG3Ap0c/pc/nfbB1sOXA/gP4+Pjw1ltvcePGDa1LFOKZLl68yPDhwwHoOasnzlWcNa7IMEmYaiA2Npa2bdty5swZvLy8OHz4MN7e3lqXJQqImYUZrUe0ZtbpWbR6pxWKqcLq1avx8vZi9uzZpKamal2iELlKSkqiV69eJCUl0bBXQxoNaiTDu08gYVrI7t+/T9u2bTl16hSenp788ccflC9fXuuyRCGwtbfl9dmvM+PoDOp2rEtyUjKTJ0+mavWqbNmyRa6nCoPz7rvvcv78ecp6l+W1r16TB2o8hXxlClFcXBzt2rXjzz//pHLlyhKkxZSLlwvvrn+XsTvG4lbDjcirkfTq1YsW/i0ICQnRujwhAFi9ejXff/89FtYWDF49GKsScpve00iYFpLU1FS6devGiRMnqFSpEn/88QcVKlTQuiyhoerNqzP1wFT6zeuHraMtBw8cxMfHh+HDhxMTE6N1eaIYCwkJYdSoUQD0/rI35WqU07giwydhWgiysrLo168fBw4cwNXVlX379lGxYkWtyxIGwNTMFP8h/sw6NYuWw1uCAitWrMDL24v58+fL/ami0N25c4du3bqRkpJC436N8e3nK9dJ80DCtICpqsrIkSPZvn079vb27N69Gw8PD63LEgbG1t6WvnP6Mu3QNKr7Vyfufhw0FhaFAAAa4klEQVRjxoyhdt3aOQ+DF6Kgpaen89prrxEREYFHAw9enfeqBGkeSZgWsBkzZrBixQqsrKz4+eefqVWrltYlCQPmVt2ND7Z+wMj1IylTqQyXL16mTZs29OjZg/DwcK3LE0Zu9OjRHDx4kFLlSvHm2jexsLLQuqQiQ8K0AK1YsYIZM2ZgYmJCQEAAzZoV74fnirxRFIX6Hesz48gMuk/tjqWtJTu276Ba9WpMnTqV5ORkrUsURmjp0qUsX74cc0tzhv4wFPty9lqXVKRImBaQffv2MXLkSACWL19Ot27dNK5IFDXmVuZ0GtuJz058RqPejUhPS2fmzJlUqVaFTZs2ya00Qm+CgoIYPXo0AH3n96WCj0yOfF4SpgXgr7/+4rXXXiMrK4sJEyYwdOhQrUsSRZiDqwNDlw/lo98+okLdCkRfi6ZPnz40929OaGio1uWJIi40NJSePXuSmZlJ6/da49PbR66TvgAJUz2LjY2lc+fOxMbG0rVrV2bPnq11ScJIeDfxZkrQFAb8bwAlnEpw6MAh6tevz8iRI7l7967W5YkiKDIykg4dOhAfH49Pdx86Tu8oQfqCJEz1KCMjg969exMWFkadOnVYt24dJibyJRb6Y2JqwiuDX2HmnzPxf9sflAfXuryreLNkyRIyMzO1LlEUEbGxsXTo0IHr16/j1dSLvkv6YmpiqnVZRZb8pNejcePGsXfvXpydnfn5558pWbKk1iUJI2Vrb0u/uf2YemAqVV+pSuy9WEaNGkX9BvXZv3+/1uUJA5eWlkb37t25cOECLlVdeGvdW0Y3c/eXX34hMjKy0M4nYaonAQEBLFiwAHNzc7Zv3467u7vWJYlioHzN8ny440OGfz8cxwqOnAs9h7+/P6/1eo2rV69qXZ4wQBkZGfTp0yfnFpjhm4dja2+rdVl6d/nyZerUqUPdunXp0qULc+bMISwsrMAm7inPeWCZPpiLCxcu4OfnR1JSEosXL86ZxVsYfoj7gRRdiixA/a9hDsNYGbtS6zI0kZ6Szu6Fu9n1zS4yUjKwtLJk/LjxTJw4EVtb4/thKZ5fZmYm/fv3Z9OmTdjY2/Duzndxq+mW7+OqqHiae9K5RGc9VKk/Bw4c4O233yYsLAwABwcHXF1dqVixIr6+vvTu3ZsaNWo863Jcni4iS5jmU0JCAn5+fly6dIl+/fqxbt06FKXwLuBLmD6qOIdptntR99jyyRZObj0JQDm3cnwx5wv69esn1/CLMZ1Ox+DBg/nhhx+wtrNm5LaRVPTRz7KmhhqmANevX6dnz54cP378sX12dnY54Vq3bl169+5N/fr1MTV95NqxfsNUURTbK1euJObpxcWEqqq8//77/Prrr3h7e7Nt2zZsbGwKtYbtCdtJ1aVKmP5rUv1JzD4jM6gBwk+H8+u8X7l+8ToAderWYdrUafIQ+mJIVVUmT57Mxo0bsbC2oP+S/pSvq78nVqmouJu509K2pd6OqU9paWmMGDGCkydPkpKS8sTX2djYUL58eSpUqECNGjXo1asXr7zyirmqqs+c2fc8Ydq2UqVKu/NevvGLi4vj3r17KIqCq6srFhaFfwE/UZeIDl2hn9dQ3Y24i5O7k9ZlGAxVVUlLSiM5NhlV9+DfupWFGaWtrDGXmZvFgg6VWxmppCalgQIly5TU+2QjFRVzxRwbk8JtJp5HeHg4JUqUIDHx2T1huXLlKF26NO7u7vzyyy+lVVV95r1nMsz7goKDg2nUqBHp6els2rSJXr16aVJHki6JVDVVk3MbotJmpYnJlMeX/VdCQgJzPp/Dt/O/JTMtEytTU8b6+jLRz4+SGvwSKApHRlYW/fbuYsu5y5hbmTNxw0Te6/ye3s/Tvml7EuMTuXzpst6PnV+qqjJ9+nSWLVvGnTt3HttvYmKCq6srbm5uVKpUiQ4dOtCxY0dKly6d/RK5ZlpQkpOTadCgAZcuXWLEiBEsXbpU65LEvxRFkWX2niI8PJxx741g2697AChrY8OnLVvyVv36mMn1VKOSmplJ782b2RkWhqWtJV9s+YL32r2n9zkdWVlZ2NnZ0blzZzZu3KjXY+dXamoq/fv3Z9euXTnDu2ZmZri5ueHm5oanpyddu3aldevW2Ns/cS3iPH3BzPRUc7Hy4YcfcunSJapXr85XX32ldTlC5FmlSpXY+stujh4+zAcjRnD8/HmG//IL3xw9yhdt29LJ27tQJ9CJgpGYnk63gAD2hYfjaGHBd0tm07UAghQe3IKSnJxM3bp19X7s/IiKiqJz585cvHiR8uXLU758eby8vOjZsyctWrTQ+wx36Uyf008//UT37t2xsLDg+PHj1KtXT+uSxEOkM807VVXZtG4dH48fT/itWwC0cHfny7Zt8XV11bg68aKuJyTQdcMGTt24gYu1NYErVlDrjTcK7Hw//vgj/fv355dffqFTp04Fdp7nNXToUExMTOjTpw8vvfQSVlZWL3qoPP0GIuM6z+H69esMGTIEgDlz5mgSpPfu3WPatGk0btyYMmXKYGNjQ7Vq1Zg7dy46nUxEEnmnKAp9BgzgYkQEX8+ciYOtLfsjImi4ciW9N2/mL1nvt8gJvXWLRqtWcerGDSqXLMmhzZsLNEgBzpw5A0DFihUZPXo0bm5u2NnZ0apVK86fP1+g536aVatWsWLFClq1apWfIM0z6UzzSFVV2rVrR2BgIG3btmXXrl2a3LMXEBDAjBkz6NSpE56enqSnp7Nx40aOHj3K3Llz+eijjwq9JkMinemLi42N5fPJk1m4ahWpGRmYKgrDfHyY2rw5rrI0psH77a+/6LNlC4np6TR1cWHHzp2U8fUt8PO2adOGgwcP4u3tTd26dWnRogVhYWEsWLAAe3t7/vrrL+zs7Aq8jgIkE5D0afny5YwYMQInJyfOnj1LuXLlNKkjKSnpsbH+jIwMqlWrRrly5Th8+LAmdRkKCdP8i7p2jU8+/JDVW7agU1WsTE1518+PCS+/TOlCvo9aPJuqqiw8cYIPdu9Gp6r0rVGD73btwqqifhZkeJYyZcoQExPDqlWrckbuAObPn8+YMWMe214EyTCvvly9epVx48YBsHjxYs2CFMgJUlVViY+PJyYmhri4OJydnUlLS9OsLmE8yleowKpNmzh37hw9W7UiNSuLeUePUnn+fKb/8QdxqXIrlqFITE+n37ZtvP/77+hUlWlt2rD++PFCC9Jr164RExODn5/fY4HZvn17AK5cuVIotWhNwvQZVFVl6NChJCYm8uqrr9K7d29N69m0aVPOTLRSpUpRpkwZypQpw7Fjx/D29ta0NmFcqteowda9e/nz+HHaN2lCQno6nx48iMc33/DpgQMSqhq7cOcODVeuJODcOUqYmbFxzBhm/PYbSokShVZDcHAwAGPGjHlsX/bM4RKFWI+WJEyfYcWKFQQFBeHk5MSSJUs0vW3go48+ok+fPtja2vLVV1+xc+dOAgMDWbZsGYAsEycKRAM/P3YdOcLBP/6gef363E9LY/r+/Tmhel9CtVCpqsq60FD8Vq7kUkwMNR0cOLlhA73/9z8wK9y7HbMnH/n5+T22L3st3AYNGhRqTVqRa6ZPcfXqVWrXrk1iYiIBAQH06dNHs1qioqKoWLEiffv2Zf369Y/smzJlCrNmzWLv3r20atVKowoNg1wzLXj79+7lk/HjOfBvV2JnYcEoPz/GNG6MszydpkDFJCcz8tdf2XzhAgD9a9Vi+Y4d2Hp6alJPz5492b59O+Hh4Xh4eDyyz9fXl+joaCIiIjRZalWP5Jppfhja8O61a9dQVZVq1ao9sv3QoUPMmzcPAB8fHy1KE8VMi9at2X/mDH8EBtKifn3i09P5/PBhPL75htG7dhEZF6d1iUbp17Awai9dyuYLFyhhZsbKIUP44eRJzYIU/n9nGhQU9Mj2VatWcerUKWbMmFHUgzTPpDN9grVr1zJo0CCcnJy4cOECzs7OmtaTkJCAh4cHmZmZjB07FmdnZ06cOEFQUBCZmZlYWloSHh6uaY2GQDrTwnfk4EE+nzyZX/6dSW6qKPSpWZMPmzbFR8PJesbibnIy4wMDWf3vSEAzV1e+X7mSyh07alrX/fv3cXBwoFGjRpw7d46xY8dSqVIl/vjjD9atW8fgwYP57rvvNK1RT+TWmBd19+5dqlWrRkxMDN9//z2DBg3SuiQAjh49ygcffEBoaCj29vZ06dKFadOmUa1aNdq2bcvWrVu1LlFzEqbaCT1zhjmTJrFp926y/v078PfwYGyTJnT09sZElil8LjpV5bszZ5iwdy/3UlKwMDFhVvfufLBqFaYODlqXx/79+/H392f9+vWkpqYyc+ZMoqOjqVq1Ku+88w4jRowwlqUpJUxf1JAhQ/juu+9o0aIF+/btM5ZviGJBwlR7EVevMn/GDFZu2EDiv7drVba3Z5SfH2/Wq4eDtbXGFRq+MzduMPK33zgWFQVAy/LlWbxwIdW6dQP5eVTY5Jrpizh06BDfffcdFhYWLFu2TIJUiOfk7uHB16tXc+3mTb6YMgX30qX55/59Ptyzh/Jff83bO3dyMjpafunJxZV79+i/bRsNVqzgWFQU5ayt2fD+++y9dIlq3btLkOrJpEmTUBSFNm3aPLZPVVX69++Poih07NgRRVHM83JM6Uwfkp6eTr169bh48SLTp0/nk08+0bok8ZykMzU8WVlZ/LplC4u+/JLAU6dyttctW5ZhPj70r1MH+0JYO9WQXU9I4LMDB1h15gyZOh0WJia826wZ05cvx65qVa3LMzoJCQl4eXlx+/ZtAgMDad26dc6+d999l8WLF9OsWTN2796NtbW1DPM+r9mzZzN58mS8vb0JDQ0tlMWRhX5JmBq2S+fPs/LLL1mzZQt3k5IAsDQ1pWvVqrxRpw7tvbywMDXVuMrCczkmhq+PHmVNSAhpWVmYKAqD6tRh+rx5uLdqJZ1oAVqyZAmjRo3C19eXkydPAjBt2jQ+++wzGjRowL59+7LXFJYwfR7h4eHUqFGD1NRUgoKCaNmypdYliRcgYVo0pKWl8dP69axctIigM2dyfrA4WVvTu2ZNetWoQTN3d6N8YLmqqhyMiODrY8f4+fLlnO09q1Zl5uzZVO/eHYzwcxuazMxMateuzaVLl9i8eTPR0dGMGTOG6tWrc/DgQUqXLp39UgnT5/Hqq6+ybds2+vXr99iiCKLokDAteq5FRLBhyRJ++PFHzv074QagtLU13apV49Xq1fGvVAmrQl7dR9+uJySwNiSE1cHBhP37eDtLExMG+fjwwdSpVOvcWUK0kP38889069aN0qVLc/fuXdzd3Tl8+DBubm4Pv0zCNK/27dtHq1atsLW15fLly//9QooiRMK0aAs9dYqAZcvYunMnYf8+sBzA2syMlpUq0cHLiw7e3lQ2gFtD8uJ2UhK/hIWx5cIFdl+5gu7f781y1tYMa9WKUdOm4ezrK8O5GqpVqxbnz5/H2dmZI0eO4Pn4IhgSpnmRmZlJvXr1OH/+PLNmzWLSpElalyTyQcLUOKiqyvngYLZ++y07du4kODLykf0VS5Wiubs7zd3decXdHS9HR4OYeZ+l0xF66xZB4eH8dPky/xcZmfND09zEhK7e3rw1aBBthw/HzNFR01oFLFiwgPfffx8AOzs7/v77b8qUKfPfl0mY5sXChQsZPXo0lStX5vz58zLpqIiTMDVON6Kj+X3DBnbt3EngiROPLa7vaG1Ng3LlHvxxdaW2szOVHRwwL+DJTLcSEwm9dYvTN25wMDKSw5GRxD/0KEQLExNaubvTrUMHXh0+nNK1a0sXaiDWrFnDm2++iaurKz4+PuzcuZNRo0axaNGi/75UwvRZYmJi8Pb25v79++zYsYNu3bppXZLIJwlT45eVlcW5P//kwM8/c3D/fg6FhHD735nBDzMzMcHTwYGqpUtTyd4et5IlcbOzw61kScrY2mJnaYmdpSUlLCweW51Jp6okpqcTl5pKfFoaMcnJRMbFERkXx7X4eP6+d4+zt2/net5KJUvyStWqdOrcmfb9+1PS01MC1MBs376dXr16YW9vz6FDhyhRogRVqlQhMzOT8+fPU6VKlYdfLmH6LO+88w7Lli2jTZs27N692yCGiUT+SJgWP6qqEhUezp9BQZw6coTTZ85wMSKCiPv38/wD67+342RkZeXpvSXNzaldpgx1vLx4+ZVXeKVDByo0bAjmebrPX2hg7969dO7cGQsLC/bt24evry8A48ePZ968efTo0YNt27Y9/BYJ06cJDQ2lfv36KIpCaGgoNWrU0LokoQcSpiJbSnIyf50+zaUzZ7j2zz9ER0YSHRVF9J073E1IID4lhfj0dBIzMnJ9fwkzM+wsLLCzssLRxoaKZcpQwdWVipUr4+7lRe2GDXGvVw9FlkcsMo4dO0br1q3JzMxk165d+Pv75+y7d+8elStXJi4ujkOHDvHyyy9n75IwfZoOHTrw+++/895777FgwQKtyxF6ImEqnpcuK4v0lBQUVYV/v3fMLC0xtbCQ4VkjcvbsWZo3b05CQgJbtmzJ9bLe559/zqRJk2jUqBHHjh3L3ixh+iR79+6lTZs2T5u9JYooCVMhhJ7JQve50el0fPTRRwBMnDhRglQIIUS+FbvOdN26dQwYMIDy5csTFhaGtVzvMCrSmQoh9Ew60/9KTU1l8uTJAHz22WcSpEIIIfSiWIXpokWLiIyMpHbt2gwYMEDrcoQQQhiJYjPMe+/ePTw9Pbl//z67du2iffv2WpckCoAM8woh9EyGeR/25Zdfcv/+fVq1akW7du20LkcIIYQRKRad6a1bt6hcuTLJyckcP34cPz8/rUsSBUQ6UyGEnklnmm3OnDkkJyfTtWtXCVIhhBB6Z/SdaVRUFF5eXqSlpREcHEzdunW1LkkUIOlMhRB6Jp0pwKxZs0hLS6N3794SpEIIIQqEUXem4eHhVKlSBZ1Ox7lz56hevbrWJYkCJp2pEELPpDP99NNPyczM5I033pAgFUIIUWCMtjO9fPkyNWrUwMTEhEuXLuHp6al1SaIQSGcqhNCz4t2Zzpo1C51Ox5tvvilBKoQQokAZZWd65coVqlatCsDff/+Nh4eHtgWJQiOdqRBCz4pvZzpnzhyysrIYMGCABKkQQogCZ3SdaWRkJF5eXmRlZXHx4kWqVKmidUmiEElnKoTQs+LZmX7xxRdkZGTQp08fCVIhhBCFwqg60xs3blCpUiXS0tI4e/YstWrV0rokUcikMxVC6Fnx60znzZtHWloaPXv2lCAVQghRaIymM71z5w4eHh4kJydz6tQpfHx8tC5JaEA6UyGEnhWvznT+/PkkJyfTqVMnCVIhhBCFyig608TERCpWrEhsbCyHDx/mpZde0rokoRHpTIUQelZ8OtNvv/2W2NhYmjZtKkEqhBCi0BX5zjQjIwMvLy8iIyPZsWMH3bp107okoSHpTIUQelY8OtPNmzcTGRlJ1apV6dKli9blCCGEKIaKdJiqqsoXX3wBwLhx4zAxKdIfRwghRBFVpId59+zZQ7t27XBxcSE8PBwrKyutSxIak2FeIYSeGf8w75dffgnA+++/L0EqhBBCM0W2Mz19+jQNGjSgRIkSXLt2DXt7e61LEgZAOlMhhJ4Zd2f69ddfA/D2229LkAohhNBUkexMr1+/jru7O6qqcuXKFdzd3bUuSRgI6UyFEHpmvJ3pkiVLyMzMpEePHhKkQgghNFfkOtOUlBQqVqxITEwMhw4d4uWXX9a6JGFApDMVQuiZcXamP/74IzExMTRo0ECWDhRCCGEQilSYqqrK/PnzgQe3wyhKnn5hEEIIIQpUkRrm3bdvH61ataJs2bJERERgaWmpZTnCAMkwrxBCz4xvmDe7Kx05cqQEqRBCCINRZDrTK1eu4O3tjbm5OZGRkZQtW1arUoQBk85UCKFnxtWZLly4EFVV6du3rwSpEEIIg1IkOtOkpCRcXV2Jj4/n1KlT+Pj4aFGGKAKkMxVC6JnxdKYbNmwgPj6exo0bS5AKIYQwOAYfpqqqsnjxYuDBxCMhhBDC0Bj8MO+xY8do0qQJTk5OREVFyaPWxFPJMK8QQs+MY5h36dKlAAwZMkSCVAghhEEy6M40JiaG8uXLk56ezt9//03lypUL8/SiCJLOVAihZ0W/M129ejVpaWm0b99eglQIIYTBMtgw1el0LFu2DJCJR0IIIQybwYbpnj17+Oeff3B3d6dDhw5alyOEEEI8kcGGafbEo+HDh2NqaqpxNUIIIcSTGeQEpOjoaCpWrIipqSnXrl2T5QNFnskEJCGEnhXdCUjff/89Op2Obt26SZAKIYQweAYXpjqdjm+//RaAoUOHalyNEEII8WwGF6b79u0jPDwcd3d32rRpo3U5QgghxDMZXJiuWrUKgLfeegsTE4MrTwghhHiMQU1AiomJwc3NjczMTK5evUqFChUK8nTCCMkEJCGEnhW9CUg//PAD6enptG/fXoJUCCFEkWEwYaqqKitXrgRg2LBhGlcjhBBC5J3BDPMeOXKEl156ibJly3Lt2jXMzc0L6lTCiMkwrxBCz4rWMG92Vzp48GAJUiGEEEWKQXSmCQkJuLi4kJycTFhYGN7e3gVxGlEMSGcqhNCzotOZbtmyheTkZJo1ayZBKoQQosgxiDBds2YNAIMGDdK4EiGEEOL5aT7MGx4eTuXKlbG2tubmzZvY2dnp+xSiGJFhXiGEnhWNYd61a9cC0LNnTwlSIYQQRZKmYarT6XKGeAcPHqxlKUIIIcQL0zRMDx06RHh4OBUqVMDf31/LUoQQQogXpmmYZnelAwYMwNTUVMtShBBCiBem2QSkpKQkXFxcSExM5PLly1SpUkVfhxbFmExAEkLomWFPQNq2bRuJiYk0adJEglQIIUSRplmYfv/994BMPBJCCFH0aTLMGxERgYeHB5aWlty8eRN7e3t9HFYIGeYVQuhbnoZ5zQq6ity4urry888/c+XKFQlSIYQQRZ7mKyAJoU/SmQoh9MywJyAJIYQQxkLCVAghhMgnCVMhhBAinyRMhRBCiHySMBVCCCHyScJUCCGEyCcJUyGEECKfnnfRhjzdbyOEht4A1mldhBCieHneRRuEEEII8R8yzCuEEELkk4SpEEIIkU8SpkIIIUQ+SZgKIYQQ+SRhKoQQQuSThKkQQgiRTxKmQgghRD5JmAohhBD5JGEqhBBC5NP/AxABr1NKdvYwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8,4))\n",
    "\n",
    "xvec = np.linspace(-1.75, 3.0, 100)\n",
    "ax.plot(xvec, f_func(xvec), 'k', lw=2)\n",
    "\n",
    "xvec = np.linspace(-0.8, 0.85, 100)\n",
    "ax.fill_between(xvec, f_func(xvec), color='lightgreen', alpha=0.9)\n",
    "xvec = np.linspace(0.85, 2.31, 100)\n",
    "ax.fill_between(xvec, f_func(xvec), color='red', alpha=0.5)\n",
    "xvec = np.linspace(2.31, 2.6, 100)\n",
    "ax.fill_between(xvec, f_func(xvec), color='lightgreen', alpha=0.99)\n",
    "\n",
    "ax.text(0.6, 3.5, r\"$\\int_a^b\\!f(x){\\rm d}x$\", fontsize=22)\n",
    "ax.text(-0.88, -0.85, r\"$a$\", fontsize=18)\n",
    "ax.text(2.55, -0.85, r\"$b$\", fontsize=18)\n",
    "ax.axis('tight')\n",
    "\n",
    "arrowify(fig, ax)\n",
    "fig.savefig(\"ch8-illustration-integral.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadrature rules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import polynomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import interpolate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = 0\n",
    "b = 1.0\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(-x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = -1.0\n",
    "b = 1.0\n",
    "\n",
    "def f(x):\n",
    "    return 3 + x + x**2 + x**3 + x**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(a, b, 100)\n",
    "xx = np.linspace(a-0.2, b+0.2, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAEYCAYAAAAnPkG+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmcTfUfx/HXd2as2beQLTshkdLCJEIY2XfZ17JkDdkmFGUnUhiUXVRIFMleQojwa7NmX2cw2/f3x72jGTOYGXPnzvJ+Ph4e5p77Pd/zuaM+53s/53u+x1hrERERERERERGRxMvD3QGIiIiIiIiIiMjDUYFHRERERERERCSRU4FHRERERERERCSRU4FHRERERERERCSRU4FHRERERERERCSRU4FHRERERERERCSRU4FHkgRjzExjzND7vG+NMYVdcNyWxpj1cd1vDI5fwPnZvNwVg4hIcmeMGWqMmXmf908aY16Ko2N5OfN+gbjoT0QkoXD3uDqhMsa8ZIw56e44JHFQgUcSNGPM38aYQGNMtru27ws/wLXWdrXWvhvf8VlrP7fWVo9OW2NMW2PMVlfHJCKSmBhjboT7E2qMuRnudUt3xxcd1tp3rbVd3R2HiEhCZ4x50Riz3Rhz1RhzyRizzRhTAWI2rk4InGN7P3fHIRKeCjySGPwFNA97YYwpDaRxXzjxwzjo/1ERSdKstenC/gDHAZ9w2z6/u31Sn7FojPF0dwwiIq5gjMkArAamAlmAx4CRwG13xuVOSf2cJvFPXx4lMVgAvB7udRtgfvgGxhg/Y8yocK/7G2POGGNOG2Pa36/zsGmPxpjBxpgLzllDLcO9n9EYM98Yc94Y848x5p2wwsvds3Kcs4q6GmOOGWMuG2OmOws1JYCZwHPOq9JX7hHLD8aY0caYbUAAUNAZT7VwbUYYYz67x/4ZjTGznZ/9lDFmlL4siEhi5sxjS4wxi4wx14FWxpjnjDE7jTFXnPluijEmhbN92C1MPYwxfznz+vvhC+bGmI7GmN+defobY0xe5/bBd80oCjLGfOp8L48xZrXzivOx8OcWZ4x+4V63dZ4vLhhj3n7A5/vMea5YZ4zxByoZY7YaY9reFe8P99g/tTFmgjHmhDHmrDHmI2NM6lj8qkVEXK0ogLV2kbU2xFp701q73lq7H+45ru7uzLnXjTHvGmMKGWN2GGOuGWOWGmNSOts+aDxfyxhzyNnPKWNMv3DvdTLG/M+Z378yxuS+K4ZIY/u7P5gzF39mjLnoPDf9bIx5NKpfgjO2gcaY/YB/uPNW4XBtIny3uWv/3MaYFcbx3eQvY0zPGPwbSBKnAo8kBjuBDMaYEs5iRVMgygIHgDGmJtAPeAUoAlS7V9twcgLZcFxJaAPMMsYUc743FcgIFAS8cRSb2t2nrzpABeBJoAlQw1p7GOgK7HBelc50n/1bA52B9MA/0Yg9vHlAMFAYeAqoDnSMYR8iIglNfWAhjly8BEee64Ujb78A1AS63LXPa0A54GmgEc4LBcaYRkB/5/vZgV3OvrHWjgk3m+gJ4AKw1NnfEhwzSnPjOA+NM8Z43x2occwynQa0wHFOyY3jHHM/LXBcxU4P7HjQL+MuHwKPA2VwnPMKAENi2IeISHw4CoQYY+YZY141xmSOxj41gfJARWAAMAtoCeQFShFulj/3H8/PBrpYa9M799sIYIx5GXgPx5g9F46x9+K7Yog0tgew1vpZa9s627TBcY7KC2TFMe6/eZ/P1RyoDWSy1gZH4/eAM14P4GvgV+fnrAr0NsbUiG4fkrSpwCOJRdgsnleA34FT92nbBJhrrT1orfUHRkTzGEOttbettZuBNUCTcAWlQdba69bav4HxOIow9/K+tfaKtfY4sAkoG83jh/Gz1v5mrQ221gZFdyfnVYJXgd7WWn9r7TlgItAshscXEUlotlprv7bWhjqv+P5srd3lzJN/4hjw311sed9ae9mZt6fw35eALsAYa+0R56B6FPCMMeaxsB2NMWmBVcCH1tr1xpjHgWeAt621t6y1e4C5RH0uaAysstZus9beBgYDka723mWltXaH8/NF+1YF50C/I468f9laew3HFxXlfRFJcJw56kXAAp8A550zZqKc6eI01lp7zVr7G3AQWG+t/dNaexX4BscFzfAijeed24OAksaYDM58uce5vSUwx1q7x5l/B+GYcV8gXJ/RGdsH4SjsFHbOTvrF+XnvZYq19oS19n5FoKhUALJba32ttYHOc+AnKO+Lkwo8klgswHGFsy133Z4VhdzAiXCv78yCMcbkCz/9Plyby85iUPh9cuO4CpCSiDNp/sFRMb+Xf8P9HACke0C8dzvx4CZRyg+kAM44p4ZeAT4GcsSyPxGRhCJCXjTGFDfGrDHG/GuMuQb44sjX99onLKeDI1dOD5cnLwChQJ5w7ecCB6y1452vcwMXojhPRHUuiHAOstbeAC7F5PPFQE4gFfBruM+zGuV9EUmgrLWHrbVtrbV5cMykyQ1Mus8uZ8P9fDOK1+HH2fcazwM0BGoB/xhjNhtjnnNuz024cb4zZ18kYn6Pzth+AfAtsNg4logYZ5y3Dt/Dw4z3c4flfGfeHwzcr0gmyYgKPJIoWGv/wTE1vhbwxQOan8ExPTJMvnD9HL9rQc8wmY0xj9y1z2kcA/8gHMk0/Hv3m0F0LzaW7fyBtOFe32uq/wkci9Rls9Zmcv7JYK19IoZxiogkNHfnxY9xXMktbK3NAAwj8iyZu88Dp50/nwA6hMuTmay1aay1uwCMMe/guOWpc7j9TwPZojhPRHUuiHAOMsakw7GYaEw+X3Tz/lkgECgW7rNktNZmfMDxRETczlr7O+CHo9ATF+41nsc58/M1HAXwVfx3++1pwo3znftnJYZjfWttkLV2pLW2JPA8jtu6Xr/fLne9DiD64/2/7jqHpbfW1opJvJJ0qcAjiUkH4OW7KvNRWQq0NcaUdE6zHx7N/kcaY1IaYyrhSMrLrLUhzv5GG2PSG2PyA324zxpA93EWyBO2GFwM7AOaGWNSGGPC1pKIxFp7BlgPjDfGZDDGeDgXoou0RoSISCKXHriKY3HKEkRefwdggDEmkzEmH9ATxxo64FjwfohzP5xtGjl/9sGxbkJ9a+2tsI6stX8Bu4ExxphUxpiyONZii/SUL2AZ8JpxLASdCsctYNEt8IfZBzQ0xqQxxhQFonxYgPMc9SkwyRiT3TjkMcYkmscMi0jy4Zx92dcYk8f5Oi+O22d3xuFhIo3nna9bGmMyOpc/uAaEONsvBNoZY8o6c/YYYJfz9t5oM8ZUMcaUdi7vcA3HBeKQB+wW3j6ghTHG07me6L3G7z8B15yLNKdxti9lnI+aF1GBRxINa+0f1trd0Wj3DY6pnhuB/zn/fpB/gcs4qvifA12dVxUAeuC4mvonsBXHiWBOjD+AI47fgH+NMRdisN9QoJAzvpHO49/L6zhuKTvkbL8cx4JxIiJJSV8cC1pexzGbZ0kUbb7GMWDeC6zEcZUYa+0yYAKOQf81YD/OBTNxrLmWAzgS7nbeaeHeK4LjfLEcGGyt3XT3QZ1Pg+mF4+LAKWf7f+9u9wAf4igKncNxvrnfRYW+OG4v+AlH0Wu9M04RkYTmOvAssMs4nhq4E8dszL5x1P/9xvOtgb+deb8r0ArAWvs9jrH2ChwzMAsRu/VscuI4N1wDDgObidkF4V6AD3AFx7pAq6Jq5Czs++BYB+gvHHcbfIpjgWcRjLUxvagkkrQYY14CPnPeCywiIomYMcYLx5XTx2N6BVZERBInjedFHDSDR0REREREREQkkVOBR0REREREREQkkdMtWiIiIiIiIiIiiZxm8IiIiIiIiIiIJHJe7g4gNrJly2YLFCjg7jBERBKcX3755YK1Nru744gryvciIveWlHK+8r2IyL1FN98nygJPgQIF2L37gU/LFhFJdowx/7g7hrikfC8icm9JKecr34uI3Ft0871u0RIRERERERERSeRU4BERERERERERSeRU4BERERERERERSeRU4BERERERERERSeRU4BERERERERERSeRU4BERERERERERSeRU4BERERERERERSeRU4BERSSIuX77s7hBERCSezJ49290hiIhIPNi7d2+026rAIyKSREyfPt3dIYiISDw4efIk/fv3d3cYIiISD0aOHBnttirwiIgkAQEBAUydOtXdYYiISDyYMGEC7dq1c3cYIiLiYocOHWLHjh3Rbq8Cj4hIEjB79mxeeOEFd4chIiIudvHiRfz8/HjrrbfcHYqIiLjY2LFj6dmzZ7Tbq8AjIpLIBQUF8eGHHzJw4EB3hyIiIi42bdo0GjRoQJ48edwdioiIuNA///zD119/Tffu3aO9j5cL4xERkXiwaNEiChUqxLPPPuvuUERExIX8/f2ZPn06W7ZscXcoIiLiYhMmTKBjx45kzpw52vuowCMikoiFhoYyduxYJk2a5O5QRETExT755BMqV65MsWLF3B2KiIi40Pnz51mwYAEHDx6M0X4q8IiIJGJff/01qVOnplq1au4ORUREXCgwMJDx48ezcuVKd4ciIiIuNnXqVBo1akTu3LljtJ8KPCIiiZS1llGjRjF48GCMMe4OR0REXGj+/PmUKFGCp59+2t2hiIiIC129epWPPvqInTt3xnhfFXhERBKpb7/9lps3b1K/fn13hyIiIi4UFBTEmDFjmD9/vrtDERERF5s2bRq1atWicOHCMd5XBR4RkUTIWsvIkSN555138PDQAxFFRJKyhQsXUqBAAV588UV3hyIiIi50/fp1Jk+ezI8//hir/VXgERFJhL7//nsuX75M48aN3R2KiIi4UHBwMKNHj2bWrFnuDkVERFzso48+omrVqhQvXjxW+6vAIyKSyISfvePp6enucERExIWWLFlCzpw58fb2dncoIiLiQv7+/kyYMIGNGzfGug8VeEREEpnNmzfz77//0qxZM3eHIiIiLhQSEsKoUaOYOnWqFtMXEUniZs6cSeXKlXniiSdi3YcKPCIiiYyvry9DhgzBy0spXEQkKVu+fDmZM2ematWq7g5FRERcKCAggPHjx/PNN988VD/6diAikohs2rSJ48eP07JlS3eHIiIiLhQSEsKIESOYNGmSZu+IiCRxM2bMoGLFijz55JMP1Y8KPCIiiYS1lqFDhzJ8+HBSpEjh7nBERMSFFi5cSLZs2ahevbq7QxERERe6fv0648aN4/vvv3/ovlTgERFJJL799lsuXrxIixYt3B2KiIi4UFBQECNGjGDOnDmavSMiksRNmTKFqlWrUqpUqYfuSwUeEZFEIGz2jq+vr56cJSKSxM2bN4+CBQvqyVkiIknclStXmDRpEtu2bYuT/lTgERFJBL766isCAwNp2LChu0MREREXun37Nr6+vixdutTdoYiIiItNmDCBOnXqULRo0TjpTwUeEZEELjQ0lKFDhzJq1Cg8PDzcHY6IiLjQJ598QpkyZahYsaK7QxERERe6cOEC06dPZ/fu3XHWpwo8IiIJ3JIlS0idOjU+Pj7uDkVERFzI39+fMWPGsHr1aneHIiIiLvb+++/TuHFjHn/88TjrUwUeEZEELDAwkHfeeYdPP/1UC22KiCRxEydOxNvbm3Llyrk7FBERcaHjx48zd+5cDh48GKf9qsAjIpKAzZo1i6JFi1KlShV3hyIiIi504cIFJk2axK5du9wdioiIuNjw4cPp1q0buXLlitN+VeAREUmgrl+/zqhRo1i3bp27QxERERcbPXo0zZo1o1ChQu4ORUREXOjAgQOsWbOGY8eOxXnfKvCIiCRQEyZMoFq1apQtW9bdoYiIiAv99ddfzJ8/n0OHDrk7FBERcbHBgwczaNAgMmbMGOd9q8AjIpIAnTt3jilTpsTpqvoiIpIwDRs2jB49evDoo4+6OxQREXGhLVu2cODAAZYvX+6S/hPE83aNMW8ZY34zxhw0xiwyxqR2d0wiIu7k6+tLq1at4nRV/YRA+V5EJKK9e/eyYcMG+vbt6+5Q4pxyvojIf6y19O/fH19fX1KlSuWSY7i9wGOMeQzoCTxtrS0FeALN3BuViIj7HD58mCVLljBs2DB3hxKnlO9FRCKy1tK3b19GjBhB+vTp3R1OnFLOFxGJaMmSJQQFBdGqVSuXHSOh3KLlBaQxxgQBaYHTbo5HRMRt+vXrx6BBg8iaNau7Q3EF5XsREaevvvqKc+fO0bFjR3eH4irK+SIiwM2bN3n77beZN28eHh6um2fj9hk81tpTwIfAceAMcNVau/7udsaYzsaY3caY3efPn4/vMEVE4sX69es5evQob775prtDiXPK9yIi/wkMDKR///6MHz8eL6+Ecs017kQn5yvfi0hyMXnyZMqVK4e3t7dLj+P2Ao8xJjPwGvA4kBt4xBgTac6StXaWtfZpa+3T2bNnj+8wRURcLiQkhL59+zJu3DhSpkzp7nDinPK9iMh/ZsyYQaFChahRo4a7Q3GJ6OR85XsRSQ7Onj3Lhx9+yLhx41x+LLcXeIBqwF/W2vPW2iDgC+B5N8ckIhLv5syZQ5YsWahXr567Q3EV5XsREeDSpUuMHj2aDz/80N2huJJyvogIMHz4cNq0aUPhwoVdfqyEMB/0OFDRGJMWuAlUBfRcYBFJVq5cucKwYcNYvXo1xhh3h+MqyvciIjgG+w0bNuSJJ55wdyiupJwvIsner7/+ysqVK/n999/j5XhuL/BYa3cZY5YDe4BgYC8wy71RiYjErxEjRlCnTh3Kly/v7lBcRvleRAT279/PkiVLOHz4sLtDcSnlfBFJ7qy19OjRA19fXzJnzhwvx3R7gQfAWjscGO7uOERE3OHgwYN8/vnnHDp0yN2huJzyvYgkZ2GD/REjRiTVJyVGoJwvIsnZokWL8Pf3j9cnJSaIAo+ISHJlraVnz54MHz4cLTApIpK0LVmyhKtXr9KlSxd3hyIiIi5048YNBgwYwNKlS/H09Iy346rAIyLiRsuXL+fChQt07drV3aGIiIgL3bhxg/79+7Nw4cJ4HeyLiEj8GzVqFC+//DLPPx+/a8urwCMi4ib+/v7069eP+fPn4+WldCwikpSNGTOGypUrU6lSJXeHIiIiLnT06FE+/fRTDhw4EO/H1jcKERE38fX15cUXX8Tb29vdoYiIiAsdPnyYWbNmsX//fneHIiIiLmStpXv37gwePJhcuXLF+/FV4BERcYMDBw4wZ84ct1T2RUQk/lhr6dq1K8OHDyd37tzuDkdERFxo0aJFXLhwgZ49e7rl+CrwiIjEs9DQULp27cq7775Lzpw53R2OiIi40Pz58/H396d79+7uDkVERFzo8uXL9OvXj5UrV7pt+QUVeERE4tmcOXMICQmhc+fO7g5FRERc6OLFiwwcOJA1a9ZoYWURkSRu8ODB1KtXj2effdZtMajAIyISj86dO8fgwYNZv349Hh4e7g5HRERcaODAgTRp0oTy5cu7OxQREXGhnTt38uWXX3Lo0CG3xqECj4hIPOrTpw+tW7embNmy7g5FRERcaPPmzaxbt87tg30REXGtwMBAOnfuzPjx48mUKZNbY1GBR0QknqxZs4bt27drYWURkSTu5s2bdOzYkY8++ogMGTK4OxwREXGhsWPHki9fPpo1a+buUFTgERGJD9euXaNbt27MnTuXRx55xN3hiIiICw0fPpzy5ctTt25dd4ciIiIudOjQIaZMmcKePXswxrg7HBV4RETiw9tvv02NGjWoWrWqu0MREREX2r17N/PmzdNsTRGRJC4kJIQOHTrw7rvvkjdvXneHA4BW+BQRcZXPP4cCBbAeHgyeNYtJzzzj7ohERMQVwuX73M8/z4oGDciRI4e7oxIREVdw5nwPLy9W7t1L53Tp3B3RHZrBIyLiCp9/Dp07Q0AABsgTEgK9e0PatNCypbujExGRuHJXvs8dFESu+fPhxReV70VEkpq7cn7O27ehSxcwJkHkfM3gERFxhSFDICAg4raAAMd2ERFJOqLI90b5XkQkaUrgY3wVeEREXOH48ZhtFxGRRMkq34uIJBsJPeerwCMi4gKhefJE/Ua+fPEbiIiIuMz163Au1T3yuvK9iEiSc+mRx6J+I4HkfBV4RERc4JMCBbjleVeKTZsWRo92T0AiIhKnzp2DKlWg7+3RBNy9rKXyvYhIkvP998H0uDGGAI9UEd9IQDlfBR4RkTi2cuVKRp38ix8mtuR6nixYA+TPD7NmJYjF10RE5OH89ZdjDeVDhywbHv2exd28OZXqMUIxBOd7TPleRCSJuXAB6tf354u0ldk8sfWdnJ/Qxvh6ipaISBw6deoUXbp2of1n7Tn+dEFOtKqEh/GgQ6YO7g5NRETiwK+/Qs2acPs2VK78LiG5ThA4phEv738NgPWbAsmfIr+boxQRkbhiLdSte5Hr19PR84sfOPnSM/RaOpKUJiUHt2Rzd3gRaAaPiEgcCQ0NpfXrrancsTIFni6Ah1GKTc4qVqzIggUL3B2GiMShH38Eb2/w9IRBg9bw+19z8Bnlg6f1dHdo4mbK+SJJ1+jRAezYkZXKbb/miZduY4xxd0j3pG8fIiJx5MMPP+TCzQtU7VMVL6MJkslZaGgoBw4coHTp0u4ORUTiyJdfQvXqkCsXLFlygrEftKHVrFakTps6QQ/2xfWU80WSrp07LcOGpSRHge00G38hwV/ATdjRiYgkErt372bsB2NpNrMZXh4q7iR3x44d4/bt25QoUcLdoYhIHJg9Gxo0gLJlYdOmYAYMbEa1HtXIUyZPgh/si+sp54skTZcuQZ06/nh6/Uvfb37DyyT82Zo6I4mIPKQrV67QuGljGo9rTNa8WXUlV9i3bx9FihRh5syZPPbYY2TLlo333nvP3WGJSAxZC++9Bx07wiuvwHffwZQpw7iV6hbe3b01W1MA5XyRpCg0FOrXv87Fiylp89FaMj9Kohjj66wkIvIQrLV06NiBoi8XpWz9sniS8Cv74nr79u3j77//xtPTkz/++IOff/6Zl19+mdq1a1OmTBl3hyci0RAaCn36wOTJjoejzJkDGzeuY878OfTb1E+zNeUO5XyRpMfXN5Aff0zPs40/o0JDi0kkszUTR5QiIgnU9OnTOfi/g/i8q0U25T+//vorTZs25c033yR16tRUqlSJUqVKceDAAXeHJiLREBgIrVo5ijtvvQXz58O5cydp07YNrT5uRYZsGRLFlVyJH8r5IknL99+Dr68XOfJ9T5uZN/BMBLdmhVGBR0Qklnbv3s0I3xG0ntuaVKlSabCfxH388cc8/vjj5MyZk2nTpt237b59+2jatGmEbefOnSNHjhyuDFFE4sCNG+DjA4sWOW7PGj8eQkODada8Gd6dvSn8fGGtu5MMKOeLJE8nTkC9erfwSvEH/TYcTRTr7oSnuaUiIrFw8eJFGjVpROMPGpO9QHYVd5K4WbNm0bVrV7y8vEiZMiU9evSgWLFivPLKK5Hanj9/njNnzpAzZ84727Zu3cqtW7d48cUX4zNsEYmhCxegVi345RfHwsrt2zu2v/3229xKfYsqvapo3Z1kQDlfJHm6dQtq1fLnhn8obyxcR8bsie8Cri4/iIjEUEhICM1bNKd0ndI8WffJRDVtU2Jn1qxZAMyYMYOvv/4aAD8/vyjb7tu3Dy8vLz777DNCQ0P57bff6NChA2PGjCFNmjTxFbKIxNA//8ALL8CBA7By5X/FnaVLl7J4xWKaf9xc6+4kE8r5IsmPtdCx420OHnyEal0XULpGykQ5W1NnKRGRGHrnnXc4f+s8nYZ10qLKycSRI0cA8Pb25vHHH2fu3LkULlw4yra//vorPj4+3Lx5kyxZspAzZ04GDhxI+7BviyKS4Bw8CDVqgL8/rF8PlSqFbT9Itze60W1FN9JnTp/oruRK7CjniyQ/06aF8vnnqSj89GIajvFItGN8FXhERGLgiy++YN7n8+izsQ8pPFNosJ8MhISEcOPGDQCyZcuGl5cXbdu2vWf7fv363fn5o48+cnV4IvKQtm2DOnUgTRrYsgVKl3Zsv3LlCvXq16PBqAbkLZ03UV7JlZhTzhdJfjZvht69LRkyb6bX6kt42hSQSIf4OlOJiETTb7/9RqcunWjr15YMWfUEleTi+vXrd35Onz59nPe/ceNGKleuTJYsWTDGMGzYMA4ePIiXlxcbNmyIVZ+rVq0iZcqUHDt2LMb77ty5E2MMH3/8cayOLZKYrF4N1apBjhywfft/xZ2QkBBatGxBkZeLUL5ped2Km4wo54skL3/9BT4+tzEef9H3u32kSpm4L+AmiAKPMSaTMWa5MeZ3Y8xhY8xz7o5JRCS8CxcuUMenDg1GNyD/U/l1JTeWEmO+Dxvsp06dGi+vuJ34euTIEWrWrMnt27d5//33WbBgAW3btqVPnz688MILUS7oGR316tWjdOnSDBw4MMb77t+/H4CnnnoqVseOayNGjMAYc88/KVKkiFY/R44coWXLlpQoUYKMGTOSNm1aihcvTp8+fThz5kyk9u+99x6NGzemYMGCGGMoUKBAlP0ePXqUYcOGUbFiRbJnz0769OkpW7Yso0ePxt/fP8p97vVZ0qVLF6ltaGgoEydOpHjx4qROnZq8efPSt2/fe/Yd0/bJmZ8f1KsHpUrB1q0Q/p94wIAB/Ov/Lz6jfPC0Ku7ElnJ+RMr5DxZXOR9ing8vXbpEv379KFy4MKlTpyZ79uxUqVKFLVu2xDq+mLaPacw3btxgzJgxlC5dmvTp05MtWzaef/55/Pz8sNZG+3eVXF2/DtWr3+TGjZt0mvsVOR9Pk6iLO5BwbtGaDKyz1jYyxqQE0ro7IBGRMIGBgTRo2IAn6z1J+SblE+09uQlEosv3YVP1o/ry/bBmz55NUFAQy5YtI1++fADs2LGDDRs2sGrVqofqu1evXrRp04bffvuNJ554Itr77d+/H09PT0qHTWVwswYNGkS59sX+/fv54IMP8PHxiVY/J0+e5MyZM9SvX588efLg5eXFgQMHmDVrFosXL2bfvn0RHmk8ePBgsmTJQrly5bhy5co9+50zZw7Tp0+nbt26tGzZkhQpUrBp0ybeeecdli5dys6dO6NcaLVSpUp07tw5wraovrjfwclpAAAgAElEQVS89dZbTJkyhfr169O3b18OHz7MlClT2Lt3L9999x0eHh4P1T45shY++AAGDoRXXoEVKyD8RI25c+ey7Mtl9N7QmxReiftKbgKgnB+Ocv6DxVXOh5jlw3/++YeXXnqJGzdu0KFDB4oWLcrVq1fZv38/p06dinV8MW0fk5hDQ0N59dVX2b59O23atKFHjx4EBASwaNEi2rVrx+HDhxk7dmy0f1/JTUgINGx4i//9LwW1+82ibO1HksQFXLcXeIwxGYDKQFsAa20gEOjOmEREwlhr6d69O0Hpg6jxTg3HlVyN9WMlseb7sKu5rpiqv3XrVooUKXJnoA+ONRyyZs1KrVq1HqrvBg0a0K1bN2bOnMnUqVOjvd+vv/5KsWLFEszTX8qUKUOZMmUibe/SpQsAHTp0iFY/VatWpWrVqpG2V65cmSZNmuDn58eAAQPubP/jjz8oWLAgAKVKlbrzpe9ujRo1YtCgQWTMmPHOtq5du1KkSBFGjx7N7NmzefPNNyPtV7BgQVq1anXfmH/77TemTp1KgwYNWLFixZ3tjz/+OD179mTx4sW0aNEi1u2To9BQ6NcPJk6Epk1h/nxImfK/97du3Ur/gf3pubon6TKlU3HnISjnR6ac/2BxlfNjmg9btWpFcHAw+/fvJ1euXHEWX0zaxzTmXbt2sXXrVnr37s3EiRPvbO/evTvFixfn448/VoHnPvr0CWbDhtQ8Wf0T6gxOlWRuxU0IJaqCwHlgrjFmrzHmU2PMI3c3MsZ0NsbsNsbsPn/+fPxHKSLJ0oQJE/hh1w80m9EML+Olwf7DSZT53hWD/eHDh2OMYceOHRw7duzOVO1ly5axatUqXnnllUizOW7evEmePHnIly8ft2/fjvBex44d8fT0ZPHixXe2pUuXjkqVKrFs2bIoY/j111+pV68eGTNmJHPmzHTu3JkbN25w4MABypYtG+vjxoeAgAAWL17MY489Rs2aNR+qr/z58wNw+fLlCNvDijsP8vTTT0co7oRp2rQp4HgK070EBgbes3AEsGjRIqy19O7dO8L2Tp06kTZtWj777LOHap/cBAVBmzaO4k6PHrBwYcTizh9//EHDxg1p+VFLHi3yaJK4kutmD8z5CS3fg3J+Usn5McmHP/74I1u3bmXAgAHkypWLoKAgAgICXBbfvdrHNIdfu3YNgNy5c0fYnjJlSrJly8Yjj0QaYonTzJmWKVO8yFVwOV0WBiep2fkJ4czlBZQDZlhrnwL8gbfvbmStnWWtfdpa+3T27NnjO0YRSYaWLVvGhxM/pMOiDqRNl1bFnYeXKPN92BfwuBzsv/rqq4wbNw6A5s2bs2DBAhYsWEC+fPm4ceMGzzzzTKR90qRJw8iRIzlx4kSEJ7UMGjSI2bNnM3XqVJo1axZhn+eee46zZ8/y+++/R9j+/fffU7FiRQ4fPsyQIUN49913+fnnn6lVqxZXr16NsBZDbI4b3saNG+nbty+dO3dm0qRJ/P3335HazJ8/nz///PP+v7Rwli5dyrVr12jXrh2enjEblN26dYsLFy5w8uRJ1q9ff+cq6sNePb/byZMnAXj00UejfH/58uWkTZuW9OnTkyNHDnr06MHVq1cjtPn555/x8PCI9N9D6tSpKVu2LD///PNDtU9O/P3htdfgs89g1CiYPBnC36124cIFar5ak5oDalKiaokkcyXXzR6Y8xNavgfl/KSS82OSD9euXQtAvnz58PHxIU2aNDzyyCMULVo0WoXxmMZ3r/YxzeHPPPMMmTJlYty4cSxbtozjx49z5MgRBg0axC+//MKIESMeGEtytG4ddO8eSrqMmxnww794eSStC7huv0ULOAmctNbucr5eThQDfhGR+LRt2za6vdGNbiu6kfWxrEkq8btRosz3YVdz43I9hooVK3L69GkAWrZsSe3atQHH2h8AhQoVinK/tm3bMnHiRN577z06derEp59+yvvvv8/IkSPp3r17pPZh/fz2228UL14cgLNnz9KkSRPKli3Lxo0b70zLb926NY8//jgQebHNmB4XHIUUHx8fvvvuO/LmzYuHhweffvopffr04dlnn+Wll14iRYoUrF27lgMHDkT6QnI/s2fPxhhD+/bto71PmE8//ZQePXrceV2gQAE+++wzKlWqFOO+7iUkJARfX1+8vLyivCXqmWeeoXHjxhQuXJhr166xdu1apk2bxubNm9m+ffud/9ZOnz5NtmzZSJUqVaQ+HnvsMbZv305gYCApndNQYto+ubh4EWrXhp9/hlmzoFOniO/fvHmTuq/VpWTtkjzX7jm8EsTwOElQzndSzo//nB+TfHjkyBHAMVOmSJEizJs3j9u3bzNhwgRat25NUFAQ7dq1i7P47tU+pjk8c+bMfPXVV3Ts2JEmTZrcaZs+fXpWrFhBvXr1ohVPcrJvH9SrF4hXiqMM/GEfjyTBC7huP4NZa/81xpwwxhSz1h4BqgKH3B2XiCRfR48epUHDBrSa0Yo8pfJomn4cSaz53lXrMezZsweAcuXK3dkWdotClixZotzH09OT999/Hx8fH+rVq8fGjRvp0aMHw4YNi7J91qxZATh37tydbWPHjuXSpUtMnjw5wpoLGTNm5JlnnuHbb7+NMF0/NscFuHLlCoGBgezdu/dOf6dOnWLJkiUsX76cqVOnkjJlSqpVq8a8efPufNF4kCNHjrB161aqVq0a7X3Cq1evHsWLF+fGjRvs3buXr776iri+NaR3797s3LmTMWPGUKxYsUjv79q1K8Lr119/nTJlyjBkyBAmT57MkCFDAMc0/qgG+uC4ohvWJmywH9P2ycHx41CjhuMxuCtWOJ6aFV5oaCitX2+NVy4var5TU+usxSHl/IiU8+M358ckH4b/N9+0adOdHFm/fn0KFizI4MGDadOmTZSL1Mc0vvu1j00OT5cuHaVKlaJu3bo8//zzXLp0ienTp9OiRQu+/PLLWD+ZLSk6fhyqVr1FUNAF+qz5jhz5k15xBxJAgcepB/C5c3X9P4F7l0hFRFzo9OnT1KhZg1pDalH85eKaph/3El2+d8V0fXAM9h999NEIizmGDTTu92jTOnXqUK5cOb7//nuaNWvG5MmT79k2rJ/wA5glS5bw4osvRnlLQHBwMHnz5r3zJSG2xwXIli0bGzduxNPTk6NHj3Ljxg2KFClCnz596NOnT4S2Fy9e5PLly2TOnPm+fYLjyic41oKIjTx58pAnTx7AUexp2LAhFSpU4ObNmwwaNChWfYY3dOhQpk2bRufOnWPUX//+/Rk5ciRr1qy5U+BJmzZthC9q4d26detOmzAxbZ/UHTrkKO5cvw7r10PlyhHft9bSq1cvjp05RqcVnZLcNP0EQjnfSTn/P/GR82OSD8MKX82bN49QPMmcOTN169Zl/vz5HDlyhBIlSjx0fPdrH9McfuDAAZ5//nkmTpxI165d72xv3rw5pUqVolOnTvzxxx8xvpU5Kbp0Cby9b3L5ciDtPl5KkWfTJdkLuAniU1lr9znvvy1jra1nrb384L1EROLWpUuXqF6jOhVaV+DZ1s/iZRJKDTzpSIz53hXT9QH27t0b4UouQNgaFJcuXbrnfkuXLmXfvn2A4wvI/b6QhvUT1u+///7L6dOnoxzoBwUFsWfPnkhXcmNzXAAvLy+WL19OoUKFKFasGOXLlydr1qzUrFmT2bNn88cff3D8+HFmzpxJ6dKlIzyG9l6Cg4OZP38+WbJkoX79+g9sHx1lypThqaeeirDWRGyNGDGCUaNG0a5dO2bOnBmjfVOkSEHu3Lm5cOHCnW1hr+9e6BQcV8azZcsW4ctITNsnZTt2wIsvQnAwbN4cubgDMGrUKL7d/C1tP29LqpSpVNxxAeX8/yjnx2/Oj0k+DCv658yZM1LbsILc3Qvxxya+B7WPaQ6fOHEit27donHjxhHapk2bltq1a/PPP/9EuQZSchMQANWq+fP33x7UHzaPZxoljceh30vS/WQiIjHg7+9P7Tq1KfBSAar0rqI1GOQOV0zXP336NP/++2+kdQ9KlSoFwLFjx6Lcb/369bRu3Zr69evTrFkz5syZw+HDh+95nP/9738R+vX39weivlrs5+fH5cuXI8UUm+OCY92H119/nVdffZVvv/2WDRs28M4773D8+HE6duxI4cKFyZ8/P71796Z169YULVr0vv0BfP3115w9e5bWrVvfcxp7bNy8efO+X7CiY+TIkYwcOZLXX3+dTz/9NMbFglu3bnHy5MkIizJXqFCB0NBQfvrpp0ht9+3bx9NPPx1he0zbJ1Vr10LVqpA1K2zbBk8+GbnNRx99xCd+n9BxWUfSZdDj0OU/yvlJI+fHJB+GFcDCFscPL2xbjhw5Hjq+B7WPaQ4PK5KFhIRE6is4ODjC38lVUBD4+ASwd28aqnaZyyu9Uib52fkq8IhIshcYGEjDRg1J83gaao+oreKOROCK6fpRrcUAjoUuM2TIwM6dOyPts2vXLho0aMALL7zA559/zqhRo/Dw8LjvbUA7d+7k0UcfvbMOTL58+fDy8uK7774jNDT0TrtTp04xcuRIgEhXc2NzXHBcQdy1axfTpk2jevXqVKtWjWHDhnHo0CF+++03Fi1axIoVKzhx4gRjx46N1sySsKntHTp0uGeboKAgfv/9d44fPx5h+7///htl+02bNnHw4EEqVqz4wOPfi6+vLyNGjKB169bMnTs3ynUawly8eDHK7UOHDiU4OBgfH58725o2bYoxhkmTJkVo+8knnxAQEEDLli0jbI9p+6Ro/nyoWxdKlHAUd6J62v3ixYvxHeNL5xWdyZQjk4o7EoFyftLI+THJh/Xq1SN9+vR89tlnd/79Ac6cOcOqVasoUqQIhQsXjlV8MWkf0xxesmRJwFGsC+/KlSt8+eWXZM6c+Z4LeCcHoaHQosVNNm5MyzP1/Wj0Hslidn7S/4QiIvcRFBRE02ZNuep1ldaTW2sNBonEFdP1wwb7d1859fT0pEGDBnz55Zfcvn37zhW+w4cPU7t2bYoWLcqqVatIlSoVhQoVokOHDsycOZNt27bxwgsvROjrxo0bbNmyJcJTOlKkSEGbNm2YPXs2NWrUoFGjRpw5c4YZM2bcaRM+ptgcN0z69OnvOfW/ZMmSdwam0XX69GnWrVvHM888Q+nSpe/Z7tSpU5QoUQJvb29++OGHO9u7devGmTNnePnll8mfPz+3bt3il19+YfHixaRPn57x48dH6GfBggX8888/gGMh1MDAQEaNGgVA/vz5ad26NQDTp09n+PDh5MuXj2rVqrFw4cII/Tz66KMRFrkcNWoUO3fupEqVKncekbx27Vo2bdrEs88+G+EJX6VLl+aNN95g2rRpNGjQgFq1anH48GGmTJmCt7d3pCd0xbR9UjNhAvTtCy+/DCtXQoYMkdusXLmSHr170G1FN3Lkz5Gkp+lL7CjnJ42cH5N8mDlzZj788EO6dOlCxYoVad++PYGBgcyYMYPAwECmTZsW6/hi0j6mObx3797Mnz+ft99+mwMHDvDCCy9w6dIlPvnkE86cOcP06dPx8kqeX/ethU6dbrJ8eRpKVv6c9nOC8EwmpY/k8SlFRKIQHBxMq9atOH3zNK/Pe50UnilU3JFIXDFdf+/evWTKlImCUUwv6NatG35+fqxevZqGDRty/PhxqlevTsaMGfnmm2/IEO5b67Bhw5g3bx4DBgxg27ZtEfpZsWIFAQEBdOnSJcL2yZMnkyJFClauXMm2bdsoXrw4vr6+rF+/no0bN1KgQAGAWB/XVfz8/AgJCYn14srNmzdn3rx5LFiwgPPnz2OMIX/+/HTp0oX+/fuTL1++CO1nz57N5s2bI2wbOnQoAN7e3ncKPD///DPg+H21adMm0nG9vb0jFHheeuklDh06xLx587h48SKenp4UKVKE0aNH06dPnztPSgkzadIkChQowKxZs1izZg3ZsmWjR48e+Pr6RjlTKKbtkwJrYeBA+OADaNwYFiyAqO6WWLNmDZ27dqbL0i7kKaknJErUlPOTRs6HmOXDzp07ky1bNsaNG8fQoUPx8PDgueeeY+HChVEWtWIaX3TbxyTm/Pnz89NPP+Hr68v333/P4sWLSZMmDWXLlmX8+PE0aNAgWrElRX373mTOnDQUKr+SN1deTzbFHQBzv1XbE6qnn37a7t69291hiEgiFhoaSpu2bTh06hDtPm9HqlSuWWDTWouH8aBDpuhN331YxphfrLVJZqGNhJDvvb29+fHHH1mxYkW8DZZq1qyJv78/W7ZsiXUf5cuXJ3/+/HzxxRdxGJlIwhIUBJ06wbx58MYbMHkyRPXAmPXr19OiVQs6LepEgXIFXLIGwwd1ajuOtSmQ/Cnyx3n/UUlKOT8h5HtQzhdJ7IYOvcWoUanJW3Itg7Ycx8t4umSMP65OLVKalBzcki3O+45KdPO9Ll2ISLITEhJC+w7tOfD3AdouaOuy4o4kTn/++SeLFy/m6NGjAHeeapQpU6Z4i2H8+PHs2LGD9evXx2r/VatWceDAAcaOHRvHkYkkHAEBUL++o7jj6wtTp96/uNNhfgeXFXck8VLOF0k63n33JqNGpSZX4e94e7PrijsJmQo8IpKsBAcH0/r11uz/ez/tF7UndZrUyS7xy/3t2bOH5s2bM3XqVM6cOXNn0F+kSJEo2wcEBNCsWTMee+wx0qVLR7ly5di1a9dDxfDEE08QHBxM9erVY7V/vXr1CAwMvGfMIondpUvwyiuOJ2bNmAFDh0JUqXzNmjW0aNWC9vPbU/DZgiruSCTK+SJJw6hRAQwbloachTYzZNsfpPBMfsUdiGWBxxhT1BhT3xjTxRjT2fmzMoqIJGhBQUE0b9GcY+eO0W5hO9KkTZMsE39MJbecX716dXLmzMn06dMpVqwYwcHB1KhRg7x580bZ/urVqzRt2pQjR45w6dIlKlWqRK9eveI5apHk4+RJqFQJdu+GpUuha9eo261atYq27dvSaWEnCj1bSMWdaEhu+R6U80WSAl9ff4YOTUvOQlt5Z/vvpEzhkWzH+NFebcgYUwLoCjQGHg3b7PzbOtucBZYCH1trD8dhnCIiD+XWrVs0adqEs4FnabOgjW7LeoDknPMzZMjAypUr6dKlC8ePH6dTp06MGzfunu1z5cpF/fr177xu1KgR69ati49QRZKdw4ehRg24cgXWrYMqVaJut3jxYt7s9SZdlnQh35P5VNy5j+Sc70E5XySxe+eda4wenYHcRbYyeOuhZF3cgWgUeIwxhYCxQH3gJrAF2AH8AVzEcQLIAhQGKgIdgR7GmC+AgdbaP10TuohI9Fy7dg2fuj6EZg3l9U9eJ2WKlMk68d+Pcr5DxYoV+fXXX6PVdvXq1UycOJHDhw/j7+9PUFAQNWrUcHGEIsnPTz9BrVrg5QWbN8NdT5y+Y8aMGYwYNYLuK7qT5wk9LetelO//o5wvkvhYC2+9dYnJk7OQt+RWBvxwkJReyfO2rPCiM4PnEHAAaAt8Ya31v19jY8wjQCOgp3Pf1PdrLyISl8aNG0eFChWo4ryse+7cOV548QU8MnnQe0VvvDy8kn3ifwDl/BjYtm0bHTt2ZNGiRTz//POkSpWK5s2bU7JkSXeHJpKkrFsHDRtCzpywfj0UKhS5jbWW0aNH8/Gcj+mxpgc58udQcef+lO9jSDlfJGGwFtq3P4efXw4KVdhKn29+I6Vn8nkU+v1E56zXxFr7tLV2wYMSP4C11t9aO89aWx5o+vAhiohEX4UKFWjSpAmbNm3izz//pFz5chw/eZxXh7yq4k70KOfHwN69e8mePTtlypQhODiY999/n6VLl1KuXDl3hyaSZHz+Ofj4QNGisG1b1MWdkJAQevXqhd8SP3qsVXEnmpTvY0g5X8T9QkKgXr3T+PnloGSVH+jzzUFSRvUIxWTqgWUua+2Xse38YfYVkcTnqafg/HkoXNidUVQhd+6lVK/emNBQHywB5C7xGd988ALffBD/xZ1zf6YnfbZbdNgf74eOFeX8mGnRogVffPEFefPmJX/+/PTp0wdrrQb7InFk8mTo3RteeglWrYKMGSO3uXnzJi1atuCvC3/xxuo3SJchnYr50aB8H3PK+SLudfs2eHsfZ9eufFRouIH2s/7Ey0Mzd8KL8W/DGNPEWrvUFcGISOJ2/jzcuOHuKCA0tAwhoT7YUD+y5OnNI5ledFsst/1T8N9alYmPcv79ZcmShY0bN0bY1qlTJzdFI5J0WAtDhsB770GDBo5ZPKmjuCHo/Pnz+NT1IXW+1HRc1pFUKbWAfmwp3z+Ycr6I+1y5YqlQ4Tj/+19+qnX/lkaj/tYC+lGITblroTEms7X24ziPRkQStbCZOz/84J7jW2uZNm0aQ0cMJXX6QKp0qsWPcz/FZ6AnxSsVd0tMH9SphaPA84hbjh8HlPNFJF4FB0OXLjBnjuPv6dMhqtn3v//+O3V86lDqtVLUGFwDL6PbcB+S8r2IJEj//BNEuXL/culSbhqPXk21bqfxUHEnSrG5OdkP+MgYMySqN40xzxljtjxUVCIiMRQYGEiXLl14f8L7BNkgus3vRv0h9ekytwsft/uY37f87u4QEys/lPNFJJ7cvOlYTHnOHBg2DGbMiLq4s27dOl6s/CKVe1em5js1VdyJG34o34tIArN16xWKFbvElatZ6TJ/Na90O6M11u4jxjN4rLUdjTHngXeNMdmstW8BGGOKAe8Br+F41KKISLy4ePEiDRo24Hba25RvXp4iFYtQonIJAIpXKk6XuV34e8/fbpvFk5gp54tIfLl8GerWdSykPH06dO8euY21lsmTJzNm7Bjaz2tPwYoF8cIrMd8Jm2Ao34tIQvPxx//QrVsWUqX1YuA3a3i87BUV8x8gVisSWWsHGWPOAuONMdmBG0B7wAIfA75xF6KIyL3t2bOH+g3r8+RrT/Lq0KiflFW8UnEVdx6Ccr6IuNrp01CjBhw9CosXQ5MmkdsEBATQtVtXtu/ZTq9ve5E9b3ZdxY1jyvcikhBYC5067WP27NJkyXOKfus2ky33LRV3ouFhlpz+BKgNtMCR9BcDQ621f8ZFYCIiD+Ln50ff/n1pPK4xZeuXxdN6KvG7jnK+iLjE0aNQvTpcvAhr10LVqpHb/Pnnn9RvUJ+MRTLSc11P0qRNo3zvOsr3IuI2AQHBPPfcL+zf/yxFnj/EG4t28kiGUDRVM3pifNnDGJPCGNML+AN4GdiDI/mnAE7EbXgiIpHdunWLrl27MnzMcN786k2eqvcUXmj9BVdQzhcRV9q9G154AQICHAv0R1XcWbt2LRWfq0jp5qVpMauFijsuonwvIu524MAFcuf+nf37n+Xlrrvo89VWZ3FHois281qPAROA80Bda+3TOCr8dYG1xph0cRifiEgER48epeJzFTnw7wHe2vAWjxV/TI9IdC3lfBFxiQ0b4KWXIF06x7o75ctHfD8oKIj+/fvTvnN72vq15cXOL2oxZddSvhcRt5k8eQ9ly4Zww78wHT5dT7Mx+/Hy1G24MRWbW7Q8gU6An7U2FMBau9QYcwVYAWwyxrxqrb0Qh3GKiLBw4UJ69OpB7UG1eb7d83iiW7LigXK+iMS5xYvh9deheHH49lvIlSvi+3///TfNmjfDZrL039yf9FnSa70d11O+F5F4FxQUQu3a37Nhw8tkynWJN5d9Sb6S1zTGj6XYFHiKWGtv3b3RWrveGFMNWANsB4o+bHAiIgBXr16lR88ebN6xme5fdCdPqTyatRN/lPNFJE5NnQo9e0KlSvDVV5ApU8T3Fy1aRI9ePXil5ytU7lY5ysXzxSWU70UkXu3Zc5xqr5zg8qXqlKp+hI6ztvNIxhC03k7sxfhSSFSJP9x7u4BKQOqHCUpEJMyWLVt4suyTnPU6S99NfclXKp+KO/FIOV9E4oq1MHSoo7jz2muOmTvhizuXL1+mWfNmDPEdQtdlXfF+w5sUnilU3IknyvciEl+stbz99rc8XSEFV68+Q+P3ttBj0Y/O4o48jId5ilaUrLWHjTEvxHW/IpK83Lp1i2HDhuH3mR9NJzTliZpP6ClZCZByvohER3AwdO8On3wCHTvCjBngFW4Uun79ejp07ECpWqV4a+NbpEqdSrdkJTDK9yISF06cOEfVats5drQuWfJeouuCVTz+5FU0ayduPPDMaYyJ4nkG92etPeHct1psghKR5G3Hjh08+dSTbD+2nQGbB1C6Rmk9JSueKOeLSFy7dQsaN3YUdwYPhlmz/ivuXL16lQ4dOtC2Y1saTmrIa++9RurUqVXciQfK9yISn6y1vP/+Oh4veIFjR+vxQutDjNz+tbO4I3ElOmfPdcaYjcaYOsY8+L4I5yMW6xtjNgNrHz5EEUku/P396dOnD3Xr16XKwCq0mtuKjNkyaqAfv5TzRSTOXL0KNWvCqlUwZQqMHg1htfo1a9ZQqnQpTnCCAVsHULJKST0lK34p34tIvPjrr9OUKLGMQYOqkjpDHt5YspY2k3eQOp1uyYpr0blFqywwEfgKuGCM2QD8BPwBXMIxlyoLUASoCLzsfP2tc18RkQdau3Yt3bp3o8AzBRi4dSAZsmbAAw/N1ox/yvkiEifOnIFXX4VDh2DhQmje3LH99OnT9Ordi52/7KTRlEYU8y6mpyK6h/K9iLhUSEgI/ft/yeQpTxAa0oQKDX+nxbifSJclEA3yXSM6BZ7GQGcgF9AdeA1oDti72hngGvAFMMNa+3McxikiSdSpU6fo/VZvdu3ZRYOJDShepbjW2nEv5XwReWj/+x9Urw7nzsGaNfDKK46B/syZMxk2YhiV2lZi4OSBpEydUrM03Uf5XkRcZtOmAzRt+ifnzzcgY84rtJq8hidfOa0xvotFp8AzDDhqrV0I7DDGNAKOAyWB7DhOAueBg8Bea22oq4IVkaTj9u3bTJgwgQ/Gf4B3e28GTB7gGOhr1o67KeeLyEPZsxmFi7EAACAASURBVMdxW5a1sGkTVKjgeCLiGz3ewKQ3vPn1m+QqlkuzdtxP+V5E4tzZsxdp1Gg927bVBFOCam/s47W39zpvx1LOd7XoFHguA2nDvV4CtLbW+rkkIhFJ0qy1rF69mt59epOjaA76fteXrPmzaqCfcCjni0isbdwI9epBliyOx6CnTXuCFi0HsnHzRl7zfY2n6j+lfJ9wKN+LSJwJCgrmrbe+YebMwoSENKfI88dp8cF35CmpRZTjU3QKPAeBDsaYb6y1p1DZTURiac+ePfTt15e/z/yNz3s+lKxWUrdjJTzK+SISK8uWQatWUKQILF9+nXnz3mPGxzPwbu/NkF1DSJVWjz5PYJTvReShWWuZOHErQ4d6EBDgQ+bcl2j83reUr3McDw+llfgWnQLPSOAb4Lgx5hiO6Zq1jDGXgT3W2rNxEYhz9f7dwClrbZ246FNEEoa//vqLYcOHsW7DOl4d+CqNWzXGy/P/7d15fNTVvf/x15k1kz0hhIQsQADZJSK4oFhRQUQtWkUFqq3YuhW0XoWi1LW/2l659mpbe9XrUhcUvdXea8W6gmIRFURkVQuBhLAlQPZ1MnN+f0zCYoKADZmZ5P3kMY9kJjPz/WTyzTuHz5zv+bp0OFZkOuaZr7wX6Xz+9CeYMQNOPTXIpElPcsaZcxly1hDmLJlDcs9kHDjUzI88GuOLyL/k+efX8PNbytmzewwxCdVMvn8pZ179JZ4Yiwb54XHYt1GstYuAYcCDwB5CP6mpwOvAdmPMdmPM68aYXzWfOrH3d6zlZmDDd3ysiITb/Pks+Lg3iz5wQO/eMH8+O3bs4MYbb2TEyBFUZVTxy09+yek/Oh2PS4tqRqoOynzlvUi0mz8fevfGOhyUJ/fmHz+bz/DhWykqHsj/vPs41758LZc9chmpPVNxGs3UjEQa44vIEWnOexz7x/jPP7+JzMxPufLKYVRWj2Di7E/57Zr/YfwNG5qbOxIuRzKDB2vt18BsAGNMEJgFFAAnAicAI4CJNK+6b4wps9amHWkRxphs4Hzg18C/HUX9IhIJ5s+Ha68lo6E2dL2wkMYf/5g7vW52XD2WX37yS+K7xeM0zvDWKUfkWGa+8l6kE2jOfGprMUByRSFPMJ17qnoQ/NNl9Dutnw6/jRIa44vItzog7wEoLKT2hz/hDZ5gr+/7nD/rY86+4UsSUprCW6fsc0QNnm+4D/jQWvsp8NeWG40xmYT+CJwI5B/lcz5E6I9LwqHuYIy5ltCpHMnNzT3KpxeRY2ru3P3B38zT1MSD3eN58b4LtaBmdGvvzFfei0S5pjlzcX0j82Np5M6GGhaMHoDBaGZ+dOrwMb7yXiTCtTHGj6Weh5Ju5ZS19cTGa7ZOpDnqBo+19p5D3L4DWNh8OWLGmAuAEmvtZ8aYM79lu48DjwOMHDlSe5JIBLFFRW2O5RN3luMy36WPLJGiPTNfeS8SnRobYelSePNN+PvfA6wqLmrzfvHb9qqZH8XCMcZX3otErrVrLYMLi9pc06V75S41dyJUJCyCcRrwfWPMFmABcJYx5vnwliQih2OtZcmSJVx08UUUH2JAX52V2sFVSYRT3otEiY0b4Y9/hAsvDJ3y/Kyz4D/+I8CGLz9muyulzcco8+UblPkiUcbvh5deCjB0aAnDhhmKyG7zfsr7yBX2Bo+19nZrbba1tjdwBbDIWvvDMJclIofQ0NDAc889x4iRI7jyJ1cSNzqOdQ9Pw+/zHHQ/v8/D8rsuDlOVEomU9yKRq7oa/vY3+NnPoG/f0KnOZ86EVasayMp+D1/cFL437Qfc+Y/32PDI95X5cljKfJHosXUrzJ5dR1paNVdc4eTrjX5OuvRV1syboLyPMjp2QkSOyObNm3n00Ud56s9PkTM0h9N+cRpDzh6C0+Gk2Dj40OthwMw3yWzYTk12CsvvupiCyaeEu2wREWmDtfDFF6HDrt56K3QIlt8PcXHwve8FGD16BavXPsCu3Us5fsoYrr9qNL5kH06cFB7Xkw+Nk37XvUO23arMFxGJQn4/vPEG/Md/lLN0aSLWeknJWM4ls77krBvqcTsc7DIn8WFSUHkfRSKqwWOtfR94P8xliEizxsZGXnvtNf77yf/m0+WfcuqUU7npjZvonte91cLJBZNP4YZnfgXArNeP6jB96YKU9yIdr7QU3n471NB5+23YtSt0+/DhcMstMGDAZlau/CMvvvRn+uT34fRZJzNk3F24XC4c5uBJ3wWTT+GC2x4B4Pdrnuvob0WijDJfJHJ8+SU8+mgDTz/dRGVlHA5HAwNOf5WL7ykj94RA8xh//5lvlffRJaIaPCISGb744gueeeYZnpv/HD2P68mJPzyRXz35K9w+Nw4cWkRTRCQK+P2wbFmoofPWW7ByZWjmTrduMH48nHsunHRSOUuWvMSTTz/Jn58v5OSpJ/Nvb/8bqb1TceFS3ouIdAJ798KLLwZ55JEqNmxIApykZHzK+devZ8ItFo/X0dzIV3sg2kXlT9Bardgt0t6Ki4tZsGABzzz3DLvLdjPqslGHnK0j0lEaGxvDXYJIVNmyZX9D5913oaoKnE449VS4916YMAGGDGngnXfe5NnnnmXmze8w7KxhnHTLSfz4rB/jdDnVyJewqKqqwlqrfU+kndTVwcKF8F//Vcn778cSDLrwxOzk+AmvcdGd1WQOMhrjd0JR2eBZv34977zzDuPGjQt3KSJRraSkhFdeeYX5L85n7bq1jJg4gvH3j6ffKf1wOBwa5EvYrV+/nnvvvZdZs2YRGxsb7nJEIk5tLbz/fqih8+ab8PXXodtzc2HKlNAsnbPPhthYP4sWLeJPf1rAX//vr+QOzmX45OHc9+B9xCTGaJAvYVdUVMS5557L7373O4YOHRruckSiUmNjqLn/xBOVvPGGl4YGLw5HHbnHv8f4m3eSP8ngxIHDuMNdqhwjUdngycrK4oYbbmDw4ME88MADDBw4MNwliUSNrVu38te//pWXX3mZ1V+sZvi44Qy/YThTzpqCy+NSU0ciyuDBg9mwYQMDBgzgN7/5DVOnTsXhCPsJIEXCxlpYu3b/LJ0lS0ID+pgYOPNMuOGG0CydAQOgvr6Od999l5tu/guv/e01MvpmMPzi4dzxjztIyEzQIVgSUQYPHsyFF17IWWedxeTJk7n77rtJT08Pd1kiEa+hAd55x/LEE+W89VYM9fU+jLFk9l/MaT8u4oxrArjdDhzK/C4hKhs8ycnJLF26lD/84Q+MGTOGSZMmcffdd5OTkxPu0kQ61AMPPMCoUaMYO3bsvtsWL17M8uXLmT17NhA6pHHVqlW89tprvPp/r1JYWEj+hHyGXz+cKWOn4PQ69c6tRCyPx8OCBQtYunQpt956K/PmzeP+++9n4sSJ2mely9i7F955Z//iyNu2hW4fMgRmzAjN0jnjjFCTp7S0lIULFzLnjv/lvffeo8+wPgy+cDC3z7qdhJ4JOHG2WjBZJBIYY5g5cybTpk3j3nvvZdCgQcyYMYNbb72VxMTEcJcnElEqK+Fvf2viySf3sHRpEo2NMRjjJLP/Uk6espUzf9qIN5bmN241W6cricoGD4DX6+W2227jmmuuYd68eeTn5/OjH/2IWbNmkZmZGe7yRDrEqFGjuOyyy3j55ZcZO3Ysixcv5rLLLuOpp57ilVdeYeHChbzx5hu4Y90MnTCUcb8aR95JeVpjQaLOaaedxrJly3jttdeYPXs2v/3tb7n33nsZO3as9mPpdJqa4NNP98/SWb4cgkFISYFzzgk1dM49F7KzIRAIsGLFCn7727/z+huv89VXXzFs7DAGTBjAvfPuxZfq00wdiSqpqak8/PDD3HLLLdx9993079+fWbNmcf311xMfHx/u8kTCZtMmeOGFChYsqGbDhnSsdeN0ucgesoRTp+5k9I/q8HhamjrOwz+hdEpR2+BpkZKSwv3338/MmTP593//d4YMGcIPf/hDZs+eTXZ2drjLk07uhBNCp53t1y9cFYylZ8+XOffcy0hNvZrS0ifxeu/noouSiE2MxZdyG3GZd+P2uSlYBQWrAI7dIH/T8nRc7uAxe37p2owxTJo0iQsuuID58+dz/fXXk56ezp133sn48eP1H1iJasXF+9fRefddKC8HhwNGjYI772w54xU4HJaCggL+/vdFvPn2myxatIjkHskMPHsgZ/zyDKafMh2H26GZmRL1evfuzTPPPMPatWv51a9+xbx58/j5z3/Oz372M83okS6hthbefLOWZ5/dxZIlsZSV9QCS8MWXMfjsd/ne1eUMHl+B29lyBiw1daQTNHhaZGZm8tBDDzFnzhwefPBBjj/+eCZNmsStt96qhdrkmCktherqjt+utUGqqqooLy+nrCwFf9P57No1j5iEa0nNPRtfog/jMM2tnI4b4NuAoQlN/Zdjy+l0ctVVVzFt2jRefvllbrnlFnw+H7fddhuXXnopbremIkvkq68PrZ/T0tRZvz50e8+e8IMfhBo655wDKSmWzZs3s2TJEv7rsUUsWryIhsYGBp4xkL5n9eX2e28nITN06JXBqKkjnc7QoUN56aWXWL9+Pffffz95eXn89Kc/ZebMmfTs2TPc5Ym0m0AAliyp4dlnt7F4sYOiohysjcWYnnTvs45xU5dz2lUVZPSvPmAmfqf577y0k063R2RkZDBv3jzmzJnDo48+yrhx48jPz+fnP/8548aN0+Kc0q5aZu68//6x3U55eTnLli1j6dKlLFm6hJWfrSSzbyb9z83Dm+pl0aOLOGP6eXz49Atcfn88g84YdGwLOoSbel0Zlu1K1+R0OpkyZQqXX345Cxcu5MEHH2TOnDncdNNNTJ8+nZSUlHCXKLKPtfDll/sbOh98EGryeDyh9XOmTw81dQYMaGLt2jV89NFH3DjjQz5Y8gFNgSaOG30cuafmct3PriOtXxoOdKZD6VoGDx7M888/T0FBAQ899BBDhgxh0qRJ3HTTTYwYMSLc5YkcNb8f3nqrlJde2sHSpU4Ki3IJBhKA40hI28wJFy7j5MsrGXjmLmJ89oD10zRTRw6t0zV4WnTr1o25c+dy6623Mn/+fObMmcOMGTO44YYbuPrqqzXwl4jl9/tZt24dn3zyCcs+XsayT5dRXFRMvxH9yDkph/wb87nkpEvwJHrY+OFGHp/+ONf/+XoGjhnI4DMG89jVj3Hd09cxcIzOLiddg8Ph4MILL+TCCy9k+fLlPPzww+Tl5XHJJZdw4403auAvYVNeDu+9t38tnaKi0O0DBsC118L48ZbjjtvB2rWf8PHHH3PjjGWsXLmS7tnd6TOqDzmn5TDjthmk9EkJLY6sho4IeXl5/P73v+eee+7hscce46KLLiIrK4sbb7yRyZMnExMTE+4SRdpUVFTLggVbePvtKlavjmP37r5Y2x3oTny3YoaNX8sJ369k0NidJKc3HDArs+UicnidtsHTIiYmhmuuuYbp06ezbNkyHnnkEe677z4mTpzI9OnTOeusszSrR8KmsbGRdevW8fnnn/PZZ5/xyYpPWL92Pem56fQa0YusEVlM/vFkeg7qicPd/G7tAVPwCz8vPKiZM3DMQK57+jq2rNyiBo90SaNGjeL555+npKSEJ554gosvvpi0tDSmT5/OlClTSE1NDXeJ0okFArBy5f5ZOh9/HLotMRHOPtty/fVlpKV9xrZtS/l0xae8+JMV+Jv85J2QR9aJWeTPyOfiEy8mJjlGM3REDiM1NZXbb7+dWbNmsXDhQh555BFuueUWpk6dyvTp08nPzw93idKFbd1azf/+7xYWLapg1Son27b1xO/PBQZjTIDUnCJOuXw1Q8aX0+/UnaT2aPxG3uv/p/LddPoGTwtjDKNHj2b06NHs2bOHF154gdtuu42ysjKmTp3K1KlTGTZs2P4HzJ8Pc+eG3m7LzYVf/xqmTQvfNyCRZ/58Fnw8l/SGIuh9+H1k165drF27ltWrV/P5F5+z6otV/POrf9KjVw9yjs+hx7AejL1rLFcdfxXuePcRDe4n3Dyh1W0DxwxUc0e6vPT0dO644w5+8Ytf8N577/HUU08xd+5cxo0bx9SpUznvvPP2v8urvJfD+ZZ9ZMeO0KnL33wzdCrzPXvAGMvAgXVMnLiRmJgP2Fn6v3zw4SreX2rpPaw3mcdnknNZDmfcfwaJOYk40KLIIt+Vy+Vi0qRJTJo0ic2bN/PMM88wadIkUlJSmDZtGldccQU5OTn7H6DMl8M5in2ksdHy4YfbeeednXzySS1ff+2lpDSTJn8OEFoHNjallL6n7mLAGUvpe8pu+uTvxhsb/EbmK/+lfXSZBs+BunXrxsyZM5k5cyarVq3ixRdf5IILLiAxMZHJkydztddL9n33YWprQw8oLAzNpQb9AZCQ+fPh2mvJaDh4H7FA6bhxbNiwgfXr17Nu3TpWr1vNunXraGpqImdQDplDMukxogcXXHkBPQf1xOlz6p1akWPE6XQyfvx4xo8fz969e3nllVf4wx/+wDXXXMOkSZOYkZrKiEcfVd7LoTXnPQfsI4FrruWlF+G3RVNYsyb0LmtcXBUpqStId/8f5WV/od542OPNpMfgHuQPHcKFQyYQ1yMOh2k9G1NE2kefPn245557uOuuu3j//fd58cUXyc/PZ8iQIVx66aVMA7rdfvtBv8/KfDlIG5lvr72WvXssizLH8dFHJaxeXc+mTQ527Uqhvj4byApdTICkHrsYeMYe+p5STG7+HnoN30tiWkMbea/8l2PDWGvDXcNRGzlypF2xYkW7PmcwGOSjjz7i1Vdf5d9+/3uyA4HWd+rVC7ZsadftSnQK5ubi2Lq11e1FDsOwpASyB2STPiCdtAFpZA3IImtwVmhgT+ce2Lcssvz7wufCXEnkmHfBRMDw1T8yO2R7xpjPrLUjO2RjHeBY5D3Atm3b+Mtf/sLlv/gFGQ0Nrb5ek9aLFX/ZQkYGZGSEDrHphL+y0oZAAHbvhl27YOdOGD21F/F7ilrdbwu55JlnSOr2MTlDN5A7qoyMAenkDM4hLS8Nh8fR6Zv3yvyDzbvgfADeXtxIL3evDtlmZ8r8Y5X3DQ0NvP3227z66qvc9+yz5ASDre+kMX6XV1MDBQXQZ2wu8Xtaj/G30Is+bAHAOPwkpJWQnldGr/w6svMryBpYRs/jKvHEfHNmTuegvG/tgQsm4jEe1n6Y1iHbO9K875IzeNricDg4/fTTOf3007EPPdTmfWxhIRv/+U/69evXKX9xZT9rLXv37mXz5s0UFBSwadMmNm7ayNcbv2bTpk0UF+9o83E51vKbjb/RO7QiES4rK4ubb74Zbrmlza/7dhdx5pn7r8fEhBo9PXoc/PHAS8ttsbEd8z3IkbM2tOBxS9PmwMuuXbB9e5CtWxvZuRPKyz1Yu3/tgwCtB/oAvUwRf9y1HuNKwsFoHM3rJSjzRSKL1+vdtxC/feaZNu8TLCxizuwqzj03gZNPhvj4Di5SjrlAIHRIbUGBZdWqClavruTLL/0UFbkp3Z1AfV3oBDwBitt8fC5FzPif10nPqyQtpwaX61B5r78BEl5q8LTB5OaGpmx+w85YD6eceQoep4fRY0Zzzhnn8L0x32PQoEEa0HWQBx54gFGjRjF27Nh9ty1evJjly5cze/bsI34eay27d++mqKiIwsJCCgsL2bJlC5u2bGLzls0UbSnCYunRqwdpfdJIyU0hdXAqI88fycS+E6k+fx6JxWWtnrc6KxWPw9Mu36uIdIBD5H1JXBwZPa9kT5GLnjkjyEzPJz62P4FAdwoKnHz0UWiWR1uTYBMSWjd9Wj7PzNx/W3p66BTZ8t3V1LTdsGn5uGNHy+eWxsbWf6eN8WMcJdjgDjwxZcQlV5E7rJakrDp69LX0HOCk4v+lkrJrT6vHVmel4nHrBygSTQ41xt9qMpk3L45588CYIDk55Zx+uoezz45n1CgYNAhc+l9TRGtogOJi2Lw5wNq1laxbV83GjX6Kix2UlPioqkrFWjehBkwykIjHV0p8t11kDd9CZv9acvMbqfj3bqSU7G71/DXZKQw/Z2fzNf2/TyKXoqotv/71wcdeAn6fhw3/eSW/vnQkOzftpOCjAp59/1nu+c091FXWkX9iPqeefCqjR43mxBNPJCsrS02fY2DUqFFcdtllvPzyy4wdO5bFixfvu94iGAxSUlLC9u3b2bZtG8XFxRQXF1O4tZDCrYUUFxezo3gHMXExpGWnkZKdQlJWEsm5yWSPzOaEnBNI7ZWKN8mLwzTPwmn51/wzXXHXDxhz83O46xr3bdfv87D8ros7/DURkX/BIfJ+/e8u597JJ1FVUUXBsi/Z+PHrFHxcwNa1W8k7Lo+JF5/EqSeOoU/vk0lOHsCePe42mgqwZk1o4d2KirY3361b62ZQW7ODunUDp7ODXpMwa2g4uElz4MdvNnKqq1s/3hhLTEwVLtdugnYHjf4igoHtJHWrITG9huTMBlL7+MkYAFkDfXTvlUZCZgIO1wEzLzEY4wAsn3svUt6LdBaHyPxND53HvLOfYvn/OFnzdgLFa7J44cVBvPBC6D5ut59+/Wo49dQYTjklhuHDYehQzdjsCNZCZSVs3w6FhU1s2FDBP/9Zy+bNjWzbBiUlXsrL42loSG5+hBNIAVJwe3cTm1xKUu/t9Ou9mszjGsgZ3kSPAbV0y6nG47WtZl9+njBJmS9RTQ2etjQvsmbvuAO2FlGdlcryuy6mYPIpuIDs/tlk989mzFVjCBKkbFcZxZ8Xs3blWhb+YSHFq4shCEOOH8LwYcPJH5LP0KFDGTRoECkpKUdcRnvNVmkX8+ez/eq5ZPiLKPHm8t99fs17PTp2MbpgMEBj46kkJz/NuHGTSUq6kvLyZ0lO/h2XTk6jsXEFDQ2N+P1+XC4nLo8blycTpycHl9eJy+PC5XXjTnSRk+/COA1gqCiHynLYug6OvCN/Phdkn8BN/3yAHLZSk52ybx8RkSjSnPc7rwmdEe/A32UHkJScxAnnnUD+hHwslob6BorXFlO4spCn3n2c7WvuobSwlLzj8hh2/DDyh+QzfvzxDB48mNzcXJzNXZn6+tZNipZm0K5docuyZaHb6+pal+l0hmb89Ojx7YeHZWRAcnI7rBfUzmeZCQSgtPTbmzUtn5e1nhwJgM9Xh89XidNVimUHDtdW3J7NWHaS1K2WpB71pOb4SenVRHJWPGnZaaRmpZKSnYIvzYchFgc9WjXsD6cl1/td9w7ZVnkvEtVazn53zR30aDj49zkZGHd9gHHXl2NtGY2BVXz9D8sXr3soWJ7Cls25bHhqAE891XwGRoL06FHN4MFBTjwxluOP9zBwIAwYEFq3LWqE6axiNTVQUhKaYbllSx0bN1ZRWFjPtm1+du6E0lI3FRU+amoSCQZbXnMX0A3ohtNdQUz8HuJSS+jVfyPdezWQcVwjmYMCpOXVkJpVg8drvyXz2z4NuTJfop0aPIcybRqBqZfzRPkTuEzbL5MxBidO0jLSSDsvjeEThmOxBGyAipIKtq/ZzrYN23j6/acpebSEnRt34va4yeufR7++/eif159+ffrRp08fcnJyyMrK2n/aXo5stkqHaF5Nvqc/9G5HRkMhs74OnXHguzZ5rLU0NTXR1OTH72/C7/cfcGmkseVj4/7bbTCIy+3C5emOw30Re/c+REzCtTiTR+H2uIjzxOHyuHC6nRgHtDRr9sd5+86oej39BzxUcAPGaXl0zdPt+twi0oGmTeOK/55GWWAvN7z2cpuZb0xogOjz+eg/qj/9R/XHWkuQIPW19ez8cifb1m3j3S/fZf7b89n51U6q9lSR0yeHfv370S+vH8f1OY4+ffrQp08vTj89i9TU1FYDTmtDs1JaGh4tDaBvNoXWrQt97ve3/nY8nsOvE9RyiYtr4/Vo4wwibZ1lxlrYu/fws2x27rSUlkIw2DqDPZ4GYmIqcLn3YsxOghTji9tKQ10hvrgKErrVktSjjpSsRhJ7eolPjye5RzLJGckk90gmMSMNd0Lmvp9Pywyclp9ZeymYfAoX3PYIAL9fowUmRaLatGlc8vhk9gT2MOv1hW3exRiD1+Vm2Jkw7Eywdi+WPTQ2LWfTxwHWveOmcFUipZt68OHSXBYv7nPQ4+PiqsnKqiUvDwYO9DJ0aBx5eS5ycyE7G7zeY/9tHpEjzPtv09QUWuOsrCx0KS1tYuvWGrZtq2P79gZKSpooLbWUlTkoL3dTXeOjvi7+gKaNAWKbL+B0V+KNLceXVE5q/230zainW24j6XlN9BjQRHLPWpIzaomJCxymYd92A+dIKPMlmqnB8y1GjnBSuOuHpOdV/YvPZHHGQtbxliZ/gN01frZ/4ue9JY3465toamiiqWEXTY3FOJxOPB43Ho8Htycd+B3nnHMpyclXUVHxHL17P8Vtt43E5arD5XLidDpxOBwcy2NBF3w8d//pwJvFBGuZUzmb7y8YSnV1NVVVVVRXV1NZWbnvsrdsL+UV5ZRVlFFWXkZ5WTkV5RVUlFdQXVlNbGIs8cnxxKfGE58ZT1JqLL5UH7HdYknolkBiWiKJaYnEd48nPi0eT4IHYwxff/g1T0x/ifEzzueDp1/gyv+MZ9AZg47Z9/9tWlaUF5Gup6XJHxcXR98T+9L3xL5Ya7GEGj8NtQ3s2byH0oJSdhbuZM3qNZQvLKdsaxllO8rw1/tJ75lOeo90MjIyyOiRQUb3DNLT0klLSyMlJYXs7GSGDEkiKSmJ+Ph44uLicDUvBGFtaDDd1oygluuFhfDppzQ3WFp/D3FxrRtAv54/l+TagzOf2lpKfjqHSx6dSEmJkz173JSXewgEWh83Zkwjbvee0No2dgd+fzEORynxCRX4EitISK0iIb2OxIw6EjIMcalxJHZLJLF7IglpCSR0T8DXbRgO5zcOkT0GjRsRkSPR0kSOcXsYMgaGjAGoxdoCgmykvj5IwSeGgo89FK+NoUaR3AAAD1lJREFUpXRLEjt2pVOwJZM330zjm80Gr7eSxMQqUlIa6N69ifR0Q0aGg4wMFz17eunZM4bMTB/du3tITjbExbXvmRytDc0W9cyZi6uNvC//2RzuWX4mZWVNVFQEqay0VFVZqqoc1NQ4qa11U1fvpbHBRyDg+8azu4Ck5gs4XVW4YyrxxFbjSywjo3cdid0bSM7wk5oToFuvJlJy/CSm15HYve4ws24O1EWOWxY5SmrwfIvdpdBQ426HZ9o3LMXtduB2u/ElAoRW52xZo9NaS7ApQFNjgEBjgIC/iSb/CXgaLmXv3odwx/yEHbsz2LZzHYGmIMFAkGAggLUWh8Ox/+IMnZLVGAcOR0s4mtAfBmP2tYJsaKP7PzZfgtZigxZrgwSDlvTG1qeHBXBu385F0y4iJj4Gb7wXT5yHmIQYvAlePPEevKleYnvHkpGUQV5iHrHJscQmx+JL9hGTFINxHDBgP8IB/JcffskT05/guqevY+CYgQw8YyCPXf3YvusiIuF04EwSV5yLuKFx5A7NBdjX/Gn5V19TT+XOSqpLqqksqaSipIJVe1dRs7aG2r211JXVUV9ZT11lHbUVtdTX1lNfU4/b48br8xLji8Hr9eKN8eJ2h94YcLlcOJ3OfZekNENSGvQLOGhsTKKhIZn6+iQaG1Kor0+lsSGFkt2pbNvejSZ/NwKB7vzBtp35aXXbWP7ZDtwxZXhi95LRvxxfUhVxKZXEpVWFDpHKaiQ5A+JT4vAl+fAl+/Al+XD5koHkg3NeTRsRiXL7Gv0+Z/NsnyBQjbVVQDFBgvgbDaVbfBSvcbLray+7t8RQtt1HZUk8O8uT2LI1BX99SvMCwIcSxJhajKMeh6MBh6MRp9OPMQEcjiaMsRgTDL3fa0Pj/mDQgbUOgkEXwaC7+eIlGPRhgz7ASYC28z6xYhsPP5wFgMNZjdNdi8tThzumDrevHl9GPSmJjcQmNRKfGiAxrYnEHkESMwMkZgaIT2kkPrWBuORGXG57FHn/3WfdiEiIGjzfom8/iG3aw6zX3whbDV9++CWPXf0C5992Ph88vYCfPpXAgDEDsPuaQ5ZgIPROcaA+QFN9E/56f6hB1BRqFAUDQYLBUEPIBg8+5YtxmH0NIafLicPlwOl24nQ5cXldOD1Oqs5NIml7eavaarJT+cWHvzhosH7gR2j/QfuWlVsOauYMHDOQ656+ji0rt6jBIyIRraX50yI+Pp74fvHQb/99rN2f7S0fD/rcWvx1fvz1fppaZoDWNxFoChD0Bwn4D878g7bvMDgcobx3OBtwuktwuHbjcDpweV2hi8dF1dnJJG1vvRBOTXYKD61e8o2sd2Ga10NQk0ZEJKQlD504cXohZ0ADOQMADp4t05L5QWuprXRTszeG6r1easo91JZ7qSlzUVPmpK7cRX21k/oqJ/46J/56F02NDgJ+B4EmL8GgwQYPHH8HcTgtDmcQp8uP09OA2xvA7WvCG9uEN74JX0ITu5/uTnpFSav6q3qm8p8fP4E3NoDD0XqMf+R5r4aNSEdTgyeChZo7jx1+tooLvIleOEYLun12zyWHXE3+UOsTHSsTbp7Q6raBYwaquSMinULLoPnARtDBdwBPnAfaWjunnXx2z6HPEtjRmS8i0pntawQZQ0JygITkGsir6bDtrx18YZt5/9k9FxGfADoMSiT6aKQWwSJltopWkxcR6TqU+SIiXYPyXqTzUYMngkXSbBWtJi8i0nUo80VEugblvUjnogMjRURERERERESinBo8IiIiIiIiIiJRTg0eEREREREREZEopwaPiIiIiIiIiEiUU4NHRERERERERCTKqcEjIiIiIiIiIhLlwt7gMcbkGGMWG2M2GGPWGWNuDndNIiLS/pT3IiJdhzJfRKTjucJdANAE3GqtXWmMSQA+M8a8Y61dH+7CRESkXSnvRUS6DmW+iEgHC/sMHmvtDmvtyubPq4ANQFZ4qxIRkfamvBcR6TqU+SIiHS/sDZ4DGWN6AycAn7TxtWuNMSuMMStKS0s7ujQREWlHynsRka7jUJmvvBcRaV8R0+AxxsQDrwA/t9ZWfvPr1trHrbUjrbUju3fv3vEFiohIu1Dei4h0Hd+W+cp7EZH2FRENHmOMm1Dwz7fWvhruekRE5NhQ3ouIdB3KfBGRjhX2Bo8xxgBPAhustb8Ldz0iInJsKO9FRLoOZb6ISMcLe4MHOA24EjjLGLOq+TIx3EWJiEi7U96LiHQdynwRkQ4W9tOkW2v/AZhw1yEiIseW8l5EpOtQ5ouIdLxImMEjIiIiIiIiIiL/AjV4RERERERERESinBo8IiIiIiIiIiJRTg0eEREREREREZEopwaPiIiIiIiIiEiUU4NHRERERERERCTKqcEjIiIiIiIiIhLl1OAREREREREREYlyavCIiIiIiIiIiEQ5NXhERERERERERKKcGjwiIiIiIiIiIlFODR4RERERERERkSinBo+IiIiIiIiISJRTg0dEREREREREJMqpwSMiIiIiIiIiEuXU4BERERERERERiXJq8IiIiIiIiIiIRDk1eEREREREREREopwaPCIiIiIiIiIiUU4NHhERERERERGRKKcGj4iIiIiIiIhIlFODR0REREREREQkyqnBIyIiIiIiIiIS5dTgERERERERERGJcmrwiIiIiIiIiIhEOTV4RERERERERESinBo8IiIiIiIiIiJRTg0eEREREREREZEopwaPiIiIiIiIiEiUU4NHRERERERERCTKqcEjIiIiIiIiIhLl1OAREREREREREYlyEdHgMcZMMMZ8ZYzZaIyZE+56RETk2FDei4h0Hcp8EZGOFfYGjzHGCTwCnAcMBqYYYwaHtyoREWlvynsRka5DmS8i0vHC3uABTgI2WmsLrLWNwAJgUphrEhGR9qe8FxHpOpT5IiIdzBXuAoAsYOsB14uBk795J2PMtcC1ALm5uR1SWP5wi6txT4dsKxqk9aoKdwkRR6/JwfR6tJYzbA8OnEBmuEuJBJGb9/mw1V+NMaZDthcN9Pt8ML0erek1OVjOsD0ECeIgJdylRIrDZn448h5g6PAmNvl3d9j2Ip1+l1vTa3IwvR6t5QzbQ7wjHkgLdykHMdba8BZgzGTgXGvtT5qvXwmcZK2deajHjBw50q5YsaJD6qsIVODH3yHbEpHOyWu8JDgSOmRbxpjPrLUjO2RjRynS874h2ECV1QBGRL47gyHFkYLDdMwk+c6U+R2Z99ZayoJlBAl2yPZEpHOKM3H4HL4O2daR5n0kzOApBnIOuJ4NbA9TLa0kOZPCXYKISGcR0XnvdXjx4g13GSIinUXEZr4xhlRnarjLEBFpd5GwBs9yoL8xpo8xxgNcAbwW5ppERKT9Ke9FRLoOZb6ISAcL+wwea22TMWYG8BbgBJ6y1q4Lc1kiItLOlPciIl2HMl9EpOOFvcEDYK19A3gj3HWIiMixpbwXEek6lPkiIh0rEg7REhERERERERGRf4EaPCIiIiIiIiIiUU4NHhERERERERGRKKcGj4iIiIiIiIhIlDPW2nDXcNSMMaVAYQdtLg3Y3UHbkuikfUQOpyP3kV7W2u4dtK1jTnkvEUb7iBxOR+8jnSbzOzjvQb/PcnjaR+RwIm6MH5UNno5kjFlhrR0Z7jokcmkfkcPRPhId9HOSw9E+IoejfSR66Gclh6N9RA4nEvcRHaIlIiIiIiIiIhLl1OAREREREREREYlyavAc3uPhLkAinvYRORztI9FBPyc5HO0jcjjaR6KHflZyONpH5HAibh/RGjwiIiIiIiIiIlFOM3hERERERERERKKcGjwiIiIiIiIiIlFODZ4jYIyZbIxZZ4wJGmMi6jRoEj7GmAnGmK+MMRuNMXPCXY9EHmPMU8aYEmPM2nDXIkdGeS+HosyXb6O8j07KfGmL8l4OJ5IzXw2eI7MW+AGwJNyFSGQwxjiBR4DzgMHAFGPM4PBWJRHoz8CEcBchR0V5L60o8+UI/BnlfTRS5stBlPdyhP5MhGa+GjxHwFq7wVr7VbjrkIhyErDRWltgrW0EFgCTwlyTRBhr7RJgb7jrkCOnvJdDUObLt1LeRydlvrRBeS+HFcmZrwaPyHeTBWw94Hpx820iItL5KPNFRLoG5b1ENVe4C4gUxph3gYw2vjTXWvt/HV2PRDzTxm22w6sQkaOmvJfvQJkvEqWU+XKUlPcS1dTgaWatPSfcNUhUKQZyDrieDWwPUy0ichSU9/IdKPNFopQyX46S8l6img7REvlulgP9jTF9jDEe4ArgtTDXJCIix4YyX0Ska1DeS1RTg+cIGGMuNsYUA6cCC40xb4W7Jgkva20TMAN4C9gAvGytXRfeqiTSGGNeBJYBA4wxxcaYa8Jdk3w75b20RZkvh6O8j07KfPkm5b0ciUjOfGOtDikUEREREREREYlmmsEjIiIiIiIiIhLl1OAREREREREREYlyavCIiIiIiIiIiEQ5NXhERERERERERKKcGjwiIiIiIiIiIlFODR4RERERERERkSinBo+IiIiIiIiISJRTg0dEREREREREJMqpwSMiIiIiIiIiEuXU4BE5CsYYnzGm2BhTZIzxfuNrTxhjAsaYK8JVn4iItB9lvohI16C8l85CDR6Ro2CtrQPuBnKAG1tuN8b8BrgGmGmtXRCm8kREpB0p80VEugblvXQWxlob7hpEoooxxgl8AaQDecBPgP8E7rbW3hfO2kREpH0p80VEugblvXQGavCIfAfGmAuAvwHvAWcBf7TW3hTeqkRE5FhQ5ouIdA3Ke4l2avCIfEfGmM+AEcACYKrVL5OISKelzBcR6RqU9xLNtAaPyHdgjLkMyG++WqXgFxHpvJT5IiJdg/Jeop1m8IgcJWPMeEJTN/8G+IHJwDBr7YawFiYiIu1OmS8i0jUo76UzUINH5CgYY04mdEzup8B5QDawAXjDWntROGsTEZH2pcwXEekalPfSWegQLZEjZIwZBCwEvgYustY2WGs3AU8Ck4wxp4W1QBERaTfKfBGRrkF5L52JZvCIHAFjTC6wFGgERltrdx3wtUxgE/C5tVZ/AEREopwyX0Ska1DeS2ejBo+IiIiIiIiISJTTIVoiIiIiIiIiIlFODR4RERERERERkSinBo+IiIiIiIiISJRTg0dEREREREREJMqpwSMiIiIiIiIiEuXU4BERERERERERiXJq8IiIiIiIiIiIRDk1eEREREREREREopwaPCIiIiIiIiIiUe7/A0X0k/w3kx1/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 4))\n",
    "\n",
    "npoints = 2\n",
    "npoints = 5\n",
    "\n",
    "X = np.linspace(a, b, npoints)\n",
    "\n",
    "ax1.plot(xx, f(xx), lw=1, color='k')\n",
    "ax1.fill_between(x, f(x), color='lightgreen', alpha=0.9)\n",
    "\n",
    "for n in range(len(X) - 1):\n",
    "    f_mid = f(X[n:n+2].mean())\n",
    "    ax1.plot([X[n], X[n]], [0, f_mid], 'b')\n",
    "    ax1.plot([X[n+1], X[n+1]], [0, f_mid], 'b')\n",
    "    ax1.plot(X[n:n+2], [f_mid] * 2, 'b')\n",
    "    ax1.plot(X[n:n+2].mean(), f_mid, 'xk')\n",
    "    \n",
    "i = (b-a)*f_mid\n",
    "#ax1.text(-1, 7, r'$\\int_{-1}^{\\,1} f(x)dx \\approx %f$' % i, fontsize=18)\n",
    "ax1.plot(X, f(X), 'ro')\n",
    "ax1.set_xlim(xx.min(), xx.max())\n",
    "ax1.set_title('Mid-point rule')\n",
    "ax1.set_xticks([-1, 0, 1])\n",
    "ax1.set_xlabel(r'$x$', fontsize=18)\n",
    "ax1.set_ylabel(r'$f(x)$', fontsize=18)\n",
    "\n",
    "names = [\"Trapezoid rule\", \"Simpson's rule\"]\n",
    "for idx, ax in enumerate([ax2, ax3]):\n",
    "    ax.plot(xx, f(xx), lw=1, color='k')\n",
    "    ax.fill_between(x, f(x), color='lightgreen', alpha=0.9)\n",
    "\n",
    "    i = 0\n",
    "    for n in range(len(X) - 1):\n",
    "        XX = np.linspace(X[n], X[n+1], idx+2)\n",
    "\n",
    "        f_interp = polynomial.Polynomial.fit(XX, f(XX), len(XX)-1)    \n",
    "        ax.plot([X[n], X[n]], [0, f(X[n])], 'b')\n",
    "        ax.plot([X[n+1], X[n+1]], [0, f(X[n+1])], 'b')\n",
    "        XXX = np.linspace(X[n], X[n+1], 50)\n",
    "        ax.plot(XXX, f_interp(XXX), 'b')\n",
    "        \n",
    "        F = f_interp.integ()\n",
    "        i += F(X[n+1]) - F(X[n])\n",
    "    ax.text(-1, 5.5, r'$\\int_a^{\\,b} f(x)dx \\approx %f$' % i, fontsize=18)\n",
    "    ax.plot(X, f(X), 'ro')\n",
    "    ax.set_xlabel(r'$x$', fontsize=18)\n",
    "    ax.set_ylabel(r'$f(x)$', fontsize=18)\n",
    "    ax.set_xlim(xx.min(), xx.max())\n",
    "    ax.set_title(names[idx])\n",
    "    ax.set_xticks([-1, 0, 1])\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch8-quadrature-rules-%d.pdf' % npoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# mid-point rule\n",
    "(b-a) * f((b+a)/2.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# trapezoid rule\n",
    "(b-a)/2.0 * (f(a) + f(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.333333333333333"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# simpsons rule\n",
    "(b-a)/6.0 * (f(a) + 4 * f((a+b)/2.0) + f(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.066666666666667"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# exact result\n",
    "integrate.quad(f, a, b)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.3125"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate.trapz(f(X), X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.083333333333332"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate.simps(f(X), X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.066666666666666"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate.quadrature(f, a, b)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.066666666666666"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate.romberg(f, a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.33333333, 1.33333333, 0.33333333])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate.newton_cotes(2)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGGBJREFUeJzt3XmQHOdZx/Hvs6vbB7LjxchW5LXjYEhBYYWtEHBxxAnOAaWYQoCCAwkERMINcYFt6Q+gcEEoiqtIASKEJLCJIQYXBhKCwVFSVMWGdWzHsuVDh7U6VtrZlbQ7u3N3P/wxPfZoNavp2Z2emZ75faq2Zqa7p/vZd1o/vftOH+buiIhIegx1uwAREWmNgltEJGUU3CIiKaPgFhFJGQW3iEjKKLhFRFJGwS0ikjIKbhGRlFFwi4ikzJokVnrNNdf46OhoEqsWEelLTzzxxIy7j8RZNpHgHh0dZWJiIolVi4j0JTM7FndZDZWIiKSMgltEJGUU3CIiKaPgFhFJGQW3iEjKNA1uM7vFzJ6q+5k3s1/tRHEiIqkwPg6jozA0VH0cH090c00PB3T3F4BbAcxsGDgJPJRoVSIiaTE+Drt3Qy5XfX3sWPU1wF13JbLJVodK3gocdvfYxxuKiPS1PXteDe2aXK46PSGtBvcu4DONZpjZbjObMLOJTCaz+spERNJgcrK16W0QO7jNbB2wA/hso/nuvs/dx9x9bGQk1lmbIiLpt21ba9PboJUe9zuBr7r7maSKERFJnfvvxzdtunDapk1w//2JbbKV4H4PywyTiIgMrLvu4swf/hmF67biZnDDDbBvX2JfTELMi0yZ2Sbg+4GfS6wSEZEUWihWOHrHnXDHndz62s1sXDec+DZjBbe754DXJFyLiEjqTM7mmi/UZjpzUkRkhc4tlpjLlzu+XQW3iMgKhKHz8uxiV7at4BYRWYHT8wUK5bAr21Zwi4i0qByEnDyf79r2FdwiIi06fjZHJfCubV/BLSLSgsVihelssas1KLhFRFpwdGYR715nG1Bwi4jENp0tkC1Uul2GgltEJI5KEHL8bOdPtmlEwS0iEsOJc3lKlS6PkUQU3CIiTSwWK5yeL3S7jFcouEVELsHde+ILyXoKbhGRS5jOFnviC8l6Cm4RkWWUKiGTPfKFZD0Ft4jIMo7NLnb1DMnlKLhFRBo4nysxs1DqdhkNKbhFRJYIQufITHcu2RqHgltEZInJszmKXbpkaxwKbhGROvOFMqfneueY7UZiBbeZbTazB83seTM7aGbfmXRhIiKdFobOkUzvDpHUxLpZMPCnwH+4+04zWwdsSrAmEZGuOH4uR74UdLuMppoGt5ldCXwP8H4Ady8BvflVq4jICs0Xypw639tDJDVxhkpuAjLA35rZk2b2MTO7bOlCZrbbzCbMbCKTybS9UBGRpAShc3h6odtlxBYnuNcAbwT+wt23A4vAPUsXcvd97j7m7mMjIyNtLlNEJDnHz+a6duPflYgT3CeAE+7+ePT6QapBLiKSenP5MlM9fhTJUk2D291PA8fN7JZo0luB5xKtSkSkAypByKEUDZHUxD2q5JeA8eiIkiPATyVXkohIZxydWaRUSc8QSU2s4Hb3p4CxhGsREemYTLbYs9ciaUZnTorIwCmUA4728LVImlFwi8hAcXcOTS8QhL13uda4FNwiMlCOn8333B1tWqXgFpGBcT5X4uT5fLfLWDUFt4gMhGIlSOWhf40ouEWk79XGtcs9eBuylVBwi0jfO342z3w+3ePa9RTcItLXzi72x7h2PQW3iPStQjngcKY/xrXrKbhFpC8FofPC6SyVPhnXrqfgFpG+dDizQC4Fd7NZCQW3iPSdk+fzzKb0OiRxKLhFpK+cWywxOZvrdhmJUnCLSN/IlSq81Ccn2VyKgltE+kI5CHn+dDbVF4+KS8EtIqkXRkeQFFN038jVUHCLSOodmVlI/RX/WqHgFpFUO342Rybbv0eQNKLgFpHUOjNf4MS5/jqdPY5Y95w0s5eBLBAAFXfX/SdFpKvOLZZSffux1Yh7l3eAt7j7TGKViEhP2b59O5lMhptvvrnbpVwkCJ1cKcC9d44gOXHsKFu+4VqefurJxLeloRIRaSiTybCw0HvHRIfee6ENkM8tMjOT6ci24va4HfhPM3Pgr9x9X4I1iUgPqPW09+/f391C6hTKAc+emqNU6a3QBviFH9/B5RtaGcRYubhbuc3dT5nZ1wOPmNnz7v7l+gXMbDewG2Dbtm1tLlNEBl2xEnBwar4nQ7vTYg2VuPup6HEaeAh4U4Nl9rn7mLuPjYyMtLdKERlopUrIwakshQE5waaZpsFtZpeZ2RW158AdwIGkCxMRgdqp7PPk+/QSrSsRZ6jkWuAhM6st/2l3/49EqxIRASpByPNTWRaLCu16TYPb3Y8A39aBWkREXlGOQnuhODinssfVma9ARURaUA5CDk7Nq6e9DB3HLSI9RaHdnHrcItIzqof8ZfVFZBMKbhHpCYVywHNT8wNzTe3VUHCLSNflShWdXNMCBbeIdNVcvsyLZ7JUAoV2XApuEema2YUih6YXGIDbRLaVgltEuuL0XGFgr6e9WgpuEekod+fYbI6puUK3S0ktBbeIdEwQOi9NZzm3WO52Kamm4BaRjiiUA148o+uOtIOCW0QSN5cv89KZLGUdOdIWCm4RSdTpuQIvzy7SY3caSzUFt4gkIgidozMLZLKlbpfSdxTcItJ2Gs9OloJbRNpqZqHIkcwigc6qSYyCW0TaIgydl2cXOTNf7HYpfU/BLSKrlitVeOnMAjldjrUjFNwismLuzun5ApOzOV1vpINiB7eZDQMTwEl3/8HkShKRNChWAg5PLzKX11mQndZKj/tXgIPAlQnVIiIpMT1f4NjZnC7F2iWx7jlpZluBHwA+lmw5ItLLipWA50/PczizqNDuorg3C/4T4DcA3VNIZBCMj8Njj8GXvgSjo/j4OKfnCjx9fE4XiOoBTYPbzH4QmHb3J5ost9vMJsxsIpPJtK1AEemw8XHYvRuK0WF9x47hP/uzzH/8kzo2u0fE6XHfBuwws5eBB4Dbzezvly7k7vvcfczdx0ZGRtpcpoh0zJ49kMtdMGkon2fbH/5ulwqSpZoGt7vf6+5b3X0U2AU86u7vTbwyEemOycmGk9dPnexwIbKcuGPcIjIAFooVStdtbTivuOX6Dlcjy2kpuN19v47hFuk/xUrAoeksz5yY4+UP7yHYsPGC+cGGjUzevbdL1clSOnNSZIBVgpBT5wtMzeVfOfNxdsdOAMJ7fpmhUonCdVuZvHvvK9Ol+xTcIgMoCKunqp86n294PPbsjp0sPPApAJ789MOdLk+aUHCLDJAwdKazRU6ez1Gq6NC+tFJwiwwABXZ/UXCL9LFKEHImW+T0XF6B3UcU3CJ9qFQJOT1X4Ey2oGuK9CEFt0gfWSxWmJorMLNQ1F3V+5iCWyTl3J1zuTJTc3nm85VulyMdoOAWSalSJWQ6W+DMfJFSRRfuHCQKbpEUcXfm8mWms0XOLpY0HDKgFNwiKVAoB2SyRaaz6l2LglukZ5WDkLOLJTLZItmCxq7lVQpukR4ShM65XImZhSLnc2UNhUhDCm4ZGNu3byeTyXDzzTd3u5QLONUTZSqBUw7D6oQecODJCdasVUT0In0qMjAymQwLCwvdLgMAd6iEIeXAqfRQWNcLwgB0e8mepOCWgVHrae/fv78r28+VKpzLlTm3WErFmPUd22/qdgmyDAW3SEIqQch8ocL5XIlzubKOBpG2UXCLtIm7ky1WmMuVmcuXWShW9OWiJELBLbJC7s5iKWA+X2a+UGY+XyEIldSSPAW3SExB6CwUK2QLZbKFCtmCglq6o+nNgs1sg5n9r5k9bWbPmtlvJ1LJ+DiMjsLQUPVxfDyRzciAGh+Hxx6DL30p9v5VKAfMLBR5eWaRAyfn+L+Xz/LcqXmOn81zPldWaEvXxOlxF4Hb3X3BzNYC/2Nmn3f3x9pWxfg47N4NuVz19bFj1dcAd93Vts3IgKrtX8Vi9XWD/aschCwWKywUKywWAxaKZd14QHpW0+B2dwdqB7+ujX7au0fv2fNqaNfkctXpCm5ZrWX2r8o993LkbTtYKFYolnXEh6RH06ESADMbNrOngGngEXd/vMEyu81swswmMplMa1VMTrY2XSSGMHRypQq+zH40fPIEswslhbakTqwvJ909AG41s83AQ2b2Le5+YMky+4B9AGNjY631yLdtq/75ukTl+q2UShU2rdN3qLI8d6dQDsmXA3KlCvlSQK4UkC8HuMP2Ldez4dSJi95X3HJ9F6oVWb2WEtHdz5vZfuAdwIEmi8d3//0XjnEDwYaNHPn1Pcwen2PjumGu3rSOzZet5Yr1azCztm1a0iMInXw5IF8KKJSDC55f6nvCybv38rr7fg0K+VfXtWEjk3fv7UDVIu3XNLjNbAQoR6G9EXgb8JG2VhGNY4f33YcdP05xy/VM3r2X2R07AciXAk6W8pw8n2ftsLF501qu3LiWzRvXsW5NrNEeSYkwdAqVgEI5fCWcC+Xq65WeeVjbj8J7fpmhUonCdVsv2L9E0iZOj3sL8EkzG6Y6Jv6P7v5vba/krrvI7fwxnjkxd8nFyoGTyZbIZEvAIpvWDXPlxrVcuWENV2xYqyDvce5OOaiGc7EcUqwL6WJl5eHczOyOnSw88CkAnvz0w4lsQ6RT4hxV8jVgewdqWZFcNJ55Osr7DWuHuGLDGi5fv5bL1g9z2bo1DA3159BKr16mNHQn9GpIh1597V6b7l27Ep4uUyr9ou/24mrvrdYjBzPYtG44+lnDZevWsGHdEOvXDHe50lUaHyfzzDMsBAHMzMCNN8K11ya6yWrmXhjKjR57lS5TKv2ip4L7tu8YY+r0NFtvuDHxbZkZQwZDQ8ZQ7bkZFj32tDNn4MUXyYQha4H9xWL10Mm9e1s+7r0ShFRCpxxE14YOQkrR83JQHbqozUv7mYK6TKn0i54K7pnMDPncYke25e4ETuMwqoU4r4Z5LeiN2utXn3fc0aMQhgT103I5wvvuY/GHf5QgdCqhv/oY3VkliAK6ElSnV4LwkkdjiEhv6qngvul1r2OxWOGjKfvyaMhgzXC15z4c9eBrPfdab74a8FHoN1hHLT+rIw21MeELx4fdIXDn2264GgM2L1mHHT/OgZPzyf2iItITeiq40yp0outadKb7WtQJJSIDTcfOpdDk3XsJNmy8YJpOKBEZHOpxp9ArJ458+EPgrhNKRAaMetwpNbtjJ8HlVxBccSVPfvkphbbIAFFwi4ikjIJbRCRlFNwiIimj4BYRSRkFt4hIyii4RURSRsEtIpIyCm4RkZRRcIuIpIyCW0QkZRTcIiIp0zS4zey1ZvZFMztoZs+a2a90ojAREWksztUBK8CH3f2rZnYF8ISZPeLuzyVcm4iINNC0x+3uU+7+1eh5FjgI6Ir9IiJd0tIYt5mNAtuBxxvM221mE2Y2kclk2lOdiIhcJHZwm9nlwD8Bv+ruF93Y0N33ufuYu4+NjIy0s0YREakTK7jNbC3V0B53939OtiQREbmUOEeVGPA3wEF3/6PkSxIRkUuJ0+O+DfgJ4HYzeyr6eVfCdYmIyDKaHg7o7v8DWAdqERGRGHTmpIhIyii4RURSRsEtIpIyCm4RkZRRcIuIpIyCW0QkZRTcIiIpo+AWEUkZBbeISMoouEVEUkbBLSKSMgpuEZGUUXCLiKSMgltEJGUU3CIiKaPgFhFJGQW3iEjKKLhFRFJGwS0ikjJx7vL+cTObNrMDnShIREQuLU6P+xPAOxKuQ0REYmoa3O7+ZeBsB2oREZEYNMYtIpIybQtuM9ttZhNmNpHJZNq1WhERWaJtwe3u+9x9zN3HRkZG2rVaERFZQkMlIiIpE+dwwM8AXwFuMbMTZvaB5MsSEZHlrGm2gLu/pxOFiIhIPBoqERFJGQW3iEjKKLhFRFJGwS0ikjIKbhGRlFFwi4ikjIJbRCRlFNwiIimj4BYRSRkFt4hIyii4RURSRsEtIpIyCm4RkZRRcIuIpIyCW0QkZRTcIiIpo+AWEUkZBbeISMoouEVEUiZWcJvZO8zsBTM7ZGb3JF2UiIgsL85d3oeBjwLvBN4AvMfM3pB0YSIi0licHvebgEPufsTdS8ADwLuTLUtERJazJsYy1wPH616fAL4jmXLAzFg7bEmtvi+pvVqj9mqN2ises861U5zgblSNX7SQ2W5gN8C2bdtWVMy3v3E7AGOjV6/o/YPm5tfdBKi94lJ7tUbt1ZrvfvNYx7Zl7hdl8IULmH0n8Fvu/vbo9b0A7v57y71nbGzMJyYm2lmniEhfM7Mn3D1W+scZ4/4/4PVmdqOZrQN2AQ+vpkAREVm5pkMl7l4xs18EvgAMAx9392cTr0xERBqKM8aNu38O+FzCtYiISAw6c1JEJGUU3CIiKaPgFhFJGQW3iEjKKLhFRFKm6Qk4K1qpWQY4tsK3XwPMtLGcdlFdrVFdrVFdrenHum5w95E4CyYS3KthZhNxzx7qJNXVGtXVGtXVmkGvS0MlIiIpo+AWEUmZXgzufd0uYBmqqzWqqzWqqzUDXVfPjXGLiMil9WKPW0RELqErwW1mP2Jmz5pZaGbLfgO73E2Ko0vMPm5mL5nZP0SXm21HXVeb2SPReh8xs6saLPMWM3uq7qdgZndG8z5hZkfr5t3aqbqi5YK6bT9cN72b7XWrmX0l+ry/ZmY/Vjevre3V7KbWZrY++v0PRe0xWjfv3mj6C2b29tXUsYK6ft3Mnova57/N7Ia6eQ0/0w7V9X4zy9Rt/2fq5r0v+txfMrP3dbiuP66r6UUzO183L5H2MrOPm9m0mR1YZr6Z2Z9FNX/NzN5YN6/9beXuHf8Bvhm4BdgPjC2zzDBwGLgJWAc8DbwhmvePwK7o+V8CH2pTXX8A3BM9vwf4SJPlrwbOApui158AdibQXrHqAhaWmd619gK+EXh99Pw6YArY3O72utT+UrfMzwN/GT3fBfxD9PwN0fLrgRuj9Qx3sK631O1DH6rVdanPtEN1vR/48wbvvRo4Ej1eFT2/qlN1LVn+l6heajrp9voe4I3AgWXmvwv4PNU7hr0ZeDzJtupKj9vdD7r7C00Wa3iTYjMz4HbgwWi5TwJ3tqm0d0fri7vencDn3T3Xpu0vp9W6XtHt9nL3F939pej5KWAaiHWSQYvi3NS6vt4HgbdG7fNu4AF3L7r7UeBQtL6O1OXuX6zbhx4DtrZp26uq6xLeDjzi7mfd/RzwCPCOLtX1HuAzbdr2stz9y1Q7act5N/Apr3oM2GxmW0iorXp5jLvRTYqvB14DnHf3ypLp7XCtu08BRI9f32T5XVy809wf/an0x2a2vsN1bTCzCTN7rDZ8Qw+1l5m9iWov6nDd5Ha113L7S8NlovaYo9o+cd6bZF31PkC151bT6DPtZF0/HH0+D5rZa1t8b5J1EQ0p3Qg8Wjc5qfZqZrm6E2mrWDdSWAkz+y/gGxrM2uPu/xJnFQ2m+SWmr7quuOuI1rMF+FaqdwaquRc4TTWc9gG/CfxOB+va5u6nzOwm4FEzewaYb7Bct9rr74D3uXsYTV5xezXaRINpS3/PRPapJmKv28zeC4wB31s3+aLP1N0PN3p/AnX9K/AZdy+a2Qep/rVye8z3JllXzS7gQXcP6qYl1V7NdHTfSiy43f1tq1zFCeC1da+3AqeoXgdgs5mtiXpNtemrrsvMzpjZFnefioJm+hKr+lHgIXcv1617KnpaNLO/Be7uZF3RUATufsTM9gPbgX+iy+1lZlcC/w7sjf6MrK17xe3VwHL7S6NlTpjZGuDrqP75G+e9SdaFmb2N6n+G3+vuxdr0ZT7TdgRR07rcfbbu5V8DH6l77/ctee/+NtQUq646u4BfqJ+QYHs1s1zdibRVLw+VNLxJsVdH/L9IdXwZ4H1AnB58HA9H64uz3ovG1qLwqo0r3wk0/AY6ibrM7KraUIOZXQPcBjzX7faKPruHqI7/fXbJvHa2V5ybWtfXuxN4NGqfh4FdVj3q5Ebg9cD/rqKWluoys+3AXwE73H26bnrDz7SDdW2pe7kDOBg9/wJwR1TfVcAdXPiXZ6J1RbXdQvXLvq/UTUuyvZp5GPjJ6OiSNwNzUcckmbZK4hvYZj/AD1H9n6gInAG+EE2/Dvhc3XLvAl6k+j/mnrrpN1H9h3UI+Cywvk11vQb4b+Cl6PHqaPoY8LG65UaBk8DQkvc/CjxDNYD+Hri8U3UB3xVt++no8QO90F7Ae4Ey8FTdz61JtFej/YXq0MuO6PmG6Pc/FLXHTXXv3RO97wXgnW3e35vV9V/Rv4Na+zzc7DPtUF2/Bzwbbf+LwDfVvfeno3Y8BPxUJ+uKXv8W8PtL3pdYe1HtpE1F+/IJqt9FfBD4YDTfgI9GNT9D3dFySbSVzpwUEUmZXh4qERGRBhTcIiIpo+AWEUkZBbeISMoouEVEUkbBLSKSMgpuEZGUUXCLiKTM/wOUfR12soFRnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.fill_between(x, f(x), alpha=0.25)\n",
    "ax.plot(X, f(X), 'ro')\n",
    "\n",
    "for n in range(len(X) - 1):\n",
    "    f_mid = f(X[n:n+2].mean())\n",
    "    ax.plot([X[n], X[n]], [0, f_mid], 'k')\n",
    "    ax.plot([X[n+1], X[n+1]], [0, f_mid], 'k')\n",
    "    ax.plot(X[n:n+2], [f_mid] * 2, 'k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8VOW97/HPLyEQbnK/5iJW0aqtgocjrVq11aqAAr4KW4L7FKxH6q3q3lXE0nppT7eX2r1bL8dqbbU9G4KIsoUiWqxc2tcWK3ipiPXCnSSQkITcJjOZy3P+mEk6kJlkApkkK/m+X6+8MrPWs9b8eDJ8s/LMs9Yy5xwiIuIdGZ1dgIiItI2CW0TEYxTcIiIeo+AWEfEYBbeIiMcouEVEPEbBLSLiMQpuERGPUXCLiHhMr3TsdPjw4W7cuHHp2LWISLe0devWQ865Eam0TUtwjxs3ji1btqRj1yIi3ZKZ7Um1rYZKREQ8RsEtIuIxCm4REY9RcIuIeIyCW0TEY1oNbjM7zczej/uqNrM7OqI4ERFPWLIExo2DjIzo9yVL0vpyrU4HdM59AkwAMLNMoAhYmdaqRES8YskSWLAAfL7o8z17os8Brr02LS/Z1qGSS4AdzrmU5xuKiHRrixeDz8c/ARcCDqIhvnhx2l6yrSfgzAEKE60wswXAAoD8/PzjLEtExCP27uUTosMQWUAdMCC2PF1SPuI2s97AdODFROudc8845yY55yaNGJHSWZsiIp53MCeHKYABXyYW2gBpPIBty1DJFOBd59zBdBUjIuIldXV1XJmVxQHgS0DfxhX9+sFPf5q2121LcBeQZJhERKSnCYVCzJkzh3f37OFXN97KgN69oytOPBGeeSZtH0xCimPcZtYP+Cbw3bRVIiLiEc45br31Vv7whz/wH489zvipBdS++y4DsnvBxo1pf/2Ujridcz7n3DDnXFW6CxIR6eoeeughnn76ae6++24u+9a3O/z1deakiEgbLFmyhB/84AfMnTuXuxbfT1V9sMNrUHCLiKTozTff5LrrruPiiy/m2Wd/w97K+k6pQ8EtIpKCDz/8kKuvvppTTz2VlStXUhlw+IORTqlFwS0i0or9+/czdepUBgwYwKuvvkr/gSdQdLhzjrYhTbcuExHpLqqrq5k2bRpVVVVs2rSJ/Px8dpbVEgq7TqtJwS0ikkRDQwPf+ta32L59O6+++ioTJkygLhCitCbQqXUpuEVEEnDOccMNN/DGG2/w3HPP8c1vfhOAXYfqcJ13sA1ojFtEJKH77ruP3//+9zzwwAPMnz8fgNIaPzX+UOcWhoJbRKSZX//61/zkJz/h+uuv50c/+hEAoXCEfRW+Tq4sSsEtIhJn7dq13HTTTVxxxRU89dRTmBkA+yvraQh18hhJjIJbRCRm69atzJ49m7POOovly5eTlZUFQF0gxIFqfydX9w8KbhERYPfu3UybNo3hw4ezZs0aBg4cCEQ/pOwKH0jG06wSEenxKioquOKKKwgEAqxfv54xY8Y0rSutCXSJDyTjKbhFpEfz+/3MmDGDXbt2sW7dOk4//fSmdQ2hCHu7yAeS8RTcItJjRSIR5s2bx1/+8heWLVvGhRdeeMT6PeV1nXqGZDIa4xaRHmvhwoUsX76cn/3sZ1xzzTVHrDvsa+BQbUMnVdYyBbeI9EiPPfYYP//5z7n11lv5/ve/f8S6cMSx81BdJ1XWOgW3iPQ4L7/8MnfccQczZ87kF7/4RdNc7UZ7K3wEOumSralQcItIj/Lf//3fXHvttUyePJmlS5eSmZl5xPpqf5ADVV1nznYiKQW3mQ02sxVm9ncz+9jMvpruwkRE2tunn37K9OnTyc3NZdWqVfTt2/eI9ZGIY2dZ1x0iaZTqrJJfAq8552aZWW+gXxprEhFpd6WlpUyZMgUzY+3atYwYMaJZm32VPuobwp1QXdu0GtxmdgJwITAfwDnXAHTNj1pFRBKoq6vjyiuvpKSkhPXr13PKKac0a1PtD1J8uGsPkTRKZajkC0AZ8JyZvWdmz5pZ/6MbmdkCM9tiZlvKysravVARkWMRCoUoKChg69atLFu2jMmTJzdrE444dpTWdkJ1xyaV4O4FnAM85ZybCNQBi45u5Jx7xjk3yTk3KdGfICIiHc05x2233cbq1at5/PHHmT59esJ2+yp8nXbj32ORSnDvB/Y7596OPV9BNMhFRLq0hx9+mKeeeoqFCxdy8803J2xTVR+kpIvPIjlaq8HtnDsA7DOz02KLLgG2p7UqEZHjtHTpUu655x4KCgp48MEHE7YJhSN87qEhkkapzir5HrAkNqNkJ3Bd+koSETk+69evZ/78+Vx88cU899xzZGQkPkbddaiOhpB3hkgapRTczrn3gUlprkVE5Lht27aNq6++mvHjx7Ny5Ur69OmTsF1ZTaDLXoukNTpzUkS6jaKiIqZOnUq/fv1Yu3YtgwcPTtjOHwyzqwtfi6Q1uqyriHQL1dXVTJ06lcrKSv785z+Tn5+fsJ1zjs9LawlHut7lWlOl4BYRzwsGg8yaNYvt27ezZs0aJkyYkLTtvor6LndHm7ZScIuIpznnuOGGG1i3bh3PPfccl112WdK2h30NFB2u78Dq0kNj3CLiaffffz+/+93vuP/++5k/f37SdoFQ2JNT/xJRcIuIZz377LP8+Mc/5jvf+Q733ntv0naN49rBLngbsmOh4BYRT1q7di033ngjl19+Ob/61a+a3Qwh3r6KeqrrvT2uHU/BLSKe8+677zJ79mzOOussXnzxRbKyspK2rajrHuPa8RTcIuIpu3fvZtq0aQwbNow1a9YwcODApG39wTA7yrrHuHY8zSoREc+oqKhgypQp+P1+/vSnPzFmzJikbcMRxycHagh1k3HteApuEfEEv9/PzJkz2blzJ3/84x8544wzWmy/o6wWnwfuZnMsFNwi0uVFIhHmzZvHn//8ZwoLC7noootabF90uJ5yj16HJBUa4xaRLu/uu+9m+fLlPPLII8yZM6fFtpV1Dewt93VQZZ1DwS0iXdrjjz/Oo48+yi233MKdd97ZYltfQ4jPuslJNi1RcItIl7Vy5Upuv/12ZsyYwS9/+csW52oHwxH+fqDG0xePSpWCW0S6pLfeeou5c+cyefJkli5dSmZmZtK2kdgMkoCH7ht5PBTcItLlfPrpp1x11VXk5uayatUq+vXr12L7nYdqPX/Fv7ZQcItIl1JaWsqUKVMwM9auXcuIESNabL+vwkdZTfedQZKIpgOKSJdRV1fHVVddRUlJCevXr+eUU05psf3Baj/7K7vX6eypSCm4zWw3UAOEgZBzTvefFJF2FQ6HmTt3Llu2bOHll19m8uTJLbavrGvw9O3Hjkdbjri/7pw7lLZKRKTHcs7xve99j1WrVvHEE08wY8aMFttX+4N8VlqL6/4TSBLSGLeIdLpHHnmEp556irvuuotbbrmlxba+hhCf9JBpf8mkGtwO+KOZbTWzBeksSER6lqVLl7Jo0SLmzJnDQw891GJbfzDMxyXV3fLCUW2R6lDJ+c65YjMbCawzs7875zbFN4gF+gIg6d2VRUTibdiwgfnz53PRRRfx/PPPk5GR/FgyEIqGdkOoZ4c2pHjE7Zwrjn0vBVYC5yZo84xzbpJzblJr03dERD766CNmzpzJ+PHjWblyJX369EnatiEU4eOSGvw95ASb1rQa3GbW38wGNj4GLgO2pbswEem+iouLmTJlCn379uXVV19lyJAhSdtGT2Wvpr6bXqL1WKQyVDIKWBm7RkAvYKlz7rW0ViUi3VZ1dTVTp06lsrKSTZs2ceKJJyZtGwpH+HtJDXUBhXa8VoPbObcTOLsDahGRbi4YDDJr1iy2bdvGmjVrmDhxYvK2sdCuDfScU9lTpTMnRaRDOOdYsGAB69at47e//S2XX3550rbBcISPS6p1pJ2E5nGLSId44IEHeP7557nvvvu47rrrkrZTaLdOR9wikna/+c1veOCBB7juuuu47777kraLTvmr0QeRrVBwi0havfbaa3z3u9/lsssu4+mnn056MwR/MMz2kuoec03t46HgFpG0ee+995g9ezZf/vKXWbFiBVlZWQnb+RpCOrmmDRTcIpIWe/bsYerUqQwdOpQ1a9YwcODAhO2q6oN8erCmx5/G3hYKbhFpd5WVlUyZMoX6+nreeOMNxo4dm7BdeW2Az0tr6cHXizomCm4RaVeBQICZM2eyY8cOXn/9dc4888yE7Q5U+Xvs9bSPl4JbRNpNJBJh3rx5bNq0icLCQi6++OJmbZxz7Cn3UVLl7/gCuwkFt4i0m0WLFvHCCy/w8MMPM2fOnGbrwxHHZ6U1VNYFO6G67kPBLSLt4oknnuBnP/sZN998M3fddVez9f5gmE8P6roj7UHBLSLH7b/+67+47bbbmDFjBo899lizudpV9UE+O1hDUDNH2oWCW0SOy+bNmykoKODcc89l6dKlZGZmHrH+QJWf3eV1Pfb+kOmg4BaRY/bZZ59x1VVXkZOTw+rVq+nXr1/TunDEsetQLWU1DZ1YYfek4BaRY1JaWsqUKVMAWLt2LfF3vtJ4dnopuEWkzXw+H1dddRVFRUWsX7+e8ePHN607VBtgZ1ldj74Le7opuEWkTcLhMHPnzuWdd97h5Zdf5itf+QoAkYhjd3kdB6sDnVxh96fgFpGUOee4/fbbeeWVV3j88ceZOXMmEL1I1GcHa/HpcqwdQsEtIil79NFHefLJJ7nzzju59dZbcc5xoNrP3nKfrjfSgVIObjPLBLYARc65K9NXkoh0RYWFhSxcuJBrrrmGhx9+mEAozI7SOqrqdRZkR2vLEfftwMfACWmqRUS6qI0bNzJ//nwuvPBCnn/+eQ7VNrCnwqdLsXaSlO45aWa5wDTg2fSWIyJdzfbt25k5cyYnn3wyL6x4id2HG9hRVqfQ7kSp3iz4F8BCQPcUEukJliyBceMoNmPKWWeR7Ry/e2Ele2szdIGoLqDV4DazK4FS59zWVtotMLMtZralrKys3QoUkQ62ZAksWED1nj1cDpSHw6z2+xn0xgbNze4iUjniPh+Ybma7gWXAN8zsP49u5Jx7xjk3yTk3Kf4MKhHxDuccH9x1F4t8PgYD24AVwKRAgPxH/08nVyeNWv1w0jl3D3APgJldDNzpnPvnNNclIh3o888/Z9myZRQWFrK9pIRMwIgGxOWxNn1KijqvQDmC5nGL9FDFxcW88MILFBYW8s477wBw3vkX8NigIcypqqTxJPbGC7QGxuR0Sp3SXJuC2zm3AdiQlkpEJO3Ky8t56aWXKCwsZOPGjTjnmDhxIv/20EOc/82ryBo0imGrVjD0B/8C/vqm7cLZfdl75w87sXKJpyNukW6utraWV155hcLCQl5//XVCoRCnnnoq9957L7P/6RoGjjqRkqr6pjMfy6fPij74/k3gHP6xuey984f/WC6dTsEt0g0FAgHWrl1LYWEhq1evpr6+ntzcXO644w4KCgo46+wJHKwJUHy4nurD9c22L58+i/B9CwF4b9P7HV2+tELBLdJNhMNh1q9fT2FhIS+99BJVVVUMHz6c+fPnU1BQwPnnnw8YpTUBPth/mIaQpvZ5lYJbxMOcc2zevJnCwkKWL1/OwYMHGThwIFdffTUFBQVccsklZGVlEYk4SmsCFB32KbC7AQW3iMc45/jwww8pLCxk2bJl7N69mz59+jBt2jTmzp3L1KlT6du3LwChcISiw/UcqKpXYHcjCm4Rj9ixYweFhYXRudbbt5OZmcmll17K/fffz8yZMxk0aFBT24ZQhANVfg7W+HVNkW5IwS3ShRUXF7N8+XIKCwv561//CsAFF1zAk08+yezZszn6LOW6QIiSKj+HagO6q3o3puAW6WIqKiqa5lpv2LChaa71I488wjXXXEN+fv4R7Z1zVPqClFTVU10f6qSqpSMpuEW6gNraWlatWtU01zoYDDbNtS4oKOC0005rtk1DKEJpjZ+D1QEaQrpwZ0+i4BbpJIFAgNdee61prrXP5yM3N5fbb7+dgoICJk6ciJkdsY1zjqr6IKU1ASrqGjQc0kMpuEU6UDgcZsOGDU1zrQ8fPsywYcOYN29e01zrjIzmF+30B8OU1QQordHRtSi4RdLOOcfbb7/dNNf6wIEDDBgwoGmu9aWXXkpWVlaz7YLhCBV1DZTVBKjxa+xa/kHBLZIm8XOtd+3a1TTXuqCggGnTpjXNtY4XjjgqfQ0cqg1w2BfUUIgkpOAWaUc7d+5smmv90UcfNc21vu+++5rNtW7UGNYVdQ1U1jWgm8xIaxTcIseppKSk6brWR8+1njVrFiNHjmy2TTAcaQrrKl9QYS1touAWOQaVlZVNc63Xr1/f6lxrAF9DiEpfkMq6Bo1Zy3FRcIukqK6urmmu9WuvvUYwGGT8+PHce++9zJkzhy9+8YtHtA+FI1T7Qxz2NVDpC2o2iLQbBbdICxoaGprmWq9atQqfz0dOTg633XYbBQUFnHPOOU1zrZ1z1ARCVPmCVNUHqQ2E9OGipIWCW+QoyeZaf/vb36agoIALLriAjIwMnHPUNYSprg9S7Q9SXR8irMFq6QAKbhFSm2udkdmL2kCI4io/Nf4QNX4FtXSO5qdoHcXMss3sr2b2gZl9ZGYPpKWSJUtg3DjIyIh+X7IkLS8jPVSS99e2bdtYvHgxJ598Ml/96ld5+umnOe+883jxxRfZW1TCv//fX3P6/7yQT0p9vLO7gu3F1eyrqOewL6jQlk6TyhF3APiGc67WzLKAv5jZWufc5narYskSWLAAfL7o8z17os8Brr223V5Geqij3l879+xh2XXXUbhoEdv27yczM5NLLrmEuxYt5qLLp5LZZwC1gSB/LwsQffuLdC2tBrdzzgG1sadZsa/2PdRYvBh8PiYBvYGbgTyfj/yFC8mZPZvevXu368tJz1BXV8e+ffvY+6//yj6fj38hGsMNAMEg55eWcu+DP+drl11J/0HDAKgKAaGGzitaJAUpjXGbWSawFTgFeNI593aCNguABUDCOawt2ruXCPA3IAi81bi8uBjLzmb06NHk5+c3feXl5R3xeMSIEc2uoibdWygUoqSkJBrMe/c2fcU/r6ioSLjtw8AcIC8YZPPseR1at0h7SCm4nXNhYIKZDQZWmtmXnHPbjmrzDPAMwKRJk9p2RJ6fT8aePZwHhIFngb3A7iFD2H3LrZQU7Wffvn188MEHrF69Gr/ff8Tm2dnZ5OXlHRHo8QGfl5dH//7921SSdB7nHIcPH04ayPv27aOoqIhwOHzEdicMGsTYnFxGj81l/JcmMmx0DqPG5HDBv93LKeVlnA0YsDDW3j8mp6P/aSLtok2zSpxzh81sA3AFsK2V5qn76U+bxiAzgdOAU7L7suNHD/Kl6bPo2zuTof16M7h/FgN6Z1JeXt7sP3Lj43Xr1lFcXIw7agLt0KFDEx6xNz4fM2YMvXppkk1H8Pv97N+/P2EgNz6uq6s7YpusrCzy8vIYm5PL5PMuYPTYXIaPHsuwUWMZOnIsI0bn0H/gwISvN8yMvB/8C+avb1oWzu7L3jt/mNZ/p0i6tJpUZjYCCMZCuy9wKdG/NttP4weQ118PgQD+sbnsvfOHlE+fBUB9Q5iihnqKDteTlWkM7pdN7vgz+NJZE+jdq/nEmGAwSFFRUcI/o3ft2sXGjRupqqo6YpvMzExycnKaDcPEB/zgwYM1JNOKSCTCwYMHkwbyvn37OHjwYLPtRo0aRV5eHl/84ul849JvMnpsDiPHRMN56Mgx9B88nGM98bDxfcT3bwLnmr2/RLwmlUPMMcDvYuPcGcBy59wf2r2Sa68l/PQz1AVCvLd0VdJmwbCjrKaBspoGoI5+vTM5oW8WJ2T3YmB2Fr17ZZCVlcW4ceMYN25c0v1UV1c3BcrRwfL222+zYsUKgsHgEdv0798/6Th7fn4+ubm5ZGdnt1OHdE01NTVJA7nxe0v9dvbZZ5Obl8/onBxGjs5hxJgcho4Yg8vMIhCKJD0t/HjPFi+fPovwfdFBkvc2vX98OxPpZKnMKvkbMLEDajkmvoYwvoYwB2IH0NlZGQzM7sWAPln075NJ/969yMhofpR8wgkncOaZZ3LmmWcm3G8kEqG0tDRpSL3//vsJjxxHjhzZ4gepo0aNSniHk66gpb9UGh8n+0slPz+fyZMnM2vWLMbm5jFqTA6jxuYyYsxYsvufQEPI0RCOEAiGm10JrzoIBHXRJZFUdbtBXX8wgj/YeEQOZtCvd2bsqxf9e/ciu3cGfXpltrifjIwMRo8ezejRozn33HMTv5bfT1FRUcKQ+/jjj3n99deTjtUm+yA1Pz+fgUnGao+wZEl0GuXevZCfH/2coIU57845ysvLE4Zx4/NEnw0MGzaMvLw8TjrpJL524YWMHZvLmJxcRuXkMnJ0DoOHjyRCBg2hCA3hMA2hI7eviUBNjabXibSnbhfcR3MO6gJh6gJhYjN4AcjMMPpmZdI3FuLZWZn0ycqgT68MemdmpDSWnZ2dzcknn8zJJ5+c5LX/MTsi0bDCxo0bE86OGDRoUIsfpOZs3EjWzTcfccKS74Yb2FdczN4JE5IOYdTX1x/xOn369IkO8eTl8fVLLmVMYyiPzWHkmByGjxpLrz59CYYjBMOu2ZmCIeBQnY6URTpatw/uZMIRR20gRG2CE+PMoHcswPv0yiArM4NemUbvzAx6xR5nZWSQmWH0yrCEQzHR/RhDhgxhyJAhnH322YnrCIcpKSlJOiSzefNmysvLj9wvMBY4EHs8GDhUXw8LF/6jjRmjRo9mbG4e40//El+75HJG5+QyYvRYRo7OZcSYHAYMGoIj+S+oegfoutEiXU6PDe6WOAeBYIRAMEJNCu0zDHplGhlmZGZEv5tBhlnsK/rLAKLLm0Vln8GMGj+YUePPIjoD3uEcRBxEnKOuro6S4iJK9u+nqGgfLLyNfcBzRD8tngXkA7lAzdJV0alyI0eT1coZp7rShog3KbjbQcQRG9tNUxRab4bnnMTwnJP4MjDxF4+QXbyfl2Orn4p994/N5b1zz0tPDSLSZXTN6Q3Sor13/pBw9pF3CNcJJSI9h464PUgnlIj0bDri9qjy6bMIDxhIeOAJvLfpfYW2SA+i4BYR8RgFt4iIxyi4RUQ8RsEtIuIxCm4REY9RcIuIeIyCW0TEYxTcIiIeo+AWEfEYBbeIiMcouEVEPKbV4DazPDNbb2Yfm9lHZnZ7RxQmIiKJpXJ1wBDwfefcu2Y2ENhqZuucc9vTXJuIiCTQ6hG3c67EOfdu7HEN8DGQk+7CREQksTaNcZvZOGAi8HaCdQvMbIuZbSkrK2uf6kREpJmUg9vMBgAvAXc456qPXu+ce8Y5N8k5N2nEiBHtWaOIiMRJKbjNLItoaC9xzr3cWnsREUmfVGaVGPAb4GPn3L+nvyQREWlJKkfc5wP/C/iGmb0f+5qa5rpERCSJVqcDOuf+AlgH1CIiIinQmZMiIh6j4BYR8RgFt4iIxyi4RUQ8RsEtIuIxCm4REY9RcIuIeIyCW0TEYxTcIiIeo+AWEfEYBbeIiMcouEVEPEbBLSLiMQpuERGPUXCLiHiMgltExGMU3CIiHqPgFhHxGAW3iIjHpHKX99+aWamZbeuIgkREpGWpHHE/D1yR5jpERCRFrQa3c24TUNEBtYiISAo0xi0i4jHtFtxmtsDMtpjZlrKysvbarYiIHKXdgts594xzbpJzbtKIESPaa7ciInIUDZWIiHhMKtMBC4G3gNPMbL+ZXZ/+skREJJlerTVwzhV0RCEiIpIaDZWIiHiMgltExGMU3CIiHqPgFhHxGAW3iIjHKLhFRDxGwS0i4jEKbhERj1Fwi4h4jIJbRMRjFNwiIh6j4BYR8RgFt4iIxyi4RUQ8RsEtIuIxCm4REY9RcIuIeIyCW0TEYxTcIiIek1Jwm9kVZvaJmX1uZovSXZSIiCSXyl3eM4EngSnAGUCBmZ2R7sJERCSxVI64zwU+d87tdM41AMuAGektS0REkumVQpscYF/c8/3A5PSUA2ZGVqala/fdkvqrbdRfbaP+So1Zx/VTKsGdqBrXrJHZAmABQH5+/jEV8z/OmQjApHFDj2n7nuaUk78AqL9Spf5qG/VX23ztK5M67LXMuWYZfGQDs68C9zvnLo89vwfAOfdgsm0mTZrktmzZ0p51ioh0a2a21TmXUvqnMsb9DjDezE4ys97AHGDV8RQoIiLHrtWhEudcyMxuBV4HMoHfOuc+SntlIiKSUCpj3DjnXgVeTXMtIiKSAp05KSLiMQpuERGPUXCLiHiMgltExGMU3CIiHtPqCTjHtFOzMmDPMW4+HDjUjuW0F9XVNqqrbVRX23THuk50zo1IpWFagvt4mNmWVM8e6kiqq21UV9uorrbp6XVpqERExGMU3CIiHtMVg/uZzi4gCdXVNqqrbVRX2/TourrcGLeIiLSsKx5xi4hICzoluM1stpl9ZGYRM0v6CWyymxTHLjH7tpl9ZmYvxC432x51DTWzdbH9rjOzIQnafN3M3o/78pvZzNi6581sV9y6CR1VV6xdOO61V8Ut78z+mmBmb8V+3n8zs2vi1rVrf7V2U2sz6xP7938e649xcevuiS3/xMwuP546jqGufzWz7bH++ZOZnRi3LuHPtIPqmm9mZXGv/7/j1s2L/dw/M7N5HVzXf8TV9KmZHY5bl5b+MrPfmlmpmW1Lst7M7LFYzX8zs3Pi1rV/XznnOvwLOB04DdgATErSJhPYAXwB6A18AJwRW7ccmBN7/Cvgpnaq6xFgUezxIuDhVtoPBSqAfrHnzwOz0tBfKdUF1CZZ3mn9BZwKjI89HguUAIPbu79aer/EtbkZ+FXs8RzghdjjM2Lt+wAnxfaT2YF1fT3uPXRTY10t/Uw7qK75wBMJth0K7Ix9HxJ7PKSj6jqq/feIXmo63f11IXAOsC3J+qnAWqJ3DPsK8HY6+6pTjridcx875z5ppVnCmxSbmQHfAFbE2v0OmNlOpc2I7S/V/c4C1jrnfO30+sm0ta4mnd1fzrlPnXOfxR4XA6VASicZtFEqN7WOr3cFcEmsf2YAy5xzAefcLuDz2P46pC7n3Pq499BmILedXvu46mp8Wo1XAAADrklEQVTB5cA651yFc64SWAdc0Ul1FQCF7fTaSTnnNhE9SEtmBvB7F7UZGGxmY0hTX3XlMe5ENynOAYYBh51zoaOWt4dRzrkSgNj3ka20n0PzN81PY38q/YeZ9engurLNbIuZbW4cvqEL9ZeZnUv0KGpH3OL26q9k75eEbWL9UUW0f1LZNp11xbue6JFbo0Q/046s61uxn88KM8tr47bprIvYkNJJwJtxi9PVX61JVnda+iqlGykcCzN7AxidYNVi59wrqewiwTLXwvLjrivVfcT2Mwb4MtE7AzW6BzhANJyeAe4GftyBdeU754rN7AvAm2b2IVCdoF1n9df/A+Y55yKxxcfcX4leIsGyo/+daXlPtSLlfZvZPwOTgIviFjf7mTrndiTaPg11rQYKnXMBM7uR6F8r30hx23TW1WgOsMI5F45blq7+ak2HvrfSFtzOuUuPcxf7gby457lAMdHrAAw2s16xo6bG5cddl5kdNLMxzrmSWNCUtrCrfwJWOueCcfsuiT0MmNlzwJ0dWVdsKALn3E4z2wBMBF6ik/vLzE4A1gA/jP0Z2bjvY+6vBJK9XxK12W9mvYBBRP/8TWXbdNaFmV1K9JfhRc65QOPyJD/T9giiVutyzpXHPf018HDcthcfte2GdqgppbrizAFuiV+Qxv5qTbK609JXXXmoJOFNil10xH890fFlgHlAKkfwqVgV218q+202thYLr8Zx5ZlAwk+g01GXmQ1pHGows+HA+cD2zu6v2M9uJdHxvxePWtee/ZXKTa3j650FvBnrn1XAHIvOOjkJGA/89ThqaVNdZjYReBqY7pwrjVue8GfagXWNiXs6Hfg49vh14LJYfUOAyzjyL8+01hWr7TSiH/a9Fbcsnf3VmlXAt2OzS74CVMUOTNLTV+n4BLa1L+Bqor+JAsBB4PXY8rHAq3HtpgKfEv2NuThu+ReI/sf6HHgR6NNOdQ0D/gR8Fvs+NLZ8EvBsXLtxQBGQcdT2bwIfEg2g/wQGdFRdwHmx1/4g9v36rtBfwD8DQeD9uK8J6eivRO8XokMv02OPs2P//s9j/fGFuG0Xx7b7BJjSzu/31up6I/b/oLF/VrX2M+2guh4EPoq9/nrgi3HbfifWj58D13VkXbHn9wMPHbVd2vqL6EFaSey9vJ/oZxE3AjfG1hvwZKzmD4mbLZeOvtKZkyIiHtOVh0pERCQBBbeIiMcouEVEPEbBLSLiMQpuERGPUXCLiHiMgltExGMU3CIiHvP/AYmleiyW/w0WAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.fill_between(x, f(x), alpha=0.25)\n",
    "ax.plot(X, f(X), 'ro')\n",
    "\n",
    "for n in range(len(X) - 1):\n",
    "    f_mid = f(X[n:n+2].mean())\n",
    "    ax.plot([X[n], X[n]], [0, f(X[n])], 'k')\n",
    "    ax.plot([X[n+1], X[n+1]], [0, f(X[n+1])], 'k')\n",
    "    ax.plot(X[n:n+2], [f(X[n]), f(X[n+1])], 'k')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py3.6",
   "language": "python",
   "name": "py3.6"
  },
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}