riemann-sums.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Riemann Sums"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition\n",
    "\n",
    "A [Riemann sum](https://en.wikipedia.org/wiki/Riemann_sum) of a function $f(x)$ over a partition\n",
    "\n",
    "$$\n",
    "x_0 = a < x_1 < \\cdots < x_{N-1} < x_N = b\n",
    "$$\n",
    "\n",
    "is a sum of the form\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^N f(x_i^ * ) (x_i - x_{i-1}) \\ \\ , \\ x_i^* \\in [x_{i-1},x_i]\n",
    "$$\n",
    "\n",
    "where each value $x_i^* \\in [x_{i-1},x_i]$ in each subinterval is arbitrary.\n",
    "\n",
    "Riemann sums are important because they provide an easy way to approximate a definite integral\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x) \\, dx \\approx \\sum_{i=1}^N f(x_i^ * ) (x_i - x_{i-1}) \\ \\ , \\ x_i^* \\in [x_{i-1},x_i]\n",
    "$$\n",
    "\n",
    "Notice that the product $f(x_i^ * ) (x_i - x_{i-1})$ for each $i$ is the area of a rectangle of height $f(x_i^ * )$ and width $x_i - x_{i-1}$. We can think of a Riemann sum as the area of $N$ rectangles with heights determined by the graph of $y=f(x)$.\n",
    "\n",
    "The value $x_i^*$ chosen in each subinterval is arbitrary however there are certain obvious choices:\n",
    "\n",
    "* A *left* Riemann sum is when each $x_i^* = x_{i-1}$ is the left endpoint of the subinterval $[x_{i-1},x_i]$\n",
    "* A *right* Riemann sum is when each $x_i^* = x_i$ is the right endpoint of the subinterval $[x_{i-1},x_i]$\n",
    "* A *midpoint* Riemann sum is when each $x_i^* = (x_{i-1} + x_i)/2$ is the midpoint of the subinterval $[x_{i-1},x_i]$\n",
    "\n",
    "Let's visualize rectangles in the left, right and midpoint Riemann sums for the function\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{1}{1 + x^2}\n",
    "$$\n",
    "\n",
    "over the interval $[0,5]$ with a partition of size $N=10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAE/CAYAAADVKysfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzs3Xl8VNX9//HXJ5MdEFFcWMQFt4LKIpvsSNgK1AU3tFKXui+1rbW27vKrtf221q8Wv2hbi1tLFcUFQRAFAQUUd7QuoCiILW6A2TOT8/vj3JAQE5LAJPdm8n4+HnmEzL259zND5jPnc86555pzDhEREREREYmmtLADEBERERERkdqpaBMREREREYkwFW0iIiIiIiIRpqJNREREREQkwlS0iYiIiIiIRJiKNhERERERkQhT0dbIzOwwM3vdzL41s8uTdMx5ZvajZBxLRJoXM5tuZtftYLszs4OTcJ4mzTNm9o6ZDW+q84lI06ord1Xbd4aZ/b9GjKWLmeWbWayxziGSbCra6snM1plZ3k786lXAYudcG+fcHfVJREGjqyBIKJ+Z2W1VE4tzbpxz7r6diCUyzKyzmT1qZl+a2RYze9vMzgohjgOC1/vpao8/aGY3JvlcR5jZ/OA5f+cGiWa2h5nNDv7vPzGz05N5fom2IMeUmln7ao+/EfyNHgDgnLvQOTe1seNpSJ4xs8Vm9uMdbK94n+UHX+vM7Opq5+vunFu8i2GHysy6m9kCM/vGzDab2atm9v0Q4hgevN7Tqj2+LNl51sxGmNmiII+vq2H7AcH2QjN7byc/R6UZCN7XRcF7/D9Be6d1xfZk5q66OqfM7CwzSwSxbDWzN81sQpVYPnXOtXbOJZIRT1jM7LjgM2Jr0LZ4ruKzoonjuDH4Pzm5ymPpVT+7kniuS81slZmVmNmMGraPDHJNYZB79k/m+cOkoq3x7Q+8sxO/18M51xoYBpwKnJPUqML3ALAe//rsCUwB/htiPAPMbFAjn6MMeBg4t5bt04BSYB/gDOD/zKx7I8ck0fIxMLniBzM7EsgJL5yk2z3IaycB15nZqLADSrKngGfx7+G9gcuBrSHFUgBMaYIGXAFwL/CLWrb/E3gdn+evAWaZ2V6NHJOEZ2LwHu8J9AJ+FWIsy4NYdgfuAmaa2e4hxpNUQdF6P/BzoC1wIP55locU0tfAzU0werkR+H/4vLOdoNPzMeA6YA9gFfCvRo6nyahoSwIzmxD0dGw2s5fM7Kjg8eeBEcCfg96e8/GN8auCn5+q69jOuTXAi/gEWHG+7Xq1zewcM/t30Ls7v2qvQtDLcbGZfWh+iuZUM+tqZsuDnpmHzSwz2Ledmc0xsy+CY80xs87VzjvVzF4MjrWgYlSgSk/6j8zs06DH55odPLW+wAznXIFzLu6ce905Ny841nAz21DtNd420hn06DxifjTsW/OjdIea2a/MbJOZrTez0XW9ttX8Hp8EGo1z7n3n3N+ooYg3s1bAJOA651y+c24Z8CRwZmPGJJHzAL4Do8KP8B/K21i10Xoz+4WZfW5mG83snBr2nW5mzwbvlReq5YeBZvZKMEryipkNrLJtW54Jeq2XmdkfgtzwsZmNC7b9BhhCZZ77c11P0jm3Cv8+qJrXqr7H08zsajNba2ZfBXlqj2BbRa45O3ivf2NmF5pZXzN7K8jDf65y3K5m9nxwnC/N7KGqDbfgvFcGv7vFzP5lZtnBtuFmtsHMfh7kls/N7OyanlOQCw8E/uKcKw2+Xgzey9tew2q/s220IPi/usv8tNT8IM/ua2a3B8/xPTPrVddrW8VmYAZwQwN+p8Gccy875x4APqq+zcwOBXoDNzjnipxzjwJv43OdpDDn3H+A+Wz/Hq+eu66qkrt+bN8dPWtnZk8HuWulmXUNfm9JsP3N4L1yah2xlONzayvgkOAYFXkkPfi5rZn9LYjnMzP7fxYUH8F790Uz+1OQXz4KcudZQQ7aZFWmkpvZePOXxWwNtt9YZdsO20rm2zcPm9n9wfN+x8z61PLUegIfO+eec963zrlHnXOf1vJ6b9e2CnLfL4LcVxA8/32CHPStmS00s3Y7em2reQbf8fzDBvxOgznnHnPOPQ58VcPmE4F3nHOPOOeKgRuBHmZ2eGPG1FRUtO0iM+uNr/YvwPck3g08aWZZzrljgaXApcEw/D3AQ8Dvg58n1uP4h+MbRGtq2X488Gv8H+pewfn+WW23scDRwAD8dM178MXjfsARVPbspwF/x49+dQGKgOoNsNOBs/G9yJnAldW2DwYOA0YC15vZ92p5aiuAaWZ2mpl1qWWfHZmIT8Lt8L2484P4OwE34/8fGmIacKjVY+qOmQ0OEndtX4MbeG6AQ4GEc+6DKo+9CWikrWVZAexmZt8LGgynAg/WtrOZjcW/B0fhGyM1/f2eAUwF2gNv4HMQ5ougp4E78LnrNuBpM9uzltP1B94PjvN74G9mZs65a9g+z11a15M0swH43FNjXsOPUB2Pn2nQEfgG/x6tHs8h+NfodvwoTh7+PXOKmQ2rOB3w2+A438PnvRurHesUfJ48EDgKOKvKtn3xvdid8KPk02ppyHwVPJ8Hzex4M9unlue2I6cA1+Jf4xJgOfBa8PMs/P9RQ/wGmGRmh9W1o5mdXkde25k83R34yDn3bZXHlNdaAPMdvuOove0yFvgZ/j17MP69Xt1k4Cb85/wa/N8zzrmhwfYeQc7Z4UhKkEvPxs92+aSW3e4D4kEsvYDRQNUp3/2Bt/C58h/ATHzn88H4IuXPVjkVtADf+bY7MB64KGirVbWjttIPguPvju+8ra0j7DXg8KCYHFHl/A0xCf/5cSi+XTUP36Zsj29TNWQtBocf4brBzDLq2jnopKot37zV4GfidcfnGB+QcwXAWlIk56ho23XnAXc751Y65xLBNSAl+AJpV7xmZgXAv4HF+CHvmlwA/NY592/nXBy4Behp28/h/Z1zbqtz7h1gNbDAOfeRc24L/g3aC8A591XQS1MYfMj+hu8m0r875z5wzhXhp/r1rLb9pqBH9U38G6dHLXGfjG/oXQd8bH6ksm8dr0lVS51z84Pn/Ai+YL3VOVeGT3YHWMOmQRTjn2+do23OuWXOud138LWsrmPUoDWwpdpjW4A2O3Esad4qRttGAe8Bn+1g31Pw78nVwYfTjTXs87RzbolzrgRf2BxjZvvhGxMfOuceCEa7/xmcr7bOpE+cc38JrgG5D+iAnwbYEF+aWRG+GLkLeLyW/S4ArnHObQjivhE4qaJXPDDVOVfsnFuAbyT90zm3yTn3GT63VOS1Nc65Z51zJc65L/CFT/W8dodzbqNz7mv8FMeqea0MuNk5V+acmwvk4xtb23HOOfzMinXAH4HPzWyJmR1Sz9cGYLZz7tWgh3g2UOycuz94zf9V8ZzqKxjtmI7vyKpr33/Ukdc+bci5A8prLc/jZvYt/vKHTdQ+0luRu95xzhXii7PqHgtGcuP4zqbq7Y26DDCzzfjP9z8AP3TObaq+U9DBMg64Ipj9swn4E3Bald0+ds79vcp7cT98XigJclApvoDDObfYOfe2c67cOfcWviO9es7ZUVtpmXNubnCuB6ilHeWc+wgYju9QehifX2c0sHi70zn33yp5c6XzM59K8DmooTnnSeALti94a9v34h3km6Mact4qUjrnqGjbdfsDP6/aQ4B/M3fcxeP2xv/xnYrv4Wm1g/P/b5Vzf43vWe5UZZ+q14oV1fBzawAzyzWzu80vgrEVWALsbtvPT/5PlX8XVvxuA7YD4Jz7xjl3tXOuO77h9wY+2Vstz7O66s/hS1d5QXFR8L2hvU5/AfYxszpHQBtBPrBbtcd2A76tYV9JbQ/gR7TPotrUyBp0xDeOKtTUi7xtu3MuH58jOgZf1ff/hO1zR1Xb3ttBIwsa/h5rH/zOlfjGRm29sfsDs6vktX8DCbYvEuub1/Y2s5nmpzxtxY9cbrfYCzvOW18Fjcbatm8TFJmXOue6Bs+hgLr/D6uq13NqoN8BY8ystg60xqS81vIc75xrg39/H85332sVqueu9TXsU6/2xA6scM7tjh+pexI/a6km++Nz0edVcs7d+BlFFaq/F3HO1ZZz+ptfAOMLM9sCXEjDck71bdnVOqy2cc6tcM6d4pzbK3h+Q/Gdc/XVGDnn2iCG7J343V2V0jlHRduuWw/8ploPQW7Qa12T76waWBvnPYzvlb5+B+e/oNr5c5xzLzXsaQD+YtbDgP7Oud3wb37wRWCjcc59ie8F64i/cLQAyK3YHhSNjX7hejBKdxN+Klmtz9nMhljlKng1fdX2wbAjHwDp1Xrle7Bzi9hIM+ac+wS/IMn38RdU78jn+E6iCjVNYdu2PeiB3QN/IfdGfGOlqi7seGSvNg3Jawnn3B/xvd8X17LbemBctbyWHfQGN9Rvg/iOCvLaD2nknAbgnFuPn9J5RPBQ9by2b2PHEMTxFX766A5X7TOzM+rIazszPfId4CAzq9rLrbzWAjjnXsBfU/mHWnb5HOhc5ef9atkvGbHk43PNmVbzdaHr8TOk2lfJN7sFnco74x/4InE/51xb/Gh3U+ScV/CfGTXmHPxU70bnnHsWP521tvwObLsFRG35ZmdzxDtUGZk0v15AV1Ik56hoa5gMM8uu8pWOH525MOhZMTNrZf4i1NqGYv8LHNTA894KnF/Lh/x04FcWrDJo/mLak2vYrz7a4HtWNpu/3qXRLmA3s9+ZXwI/PXitLgLWBA2MD/A9S+ODedHXAlm7cK4bzWxxPXd/IDjX2Np2cM4tdX4efW1fS2uJw8wvcFCx8Eu2mWUFxyzAJ9ubg7+hQcBxQTzS8pwLHBv8XezIw8BZZtbNzHKp+T37ffPXYWbiG+4rg4JiLv46ztOD9+GpQDdgzk7Eu7N57argPVHddOA3FdO8zWwvMztuJ+ICn9fy8XmtE7WvcrhLzC/kdJOZHWx+IZX2+FV/VwS7vAl0N7OewXO+cRfPN8NqWO66FrcBA/HX9NXIOfdQHXmtxumRwXPNxo9UWJDXMoNjfoCfRXFD8PgJ+GsGH63/M5Vm7HZglJnVNK3xYeBs89fv5lJ7x3RtGpRzgrbFX2s6j3Puc2AB8Ecz2y34m+5qldfFNlQb4GvnXLGZ9cPPnEi6IK+fZ2Z7Bz8fjr8eriLnvIHP/3sE7ccrdvF866z+twy5Br+GQq2cvwVEbfmm1oI5+LzKBmJArEp7HPyUziPMbFKwz/XAW8659+oZd6SpaGuYufiipuLrRudXQTsPf6HoN/jehbN2cIy/Ad2CIfjarufYjnPubeAFamhsOOdm46e/zDQ/9Wc1fm72zrgdv7z4l/g3/TM7eZz6yMW/uTbjVx3bH59scP5au4vxCfYzfG/RhpoPUy/74VfgrFMwxfIG/GhEsu2P/7up6PEpwi/sUOFi/Ou/CT8H/iLnr0OUFsY5tzbILXXtNw//vn0en3uer2G3f+D/pr/GL0h0RvC7XwET8CPsX+E/YCcEI98N9b/4a86+MbM76vk7T+Nz5nm1HO9JYIH562NW4KeJ74yb8NPNtwTnrGv0cmeVAgcAC/HL/K/G996fBdsKmJuD7R8CO3Pta1UNyWtb8YvHNEZeG4rPZXOpXMBqQZXtpwF98P/XtwInOX9toaS44P/5fvy169W3zcMvgrQIn7uWB5tK6nn4G4H7grbUKfX8ndvxRUxN10tNwXeovov/W52Fv253Z1yM74D9Fl80PLyTx6nLZny76W0zy8e32Wbj3+vgO33fxF9nu4BdWPo+6IjZk8qCcIeccy8CL+/s+epwLT7PXI2fOVEUPFbxNzcJv0bBN/jPjdNqPkzzY87Ve1aLSLNkZm8AI4NGqkiLEYzEbHDOXRt2LJI8QQPqTfyUz7Kw4xHZVeZXT1wNZFW7hlQiwPyq2Jc45ybXubM0mhovbBRJJc65hq44JSISWc65UnYw1VGkOQimyz6NX2jtd8BTKtiiyflVsXd1doDsIk2PFBEREZGmdgF+efi1+NVhLwo3HJFo0/RIERERERGRCNNIm4iIiIiISISpaBMREREREYmw0BYiad++vTvggAPCOr2INIJXX331S+dco98IvbEpP4mknlTIT8pNIqmnvrkptKLtgAMOYNWqOm9DJCLNiJl9EnYMyaD8JJJ6UiE/KTeJpJ765iZNjxQREREREYkwFW0iIiIiIiIRpqJNREREREQkwlS0iYiIiIiIRJiKNhERERERkQhT0SYiIiIiIhJhKtpEREREREQirM6izczuNbNNZra6lu1mZneY2Roze8vMeic/TBGR71J+EpEoUm4SkWSrz0jbDGDsDraPAw4Jvs4H/m/Xw/ISCZgzB6ZO9d8TiWQdWURSxAxCyk8iIjswA+UmEUmi9Lp2cM4tMbMDdrDLccD9zjkHrDCz3c2sg3Pu810JLJGAMWNg5UooKIBWraB/f5g/H2KxXTmyiKSKsPKTiMiOKDeJSLIl45q2TsD6Kj9vCB7bJfPm+YItPx+c899fegnmzt3VI4tIC9Io+anCoEFw5pnwj3/Ali3JOqqItACNmpvuvBOGDoVbb4V33knWUUUkTHWOtNWD1fCYq3FHs/Px0wDo0qXLDg/6+ut+hK2qoiK4/IpyLKuMzKydivU72raO0f/oZLwMIhJBjZKfAIqL4aCD4Jln4MEHoWNHeOop6K0rU0Skbo2WmwB22813dv/qV/7rhhv8l9V0VhFpFpJRrWwA9qvyc2dgY007OufuAe4B6NOnT43JqUKvXn5KZH5+5WMZGY51H6Vx/bUxbv9rIe322OEh6uXztdkk52UQkQhqlPwEkJ0NDzzgp3IvWwZTpsCQITBzJkycmIzQRSSFNVpuAvjRj/zXxo3w61/DTTfBmjXwt79BVpI6vUWkaSVjeuSTwJRgJaQBwJZkzMkeN85fw9a6NZg5cnIdvfvF+f20At59K8bl5+QSj+968CKS0holP1UVi8GwYX46d7ducMIJ/t8iIjvQ6LkJ/AyAv/8dfvMbeOghuOSSZJ9BRJpKnUNMZvZPYDjQ3sw2ADcAGQDOuenAXOD7wBqgEDg7GYHFYn7RkXnzYNbsBH2GlDB4RJxYDMrLi7j6slz+cmcWF/20JBmnE5FmKKz8VJN994Vnn4WePeGMM/wU7zZtGutsIhJlUcpNZn60LT8ffvtb3yk+aVJjnU1EGkt9Vo+cXMd2BzRK300sBhMmQGabBB26Vg6rff/4MpYtLuXu27M4Zkicnn10LwCRlijM/FST3Xf317cNGwaXX+57uEWk5YlabgI/RXLhQjjvPD+TqXPnpjy7iOyqZEyPDMWvpxbRsXM51/40h7KysKMREfEGD4ZrroEZM/xMARGRKMjI8CvdlpbCxReHHY2INFSzLdpat4Grby7m03UxHv9XZtjhiIhsc911fmXJX/0KysvDjkZExDv4YLj2Wr/S7Ysvhh2NiDREsy3aAIYcG6dX3zh3355FcVHY0YiIeBkZcPPN8Oab8MgjYUcjIlLpssv8Nbi//rW/D66INA/Numgzg8uuKmbTf9OYeZ9G20QkOk47DY44wo+6aQq3iERFq1Z+tG3JEliwIOxoRKS+mnXRBtBnQIJBw8v427QsCgvq3l9EpCnEYn6Z7Q8/9PdzExGJivPOgwMO8MWbRttEmodmX7QBXHhFCVs2p/HUoxptE5HomDgRevSAP/1JDSMRiY7MTH/N7apVurZNpLlIiaLtqN4Juh0ZZ+Z9mWoYiUhkmPnrR1avhhdeCDsaEZFKZ5zhb1Ny551hRyIi9ZESRZsZTD67lLUfxFj5YizscEREtjn9dGjXzl/0P3UqzJkDCd1aUkRC1qoVnHsuPPoofPZZ2NGISF1SomgDGDuxjHZ7lDNzRlbYoYiIbJOZCW3awPLlcMMNMHkyjBmjwk1Ewnfxxf62JHffHXYkIlKXlCnasrLhxMmlLH42nY0bLOxwREQAf4PtL7/0/3YO8vNh5UrdeFtEwnfQQTB+vC/aSkrCjkZEdiQ97ACS6ZQzS7n3riyeeCSTi36q7CMijW/lq3G25Nc+bPbI7BhFRTGgsjOpoMAxa3aCzDa1/17b1jH6H51SKVpEIuiSS/y07aefhhNPDDsaEalNSrUIOnRy9B2YYM6jGVx4RQmmATcRaWRb8hN06Fpc6/a+Q9KZNTOXwsLKx7JzoM+QEjp0jdf6e5+vzSbFUrSIRFBenr/Z9v33q2gTibKUmR5ZYeKkUtZ/EuPNV7UgiYiEb/CIOEf2ipOT6wBHLOY4qlecwSNqL9hERJpKerpfSXLuXPjqq7CjEZHapFzRljeujOxsx1OPZoQdiogIsRhMf6iQ308rpPtRCdJicNs9hcTUryQiEXHmmVBWBv/6V9iRiEhtUq5oa9Uajh1bxvynMijVZW0iEgGxGAzLi/OLG4opKzUWLVCnkohER48ecOSR8MADYUciIrVJuaINYMKJZWzdksaS53U9iIhER6++CTp1KWfOYyraRCRazjwTVqyADz4IOxIRqUlKFm0DhsTZc69y5j2RGXYoIiLbmMH440tZuSydr77QSkkiEh2nn+5z1MyZYUciIjVJyaItPR2OHVPG0ufTKSoKOxoRkUqjJ5ThnPH8fI22iUh0dOoEAwfCo4+GHYmI1CQlizaAUePLKC4yXlqsKZIiEh2HHF7O/gcmWDhXuUlEouWkk+Ctt2DNmrAjEZHqUrZo6zMgwe7tylk4T73ZIhIdZpD3/TJefimdzd9oiqSIREfFfdo02iYSPSlbtKWnw4jRcV5YqFUkRSRa8sbFSSSMxQs02iYi0dGlC/Ttq6JNJIpStmgD35ud/62xYpkaRiISHd2OStCxs2YCiEj0TJoEr7wCn3wSdiQiUlVKF239B8Vps5tj4Vw1jEQkOswgb1wZLy1J59utYUcjIlJp0iT//bHHwo1DRLaX0kVbZhYMHVnGCwvTSSTCjkZEpNLIcWXEy4xli9SpJCLRcfDBcNRR8PjjYUciIlWldNEGMCwvzjdfp/HWa7GwQxER2eao3gna7VHOCws1fVtEomXiRHjxRfj667AjEZEKKV+0DRpeRnq6Y8lzahiJSHTEYjD42DjLFqUTj4cdjYhIpYkTIZGAZ54JOxIRqZDyRVub3aB3/wSLn9UUJBGJlmF5ZWzdksYbqzQTQESio29f2GcfeOqpsCMRkQopX7SBbxit/SDGhk91TyQRiY6BQ+OkZzheWKhOJRGJjrQ0GD8e5s2DsrKwoxERaCFF2/A8P/dIDSMRiZLWbaDvgLiuaxORyJk4EbZsgWXLwo5ERKCFFG37HVDOQYck1DASkcgZNirOurUxPvm4RaRjEWkm8vIgK0tTJEWiosW0EoaOjLNqRToF+WFHIiJSaehIP/doiTqVRCRCWreGESNgzpywIxERaElF27H+nkgrl6lhJCLR0bmL46BDEixbpNwkItEyfjx8+CGsXRt2JCLSYoq2Hn0StGrtWLZYDSMRiZZBw+KsWplOYWHYkYiIVBo71n+fNy/cOESkBRVtGRkwYHCcZYsycC7saEREKg0eEaes1Fi1XJ1KIhIdBx/sv1S0iYSvxRRtAINGlPGfjWl89GGLetoiEnFH94+TneM0RVJEImfcOFi0CIqLw45EpGVrUdXL4OF+6X81jEQkSjKzoN+gOC9q+raIRMzYsVBUBEuWhB2JSMvWooq2fTs6uh6aYNli3a9NRKJl8PA46z/R0v8iEi3Dh/ul/zVFUiRcLa51MHhEnNdejlFYEHYkIiKVBg/3S/9rtE1EoiQ31xduzzwTdiQiLVuLK9oGDS+jrNR4+SU1jEQkOjrv79j/IC39LyLRM3YsvPcerFsXdiQiLVeLK9p6902Qne1YsVQNIxGJlkHD4qxakU5pSdiRiIhUGjPGf3/22XDjEGnJWlzRlpkFfY6J89ISFW0iEi3HDI1TXGS8vioWdigiItscfjh07gwLFoQdiUjL1eKKNvANo3VrY2zcYGGHIiKyTd9j4qRnOJarU0lEIsQMRo+GhQshkQg7GpGWqUUWbQOH+qX/l2uKpIhESG4r6Hl0gpeWaIVbEYmW0aNh82ZYtSrsSERapnoVbWY21szeN7M1ZnZ1DdvbmtlTZvammb1jZmcnP9TkOeiQcvbet1y92SLNXKrlJvAzAd5bHWPzN2FHIiK7ItXy08iRfsRNUyRFwlFn0WZmMWAaMA7oBkw2s27VdrsEeNc51wMYDvzRzDKTHGvSmMHAIXGWPp/OAzNizJmj4X6R5iYVcxNUzgR47dUWORFCJCWkYn5q3x5694aZM2HqVNR2Emli9WkV9APWOOc+cs6VAjOB46rt44A2ZmZAa+BrIJ7USJMokYDVb8UoKkzjwRkxJk/2KyMp+Yg0KymXmwAOPyLB7u3Kee0VFW0izVjK5adEAr78Et59F264AbWdRJpYfVoFnYD1VX7eEDxW1Z+B7wEbgbeBnzjnypMSYSNYtiidz9ZXPHUjPx9WroR580INS0QaJuVyE0AsBgMGx3l1VRrOhR2NiOyklMtP8+bBpk3+386htpNIE6vPRV01LbFYvSkxBngDOBboCjxrZkudc1u3O5DZ+cD5AF26dGl4tEny3uoYxUXbP1ZQ4Jg1O0Fmm+R0GbVtHaP/0bpmTqQRJS03QXTyE8CAoXGeeSqTd9+F7t1DDUVEdk7obaeVr8bZkp+8YbBHZscoLo5R9akVFMAbb8CECUk7jYjUoj5VxQZgvyo/d8b3ClV1NnCrc84Ba8zsY+Bw4OWqOznn7gHuAejTp09ofciHH5EgJwcKCysfy86BPkNK6NA1OTMTPl+bTf1eXhHZSUnLTRCd/AR+pA38jWxVtIk0S6G3nbbkJ+jQtXgnQq9Z3yHpzJqZu13bqVUr6NkzaacQkR2oz/TIV4BDzOzA4ALZ04Anq+3zKTASwMz2AQ4DPkpmoMk0eEScI3vFycryuS8ry3FUrziDR0R2KrmIfFfK5aYKHTs7OnUuZ+HCsCMRkZ2Ucvmpou2UkeHbTrm50L8/jBsXcmAiLUSdRZtzLg5cCswH/g087Jx7x8wuNLMLg92mAgPN7G3gOeCXzrkvGyvoXRWLwfSHCvntHYWkpTmOGRpn+kOFxGJhRyaRzJRhAAAgAElEQVQi9ZWKuamq3n0cixdDaWnYkYhIQ6VifqpoO115vb++5Mc/hvnzUdtJpInUa/6ec24uMLfaY9Or/HsjMDq5oTWuWAzyvh/noEMTrP8kTUlHpBlKxdxUodfR5Tz1eIwVK2Do0LCjEZGGSsX8FIvBaT8q4+4/OTZtMrWdRJpQi19T+uDDEqz9IMYX/63pmmERkXD06FlOWhqaIikikWLmO5Weew7KI7vWpUjqafFF2yGH+ZWVVizToiEiEh2t20C/fn4xEhGRKOnVp5wvvoC33go7EpGWo8UXbft2KqfdHuUsX6qiTUSiJS8PXn4ZNm8OOxIRkUq9evshNnUqiTSdFl+0paVB/0FxVi5N141sRSRSRo3y048WLw47EhGRSu33gm7dNH1bpCm1+KINoP+QOF9sSuOjD/VyiEh0DBjg74OkhpGIRE1eHixdCsXJuxWciOyAqhTgmCH+/my6rk1EoiQzE4YNU9EmItEzahQUFcHy5WFHItIyqGjD38i2ywEJVui6NhGJmLw8eP99WL8+7EhERCoNG+ZvAaBOJZGmoaIt0H9wnFUr0ikrCzsSEZFKeXn+uxpGIhIlbdr4KdxajESkaahoCwwYEqcg31j9hu4UKSLRccQRsPfeKtpEJHry8mDVKvjmm7AjEUl9KtoC/Y5JYOZ0XZuIRIqZbxgtXIhWuBWRSMnL83lp0aKwIxFJfSraAm3bObofpevaRCR6Ro2CTZtg9eqwIxERqdS/P7RurZkAIk1BRVsV/QfHefv1GAX5YUciIlJp5Ej/XdeOiEiUZGTA8OHKTSJNQUVbFQOGxInHjVUrNNomItGx335w2GHqzRaR6MnLgzVrYN26sCMRSW0q2qroeXSC7GzHck2RFJGIGTUKXngBSkrCjkREpNKoUf67OpVEGpeKtiqysqF3vzgrtRiJiERMXh4UFsKKFWFHIiJS6Xvfgw4dVLSJNDYVbdX0HxJn7Qcx/vu5hR2KiMg2w4f7G9nq2hERiZKKFW6few7Ky8OORiR1qWir5pghcQBWvqjRNhGJjrZtoV8/9WaLSPSMGgVffglvvhl2JCKpS0VbNYd+r5x2e5Zr6X8RiZy8PHjlFd3IVkSiRSvcijQ+FW3VpKVB/0FxVixL141sRSRSRo3y048WLw47EhGRSh07Qvfumgkg0phUtNVgwJA4X25KY837enlEJDr694dWrdSbLSLRk5cHS5dCcXHYkYikJlUlNai4rk1TJEUkSjIzdSNbEYmmUaN8wbZsWdiRiKQmFW016NDJcUDXBC8tUdEmItEyerS/ke3HH4cdiYhIpWHDICMDFiwIOxKR1KSirRYDh8Z5dUU6pbqRrYhESMWNbDXaJiJR0ro1DByo3CTSWFS01WLAkDjFxcbrq2JhhyIiss3hh0PnzmoYiUj0jB4Nb7wB//1v2JGIpB4VbbXoe0yc9HSn69pEJFLM/Gjbc89BIhF2NCIilSpmAjz3XLhxiKQiFW21aNUaehyd4KUlGWGHIiKyndGj/b3aXn017EhERCr17g177KHr2kQag4q2HThmSJz3Vqfx9VcWdigiIttU3MhWDSMRiZJYzOenZ59F97oVSTIVbTtwzNA4zpmmSIpIpOy1l+/RVtEmIlEzejRs3AjvvBN2JCKpRUXbDnQ7KsFubctZrqX/RSRiRo+G5cth69awIxERqVRxXZs6lUSSS0XbDsRifhXJ5UvSNcwvIpEyejTE47BoUdiRiIhU2n9/v8rt/PlhRyKSWlS01WHgsDib/pvGh+/ppRKR6Bg0CFq1UsNIRKJnzBhYsgSKisKORCR1qBKpw8BhcQBeekFTJEUkOjIzYcQIFW0iEj1jxkBxsS/cRCQ5VLTVYd8Ojq6HJnhJ17WJSMSMHQsffQRr1oQdiYhIpWHDICtLnUoiyaSirR4GDYvz2svpGuYXkUjJy/Pfr7gC5szRzbZFJBpyc2HIEBVtIsmkoq0eBg6LU1pirFqh0TYRiYZEAi65BMzg6adh8mQ/JUmFm4hEwZgx8O67sH592JGIpAYVbfXQu1+crCyn69pEJDLmzYOVKytvYJuf73+eNy/cuEREwE/fBi39L5IsqkLqITsHjh4Q58XFerlEpGm890ECKKl1+yOzYxQUxADb9lhBgWPW7ASZbXY83Na2dYz+RyufiUjj6d4dOnWCZ56Bc88NOxqR5k+f2vU0eESc39+Yw4ZPjc5ddNM2EWlchSXldOhaVuv2vkPSmTUzl8LCyseyc6DPkBI6dI3v8Nifr81G6V9EGpOZH22bNcvfUzJdKUdkl2h6ZD0NHuEbQcsWZYQciYiIz0lH9oqTk+sAh5njqF7xbblKRCRs48bBli2wfHnYkYg0fyra6mn/A8vZb/8EyzRFUkQiIBaD6Q8V8vtphQwaFsc548b/KSIWCzsyEREvL8+PsOlaW5Fdp6Ktnsx8z/bLy9IpKQ47GhERX7gNy4tz5fU+KS3X/SRFJELatoVBg2Du3LAjEWn+VLQ1wOARcYqLjVdXqmEkItFx0CHldOxcrunbIhI548bBm2/Cxo1hRyLSvKloa4A+x8TJzHIsW6SiTUSiwwwGDS9jxbJ0ykrDjkZEpNK4cf77M8+EG4dIc1evos3MxprZ+2a2xsyurmWf4Wb2hpm9Y2YvJDfMaMjJgT4D4ixV0SYSCcpNlQaPiFNYYLz2ii5qE4kC5SfvyCP90v+aIimya+os2swsBkwDxgHdgMlm1q3aPrsDdwE/cM51B05uhFgjYejIOJ98FOPTjzVIKRIm5abt9R8UJyPTsfR5TZEUCZvyUyUzP9r27LNQqpkAIjutPpVHP2CNc+4j51wpMBM4rto+pwOPOec+BXDObUpumNEx9Fh/36Qlz2u0TSRkyk1V5LaCvsfEWfKccpNIBCg/VTFhAmzdCsuWhR2JSPNVn6KtE7C+ys8bgseqOhRoZ2aLzexVM5tS04HM7HwzW2Vmq7744oudizhknfd3HHRIQg0jkfAlLTdBauSnYXlx1q2N8YlmAoiETW2nKkaOhKwsmDMn7EhEmq/6fLJbDY+5aj+nA0cD44ExwHVmduh3fsm5e5xzfZxzffbaa68GBxsVQ0fGWbUinYL8sCMRadGSlpsgNfLTkIqZAOpUEgmb2k5VtG4NI0aoaBPZFfUp2jYA+1X5uTNQfeHWDcAzzrkC59yXwBKgR3JCjJ6hI8uIlxnLl6phJBIi5aZqOndxdD00wZLndF2bSMiUn6oZPx4+/BA++CDsSESap/oUba8Ah5jZgWaWCZwGPFltnyeAIWaWbma5QH/g38kNNTp69knQpq1jyUI1jERCpNxUg6Ej47y6Ikb+t2FHItKiKT9VM368//700+HGIdJc1Vm0OefiwKXAfHwyedg5946ZXWhmFwb7/Bt4BngLeBn4q3NudeOFHa70dBg8vIwlz6dTXh52NCItk3JTzYaOLCMeN15aopkAImFRfvquAw+E7t01RVJkZ9XrU905NxeYW+2x6dV+/h/gf5IXWrQNHRln3hOZrH4zxlG9EmGHI9IiKTd9V4+jE+zWtpwXFmYwenw87HBEWizlp+8aPx5uuw22bIG2bcOORqR50RJjO2nw8DixmGPxAvVmi0h0pKfDkGPjLH0unbhqNhGJkB/8AOJxmDcv7EhEmh8VbTupbTtH7/4JFi3QdW0iEi0jRpex+Zs03lgVCzsUEZFtBgyAvfaCJ54IOxKR5kdF2y4YMbqMtR/EWL9OL6OIRMeg4XEyMh2L1akkIhESi8HEiTB3LpSWhh2NSPOiamMXDB/l74m0SFMkRSRCWrWGfgPjLFqQjqt+ZygRkRAddxxs3QovvBB2JCLNi4q2XdC5i+OQwzVFUkSiZ8ToOOs/ibH2A6V5EYmOvDzIydEUSZGG0qf5LhoxuozXX4mx+RsLOxQRkW0qZgIsfladSiISHbm5MHo0PPkkmgkg0gAq2nbRiNFxysuNFxZqiqSIRMfe+zqO6BHn+fnKTSISLccdB+vXw+uvhx2JSPOhom0XdTsqwT4dynn+GfVmi0i0HDs2zuo30vnP55oJICLRMXEipKXB7NlhRyLSfKho20VmMHJcGS+9kE5hQdjRiIhUGjnOT5FUp5KIREn79jBsGDz6aNiRiDQfKtqSIG9sGSUlxrJFmoYkItFxYNdyuh6a4Ll5KtpEJFpOPBH+/W//JSJ1U9GWBL36JWi3ZzkL56phJCLRMnJsGa+ujPH1V5oiKSLRccIJ/rumSIrUj4q2JIjF4NgxcZY8n0FJcdjRiIhUyvt+GeXlpvtJikikdOoEAwZoiqRIfaloS5K8cWUUFhjLl6phJCLRcVi3cjp1KdcUSRGJnEmT4LXXYN26sCMRiT4VbUnSb2CcNm2dpkiKSKSY+SmSK5als3VL2NGIiFQ68UT//bHHwo1DpDlQ0ZYkGZkwYlQZixZkUFoSdjQiIpVGTygjXmYsWqBOJRGJjoMOgt694eGHw45EJPpUtCXRmIllfLtVUyRFJFqO7JmgY+dyFsxR0SYi0XLKKbBypaZIitRFRVsSDRgcZ7e25cx/Sg0jEYkOMz/atnxJOlu+0SqSIhIdp5zivz/ySLhxiESdirYkysiEY0fHeXZuBvffG2POHEgkwo5KRATGTCgjHjeen6+ZACISHQceCH36wD33wNSpqO0kUgsVbUmUSMC/30mjpNh46P4YkyfDmDFKPiISvm5HJejcJcF8TZEUkQhJJGDzZlizBm64AbWdRGqhoi2Jli1K59N1seAnIz/fz9OeNy/UsEREMPPX3a5cls6WzWFHIyLizZsHGzf6fzuH2k4itdA8mSR6b3WM4qLtHysocMyanSCzTfK6jNq2jtH/aP3XiUjDjJ5Qxt+mZbNsSRonHxd2NCLSXL33QQJIzlLZj8yOUVQUAyqvty0ogDfegAkTknIKkZSgln8SHX5EgpwcKCysfCw7B/oMKaFD13jSzvP52mz0XyciDXV493IOPDjB8wtjde8sIlKLwpJyOnQtS8qx+g5JZ9bM3O3aTq1aQc+eSTm8SMrQ9MgkGjwizpG94uTkOsCRFnMc1SvO4BHJK9hERHaWGYw/vozVb6Xx6adhRyMiUtl2ys7xbaeMDOjfH8aNCzsykWhR0ZZEsRhMf6iQ308rZP+DEuDgt3cWEVOntohExPePLwXgn/8MORARESrbTv9zVyEdOznatYNnnkFtJ5FqVLQlWSwGw/LiHHdyKeXlxnPztFKbiERH5/0d3+tezkMPhR2JiIhX0Xb64Y8SbNrkFyIRke2paGsk+3Yqp+uhCeY+rqJNRKLl2LwEb78Nb78ddiQiIpUGDi4nJwd1KonUQEVbIzGD7x9fxuuvpLPhU6v7F0REmsjQEeXEYvDgg2FHIiJSKScXjjsOHn4YSkvDjkYkWlS0NaIJJ5Zi5pjzWGbYoYiIbLP77v4i/wcf1A1sRSRapkyBr76CuXPDjkQkWlS0NaIOnRz9BiV48pEMnAs7GhGRSmed5W9ou3Bh2JGIiFQaNQo6dIAZM8KORCRaVLQ1sh+cVMqGT2O89rKWQRKR6JgwAdq1U8NIRKIlPR1++EN4+mnYtCnsaESiQ0VbIxs5rozcVo4nH9EUSRGJjqwsOP10mD0bNm8OOxoRkUo/+hHE4/CPf4QdiUh0qGhrZLm5MGZCGfPnZFBYGHY0IiKVzjoLSkr8Rf8iIlHRvTv07Qv33Rd2JCLRoaKtCfzg5FIKC4yFc7X8v4hEx9FH+8bRvfeGHYmIyPbOOgveeANefz3sSESiQUVbE+jdL8H+ByZ47J+aIiki0WEG55zjb2S7enXY0YiIVJo8GbKz4a9/DTsSkWhQ0dYEzOCEyaW89nI6H6/RSy4i0TFlCmRmwl/+EnYkIiKV2rWDk07yN9rW5SUiKtqazHEnlZGe7nhUo20iEiHt28MJJ8ADD0BxcdjRiIhUOu882LIFHnkk7EhEwqeirYnsuZdj+Og4T83KoLQk7GhERCqddx588w089ljYkYiIVBoyBA47TDMBREBFW5OaNLmUb75OY9ECLUgiItExYgQcdJAaRiISLWbw4x/Diy/Cu++GHY1IuFS0NaEBQ+J07FzOww9qiqSIREdamm8YLV4M770XdjQiIpWmTIGMDLj77rAjEQmXirYmFIvBSWeU8spL6Xz0oV56EYmOc8/1DaO77go7EhGRSnvv7RckmTEDCgrCjkYkPKocmtiJk0vJyHT8636NtolIdOy9N5x8sr+ZbX5+2NGIiFS65BLYutWvJCnSUqloa2J77OkYM6GMJ2dlUqCGkYhEiBpGIhJFAwdCjx4wbRo4F3Y0IuGoV9FmZmPN7H0zW2NmV+9gv75mljCzk5IXYuo5dUopBfnG07M12iayK5SbkuuYY6BnTzWMRJJB+Sl5zODii+Gtt+Cll8KORiQcdRZtZhYDpgHjgG7AZDPrVst+vwPmJzvIVHNU7wSHH5HgnzMy1TAS2UnKTcln5kfb3n4bliwJOxqR5kv5KfnOOAPatoU77ww7EpFw1GekrR+wxjn3kXOuFJgJHFfDfpcBjwKbkhhfSjKDM84uYe0HMVYsTQ87HJHmSrmpEZx+OuyxB9x+e9iRiDRryk9J1qqVXzBp1ixYvz7saESaXn2Ktk5A1bfHhuCxbcysE3ACMD15oaW2cceVsede5TzwV02RFNlJyk2NIDcXLrgAnngCPvoo7GhEmi3lp0Zw2WV+6va0aWFHItL06lO0WQ2PVZ/UdzvwS+dcYocHMjvfzFaZ2aovvviivjGmpMwsf23bskUZWv5fZOckLTeB8lNVl1zi7912xRUwdSrMmQOJOl9BEalCbadGcMABcOKJcM89Wv5fWp76VAsbgP2q/NwZ2Fhtnz7ATDNbB5wE3GVmx1c/kHPuHudcH+dcn7322msnQ04dJ/+wlMwsx0P3arRNZCckLTeB8lNV++4L7dvDU0/BDTfA5MkwZowKN5EGUNupkfz0p/DNN/72JCItSX2KtleAQ8zsQDPLBE4Dnqy6g3PuQOfcAc65A4BZwMXOuceTHm2K2bO9Y/wJZTw1K5Nvvq6pU05EdkC5qZHMm+eX/gc/FSk/H1au9I+LSL0oPzWSY46Bfv38dbfqSJKWpM5VMJxzcTO7FL+yUQy41zn3jpldGGzXXOxdMOW8EmbPzOQff8/kkp+XhB2OSLOh3LTz3vsgAdSebx6ZHaO4OEbVGV4FBY5ZsxNktqm9ldS2dYz+R2txJRHlp8ZjBldeCaecAo8/DpMmhR2RSNOo16erc24uMLfaYzUmHOfcWbseVsvR9dByRowuY+aMTM6+sITcVmFHJNJ8KDftnMKScjp0Lat1e98h6cyamUthYeVj2TnQZ0gJHbrGa/29z9dmU8+PFZGUp/zUeE48EQ4+GG691f/bNFlJWgCtgBEB51xSwpbNaTz6T13bJiLhGzwizpG94uTkOsBh5jiyV5zBI2ov2EREmkosBlddBatWwaJFYUcj0jRUtEVAj94Jju4f54G/ZFFWGnY0ItLSxWIw/aFCfj+tkNETynDOOG1KKbFY2JGJiHhnnukXTbr11rAjEWkaKtoi4pxLSvjPxjTmPJYRdigiIsRiMCwvzq13FtG5S4J778rCVV+wXEQkJNnZfiXJZ5+FV14JOxqRxqeiLSIGD4/T7cg4f/1zFnHNQBKRiEhPh3MvLWH1m+ksW6zr1UQkOi66CPbYA26+OexIRBqfiraIMIMLrihh/Scxnp6t0TYRiY4fTCqjY+dypv9Jo20iEh1t2vjRtjlz4LXXwo5GpHGpaIuQ4aPiHN49wV/u0GibiERHRiace0kJb7+ezksvaLRNRKLjsstg99012iapT0VbhJjBBT8p5tN1MeY9odE2EYmO408pZd+O5Uy/XaNtIhIdbdv60bYnnoDXXw87GpHGo6ItYkaMiXNYtwT/d1sWZbXfRklEpEllZMKPLy3hzVfTWbZIo20iEh2XX+5H266/PuxIRBqPiraISUuDy64qZsOnMR7/l+7bJiLRccKppXTukuCO32VTXh52NCIi3u67+/u2zZkDL74YdjQijUNFWwQNOTZOr75x7r49i+KisKMREfEyMuGSK0t4/90YC+ZoCreIRMfll8M++8Cvf42mcEtKUtEWQWZ+tG3Tf9OYeZ9G20QkOsb+oIyDD0sw7Q+awi0i0dGqFVx3HSxZAgsWhB2NSPKpaIuoPgMSDBpexl//nMXWzWFHIyLixWK+U+mTj2PMnqlOJRGJjvPOgwMPhF/+EhKJsKMRSS4VbRH2018X8+1W4547ssMORURkm+Gj4vTuF2faH7PI/zbsaEREvMxMuOUWePNNeOCBsKMRSS4VbRF26PfKOf7UMv4xI5P16/RfJSLRYAa/uL6Yb75K4967ssIOR0Rkm1NPhf794ZproKAg7GhEkkeVQMRdemUx6elw+60abROR6OjeI8H4E0p54C9ZfP6ZhR2OiAjgO5X++EfYuNF/F0kVKtoibq99HGdfVMKzT2ewanks7HBERLa5/JfFAPzpFnUqiUh0DBoEkybB734H69eHHY1IcqhoawbOurCEDp3KueW6HJ6fn85D98eYM0cX2YpIuDp0cpx9cQnPPJmpTiURiZQ//MG3k844A6ZORe0mafZUtDUDOTnws2uLWPN+jF9cnMsDf48xeTKMGaMEJCLhOueiEjp2Lue31+eQiIcdjYiIt99+0LEjLF0KN9yA2k3S7KloayaysiAtzVFWajhn5OfDypUwb17YkYlIS5adA1deX8SH78V46gmNtolINMybB5s2+X87h9pN0uylhx2A1M/778RwbvvHCgocs2YnyGyTnG6jtq1j9D9afxIi0jAjx8YZOLSM+/6WzrVX+95tEZGGeu+DBFCSlGM9MjtGYWEMqFwoqaAA3ngDJkxIyilEmpRa6M3E4UckyMmBwsLKx7JzoM+QEjp0Tc6cpM/XZqM/CRFpKDO45jfFnDCyNT/5CTzySNgRiUhzVFhSToeuZUk5Vt8h6cyambtduyknB3r2TMrhRZqcpkc2E4NHxDmyV5ycXAc4zBxH9IgzeIQuIhGR8O13QDmnT0kwa5a/4F9EJExV201mfqpSTg6MHRtyYCI7SUVbMxGLwfSHCvn9tEJ694/jnHFEzwQxXUIiIhFx0qkJuneHiy+GrVvDjkZEWrKq7aYpZyf48Y/hq69g1qywIxPZOSrampFYDIblxTn5hyUcf0op99+TxTtvqmoTkWjIyIC//hU++wyuuirsaESkpatoN50+JcH06dC/P1x6aeUCJSLNiYq2ZurK64vYs73jup/nUJqca3ZFRHbZgAHws5/B3XfDwoVhRyMi4sVi8Pe/+1UkL7kk7GhEGk5FWzO1W1u4/nf+3m3Tb88KOxwRkW1uvhkOPRTOPVfTJEUkOr73PbjpJj9F8l//CjsakYZR0daMDR0Z57iTS7n3rixef0XTJEUkGnJyYMYM2LABLr887GhERCr9/Od+muRFF/kcJdJcqGhr5n55UxEdOjt+fUUu+d+GHY2IiHfMMXDNNXDffboFgIhER3o6PPAAlJbCWWdBeXnYEYnUj4q2Zq51G7jl9kI+32Dcen1O2OGIiGxz3XXQrx9ccIF6tEUkOg45BG6/HZ57zn8XaQ5UtKWAXn0TnHdZCU/OymTOYxlhhyMiAvjVJB980Pdon346xHVbSRGJiHPPhRNOgKuvhldeCTsakbqpaEsRF1xRQu9+cab+KoeP1+i/VUSi4ZBDYPp0WLoUbrwx7GhERDwzf4uSDh3g1FNh8+awIxLZMbXuU0R6Ovzuz4VkZzuuvCiX4qKwIxIR8X74QzjnHLjlFpg/P+xoRES8Pfbwq0iuX+9H3pwLOyKR2qloSyH7dHD85vYiPnwvxs2/ylHyEZHIuPNO6NbNT5P8299g6lSYMwcSibAjE5GWbMAAuPVWeOwx+MMfwo5GpHbpYQcgyTV4RJyLflbM/92WzRE9Epx+dmnYIYmIkJvr7410xBFw/vm+R7tVK7/09vz5/sa3IiJh+NnPYMUKf31br16Qlxd2RCLfpZG2FHTBT0oYllfGH27OZtUKtYREJBrWrPGLk5SX+6ItPx9WroR588KOTERaMjO49144/HA47TRYty7siES+SyNtKSgtDW7530JOn9ian52fyz+eyqfz/porKSKN670PEkBJrdsfmR2jpCQG2LbHCgocs2YnyGxT+zzJtq1j9D9aH1ci0njatIHZs/1tSiZOhBdfhN12CzsqkUr6FExRbXaDP88o5IyJrbjsnFbcPzufNko+ItKICkvK6dC1rNbtfYekM2tmLoWFlY9lZEKfISV06Fr7/QA+X5uNPq5EpLEdeqifxj12rL/+9oknNHVbokPTI1PY/geWc9s9hXzyURq/uCiXstrbUiIijW7wiDhH9oqTk+swc1iao6wM9thDMwFEJBry8vzCSU8/7a9106JuEhUq2lJcv4EJrr+1iJeWZHDTVVpRUkTCE4vB9IcK+f20Qi75eQlT/1jEvh0cl52by4ZPrO4DiIg0gYsugiuugDvugNtuCzsaEU/zTVqA408t4z+fF3PXH7PZe99yLv9l7deciIg0plgMhuXFGZbnp0Me2TPBlBNacdGZrbhvdgF77KmeJREJ3x//CJ99BldeCR07wuTJYUckLZ1G2lqIC35SwqTTS/nrn7O5/57MsMMREQHgwIPLuePeQv6zMY2LzmzFt1vDjkhExC/qdv/9MGwYTJkCc+eGHZG0dCraWggzuPaWIkaNL+MPU3N4bGZG2CGJiADQq2+C2+4p5MN/p3HZ2a0oKgo7IhERyM72i5EcdRRMmgRLloQdkbRkKtpakFgMbr2jkEHDy7jpqhyenq3CTUSiYcixcW753yJefyXGT3+cS0lx2BGJiEDbtjB/Phx4IEyYAMuXhx2RtFT1KtrMbKyZvW9ma8zs6hq2n2FmbwVfL9C7VMMAAB4KSURBVJlZj+SHKsmQkQm33VNIn2MSXHOFCjdp3pSbUsvYH5Rx0/8UsXxpOj9R4SbNnPJT6mjfHp59FvbZB8aMgRUrwo5IWqI6FyIxsxgwDRgFbABeMbMnnXPvVtntY2CYc+4bMxsH3AP0b4yAZdfl5MCdfy/g0h+14porckgkoO3ujleWxij9FsaN031JJPqUm1LT8aeW4Rzc8ItcfvLjXP70l8K6f0kkYpSfUk+nTrBoEQwf7gu3OXNgyxZ4/XXo1UttJ2l89Vk9sh+wxjn3EYCZzQSOA7YlHufcS1X2XwF0TmaQkny5ufDn+3zhdu1Pc8jIhHgZPPYw9O/vpwIo+UjEKTelqBNOKwMKufGqHC6Z0oprrk+EHZJIQyk/paDOnX3hduyxvnjLyoLiYmjVSm0naXz1mR7ZCVhf5ecNwWO1OReYtytBSdPIzYXJZ5eQFoOyUsM5Iz8fVq6EefoflOhTbkphJ5xWxq13FvHGqhhX/zyDL78MOyKRBlF+SlH77QfXX+9vul1U5L+r7SRNoT4jbTXd8bTGG+mY2Qh84hlcy/bzgfMBunTpUs8QpTF99EEMV779YwUFjlmzE2S2SU7vdtvWMfofrVsCyv9v787jo67OPY5/zsxkD5UdIpuABBAFRDZrKKLIooBarFfcqpaqV9vaXm21datLpW6t0GqRK8WqCFVAL1ZUxCKCKKsRQQQBF7bKZi2ThExm5tw/TsbMsGSBmcwk+b5fr98rM2Fm8ojwcM75Pec5cRe33FT+GuWnFDPy/DKysi23XJ9NQYFbxe7QIdlRiVSLxk4p5JONISB+Z9QueNsLeIn+31xUBIWFrlmJSCJUZyS9DWgX9bwtsOPgFxljegJPASOttXsP90HW2im4mm369u2rE1RTQLeTQ2RlQXHUthGPB3oPLCWvczAuP2Pn5kx0jrskQNxyEyg/paozzwky4ZEy7rsrndNPd2cl9e6d7KhEqqSxUwopLg2T17ksbp/Xb5CPWTOzY8ZOaWnQS61kJIGqUx65AuhijOlojEkHLgHmRr/AGNMemANcYa3dGP8wJVEKhgQ55dQgWdkWsKSlWUIhw9+mZPD1vsMtFIqkDOWmBuLknpbFi91ekYIC1wBAJMUpP9Vj0WMnYyw+HwQCMHMmlMbvhp5IjCpvf1hrg8aYnwBv4O4F/9Vau84Yc335r08G7gKaAU8YYwCC1tq+iQtb4sXrhcnTi1my0MfcWT7GXBSkuMhw581ZXD4mh0l/LaZzfrjqDxKpZcpNDcvJJ8Py5TB6NJx/Pjz8MPziF2C0tiQpSPmpfoseO61cnMHYC3wUFrq9bl9+CbNnQ8uWyY5S6ptq1axZa+cB8w763uSox+OB8fENTWqL1wuDhwbx+y2Dh7p9bMe3C/Pz8dlcNiaXCROLGTI8PqWSIvGk3NSw5OXBokXwwx/CzTe7VttTprhjTERSjfJT/RYZO+V39DFssI/Ro6FrV7jqKujbF156CU47LdlRSn1SrcO1peHp1SfEjFf9dDoxxE3jc/jzwxmE1HVbRJIsJwdeeAHuvReee86VS372WbKjEhGBiy+GJUvc44ICmDYtufFI/aJJmxxR6zzLtFlFXPhfAaZMyuT6y3PYu1u1SCKSXB4P3HknzJ0LmzdDnz7ucSjk9rvdd5/7qoUmEaltffrAypXw3e/CNdfA+PHuaACRY6VJm1QqIxPueaSEex4upnCFl4uG5/LeO+oEKSLJN3o0rF4NnTq5fW4dO8Ill8Ddd8O4cTB8uCZuIlL7WraE+fPh9tth6lTo1w/Wrk12VFLXadIm1XLhJWVMf8XPcY0t112Wwx9+l0lAHZJEJMk6dYJ333VnI23d6s5K0mG3IpJsXi/cfz+8/jrs3u0mbk884fKTyNHQLROptvzuYWa86ufhezJ5enIGSxf5eGBiMfnd1V1SRBKnOgfjNm15uMNuLbNeCpHe6Mi3247L9TLgNP1TKCKJMXw4rFnjGijdeKMr5Z46Fdq0SXZkUtfoTpvUSFYW3PX7A/xpWhF79xguOS+XJx/LoCx+Z1aKiMRwB+MeqPTqN6iU7OzY9xkD7fMDlb7vG7/qJ0UksVq1cnf9H38c3nnHHWHy9NO66yY1o0mbHJXBQ4PMWeDnnHPLePzRTC4dlcvaQm+ywxKRBurgw27T0y3GwIN3ZTFtcroWlkQkqYyBG26ADz90k7arr4aRI9X9VqpPkzY5ak2aWh78cwkTnypi3x7DZWNymHBnJvv/k+zIRKShiRx2+9Djxdx4cymPPlnMvHf3893BQf74uywuHZXLmg+0sCQiydWliztvctIkdzxAjx7w0ENoYUmqpEmbHLMhw4O8vHA/l/wwwMy/pTPmzEbMnZVGWRksWuBj+jNetd8WkYSLHHZ77U2lDB4aJK+NZeLUYv44pYiv9xkuH5PL3b/MYt9eHV0iIsnj8cBPfwrr17s9b7feCr16wVtv6egSOTLtvpa4aPQd+PV9BxhzURkP3JnJHb/I5vd3hSkrMwRKYc4LMGAAvPGGG1iJiNQGY+DskUEGDtrP5McymT41nQXz0rjupgOMuyqQ7PBEpAFr1w5eegleeQV+/nMYOhRatHBdcEtKICdHYyepoDttElc9eoV49uUiLrvmAP79htIDBmuN2m+LSFLl5MLNdxxg1nw/p5wa5JH7srjgrFwWv+1RMwARSarRo2HdOrj8cnc8QHGxji6RQ+lOm8SdxwONm7gV7ujBkN9veeb5yttv15TadYtITXTqEmbyc8UsWejj0fszuf+3abzxKkyYAGed5fKWiEh1VedIkurypHsxxou1FYnI74fly91ZlNKwabQrCdHt5BBZWW61KNrLs7zkHBfiqutLadrs2Je3d27ORH+MRaSmCoYEOf17fp77SzYvPJ/G0KEweDDcey9873vJjk5E6gp3JEl8uoj0G+Rj1szsQ8ZOEyfCccfBf/83hxxtIg2HyiMlIaLbb4MlK9vS67QQw0aV8bcn0xlxeiMeuS+T3V9pWVtEksPrhWEjw2zc6AZFn3ziJm5nnukaAqhsUkRq08FHl+TmwmmnQf/+cMstcMIJ8OCDsH9/siOVZNCkTRIiuv32OecFeOjxYp6eXcSESSW89JafoSPLeO6pdEZ8txH33JrF51v0R1FEkiMzE372M9iyBR57DD791DUE6N8fZs2CQEDd3EQk8aLHTldeHWLGDLen7c03YfFiOPVUuO02aN8ebr8dvvoq2RFLbdJIWRIm0n77rBFlDB4a/LbzUacuYR6YWMLct/1ccHGAV2ancf6ZufzsmmxWvufV6raIJEV2Ntx0E2zeDJMnw9dfww9+4MqSxo6Fu++GceNci25N3EQkESJjp0uvDDFqVEXXyIIC10Vy2TK3/3bCBOjQAX70I1i7NrkxS+3QpE2Spn3HMHdOOMDrS/fz45+WUrjSyzUX5/KD4bnMej6NkpJkRygiDVFmJlx3HWzY4Fa1AwF3Rbq5vfeeurmJSHL07w+zZ7ty7quvhhkz4JRT3ERuzhwIBpMdoSSKOjhI0jVvafnJL0v50U9KefWlNGZMy+DeW7P54+8so8cGGHtpgC7dwskOU0Tqoao6v+3Y5cVaL1Cx/7a4GH58bZhrrg1x+hlh0tMPfZ8624pIIuXnw1/+AvffD//7v/DEE64ioE0bGD/e3YFr1y7ZUUo86U6bpIysLLjo0jJmzfczbZafQWeV8eL0dMae04hLR+cwa3oa//nGvTYUgkULfEx/xqs9JiJy1FzntwNHvPoNKiUrK/Y9aWmWQAAeuCeNKy9J49lnYX8gEPO+b/xKSiKSeM2auYqALVvcQd2nnOK64HboACNHwosvwoEDbpykvbl1m5YBJeUYA6cNCHHagBJuvecA/5idxpyZ6dx7Wza/v9ty5tAyPt/iYesXXg6UwJwXYMAAV+sdqf0WEYmHSDe3NR/4OFACmVnQ89Qgjz9TzPJ3fbz093RmPpPOc1Mz6NYjxKjvBxgxJj7tv0VEqsvngwsucNdnn8G0ae66+GK3Lzc3F/btcxO4nByNm+oi3WmTlNakqeWKHweYs8DPjH/4GTsuwLtv+9i43kdJscFag98P77+vPSYiEn/R3dxuvLmUhx4vZvL0YtLT3YTu0cnFvLVyP7++rwSv1/LIfVmc078Rt/5PGlOmwJ49yf4vEJGGpmNHd7ft889h/nzo3Ru2b4eSkoq9uUuXwrx5yY5UakJ32qROMAZ69ArRo1eI45pYJj+WAbZij0lRkeXmWyxr1gfpdaolLe3Yfp72o4hIRKSb2+Chh9/h36SpZdxVAcZdFeCzzR7mvZzGq7PTue46uOEGd+7b2LFuBTwvr+J9oZBbbPrgA9fKe+RIrXqLiFPVfttqS4cOnb2Yd7zYqHFTSQlceaXb/zZ2rGtw4tGtnJSmUanUOT16hsjOcs0AIrxe+OJzw+2/Sie3kaVgiDtmYNCQMr7TuOY/Y+fmTPTXQ0RqqmPnMDfeXMqF5xtaNc7ghRdcp7cbbnDXwIEwZoyboN18MyxfDkVFKlcSkVhuv218Sq37DfIxa2Z2zLgpIwM6dYKJE+GRR6B1a5ebRo+Gs8/mkL28knyaU0udE9ljkpVtAUtWtqXvwCBvF/6HP00r4pzzyli+1Mevf5bN4N7f4Yffz+GpP2ewfq2HsJpQikgtMMaVJD3wgGvN/dFHrstbMAi/+Y27s7ZwoStTipQrLVumMm8Rib/ocZMxltxcd+7b8uWwaxdMn+6eP/+8m7Q1bQrnnQd//jNs2oTOz00RupUgdU5kj8mShT7mzvIx5qIgBUOCMSVM4TCs/dDLojd9LHk7jUkPZjLpwUyaNg8zsCDIwIIgAwqC5LWJzUShECxZ6GPFYi+B/SpXEpGjc7jSpn4F7tq7Byb90cf773qIPkrA77fceXeYTV+G6NrV4j3Mv9Aq3RaRmooeN61cnMFFF/q+Hd80bgyXXuqu0lJ45x3XXfK11yr2vHXsCMOGuTtwZ53lOlZGU6l37VDmlzopMkHz+y2Dhx7at9bjgZ6nhuh5aoif/qqUPbsMS9/xsXSRj/eX+Jj3sjtYqV2HEP2+G6LvwCCn9gvy219m8dEHPkrUlVJEjkFlpU15neGy8T7WfBBbrmQMFK72ULjaS3aOpU//IH1PD9J3YIjuJ4dIS1Pptogcnci4Kb+jj2GDD59DMjLgnHPcNXGiu8s2f74bBz3/PDz5pHtdz54wZAgMHgxnnOEmfMuWqdQ70ZT5pUFo3tIy5qIyxlxUhrXw6Scelr3rY8VSH2++msacGW4SZ4z9dqNudFfKUaOSGb2I1DdHOkpgwp9KWL3My/KlPla852PJQrexJCvb0rNPiBM7A6VuUHTcccn9bxCR+u3EE911ww2utHv5cvjnP11p95NPuokduAWnSAlldKm3xk7xpUmbNDjGQH73MPndA1wxPkAo5CZxkx7MZMnC2L8SRUWWH19rOfPsMF27Wbp2Dx9SFlATKm0SEYgtV9qwzkvXHqFvy7yHjQoybJTrVLl3t2HVMi+rlvn4YIWP55/xMP1v5R11e7jJ24ABrvNbjx7urCZQuZKIHN6xdqXsP8hdgQB8usHw7NNePlh1cKk33HILfPihy019+0KTJscee0On0aM0eF4vdOsR5r+uDLB6mS+mXMnnA1+65cUZXkIhl5Batg5z0ikhTurpSpa6nxyiRSuLMUf4AVFU2iQiEVUdJQDQrIWNmcRtWpNJbloGS5e61ew5c2DqVPfarCzX/CTS5OSLL1xbb5UriUhEPLtSdugOabk+br0x+5CxU1ER3HFHxfdOPBFOO81dffq4PNW0aVzCaDA0ehQpF12uVFIMWdmuXGny9GICAdiwzstHhV4+XuNl3Rovixb4vi2lbNIsTLeTQuSfFCa/e4j87iE6dg6TnuE+Ww1ORCQecnLgnMFuzwm4kqRNm2DFCnetWuUmcaVRC+l+PyxaBNdfD+PGwSmnQIsWyYlfROqXg0u9c3LMt4tE+/fDypWurHLVKnjvPfj73yve27499Orlrp49XW468URVDByJJm0i5SrrSpmVBb37hujdt6LpSZEfNnzsZf1aLxs/9rJ+nZcZT/sIlJryz7O07ximU5cwn6z1sHuXh0AAZv/dndWkVW8RqakjlTY1bwMj28DIC+DZp71M/1vsQbrBoOWppwxPPeWeN2liaX+CpcMJlg4nhGnfwdLjJA8jhvoqrRrQIEpEolXVmXLoUHdF7Nnj8kdhYcXXefNcbgHIzIRu3aB7d1dNsGOHW4RSxYAmbSIxqupKGS0nF/r0D9Gnf8XrgkH4YouHTz/x8uknHjZv9LJmtYfdX1XUexcVuY28AwbAoEHQtSt06eKutm1d50sRkcOpTmnTgME+5rwQW66UlQ13/K6YFq3C5fnJy6aNHhbM91JcVDECatrU5aTIlZ/vclPnzq6z3PDh6hInIrGq05kyonnzig6VEQcOwMcfu/MsP/oI1q2DN990E7wIv9+NnQoKXNfK/PyK/NSyJdXaolLXadImEkc+H3TOD9M5P8yIMe57Tz6WwRN/yIg5nNJay5bPLGvXGkpLKzJNerol73hLXhtLXl754+OhdZ6lVWtLRsahPzMUghXLPGz/0sPo87xa+RZp4I7UmfLcC8vwemHgoIqFJmvhq52GzRu9FL6fgS3zsWEDvP46PP107Oc2awZffw3hsHvu98PSpTBzpmv53RAGTSJyZMfa5CTvBHcNGw3Tn/Hy7LTYioFIOfiqVVAWtXbVqJErq+zc2V2dOrmrY0dXgpmWdujPqotVA5q0iSRYt5NDZGVxyKr3vY+WMOisILu+MnyxxcOXn3n58nMPX37mYesXHgpXeTlwIHYU1LxlmOPbVlyt8ixzZqTxxRYvpaXw/HNa+RZp6CrrTHkwY6D18ZbWxwcJlBi65UeVgBfBju2G7VsNO7YbFr/jYe9eQ3SXuJISuPxyGD/eLSy5C1q2st9enTt6OHeYr9KcVBcHUCISK55NTvoN8jFrZmzFQG4uTJsGI0a4RksbN7pJ3Kefuq9r1sD//V/shM7jcVVMJ5wAHTq4q317eOIJ9/661KxJkzaRBDtSg5OCIUE8HncXrXVeiAFnxJZjWuvafW/70sP2re5O2vatHrZv87C20MuC19IIlsVO6iLlA/36uYFPmzYV1/HHu6tFi9ROSiJy7KrTmfJghxtwndiz4nHfBYd2iUvPsIz6foCsbNix1cOObR7Wr/ew/5vY3OTzufzTti20axebm1q3ht/8xrUHLy6uOwMoEUmcIzU4iSzoRO6qHSwUgu3bYcsW+Owzd33+ubsWLnR75CLVAhF+v/u1gQNdd8vDjZ2aNq1eNUEiF6A0aRNJsMoanFTGGHcoePOWsQ1QIkIheGxCBs9MyTiofMCydZvly62Gr/dBOBybZTweS5Om0LSppWmzisdNmrrHTZpYGjdxX7NzoHEjnS0nIkcuu7xzwoFD8pl/P+zc7uFfOzx8UphOo6w0tm51g6nVq+GVV2KrD2Le64e334azz3YDqNat3dWqVcXVvHlFhzkRqX+ix06vvZzGgAGGfgPCvLWkmh9goE0ndxVEfTsYhKlPepnzopfoqoFw2OWnL76A3bsP/bj0dJeH8vIqrkheinxt1gyuvdZ1zEzEvl+lPJFaUJMGJzX5zL4DQ7z47KGll3c/VMLgoUGCQdiz27D7Xx52/cuw+ysPu3cZ9uyq+Lppk2HfHk/MxC8iLd3SooVl+7a4hCwidVhNyi5zG0GXbmG6dAsTLDV0y49d2rbWdeDds8fw4gzXECV6ABUKWVathvffJ2bfb4QxlkaN4Lbbw/z6V7odJ1IfRY+dzrsgPmMngCHn+nj91UNLL6dMgVGjXLfKnTvdJG77dvd4xw73dedOV4q5eDHs3Vv5z/H7XeOm115zn3usNGkTqcMqK70EtxIdKb+sTCgE/95n2LPbsHe3h717DPv2GPbt9fDNHi+glpYiEr+yy4guQAAfS5cc2u1ywqRivnd2kCI/7N3jYe9uw949hr27POzb53LUdxobQJM2Eam+g8dOmVlwYn4Yb3aQ+YtiX9u4lbu69z70c8rK4Jt/w9f7DF9/bXjtHx6WLqnoFg7ujlthoSZtIg3e0ZZeHu5zmrWwNGthgdgV8Z2bM1GqEJFEOVLZZcGQIMa4u3a5jcJ06Hjoe11+EhGpvniNnQ7WuJWPwtWxC1A5OdD7MBO+o6GRmEgdl4jSSxGR2lKTsksRkXhIxNipsuYp8aBJm4iIiCTV0ZRdioikkugFqJWLM7joQp+6R4qIiIiIiKSSyAJUfkcfwwbHd5ql7gIiIiIiIiIprFqTNmPMCGPMBmPMJmPMbYf5dWOMmVT+62uMMX3iH6qISCzlJhFJVcpPIhJPVU7ajDFe4HFgJHASMM4Yc9JBLxuJ69zbBbgW+Euc4xQRiaHcJCKpSvlJROKtOnfa+gObrLVbrLUBYCZw/kGvOR94xjrvA42NMXlxjlVEJJpyk4ikKuUnEYmr6kza2gBbo55vK/9eTV8jIhJPyk0ikqqUn0Qkroy1tvIXGPMDYLi1dnz58yuA/tban0a95lVggrV2Sfnzt4BfWWtXHfRZ1+JKAAC6AhuqGWdzYE81X5sqFHPtqGsx17V4oWYxd7DWtkhkMBHxzE3lv6b8lLrqWrygmGtLvc9Pyk0pr67FXNfihfofc7VyU3V6UW4D2kU9bwvsOIrXYK2dAkypxs+MYYxZaa3tW9P3JZNirh11Lea6Fi+kdMxxy02g/JTK6lq8oJhrSwrHrLHTUVDMiVfX4gXFHFGd8sgVQBdjTEdjTDpwCTD3oNfMBa4s74Q0EPjGWrsznoGKiBxEuUlEUpXyk4jEVZV32qy1QWPMT4A3AC/wV2vtOmPM9eW/PhmYB5wLbAKKgasTF7KIiHKTiKQu5ScRibdqHdVtrZ2HSy7R35sc9dgCN8Y3tBg1LgtIAYq5dtS1mOtavJDCMadAboIU/v2pRF2Lua7FC4q5tqRszCmQn1L296YSijnx6lq8oJiBajQiERERERERkeSpzp42ERERERERSZKUn7QZY0YYYzYYYzYZY25LdjxVMcb81RizyxizNtmxVIcxpp0xZqExZr0xZp0x5qZkx1QVY0ymMWa5MebD8pjvSXZM1WWM8RpjPjDG/CPZsVSHMeZzY8xHxphCY8zKZMeTSpSbEk/5qfYoN9Uvyk+JpdxUu5Sfyj83lcsjjTFeYCNwDq417gpgnLX246QGVgljzPcAP/CMtfbkZMdTFWNMHpBnrV1tjGkErAIuSPHfYwPkWGv9xpg0YAlwk7X2/SSHViVjzP8AfYHvWGtHJTueqhhjPgf6Wmvr2vkoCaXcVDuUn2qPclP9ofyUeMpNtUv5yUn1O239gU3W2i3W2gAwEzg/yTFVylr7DrAv2XFUl7V2p7V2dfnj/cB6oE1yo6qcdfzlT9PKr9RdfShnjGkLnAc8lexY5JgpN9UC5afaodxU7yg/JZhyU+1RfqqQ6pO2NsDWqOfbSPG/FHWZMeYE4FRgWXIjqVr5rfJCYBfwprU25WMGHgN+BYSTHUgNWGC+MWaVMebaZAeTQpSbapnyU0IpN9Uvyk+1SLkp4ZSfyqX6pM0c5nspvypQFxljcoHZwM+ttf9JdjxVsdaGrLW9gbZAf2NMSpdTGGNGAbustauSHUsNnWGt7QOMBG4sL2ER5aZapfyUOMpN9ZLyUy1Rbkos5adYqT5p2wa0i3reFtiRpFjqrfLa5tnAdGvtnGTHUxPW2n8DbwMjkhxKVc4AxpTXOc8EzjLGPJfckKpmrd1R/nUX8BKu7EaUm2qN8lPCKTfVP8pPtUC5qVYoP0VJ9UnbCqCLMaajMSYduASYm+SY6pXyjalTgfXW2j8kO57qMMa0MMY0Ln+cBQwFPkluVJWz1v7aWtvWWnsC7s/xP621lyc5rEoZY3LKN1hjjMkBhgF1orNXLVBuqgXKT4mn3FQvKT8lmHJT7VB+ipXSkzZrbRD4CfAGbpPnC9badcmNqnLGmBnAe0BXY8w2Y8yPkh1TFc4ArsCtXhSWX+cmO6gq5AELjTFrcP84vWmtrRNtYOuYVsASY8yHwHLgVWvt60mOKSUoN9Ua5Sc5HOWmSig/1QrlJjmShOWnlG75LyIiIiIi0tCl9J02ERERERGRhk6TNhERERERkRSmSZuIiIiIiEgK06RNREREREQkhWnSJiIiIiIiksI0aRMREREREUlhmrSJiIiIiIikME3aREREREREUtj/AxvCG/ziovBGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = lambda x : 1/(1+x**2)\n",
    "a = 0; b = 5; N = 10\n",
    "n = 10 # Use n*N+1 points to plot the function smoothly\n",
    "\n",
    "x = np.linspace(a,b,N+1)\n",
    "y = f(x)\n",
    "\n",
    "X = np.linspace(a,b,n*N+1)\n",
    "Y = f(X)\n",
    "\n",
    "plt.figure(figsize=(15,5))\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "plt.plot(X,Y,'b')\n",
    "x_left = x[:-1] # Left endpoints\n",
    "y_left = y[:-1]\n",
    "plt.plot(x_left,y_left,'b.',markersize=10)\n",
    "plt.bar(x_left,y_left,width=(b-a)/N,alpha=0.2,align='edge',edgecolor='b')\n",
    "plt.title('Left Riemann Sum, N = {}'.format(N))\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "plt.plot(X,Y,'b')\n",
    "x_mid = (x[:-1] + x[1:])/2 # Midpoints\n",
    "y_mid = f(x_mid)\n",
    "plt.plot(x_mid,y_mid,'b.',markersize=10)\n",
    "plt.bar(x_mid,y_mid,width=(b-a)/N,alpha=0.2,edgecolor='b')\n",
    "plt.title('Midpoint Riemann Sum, N = {}'.format(N))\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "plt.plot(X,Y,'b')\n",
    "x_right = x[1:] # Left endpoints\n",
    "y_right = y[1:]\n",
    "plt.plot(x_right,y_right,'b.',markersize=10)\n",
    "plt.bar(x_right,y_right,width=-(b-a)/N,alpha=0.2,align='edge',edgecolor='b')\n",
    "plt.title('Right Riemann Sum, N = {}'.format(N))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that when the function $f(x)$ is decreasing on $[a,b]$ the left endpoints give an overestimate of the integral $\\int_a^b f(x) dx$ and right endpoints give an underestimate. The opposite is true is when the function is increasing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the value of each of the Riemann sums:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Partition with 10 subintervals.\n",
      "Left Riemann Sum: 1.61348869661\n",
      "Midpoint Riemann Sum: 1.37354342832\n",
      "Right Riemann Sum: 1.13271946585\n"
     ]
    }
   ],
   "source": [
    "dx = (b-a)/N\n",
    "x_left = np.linspace(a,b-dx,N)\n",
    "x_midpoint = np.linspace(dx/2,b - dx/2,N)\n",
    "x_right = np.linspace(dx,b,N)\n",
    "\n",
    "print(\"Partition with\",N,\"subintervals.\")\n",
    "left_riemann_sum = np.sum(f(x_left) * dx)\n",
    "print(\"Left Riemann Sum:\",left_riemann_sum)\n",
    "\n",
    "midpoint_riemann_sum = np.sum(f(x_midpoint) * dx)\n",
    "print(\"Midpoint Riemann Sum:\",midpoint_riemann_sum)\n",
    "\n",
    "right_riemann_sum = np.sum(f(x_right) * dx)\n",
    "print(\"Right Riemann Sum:\",right_riemann_sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "We know the exact value\n",
    "\n",
    "$$\n",
    "\\int_0^5 \\frac{1}{1 + x^2} dx = \\arctan(5)\n",
    "$$\n",
    "\n",
    "and we can compare the Riemann sums to the value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.37340076695\n"
     ]
    }
   ],
   "source": [
    "I = np.arctan(5)\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Left Riemann Sum Error: 0.24008792967\n",
      "Midpoint Riemann Sum: 0.000142661371648\n",
      "Right Riemann Sum: 0.2406813011\n"
     ]
    }
   ],
   "source": [
    "print(\"Left Riemann Sum Error:\",np.abs(left_riemann_sum - I))\n",
    "print(\"Midpoint Riemann Sum:\",np.abs(midpoint_riemann_sum - I))\n",
    "print(\"Right Riemann Sum:\",np.abs(right_riemann_sum - I))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "## Error Formulas\n",
    "\n",
    "A Riemann sum is an approximation of a definite integral. A natural question arises: how good of an approximation is a Riemann sum?\n",
    "\n",
    "---\n",
    "\n",
    "**Theorem.** Let $L_N(f)$ denote the left Riemann sum\n",
    "\n",
    "$$\n",
    "L_N(f) = \\sum_{i=1}^N f(x_{i-1} ) \\Delta x\n",
    "$$\n",
    "\n",
    "where $\\Delta x = (b-a)/N$ and $x_i = a + i \\Delta x$. The error bound is\n",
    "\n",
    "$$\n",
    "E_N^{L}(f) = \\left| \\ \\int_a^b f(x) \\ dx - L_N(f) \\ \\right| \\leq \\frac{(b-a)^2}{2 N} K_1\n",
    "$$\n",
    "\n",
    "where $\\left| \\, f'(x) \\, \\right| \\leq K_1$ for all $x \\in [a,b]$.\n",
    "\n",
    "---\n",
    "\n",
    "**Theorem.** Let $R_N(f)$ denote the right Riemann sum\n",
    "\n",
    "$$\n",
    "R_N(f) = \\sum_{i=1}^N f(x_{i} ) \\Delta x\n",
    "$$\n",
    "\n",
    "where $\\Delta x = (b-a)/N$ and $x_i = a + i \\Delta x$. The error bound is\n",
    "\n",
    "$$\n",
    "E_N^{R}(f) = \\left| \\ \\int_a^b f(x) \\ dx - R_N(f) \\ \\right| \\leq \\frac{(b-a)^2}{2 N} K_1\n",
    "$$\n",
    "\n",
    "where $\\left| \\, f'(x) \\, \\right| \\leq K_1$ for all $x \\in [a,b]$.\n",
    "\n",
    "---\n",
    "\n",
    "**Theorem.** Let $M_N(f)$ denote the midpoint Riemann sum\n",
    "\n",
    "$$\n",
    "M_N(f) = \\sum_{i=1}^N f(x_i^* ) \\Delta x\n",
    "$$\n",
    "\n",
    "where $\\Delta x = (b-a)/N$ and $x_i^* = (x_{i-1} + x_i)/2$ for $x_i = a + i \\Delta x$. The error bound is \n",
    "\n",
    "\n",
    "$$\n",
    "E_N^{M}(f) = \\left| \\ \\int_a^b f(x) \\ dx - M_N(f) \\ \\right| \\leq \\frac{(b-a)^3}{24 N^2} K_2\n",
    "$$\n",
    "\n",
    "where $\\left| \\, f''(x) \\, \\right| \\leq K_2$ for all $x \\in [a,b]$.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "There are several points to notice:\n",
    "\n",
    "* Left and right Riemann sums have the same error bound which depends on the first derivative $f'(x)$.\n",
    "* Midpoint Riemann sum error bound depends on the second derivative $f''(x)$.\n",
    "* We expect the midpoint Riemann sum to give a better approximation as $N \\to \\infty$ since its error bound is inversely proportional to $N^2$ but left/right Riemann sum error bound is inversely proportional only to $N$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "## Implementation\n",
    "\n",
    "Let's write a function called `riemann_sum` which takes 5 input parameters `f`, `a`, `b`, `N` and `method` and returns the Riemann sum\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^N f(x_i^*) \\Delta x\n",
    "$$\n",
    "\n",
    "where $\\Delta x = (b-a)/N$ and $x_i = a + i\\Delta x$ defines a partition with $N$ subintervals of *equal length*, and `method` determines whether we use *left* endpoints, *right* endpoints or *midpoints* (with midpoints as the default method)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def riemann_sum(f,a,b,N,method='midpoint'):\n",
    "    '''Compute the Riemann sum of f(x) over the interval [a,b].\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    f : function\n",
    "        Vectorized function of one variable\n",
    "    a , b : numbers\n",
    "        Endpoints of the interval [a,b]\n",
    "    N : integer\n",
    "        Number of subintervals of equal length in the partition of [a,b]\n",
    "    method : string\n",
    "        Determines the kind of Riemann sum:\n",
    "        right : Riemann sum using right endpoints\n",
    "        left : Riemann sum using left endpoints\n",
    "        midpoint (default) : Riemann sum using midpoints\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        Approximation of the integral given by the Riemann sum.\n",
    "    '''\n",
    "    dx = (b - a)/N\n",
    "    x = np.linspace(a,b,N+1)\n",
    "    \n",
    "    if method == 'left':\n",
    "        x_left = x[:-1]\n",
    "        return np.sum(f(x_left)*dx)\n",
    "    elif method == 'right':\n",
    "        x_right = x[1:]\n",
    "        return np.sum(f(x_right)*dx)\n",
    "    elif method == 'midpoint':\n",
    "        x_mid = (x[:-1] + x[1:])/2\n",
    "        return np.sum(f(x_mid)*dx)\n",
    "    else:\n",
    "        raise ValueError(\"Method must be 'left', 'right' or 'midpoint'.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "Let's test our function with inputs where we know exactly what the output should be. For example, we know\n",
    "\n",
    "$$\n",
    "\\int_0^{\\pi/2} \\sin(x) \\, dx = 1\n",
    "$$\n",
    "\n",
    "and, since $\\sin(x)$ is concave down on $[0,\\pi/2]$, we know that left endpoints will give an over-estimate, and right endpoints will give an over-estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0000102809119054"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "riemann_sum(np.sin,0,np.pi/2,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0078334198735821"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "riemann_sum(np.sin,0,np.pi/2,100,'right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.992125456605633"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "riemann_sum(np.sin,0,np.pi/2,100,'left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "We also know that $\\int_0^1 x \\, dx = 1/2$ and midpoint should give the result exactly for any $N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "riemann_sum(lambda x : x,0,1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "## Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "### Approximate Pi\n",
    "\n",
    "Find a value $N$ which guarantees the right Riemann sum of $f(x)=\\frac{4}{1 + x^2}$ over $[0,1]$ is within $10^{-5}$ of the exact value\n",
    "\n",
    "$$\n",
    "\\int_0^1 \\frac{4}{1 + x^2} dx = \\pi\n",
    "$$\n",
    "\n",
    "Compute\n",
    "\n",
    "$$\n",
    "f'(x) = -\\frac{8x}{(1+x^2)^2}\n",
    "$$\n",
    "\n",
    "Use brute force optimization to find a bound on $\\left| f'(x) \\right|$ on $[0,1]$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5980759093919907"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linspace(0,1,1000)\n",
    "y = np.abs(-8*x/(1 + x**2)**2)\n",
    "np.max(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "Therefore, $\\left| f'(x) \\right| \\leq 2.6$ for $x \\in [0,1]$. Use the error bound\n",
    "\n",
    "$$\n",
    "\\frac{(b-a)^2}{2 N} K_1 \\leq 10^{-5} \\ \\Rightarrow \\ \\frac{1.3}{N} \\leq 10^{-5} \\ \\Rightarrow \\ 130000 \\leq N\n",
    "$$\n",
    "\n",
    "Let's compute the right Riemann sum for $N=130000$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.14158496127\n"
     ]
    }
   ],
   "source": [
    "approximation = riemann_sum(lambda x : 4/(1 + x**2),0,1,130000,method='right')\n",
    "print(approximation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "Verify the accuracy of the approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(approximation - np.pi) < 10**(-5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "### Approximate ln(2)\n",
    "\n",
    "Find a value $N$ which guarantees the midpoint Riemann sum of $f(x)=\\frac{1}{x}$ over $[1,2]$ is within $10^{-8}$ of the exact value\n",
    "\n",
    "$$\n",
    "\\int_1^2 \\frac{1}{x} dx = \\ln(2)\n",
    "$$\n",
    "\n",
    "Compute\n",
    "\n",
    "$$\n",
    "f''(x) = \\frac{2}{x^3}\n",
    "$$\n",
    "\n",
    "Since $f''(x)$ is decreasing for all $x>0$ we have $\\left| \\, f''(x) \\, \\right| \\leq 2$ for all $x \\in [1,2]$. Use the error bound:\n",
    "\n",
    "$$\n",
    "\\frac{(b-a)^3}{24 N^2} K_2 \\leq 10^{-8} \\ \\Rightarrow \\ \\frac{1}{12 N^2} \\leq 10^{-8} \\ \\Rightarrow \\frac{10^4}{\\sqrt{12}} \\leq N\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2886.7513459481288"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "10**4 / np.sqrt(12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "Therefore a partition of size $N=2887$ guarantees the desired accuracy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.693147176811\n"
     ]
    }
   ],
   "source": [
    "approximation = riemann_sum(lambda x : 1/x,1,2,2887,method='midpoint')\n",
    "print(approximation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "Verify the accuracy of the approximation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(approximation - np.log(2)) < 10**(-8)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}