08_quadrature.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Quadrature\n",
    "\n",
    "**Goal:** Evaluate integrals\n",
    "\n",
    "$$ \\int^b_a f(x) dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Many integrals do not have closed form solutions\n",
    "$$ \n",
    "   \\int^b_a \\sqrt{1 + \\cos^2 x} dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Solution to ordinary differential equations\n",
    "   \n",
    "   $$\\frac{\\text{d}^2 u}{\\text{d}t^2} = f\\left(u, \\frac{\\text{d} u}{\\text{d}t}, t \\right)$$\n",
    "   \n",
    "   Defining $v = \\frac{\\text{d} u}{\\text{d}t}$ then leads to\n",
    "\n",
    "   $$\\begin{bmatrix}\n",
    "   \\frac{\\text{d} v}{\\text{d}t} \\\\ \\frac{\\text{d} u}{\\text{d}t} \\end{bmatrix} = \\begin{bmatrix} f(u, v, t) \\\\ v \\end{bmatrix}$$\n",
    "   \n",
    "   which can be solved by integration\n",
    "   \n",
    "   $$\\begin{bmatrix}\n",
    "   v \\\\ u \\end{bmatrix} = \\begin{bmatrix} v(t_0) + \\int^t_{t_0} f(u, v, \\hat{t}) d\\hat{t} \\\\ u(t_0) + \\int^t_{t_0} v d\\hat{t} \\end{bmatrix}$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Solving partial differential equations\n",
    "$$\n",
    "    u_t = \\nabla^2 u\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Basics of Quadrature\n",
    "\n",
    "We want to approximate an integral $I$ with some approximation $I_N$ such that\n",
    "$$\n",
    "    I = \\int^b_a f(x) dx \\approx I_N = \\sum^{N}_{i=1} w_i f(x_i)\n",
    "$$\n",
    "where the $x_i$ are the *quadrature points* or *nodes* and the $w_i$ are the *weights*.  Usually a particular quadrature rule specifies the points $x_i$ resulting in a particular set of weights $w_i$.\n",
    "\n",
    "Convergence requires that\n",
    "$$\n",
    "    \\lim_{N \\rightarrow \\infty} I_N = I.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Riemann Sums\n",
    "\n",
    "Given $f(x)$ and a partition of the interval $[a,b]$ with $\\{x_i\\}^N_{i=0}$ and $a = x_0 < x_1 < \\ldots < x_N = b$ and $x^*_i \\in [x_i, x_{i+1}]$ we define the Riemann integral as\n",
    "\n",
    "$$\\int^b_a f(x) dx = \\lim_{N\\rightarrow \\infty} \\sum^{N-1}_{i=0} f(x_i^*) (x_{i+1} - x_i)$$\n",
    "\n",
    "This is a general definition and leads to a number of quadrature approaches based on how we pick $x_i^* \\in [x_i, x_{i+1}]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Midpoint Rule\n",
    "\n",
    "Choose $x_i^*$ such that\n",
    "\n",
    "$$x_i^* = \\frac{x_{i+1} + x_i}{2}$$\n",
    "\n",
    "so that\n",
    "\n",
    "$$I[f] = \\int^b_a f(x) dx \\approx \\sum^{N-1}_{i=0} f\\left(\\frac{x_{i+1} + x_i}{2} \\right ) (x_{i+1} - x_i) = I_N[f]$$\n",
    "\n",
    "over $\\Delta x_i = x_{i+1} - x_i$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Integrate using midpoint rule\n",
    "\n",
    "Calculate and illustrate the midpoint rule.  Note that we are computing the cummulative integral here:\n",
    "\n",
    "$$\n",
    "    \\int^x_0 sin(\\hat{x}) d\\hat{x} = \\left . -\\cos \\hat{x} \\right|^x_0 = 1 - \\cos x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Note that this calculates the cummulative integral from 0.0\n",
    "\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 10\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "x_star = 0.5 * (x_hat[1:] + x_hat[:-1])\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(hspace=.5)\n",
    "axes = fig.add_subplot(2, 1, 1)\n",
    "\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "axes.plot(x, f(x), 'b')\n",
    "\n",
    "for i in range(num_partitions):\n",
    "    axes.plot([x_hat[i], x_hat[i]], [0.0, f(x_star[i])], 'k--')\n",
    "    axes.plot([x_hat[i + 1], x_hat[i + 1]], [0.0, f(x_star[i])], 'k--')\n",
    "    axes.plot([x_hat[i], x_hat[i + 1]], [f(x_star[i]), f(x_star[i])], 'k--')\n",
    "    \n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Partition and $f(x)$\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "I_hat = numpy.zeros(x_star.shape)\n",
    "I_hat[0] = f(x_star[0]) * delta_x\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + f(x_star[i]) * delta_x\n",
    "    \n",
    "axes = fig.add_subplot(2, 1, 2)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_star + delta_x / 2.0, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Newton-Cotes Quadrature\n",
    "\n",
    "Using $N+1$ equally spaced points, evaluate $f(x)$ at these points and exactly integrate the interpolating polynomial:\n",
    "\n",
    "$$I_N[f] = \\int^b_a P_N(x) dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Trapezoidal Rule\n",
    "\n",
    "Use $N = 1$ polynomial to derive the trapezoidal rule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Trapezoidal rule uses $N = 1$ order polynomials between each point (i.e. piece-wise defined linear polynomials).  The coefficients of the polynomial in each interval are\n",
    "\n",
    "$$p_0 = f(x_i) \\quad \\quad p_1 = \\frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i}$$\n",
    "\n",
    "which gives the interpolating polynomial\n",
    "\n",
    "$$p_1(x) = \\frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i} ( x- x_i) + f(x_i)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Integrating this polynomial we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    I_N[f] &= \\int^{x_{i+1}}_{x_i} (p_0 + p_1 (x - x_i)) dx = \\left . p_0 x + p_1 \\left (\\frac{x^2}{2} - x_i x\\right) \\right |^{x_{i+1}}_{x_i} \\\\\n",
    "    &= p_0 \\Delta x + p_1 \\left (\\frac{1}{2} (x_{i+1} + x_i) \\Delta x - x_i \\Delta x\\right) \\\\\n",
    "    &= f(x_i) \\Delta x + (f(x_{i+1}) - f(x_i))\\left (\\frac{1}{2} (x_{i+1} + x_i) - x_i\\right) \\\\\n",
    "    &= f(x_i) \\Delta x + (f(x_{i+1}) - f(x_i)) \\frac{\\Delta x}{2} \\\\\n",
    "    & = \\frac{\\Delta x}{2} (f(x_{i+1}) + f(x_i))\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also simplify the sum over all the intervals by noting that all but the end points will have total contribution of $\\Delta x$ to the entire sum such that\n",
    "\n",
    "$$\n",
    "    I_N[f] = \\frac{\\Delta x}{2} (f(x_0) + f(x_N) ) + \\sum^{N-1}_{j=1} \\Delta x f(x_j)\n",
    "$$\n",
    "\n",
    "This is known as the composite trapezoidal rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEYCAYAAAC3LjroAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4FGXXwOHfSagh9N6SICBVihQb+ClYAEFEBBWUJsQuvuJroaivCtjFCsQGQiwIqIhgBVRUqqL0ThKkhia9JOf7Yya6kMIm2WRLzn1deyW7M/PMmc1mzs7TRlQVY4wxxlOYvwMwxhgTeCw5GGOMSceSgzHGmHQsORhjjEnHkoMxxph0LDkYY4xJx5KDMcaYdCw5GGOMSceSgylQRGSliFyW0+WBQEQmiMjTWSyvKCLfisg+EXlXREaLyP1elr1IRBr5LloTrCw5mIAiIltE5KiIHBKRne6JMDIXZV3h+ZqqNlLVed4sD2KPAutVtSzwMNAHGO/lti8AT+ZVYCZ4WHIwgaiLqkYC5wMtgeHZ2VhECuVJVMHjCuAT9/d+wCxVPerltjOAy0WkSl4EZoKHJQcTsFT1L2A20BhARB4RkY0iclBEVolIt7R13auAh0XkT+CwiHwIRAFfuFchD3msd4WITMpquft7AxGZJyL73eqma8/Y34Mi8qeIHBCRj0WkWGbH4kXsmZYlIs1F5Dd324+BDPcjIkVE5ABwnntcy4GOwA9nrPeciHzm8fx5EfleRIqo6jFgKXB1Vn8bE/osOZiAJSI1gU7A7+5LG4G2QGngf8BkEanqscnNwDVAGVW9GUjEvQpR1ec8y1bVW7NaLiKFgS+Ab4BKwL1AvIjU81itJ9ABqAU0wfmWnpmzxZ5hWSJSBPgMmASUw7ki6J7RDlT1BHARsMs9pvNwEsXaM1Z9FufqoLmI3OHu93p3e4DVQNMsjsUUAJYcTCD6TET2A/NxvvWOAlDVT1R1m6qmqurHwHqgtcd2r6pqUjaqULJyIRAJPKOqJ1R1DjATJwF57m+bqu7FSSTNMivMy9gzKutCoDAwRlVPqupUYHEWcTcD/vB4XgY4eEYse4CXgYk47ROdVPWAxyoH3e1MAWbJwQSi61S1jKpGq+pdaSd7EekjIsvcap79ONVNFTy2S/JhDNWAJFVN9XgtAaju8XyHx+9HcJJJhryIPbOyqgF/6elz6ydkEfeZyWEfUDKD9X7Huap4VFXPfN9KAvuz2IcpACw5mKAgItHAW8A9QHlVLQOsAMRjtTNvTnK2m5VktXwbUFNEPP9HooC/vIv4X17GnpntQHUR8Vw3Kov1m3J6cvgTOPeMeM4DxuJcOQzIoIwGZ5RhCiBLDiZYlMA5me8GEJH+uA3VWdgJnJPD5QtxvsE/JCKF3bEPXYCPshFzmpzEnuZX4BRwnxvH9ZxeHXWmM5PDLOD/0p6ISHWcaqs7gLuA8zzHdbgN4S2Ab72Mz4QoSw4mKKjqKuBFnJPlTpwqkZ/PstloYLhblfNgdpa7jbNdcHr7JANvAn1UdU0+xe4Zx/U4DdR7gRuB6Rmt63Y/LQt4xvg+0ElEiotIKZxk8ZKqzlDVI8DzwEiP9bsA81R1m9cHaEKS2G1CjQltIjIKpwfTGC/WXQjcpqor8j4yE8gsORhjjEnHqpWMMcakY8nBGGNMOpYcjDHGpBO0E5RVqFBBY2Ji/B2GMcYElaVLlyarasWzrRcwyUFE3gU64/SqOGsf8JiYGJYsWZL3gRljTAgRkaxG2P8jkKqVJuBMAGaMMcbPAubKQVV/FJEYf8cRilRh3z7Ytg327oX9+5X9+1M5diycnTs3s3nzElJTU6hUKYqqVWtStWpVypcvRMWKULEiVKoEhQv7+yiMMfkpYJKDN0QkFogFiIrKanqZgunAAVi2DH777TDr18OmTSVYu3YrSUlvkZKShDNDddrPT3Bq8VYDjwCbPEoKx5nvbTbQCFhCmTILqFatJuecE8X550fRunU5mjQRatQA8WaGIGNMUAmq5KCqcUAcQMuWLQv06L3UVPjpp11MmvQDP/74A9u2JXH4cNrJfw9FioylceM7qFdvH1u2PEWpUlWpXDmKatWaUaNGF3r2jKJFCzhx4nIOHJjH4cMH2bw5kS1bEklKSiIhIZGePStx4gR8/PFXzJs3gv37YdUqmDkToDiwhfLlK1Gr1kxq1tzFNde05uabGxERYdnCmGAXUCOk3Wqlmd40SLds2VILWoP09u0wY8ZRvv++ON9/D3v31se5h0whSpU6l6pVa1KnThTnnVeTHj06cP75zUlJSSElJYUiRYpkWu5zzzn3uXnooYcyXJ6amsquXbtISkoiMTGRtWuT+OOPRNq2fZ7ffgvno48u5ejRjcA2RBoQHd2Trl170K9fI5o2tSsLYwKJiCxV1ZZnXc+SQ2DbsAEmTEhm0qRPSUycAvxJ1apbufrqwlSo8C1jx15PoULh7N+f8+n3y5Rx7uuS0zLKlClDSkoqffo8y4wZU9i69QecSUivJzp6GtdeC926waWXQnh4jsM0xvhA0CUH956/l+HcAGUn8LiqvpPZ+qGcHLZvh/h4GD/+FzZs+B/wPZBC+fJ16NbtRl544SFKly4F5P7E7osyztx+x44dTJgwnQ0bSrFr1y18881Rjh+/nMjIa+jevSf33FOPFi3sisIYfwi65JBdoZYcTp6EyZP38tJLn7FyZStUz6NBgx/ZvXsAN97Yk4EDe9K0aVPkjDNqICaHM61YsYmePfuxevV8nCuKJlSq1JO77x7A3XdXpXz5HO3WGJMD3iaHQBrnUCAtX76PLl3eIzKyEwMGVGbFitu4+OJPWLMGVq5sy65d63n99VE0a9YsXWIIFo0bn8OqVT+ydWsSo0e/Qu3aJdm1aziPP76NatWga9eNTJu2wd9hGmM8WHLwk88+28Q55wyiSZN6zJw5gPDw1XTv/gALFy7hp5/+R716ICJBmxAyUr16dR555D42bJjP1q1bWbbsfG6/HWbPfpYbbqhLxYr9GDr0W06eDM6rWWNCiqoG5aNFixYabFJTVWfNUo2JeVqhkEIRPffcG/XTTxdpampqjsqcNWuWzpo1K1dx5baM3G6/alWiXnhhfw0Pr6SAFit2mT7wwM96+HCOizTGZAJYol6cY63NIR+kpsKkScmMGRPBsmURlC//CfXr/8C77w7j3HOr+ju8gHH06HEGD45jwoSnOXlyF0WKdOY//3md4cOjiYz0d3TGhAZrcwgAqvDRR39TteoT9Ot3DomJr/HOO7BtWw/mz3+dSZPeZMSIEbnax4gRI/xehq9iGDXqaeLi7mXv3k3Exo4mLGw5zz4bSa1aMHr0SQ4fztUujDHZ4c3lRSA+Ar1aaebMI1qz5vMK5RTQVq266x9/rDxtndKlS2vp0qVztZ9AKCOvYjh16pT++qvqVVelKrTWYsUG6JNPbtHjx3O1K2MKNLysVrIrBx/780+4+mro3LkvSUn/pXHj1ixYsIRFi6bSpElDf4cXVMLDw7nwQvjss2P06HERx49P5rHHzqVSpcGMG7eTIK0RNSYoWHLwkaSkFC699H2aNt3G4sXwwAOP8u23P7B8+WwuuKCFv8MLasWLF2fKlDFs3ryeq6/uw4EDb3Dnnedw3nkLWbDA39EZE5osOeTS0aNK797TiY4+j59+6ssll7zHxo3w4ovNueKKS/0dXkiJjo7iq6/eYvXqVbRvP4A9e5pz0UXQocPvrFlzyN/hGRNSLDnkUGqq8sQTX1OmTCs++KA7JUoob7wxlZ9+GkrZsv6OLrTVr38u3333GuvXF+HRR0/y9dddaNCgNp07v8KBA8f8HZ4xISGopuwOFIsW7eDaa+9i587dFC68h//+dwKjR99CeDZnlfv6669zHUsglOGvGCIjYdSowrRqNZXbbx/Kl1+OonLlSYwa9SoPPHBxrmMypiCzcQ7ZcPy40rv3ZKZNGwwc4ZprnuTjj++nRInMp8M2+eepp2bz5JN3cOpUEnXrDmb27JHUrh3h77CMCSg2zsHHfvkF6tdfxLRpfShXrj7XX9+D2rX/ylViGDx4MIMHD85VXIFQRiDEAJCc/BX9+3eides7Wb9+DHXrNuH+++eRmpqrYo0pkOzK4SwOHFAGDVrG1KnNqVEDBg78imHDrqS8O5WoP2dDDZQyAiGGM8uYPHkesbG3cfToZVx88Tu8/TY0aJDjoo0JGXbl4AMTJiRQufLVfPJJa265ZS0rV8Jjj3XIdtuCyX+33HIZu3b9ybhxL7NmDTRpsowBA77l5El/R2ZMcLDkkIE9e1K58MKx9O/fmBMnfuHBB19hwoS6lCzp78hMdkRGluD220uxahVUqTKS9967iipVBjF//gF/h2ZMwLPkcIYvvkilWrWrWbjwLmrVuojVq1fw/PN3ERZmb1WwqlwZ1q17n27dHmbv3ndp27YRt9zypV1FGJMFO+O5DhxQbrsNrr02jLJl2zNixNts3Pg19erF+Ds04wPFixdn+vRn+OabBZQuXZb4+M7UqzedlSv9HZkxgcnGOQATJ64hNvY2Tp58jEcfvZrHH3+EokWz3mbhwoW53m+olBEIMXhbxpVXtmLnziXcddd4ZszoTIsWMHToboYPr4hdHBrzrwLdW+nQoVN07Pgi8+c/TlhYBE8//TaPPnq9jyI0gW7XLujf/wCzZjWiYsVLmD37dVq0qOjvsIzJU9Zb6SymTVtBpUoXMX/+I9Sq1Yn161dlKzEMHDiQgQMH5iqGUCkjEGLISRmVKsGnn0Zw/fV3snv3p7Rq1ZC77/6I1NTg/MJkjC8VuCuHlBSle/dXmDHjC+BPHn74DUaN6pHtezUHYt9+f5URCDHktoxvvllJ9+79OXQIoqLq8OOPY4mOLp3jWIwJVHblkIHly3dTpUoXPv/8P6j+QIkSxxk9ume2E4MJPVdd1Yjk5F8ID/+TxMQPqV27OePGLfJ3WMb4TYFJDqNG/UTTpk1JTv6Wm256ldKlIwkPLzCHb7xQtGghIiOLUbx4CURSuPPOS2jf/jmOH7f5N0zBE/Jnx6NH4e67YdiwIhQtWp5PP13Ihx/e6++wTAArUqQQmzYtIzq6K3PmPEzTprPYtMnfURmTv0I6OXzzTQK1a7/Om2/CkCEXsHfvH1x3XTN/h2WCQM2aZdm8+ROGDfuG7duvoXlzGDdup7/DMibfhGSDtCrcccc04uIGIpLK5Mlr6dWrymnrJCcnA1ChQoUc7T+324dSGYEQQ16WkZAA3bqt5/ffz6dx4zv44YeRlCtn07Sb4ORtg3S2k4OIlACOqWpKToPzhcySw7ZtR2jT5j9s3hxHqVKt+eabD7nggnP8EKEJJQcPHuXyy4ewdOlYihVrxSeffEjnzrX9HZYx2eaz3koiEiYivUTkSxHZBawBtovIKhF5XkTq+CJgX/jll1Rq1bqczZvjuOyyh9m1a36mieGmm27ipptuyvG+crt9KJURCDHkdRklSxZnyZI3eeKJaRw/vp4uXZrTv/8HBOmFtzFnddYrBxH5AfgO+BxYoaqp7uvlgMuBXsCnqjo5j2M9jeeVQ0qK8uKLMGyYUKbMh4wYUYH77rsyy+2tb7/vygiEGPKzjN9+S+SKK3qxb19Trr/+Dd5+G7tvuAka3l45eDO30hWqmm7+SlXdC0wDpolI4RzE6BPr1++jbdtB7NzZme7d+/H22zfj/n8bkyfOPz+KHTvm8dJLKYwYAQ0b/sYzzwh9+zb3d2jG+MxZq5XSEoOIvCKZjBbLKHlkl4h0EJG1IrJBRB7xZpvXXptP/fpN2bnzc3r0+JtPPsESg8kXRYoU4pFHijJ/Puzb9x/69buQrl1fISXF6plMaMhOV9aDwAy3QRoRuVpEfvZFECISDrwBdAQaAjeLSMOstlm3bjv33fd/hIcXZvLkX5gy5T5soLPJbxdcACtXTqNq1auZMeN+qlbtwurVu/0dljEZOnzY+3W9nrJbVYeLSC9gnoicAA4BXn3D90JrYIOqbgIQkY+ArsCqzDY4eHAbxYpdTrNm4bz11n956y3o2bMnd911F0eOHKFTp07ptunXrx/9+vUjOTmZQ4cOAXDZZZf9s/zOO+/kxhtvJCkpiVtvvTXd9kOGDKFLly6sXbs2w+2HDx/OFVdcwbJly7j//vvTbT9q1CguvvhifvnlF4YOHZqujDFjxtCsWTO+++47nn766XTbjx8/nnr16vHFF1/w4osvApxWxqRJk6hZsyYff/wxY8eOTbf91KlTqVChAhMmTGDChAn/vJ5WxpEjR4iIiODNN99kypQp6bafN28eAC+88AIzZ848bXvPi8qnnnqK77///rRty5cvz7Rp0wB49NFH+fXXX09bnrZvgPvvv59ly5adtvzcc88lLi4OgNjYWNatW3fa8mbN/h2/csstt7B169bTll900UWMHj0agO7du7Nnz57Tlrdv354RI0YAcPjw4dP+rgCdO3fmwQcfBEi3DGDo0B78+OOVfPLJgzRq1ItzzilJjRp7/1nu+dm74YYb0m2fnc/e7bffnm55dj97Z8rJZ89TTj97aWbNmpWjzx449+qYPXs2kLPPXo0aNZg82Wkyzelnb8yYMUDuP3sdO3bk6NGjpy0/22fP2/PeuHFzue++X9Mtz4zXyUFE2gODgMNAVWCAqq71ek9Zqw4keTzfClyQQQyxQCxAeHg5WrdOQSRnUxu0adMmR9v5avtQKqNNmzYUL148V2X07Nnzn3/QnEprRL7llltyVUZG/6BnExYmTJlyLxMntqR//9vYuHENycnX0KTJwWzfJ2L+/PlA7v4ugVCGL2PIjRkzZlCmTBmfdFTo2bNnrrZv3759rmI4dOhQtt/P1FSle/eXmT79YaDKWdf/h6p69QDmAG3c388DlgHtvN3+LGXfALzt8fxW4PWstmnRooUaE4i2bz+otWv3V0BLl75EFy3akq3tS5curaVLl85VDIFQRiDEEChl+CuGtWt3auXKnRTQypW76tq1yQosUS/Oy15/p1HVdqo63/19OU77QPrrz5z5C6jp8byG+5oxQadKlUjWr3+X22//gAMH/uSCC5rz8cfWDmHy14IF0LbtUnbunEO3bq+zbdunnHtuea+392YQXGY9lLYD7bNaJxsWA3VFpJaIFAFuAmbkskxj/EYExo27mdmzf6dKlSe56aaKPPIINsOryXPHj5/kjjt+oG1bKFGiI7NmbWb69LsJC8veadqbK4c5InKviER5vuiexC8SkYlA32zt9Qyqegq4B/gaWA1MUVW79bsJeh061GbDhnuIjYVnn/2FcuVa8PXXK/wdlglRixZtpnLlSxk/vj1XXrmJ33+Hjh2z0c7gwZsG6fVACvCpiFQF9gPFgHDgG2CMqv6eo717UNVZwKzclmNMoImIgPHjoUqV4zz11DY6dGjFoEEvM3787RneaKpDhw653mcglBEIMQRKGfkRw/DhUxg1ahCqMHDgB8TFnZOr7v3eTJ/xu6o2F5HfcLqcVgSOqmrOm/59IKe3CTXGnxYt2slVV/XhwIFvqFWrOz/99BbVq9vcGybnTpxQLrzwLn7/fRzFi1/ItGkf0LFjrUzX9+VtQr8XkV+BykAfoBqQvX5+xhgAWreuzM6ds7n88ufYvPlzWrSYyJ9/nr5OcnLyP1OH51QglBEIMQRKGXkVw6ZN0Lat8PvvNWna9BG2b/8xy8SQHV5N2S0itYG5wEScbqyNgBM4E/Hd6JNIssmuHEywe+utFYwY0YD9+8N56KF1PPZYbQoVCg+qSQgDPYZAKcPXMagq/fuPZcqUcyhSpANvvQU9enhXji8n3kNVN4rIFar6z9BAEYkEGnsXjjHmTIMGNaZrV7j11n089VQb3nijEXPm5OvkxiYIbdmyl0svHUhS0qdUrNiHxYs7EB3t+/1kZ5zDujOeH1LVBb4PyZiCo1IlmD27DL17P8fevYto3rwpx47leh5LE6JOnDhFnTpNSUqayVVXvcjWre/lSWKAEL+HtDHBICxMmDy5H59//htFi9bk+PEjHDx4lOTkbMySZkLakSMnOHToOEePHkakGO+88ytff/0ARYrk3SnckoMxAeLaa+uxY8cCWrd+nNTUVKpVO4+4uLn+Dsv42fTpS6lQoSUpKcdo2PAhtmz5jQEDWuT5fr2eeM8Yk/dKly7KwoVP8OKL7Xj44du4/fZ2vPPOHXz11XOULVvS63Iymvk1u3JbRiDEEChl5GR7Vbj++mf47LPhiFTmoYe+4NlnO+cqjuzwqrdSILLeSibUJSUd4YorRrBu3csUK1aXpUuX07BhEX+HZfLBX3/BbbfB119/SLVq3/Htty/SsKFv7mTmy3EOxhg/qFkzgjVrXuSJJ35B5GFatSrCa6/BwYNnb4tYu3Yta9fmbkb93JYRCDEEShnebn/48BE6dRpC3bqv8uOP8NprN5OU9I7PEkN22JWDMUFg61YYNAi++mo6RYrcyxtvjGfgwMyrGEKxb38wl+HN9tOn/0CfPrdx+PBGqlZ9gB9+eJG6dXO0uyzZlYMxIaRGDZg1C0aMiOHUqfIMGtSF1q1vZffuPWff2AS0gwcPceWV99C9+2UcPqzExs4hKSlvEkN2WHIwJkiIwJNPns+6dUuoXftxFi/+iOrVGzF27Bf+Ds3kUFISdOiwjO++G0fVqoP57bc/GT/+csLD/R2ZJQdjgk7t2kVYv/4JnnpqCamp1bj33r8ZORJOnPB3ZMZbe/ceYMCAT2jYEJYta8Pw4etJShpD8+Yl/B3aPyw5GBOERGD48KYkJCyiW7deDB8OMTHv8MQTHxKs7YgFxZtvzqZq1ca8915vWrTYyooV8NRTtQLiasGTjXMwJohVr16ITz6BL75QevT4gP/9bw4TJnxM5843UL587r6F9u2bq3t45Xr7UCqjb9++HDx4jEaN+rFq1UTCwxvyv/9NY8SIGrm650Jest5KxoSIv/9OoWvXMcybNxyI4MYbn+S99wZRvLiNjfCnU6dSGDx4MuPHP0FKShLNmz/KzJnDqVatqF/isd5KxhQwpUqFM3fuEGbO/JMqVXrw8cf3UKZMfR577H1SUlKyXd7ChQtZuHBhjuPJ7fahUMa4cd9QunQT3nyzH6VK/R8TJizmt9+e8ltiyBZVDcpHixYt1BiTsZSUVH344dlaqND5CmipUg117Nipmpqa6nUZpUuX1tKlS+c4htxuH6xlpKam6o4dp/Suu1RF3tewsHo6cOAUPXEiJVcx+AqwRL04x9qVgzEhKCxMeOaZDuzdu4Tu3ady8GAqd931P0aMUP7+29/Rha45c36mdu3LiY5+mXHj4M47e7F9+wreeqsHhQsH1+k2uKI1xmRLyZLC1Knd2bhxBV26fMnIkWHUqrWf8867ljlz5vs7vJCxZMkymjXrTPv2bdi8eQ3165dn+XJ4441wKlUKzn4/lhyMKQBq1Qrn889rsngx1Kq1lhUrFtG+fVuaNOnE4sW/+zu8oJWaCjfe+BStWjXnjz9+pnr10cyatZFly/rTsKG/o8sdSw7GFCAtW8LixRfw+ecbqVbtGZYvX0Dr1ufTunVPDh8+7u/wgsbmzQm8++5uWrSAKVMuo0KF4bz//maSkh6hY8fAGciWG8F5vWOMyTERuPbaEnTp8jAffHA7Q4a8xOLFa2natCgPPgjduh2gcuXS3HvvvbnaT263D8Qy/vprJwMGjOTbb8ejegd1677C+++3pVevtgE3iC23bJyDMQWcKnz2mTJ6tLB48SZEmnDBBf15991hNGhQxd/hBYTExH3cdtvzfP/9K6gep2zZATz55AjuvLNm0CUFG+dgjPGKCHTrJixcCFOnRlC9em8WLBhLw4ZtqVbtMl555TOOHTuW7XJnz57N7NmzcxWbv8uYM2czjRsPICbmUr777hkqVLiON95YTXJyHPfcE3yJITvsysEYk86XX27g8cc/YOnSV4E9hIeXpHXra7n//h7ccEMXwsLO/r0yEO6jkN0yTp48yddfz+HPPyOZN+8Svv12K3Ae9erFMmJEb3r3bpLjOAKFt1cOlhyMMZnaseMkw4fP45NPpvD339MJCyvDoEEb6N1bKFnyTxo0qEfRohmP9g2W5HDq1CnmzJnL2LFT+Oqr6Rw7the4nho1phEbC337niQqqnCO9x9oLDkYY3xGFebOPcmrrybwzTd1OHr0JGFhVShU6BTt2nXlzjt7cvXVV56WKAI5OTjnPeGPP+C669qRkDAXiCQsrCsXXdSTIUOuokuXYhQKwS471uZgjPEZEWjXrjCffVaHXbtg0qQwWrSI5+TJ7nz11Rd07dqFUqUqc//9k8jFeTxPpaSkMHv2XDp1upMyZWpTq9YRmjeHxMTBNGnyKa++uovk5MnMn38t3bqFZmLIDr9fOYhID+AJoAHQWlW9uhywKwdj/G/fPpgx4wRvv/0dCxZM4dSp2wkLu4gGDRaxenUbChUqxHffLaJp0yhKlSqV7fJze+Xw118HqV27CidPniI8vCwnT+4EIggP78Lll79Ez57VuO46qFgxR8UHJW+vHAIhN64ArgfG+zsQY0z2lC0LffsWoW/fTpw40YkFC+C77+DDD9eRmlqcEyfqc+ml5wFQqFApypWLolGjFrz55gRiYuDXX+cSFhZGVFQU1atXp0iR06cXHzp0aJb7P3LkCJ999hlbtiSxYkUi69cnsnVrIlFRg9m3bwDr1ycB0cAWSpe+lDZtejJoUCeuvDKCYsXy5j0JFX6/ckgjIvOAB+3KwZjQsHfvMaZPX8ncuRv5889EkpISOXAgCSgDvAdAoULNOHXqD3cLITKyCvXrd6J377eJiIDFi99n376d/P33LvbsSSQ5OYl9+xKpU6cPtWqNIiHhb5YuLe1uXxaIAmpSocIA2rTpRtOmxylbdhk33tiYKlVCY+RybgVdg7Q3yUFEYoFYgKioqBYJCQn5FJ0xxhcOH4ZVq2DdOti8Gf74YyMbN25h165EDhxI5PDhJFRrAcPcLSoAp4BjQE0gCpGaREZeQ40aPaheHSIj11C3bg0aN46kXj2oVw/c2iiTgYBKDiLyHZDRUMthqvq5u8487MrBmAJNFY4dg6NHnceuXXuWQeP9AAAgAElEQVQoVEgoX74shQoJJUtCsWIE7K01g0FAtTmo6hX5sR9jTHATgeLFnQdA9erl/RtQAWZdWY0xxqTj9zYHEekGvAZUBPYDy1T1ai+2OwiszePw8lsFINnfQeSBUDwuO6bgEIrHBLk7rmhVPWvnXb8nh5wSkSXe1JsFk1A8JgjN47JjCg6heEyQP8dl1UrGGGPSseRgjDEmnWBODnH+DiAPhOIxQWgelx1TcAjFY4J8OK6gbXMwxhiTd4L5ysEYY0weseRgjDEmnaBLDiLSQUTWisgGEXnE3/H4goi8KyK7RGSFv2PxFRGpKSJzRWSViKwUkcH+jim3RKSYiCwSkT/cY/qfv2PyFREJF5HfRWSmv2PxFRHZIiLLRWSZiITEXDsiUkZEporIGhFZLSIX5dm+gqnNQUTCgXXAlcBWYDFws6qu8mtguSQilwKHgPdVtbG/4/EFEakKVFXV30SkJLAUuC6Y/1YiIkAJVT0kIoWB+cBgVV3g59ByTUQeAFoCpVS1s7/j8QUR2QK0VNWQGQQnIhOBn1T1bREpAkSoap7cXinYrhxaAxtUdZOqngA+Arr6OaZcU9Ufgb3+jsOXVHW7qv7m/n4QWA1U929UuaOOQ+7Twu4jeL5dZUJEagDXAG/7OxaTOREpDVwKvAOgqifyKjFA8CWH6kCSx/OtBPkJpyAQkRigObDQv5Hknlv9sgzYBXyrqkF/TMAY4CEg1d+B+JgC34jIUne6/2BXC9gNvOdWAb4tInl2k4pgSw4myIhIJDANuF9V//Z3PLmlqimq2gyoAbQWkaCuBhSRzsAuVV3q71jyQBtVPR/oCNztVt8Gs0LA+cBYVW0OHAbyrN012JLDXzh3/EhTw33NBCC3Xn4aEK+q0/0djy+5l/NzgQ7+jiWXLgGudevnPwLaichk/4bkG6r6l/tzF/ApTrV0MNsKbPW4Wp2KkyzyRLAlh8VAXRGp5TbG3ATM8HNMJgNu4+07wGpVfcnf8fiCiFQUkTLu78VxOkas8W9UuaOqj6pqDVWNwfl/mqOqt/g5rFwTkRJuRwjcqpercO5XH7RUdQeQJCL13JfaA3nWwSNfbvbjK6p6SkTuAb4GwoF3VXWln8PKNRH5ELgMqCAiW4HHVfUd/0aVa5cAtwLL3Tp6gKGqOsuPMeVWVWCi22suDJiiqiHT9TPEVAY+db6jUAj4QFW/8m9IPnEvEO9+Od4E9M+rHQVVV1ZjjDH5I9iqlYwxxuQDSw7GGGPSseRgjDEmHUsOxhhj0rHkYIwxJh1LDsYYY9Kx5GCMMSYdSw7G+IiItBKRP937PpRw7/kQ1HMvmYLLBsEZ40Mi8jRQDCiOMw/OaD+HZEyOWHIwxofcaQ0WA8eAi1U1xc8hGZMjVq1kjG+VByKBkjhXEMYEJbtyMMaHRGQGztTXtXBuk3qPn0MyJkfsysHkGRF5Ir/uDeDeTP6K/NhXFjH0AU6q6gfAM0ArEWmXzTJmi0jfPAkwl0RknogM9HccuSUiMSKiIhJUs1LnN0sOQSY7J8FQ+Wf2FTdZqYhckBflq+r7qtrd/T1FVS9Q1TnZLKOjqk70dWx5fULMzhcBEeknIvPzIg7jO5YcTI4E27cu9+ZDfYC97s+82k9QvS+hyL3fhsklSw5BLO0bmIi8ICL7RGSziHR0l40E2gKvi8ghEXndfb2+iHwrIntFZK2I9PQor7yIfCEif4vIYhF52vMbnvvN824RWQ+sd197RUSS3G2WikhbL2MvKyIzRWS3G/tMEanhsXyeiDwlIj+LyEER+UZEKngsv1VEEkRkj4gM82KXbXFu1nMfcJPbq8jzffxZRF4XkQMiskZE2p8Ry2gRWeQe5+ciUs5dlvaN/DYRSQTmuK9f645z2O9u38B9vbb73p/vPq/mvgeXeexr4BlxveyWs0lELnZfTxKRXZ5VUCJyjTg3nv/bXf6Ex/H/6P7c734eLnK3GSAiq92/wdciEu1R3pXue3HA/fyIF+9z2rYqIneIyHo39jfE0QAYB1zkxrHfXb+o+zlOFJGdIjJOnLvtpZX3kIhsF5FtIjLQLb+Ou2yCiIwVkVkichi4/CzvhfGGqtojiB7AFuAK9/d+wElgEM6d8e4EtvFvR4N5wECPbUsASTh3jyoENAeSgYbu8o/cRwTQ0F13vsf2CnwLlAOKu6/dgtNDpxAwBNgBFHOXPQFMzuQ4ygPd3X2VBD4BPvNYPg/YCJyLM2ZgHvCMu6whcAi4FCgKvAScSntfMtnfO8AUoDCwB+jusayfu/1/3OU3AgeAch6x/AU0dt/DaWnHBcS478v77rLibsyHcW4jWhh4CNgAFHG3GYRze8cInLsavnDGcQ88I67+7t/3aSAReMM97quAg0Cku/5lwHk4X/qaADuB686Is5DHvrq6cTVw/37DgV/cZRXcsm9wj+E/biwDM3l/T/tbu/uaCZQBooDdQAeP45p/xvYv49zyt5z7efgCGO0u64DzuWrkvmeT3fLruMsnuH+vS9xjL5bd98IeGfxN/R2APbL5B0ufHDZ4LItwP/RV3Of/nGjc5zcCP51R3njgcffkcxKo57HsadInh3ZniW8f0NT9/bQTxlm2awbs83g+Dxju8fwu4Cv398eAjzyWlQBOkElycN+Xvz1ODuOBzz2W98MjqbqvLQJu9YjlGY9lDd39hXucaM7xWD4C5xaiac/DcJLLZR6vzQCWA38CRc84bs/ksN5j2Xnuvip7vLYHaJbJcY8BXnZ/T3dCBGYDt50R5xEgGqfqbYHHMsG5wX12kkMbj+dTgEc8jmv+GWUfBmp7vHYRsNn9/V3cROE+r0P65PD+WT5fWb4X9kj/sGql4Lcj7RdVPeL+GpnJutHABe5l/n73kr43UAWoiPPtMclj/aQMyjjtNRF50K2WOOCWVxrnW2eWRCRCRMa7VUN/41R7lJHT64t3ePx+xOO4qnnGoaqHcU6SmemG86037f7V8UBHEanosc5f6p45XAnuftIknbGsMKcfp+fyau46afGlusure6zzFs6VyGuqejyL2Hd6/H7ULe/M1yIBROQCEZnrVlMdAO4g679FNPCKx2dhL86Jujrp32Ml489DVjL7+52pIk4CX+oRy1fu65wZSyZxnPm5zO57Yc5gySG0nTmIJQn4QVXLeDwiVfVOnMv+U0ANj/VrZlWmOO0LDwE9gbKqWgbn8t6buukhQD3gAlUthVNFhJfbbveMTUQicKqpMtMX58SUKCI7cKqwCgO9PNapLiKe+47CuZpIU/OMZSdxquTSeL7X23BOvGnxibv9X+7zSJxvsu8AT6S1X/jABzhXJDVVtTRO3X7aMWU0oCkJuP2Mz0NxVf2F9O+xkPHnISfOjCUZJ8k18oijtKqmJZPtZONz6crqvTBesOQQ2nYC53g8nwmcK05jbmH30UpEGqgzzcN0nJNVhIjU5+y9ekriJJTdQCEReQwo5WVsJXFOCPvdk+Pj2TiuqUBnEWkjTsPyk2TyWRaR6kB7oDNO1VUzoCnwLKcfXyXgPvc96YFTDz/LY/ktItLQTURPAlM186kxpgDXiEh7ESmMkwiPA7+4y18BlqjqQOBLnBOXL5QE9qrqMRFpzenJbzeQyumfh3HAoyLSCEBESrvHjhtXIxG5XpweWPfhXGH6wk6ghvu3S7uyegt4WUQqubFUF5Gr3fWnAP1FpIH7/o/wYh9ZvRfGC5YcQtsrwA1uT5RXVfUgTiPmTTjfbnfgnCSLuuvfg1MttAOYBHyIc1LLzNc4l//rcKpRjuF91cMYnMbbZGCBW45XVHUlcDfOt8PtOO0cWzNZ/VZgmap+o6o70h7Aq0AT+XfW1IVAXTeekcANqupZVTUJp257B06D531ZxLcWp6H+Nbe8LkAXVT0hIl1xGljvdFd/ADhfRHp7e/xZuAt4UkQO4rTLTPGI6Yh7XD+7VTcXquqnOH//j9yqvRVAR3f9ZKAHzmC+PTjvzc8+iBGcHl0rgR0iknb19TBO4/gCN5bvcK4sUdXZOH+vuWnruNtk9dnM9L0w3rHpM0ymRORZnMbtgByx6ysi0g+nobVNJsvn4TS2vp2fcZmMud1hV+A05J/ydzyhyq4czD/EGQPRxO2P3hq4DfjU33EZIyLd3LEQZXGudr6wxJC3LDkYTyVx2h0OAx8DLwKf+zUiYxy3A7twxr6k8G+1nMkjVq1kjDEmHbtyMMYYk07QThJWoUIFjYmJ8XcYxhgTVJYuXZqsqhXPtl7QJoeYmBiWLFni7zCMMSaoiEjC2deyaiVjjDEZsORgjDEmHUsOxhhj0rHkYIwxJh1LDsYYY9Kx5GCMMSYdSw7GGGPSseRgjDEmHUsOxhhj0rHkYIwxJh1LDsYYY9IJiOQgIjVFZK6IrBKRlSIy2N8xGWNMQRYQyQHnJvVDVLUhcCFwt4g09HNMxmRLfHw8MTExhIWFERMTQ3x8vL9DMibHAmJWVlXdjnOjeFT1oIisBqoDq/wamDFZSU2Fgwfh8GHiP/qI2GHDOHLsGAAJCQnEDhwIu3fT++abITISIiJAxM9BG+OdgLsTnIjEAD8CjVX178zWa9mypdqU3caX4uPjGTZsGImJiURFRTFy6FB6N2oEa9dCQgIkJjo/d+6E3bthzx4nQQAxQEbzIEcDW9KeFC0KFStChQpQrRpER0NUFNSqBfXrw7nnEj99+ukxjBxJ79698+PwTQEhIktVteVZ1wuk5CAikcAPwEhVnZ7B8lggFiAqKqpFQoJX05IbkzVV4l9+mdhHH+XIiRP/vBwBxAG9wfnGX60a1Kzp/Ew7yZcrByVKEHbHHWT0nyRA6htvwKFDkJzsPHbtgr/+cpLN3r3/rBsPxIpwxON/MqJ4ceLeessShPGZoEsOIlIYmAl8raovnW19u3IwOZaaCn/8AXPnwvz58PPPxOzalfE3/4oV2fLrr843/MKFMy0yJiaGjL6sREdHs2XLlsxjOXQINm2CNWuIiY0l4cCB9GUUKsSWvn2hbVto3x5q1PDiII3JWFAlBxERYCKwV1Xv92YbSw4mW5KT4csv4auv4PvvnWohgNq14ZJLCHv//Yy/+YuQ6lYdZSU+Pp7Y2FiOHDnyz2sRERHExcV5/a0/LCyMjP4fBUgtV+7fq4z69eHKK+Gaa+Cyy5zqKmO85G1yCJTeSpcAtwLtRGSZ++jk76BM8Miwp9CWLfDCC8437sqVoV8/mDcPrr4aJk50qnY2bICJE4mKjs6w3KioKK/237t3b+Li4oiOjkZEiI6OzlZiyGpfUdHRTjL74w948UWIiYG334YOHZyqrR49ID4e/v478/fCmOxS1aB8tGjRQo1RVZ08ebJGREQo8M8jIixMJ4MqqDZtqvrYY6pLlqimpnpfRkSETp482b/HkVkMR4+qzpypGhurWrWqc5xFi+rkli01okgRvx6HCWzAEvXiHOv3k3xOH5YcTJromjVPOxmmPaLLlFHduNHrciZPnqzR0dEqIhodHe2XE2qOYkhJUZ0/X/XeezU6LCzj9yI6Os9jN8HB2+QQEG0OOWFtDoZVq2D8eMJefTVX7QWhJMt2iz17nN5VpkALtjYHY7xz6hRMmQKXXgqNGsHYsURFRGS4qrftBaEk03YLcLrg9ukD9qXKeMGSgwkOe/fCM884A8ZuvNFpTH7uOfjrL0bGxRFxRoKIiIhg5MiRfgrWf0aOHJnxezF6NAwcCJ99Bq1awSWXOEn21Ck/RWoCnjd1T4H4sDaHAmLzZtX77lONiHCayNq1U/38c9VTp05bLRDaCwJFlu/FgQOqY8ao1q7tvJ8xMaqvvqp66JD/Ajb5CmtzMEFt5UoYPRo++sgZndyrFwwZAk2a+Duy0JCSAjNnwvPPw88/O20R993nPMqW9Xd0Jg9Zm4MJGqf1y69WjfjWraFxY6cKZPBgZwTxxImWGHwpPBy6dnVGiM+fD23awBNPOPM9DR1K/LhxNlaioPPm8iIQH1atFBoy7NsPOvm661STk/0dXsHyxx+qPXvqZPdvgI2VCElYtZIJBjHVqpGwfXu61886J5HJMzHVq5OwbVu61+1vEhqsWskEtsREGDCAxAwSg7M4MZ8DMmky/ZskJMDRo/kcjfEXSw4mf+3Z4zQsn3sufPABUSVLZrhaQRyjECiyHCtRty68847ToG1CmiUHkz+OHXMmwatTB8aMcXofrVvHyLFjbYxCgMl0rMSwYc79LAYOhGbNYPZsZ/YqE5q8aZgIxIc1SAeJ1FTVjz5y+tODaseOqsuXn7aKjVEIPJn+TVJTVadO/XecxBVXOA3ZJmhgDdLG75Ysgfvvd/rRN23q9Km/8kp/R2V84cQJGDsWnnwS9u+HQYOc3ytV8ndk5iz82iAtIiVEJDwvyjZBYMcO594JrVrB+vXOvQeWLrXEEEqKFHHGoKxfD/fe67RD1K3r3G/C41arJnj5JDmISJiI9BKRL0VkF7AG2C4iq0TkeRGp44v9mMBz2gC26Gjie/VyGps//BAeftg5edx2mzPoyoSecuWcNqTly52BdA8+CE2bEv/IIzaILth5U/d0tgfwAzACaAKEebxeDugOTANu8cW+0h7W5uB/mQ5ga9ZMdf16f4dn/GHmTJ1cqZINogtg5Gebg4gUVtWTuV0nO6zNwf9iYmJISEhI97oNlirYYqKjSchgnIp9LgJDvrY5pJ30ReQVEZGs1jEh4vhxZ1BUBmwAW8GWmJSU8ev2uQgqvm6QPgjMEJESACJytYj87ON9GH+bMweaNiWzYWo2gK1gy3QQnSpcf70zOt4EPJ8mB1UdDnwIzHOTwgPAI77ch/GjHTugd29o3x5OnmTkf/9rA9hMOpkOouvZE776Cho0cG7UdNIqEwKaNw0T3j6A9sBcYB6wFqjny/I9H9YgnY9OnVJ9803V0qVVixRRHTFC9cgRVbUBbCZjmX4utmxR7drVGUDXuLHq/Pn+DbQAwh+D4ERkDvCYqs4XkfOAScADqjrHZztxWYN0Plm2DG6/HRYtcq4Y3nzT6apqTG7MmAH33ANJSc50HM8+63SLNXnOL4PgVLWdqs53f18OdASe9uU+TD45dMjps96yJWzZApMnw7ffWmIwvnHttbBqFfz3v/Dee1C/vvMZC9IZG0KRrwbBZdZDaTtOVVOm65gA9MUX0LChM9r1tttgzRqnrcH+hMaXIiOdtoelS+Gcc+DWW51R9OvX+zsyg++uHOaIyL0iclo3BREpAlwkIhOBvlkVICLvisguEVnho5iMF04b4VyjBvGtWjnf6kqVcm4fOX683VPY5K2mTZ35t958ExYvhvPOg6efJn7iRBtl7U/eNEyc7QHEAXcBS4FtwCpgE5AAvAU096KMS4HzgRXe7NMapHMv0xHOPXuqHj/u7/BMQbRt27+3KhWxUdZ5gHweIf27qjYXkd+A1kBF4Kiq7s9mOTHATFVtfLZ1rUE692yEswlUMZUqkbB7d7rX7bOZe/ndIP29iPwKVAb6ANUAn99PUERiRWSJiCzZncEHx2TD4cM2wtkErMTk5Ixft89mvvHV9BkPArcAKUAtnEn4VojIShH52Bf7cPcTp6otVbVlxYoVfVVswTN7NjRubCOcTcDKdJR10aKwaVM+R1Mw+awrq6puBK5Q1RGqep2q1gUuAF721T5MLm3fDjfeCJ06QbFijBw+3EY4m4CU4SjrwoUZqQqNGsEzz9gI67zmTcNEfj2AGKxB2vdSUv4d4Vy0qOqTT6oeO6aqNsLZBK4MP5tJSardutkI61wg2G4TKiIfApcBFYCdwOOq+k5m61uDtJeWLYM77oCFC6FdO+fWjjaQzQQ7G2GdY369TWhOqOrNqlpVVQurao2sEoPxwsGDMGSIM8J50yaYNAm++84SgwkNaSOshwz5d4T1++/bCGsfCpjkYHxEFaZPd0Y4v/QSDBjgjHC+5RYb4WxCS2QkvPCCM8K6dm3o29e5Ol692t+RhQRLDqFk40bo3Bm6d3cusX/5BeLi7HLbhLa0Edbjx8MffzjPhw6Fw4f9HVlQs+QQxE6b+qJMGeLr1YMff3TmRFq6FC66yN8hGpM/wsIgNta5Sr75Zhg9Gho2JP4//yEmOtqm4MiBgGmQzq6C3iAdHx9PbGwsR44c+ee1iPBw4saMofc99/gxMmMCwE8/Ed+rF7Fbt3LE4+WIiAji4uLo3bu330LzN28bpC05BKmY6tVJ2LYt3es2vYAxjpjoaBIyGFFd0P9Hgq63kvHSoUMwdCiJGSQGsOkFjEmTmJSU8euJidaryQuWHIJFaqrTVe/cc2H0aKJKlMhwNZv6whhHplNwqELbtk67nMmUJYdg8MsvTuNy375Qowb88gsjx4+3qS+MyUKGU3BERDDytttg3Tpo1crp6r19u58iDHDeDKMOxEeBmD5j0ybVnj2dqQKqVlV97z1nKgyXTX1hTNYy/R/Zv191yBDVwoVVS5RwppQ5fNi/weYTgm36jOwK6QbpvXth1Ch47TUID4eHHnLu5xwZ6e/IjAktGzfCww/DtGlQvTo89RT06eP834Uoa5AORseOwfPPO6M9X3oJevVyLn+feMISgzF5oXZtmDrVGR9UvbpTzdSsGXz5ZYFvtLbk4CenDWCLjiZ+4ECoW9e5Srj4Ymek53vvOW0Mxpi81bYtLFgAU6Y4X9I6d4b/+z9nvITn/2pBGkjnTd1TID6Cuc0h03s316mjOneuv8MzpmA7ccKZ4r5qVede1mFhIXUva6zNIXBlOjgnKootmdy60xiTz44ccQab7t+fblEwD6Tzts2hUH4EY1ynTsFHH2U6UC2zQTvGGD+IiCDxwIEMFxWEwabW5pAfjh379yY7t95KVOHCGa5mA9iMCSxZDqRr1w6+/TZkG64tOeSl3bvhySchOhruugsqV4YZMxj57rs2gM2YIJDhQLrixRnZq5czA+xVV0GLFjB5cujd09qbholAfAR0g/Ty5aqxsarFijkD2Dp1Up0zRzU19Z9VbACbMcEh0//VY8dU4+JU69d3/s+rV1cdPVp1927/BnwWeNkg7feTfE4f/kwOGX5YTpxQnTJF9dJLnbe1aFHVQYNUV63yW5zGmHyQkqL65Zeq7dr9+7/fp4/qwoWqGnhfBC055JEMu6EWKqSTS5Vy3s5atVSfe041Odkv8Rlj/Gj5ctU771SNjFQFnRwVpRGFCwdUV1hvk4N1Zc2mmJgYEjLobhpdvDhbPvkEOnQI6aH3xhgv/P03TJ5MzP33k5BBW4Q/u8La9Bm+tn8/TJpEYibjEBKPHYNrrrHEYIyBUqXgrrtIPHUqw8WJCQnw3XcB3YhtySEriYkwbpxzNVCxIvTpQ1QmJ3/rhmqMOVOmXWFF4MoroUoV6N8fPv3UudoIIAUyOWQ6V8rBgzB7Nvz3v9C4sdMF9c47Yf16+M9/4NdfGTlhgnVDNcZ4JdN7SrzzDkyfDp06OYnh+uuhfHln7MSzz8KiRc6gWZdf5nfypmEiPx5AB2AtsAF45Gzr57RBOtMG5dq1VcPDnUblwoVV27dXfeEFp7eRRxfUtDICqfeBMSZwnfV8ceKE6rx5qg89pNq4sXMOAtWSJVU7dtTJ3btrRNGiuW7UTovDOe0HSYO0iIQD64Arga3AYuBmVV2V2TZeN0inpjp3etq4EdatI+aBB0g4eDDdatFFi7LlwQfh8sudu66dke2NMSZf7NgBP/wA8+bBDz8Qs3o1GbV0RleowJbZs6FOHShTJssi4+PjiY2N5ciRIwCoqpwtjEBJDhcBT6jq1e7zRwFUdXRm27SsU0eXjB7tNOgcPw4HDjiP/fth507Yts1JCklJcPToP9uF4aTeDGIgNTXVtwdmjDG5FBYWRkbnaQH+OWOVK+dM71+1qvOoWBFKl3YepUoRM2QICcnJ/2zrTXIIlIn3qgOes85tBS44cyURiQViAVoA9OyZvqTISGeaiqpVnZt2dO7s3NCjTh2oXZuodu0ynBHVGpSNMYEoKioqw+7zUVWrwhtvOLUiGzf++4V4xQpITna+NLtyMk1goCQHr6hqHBAH0LJhQ2XKFChcGIoWdbqOlSoFhbI+pJGjRp12eQXWoGyMCVwjR47M+Jz1/PPQrVvmG6bVqBw8SFTbtiRs356t/QZKb6W/gJoez2u4r2UuIsLpUVSvHsTEOJdVZ0kMAL179yYuLo7o6GhEhOjoaOLi4ujdu3du4jfGmDyR43NW0aJQqRLUrs3I559P12vqbAKlzaEQToN0e5yksBjopaorM9smmG/2Y4wx+S0+Pp5hw4aRkJAQPA3SACLSCRgDhAPvqmqW9TwichCn62soqQAkn3Wt4BOKx2XHFBxC8Zggd8cVraoVz7ZSwCSH7BKRJerF/CDBJBSPCULzuOyYgkMoHhPkz3EFSpuDMcaYAGLJwRhjTDrBnBzi/B1AHgjFY4LQPC47puAQiscE+XBcQdvmYIwxJu8E85WDMcaYPGLJwRhjTDpBlxxEpIOIrBWRDSLyiL/j8QUReVdEdonICn/H4isiUlNE5orIKhFZKSKD/R1TbolIMRFZJCJ/uMf0P3/H5CsiEi4iv4vITH/H4isiskVElovIMhEJiRGzIlJGRKaKyBoRWe1OWpo3+wqmNoecTO0dDETkUuAQ8L6qNvZ3PL4gIlWBqqr6m4iUBJYC1wXz30pEBCihqodEpDAwHxisqgv8HFquicgDQEuglKp29nc8viAiW4CWqhoyg+BEZCLwk6q+LSJFgAhV3Z8X+wq2K4fWwAZV3aSqJ4CPgK5+jinXVPVHYK+/4/AlVd2uqr+5vx8EVuPMvhu03HumHHKfFnYfwfPtKhMiUgO4Bnjb37GYzIlIaeBS4B0AVT2RV4kBgi85ZDS1d1CfcAoCEYkBmgML/RtJ7rnVL8uAXcC3qhr0x4Qzbc1DeNweIEQo8I2ILHWn+w92tYDdwHtuFeDbIlIir3YWbMnBBBkRiQpJwUcAAAJNSURBVASmAferamDdQT0HVDVFVZvhzBzcWkSCuhpQRDoDu1R1qb9jyQNtVPV8oCNwt1t9G8wKAecDY1W1OXAYyLN212BLDtmf2tv4jVsvPw2IV9Xp/o7Hl9zL+bk49z4PZpcA17r18x8B7URksn9D8g1V/cv9uQv4FKdaOphtBbZ6XK1OxUkWeSLYksNioK6I1HIbY24CZvg5JpMBt/H2HWC1qr7k73h8QUQqikgZ9/fiOB0j1vg3qtxR1UdVtYaqxuD8P81R1Vv8HFauiUgJtyMEbtXLVUBQ9wZU1R1AkojUc19qD+RZB49guxPcKRG5B/iaf6f2zvSeD8FCRD4ELgMqiMhW4HFVfce/UeXaJcCtwHK3jh5gqKrO8mNMuVUVmOj2mgsDpqhqyHT9DDGVgU+d7ygUAj5Q1a/8G5JP3AvEu1+ONwH982pHQdWV1RhjTP4ItmolY4wx+cCSgzHGmHQsORhjjEnHkoMxxph0LDkYY4xJx5KDMcaYdCw5GGOMSceSgzE+IvL/7d2xTcRQEATQ2QwERFRAESQXUggVUAbSdUIBFEB4yWVUcxJCS2BnP8KyZCzeq2Cz0X7Ls/VYVZ/z3Yeb+ebDrruX+L/8BAcrqqrXJFdJrjP14Bw3HgkWEQ6wornW4JzkkuTQ3d8bjwSLeFaCdd0nuU1yl2mDgF2yOcCKquo9U/X1Q6YzqS8bjwSL7KqVFf6yqnpO8tXdb3Nz66mqnrr7Y+vZ4LdsDgAMfHMAYCAcABgIBwAGwgGAgXAAYCAcABgIBwAGP9B1ceWpaAQVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Note that this calculates the cummulative integral from 0.0\n",
    "\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 20\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(hspace=.5)\n",
    "axes = fig.add_subplot(2, 1, 1)\n",
    "\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "axes.plot(x, f(x), 'b')\n",
    "\n",
    "for i in range(num_partitions):\n",
    "    axes.plot([x_hat[i], x_hat[i]], [0.0, f(x_hat[i])], 'k--')\n",
    "    axes.plot([x_hat[i + 1], x_hat[i + 1]], [0.0, f(x_hat[i+1])], 'k--')\n",
    "    axes.plot([x_hat[i], x_hat[i + 1]], [f(x_hat[i]), f(x_hat[i+1])], 'k--')\n",
    "    \n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Partition and $f(x)$\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "I_hat = numpy.zeros(x_hat.shape)\n",
    "I_hat[0] = (f(x_hat[1]) + f(x_hat[0])) * delta_x / 2.0\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + (f(x_hat[i + 1]) + f(x_hat[i])) * delta_x / 2.0\n",
    "    \n",
    "axes = fig.add_subplot(2, 1, 2)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_hat + delta_x, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Simpson's Rule\n",
    "\n",
    "Simpson's rule uses $N = 2$ order polynomials between each point (i.e. piece-wise defined quadratic polynomials).  \n",
    "\n",
    "The polynomial has the form\n",
    "\n",
    "$$P_2(x) = \\frac{2 f(x_i)}{\\Delta x^2} \\left (x - \\frac{\\Delta x}{2} \\right ) (x - \\Delta x) - \\frac{4 f\\left(x_i + \\frac{\\Delta x}{2}\\right)}{\\Delta x^2}  x (x - \\Delta x) + \\frac{2 f(x_{i+1})}{\\Delta x^2} x \\left (x - \\frac{\\Delta x}{2} \\right )$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Integrating this polynomial we have\n",
    "\n",
    "$$\n",
    "    I_N[f] = \\int^{x_{i+1}}_{x_i} P_2(x) dx = \\frac{\\Delta x}{6} f(x_i) + \\frac{2 \\Delta x}{3} f\\left(x_i + \\frac{\\Delta x}{2} \\right ) + \\frac{\\Delta x}{6} f(x_{i+1})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also show this by using the method of undetermined coefficients.  \n",
    "\n",
    "Use the general form of the quadrature rule and determine weights $w_j$ by using functions we know the solution to.  These functions can be any representation of polynomials up to order $N=2$ however the monomials $1$, $x$, $x^2$ are the easiest in this case.\n",
    "\n",
    "$$\n",
    "    I_{\\Delta x}[f] = w_0 f(0) + w_1 f(\\Delta x / 2) + w_2 f(\\Delta x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{aligned}\n",
    "    &\\text{if}~f = 1:  &I[f] =  \\int^{\\Delta x}_{0} 1 dx = \\Delta x  \\\\\n",
    "    &\\text{if}~f = x:  &I[f] =  \\int^{\\Delta x}_{0} x dx = \\frac{\\Delta x^2}{2} \\\\\n",
    "    &\\text{if}~f = x^2:  &I[f] =  \\int^{\\Delta x}_{0} x^2 dx = \\frac{\\Delta x^3}{3}\n",
    "\\end{aligned}$$\n",
    "\n",
    "What are the corresponding systems of equations?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    &\\text{if}~f = 1:  &I[f] =  \\int^{\\Delta x}_{0} 1 dx = \\Delta x & & I_N[1] &= w_0 + w_1 + w_2 \\\\\n",
    "    &\\text{if}~f = x:  &I[f] =  \\int^{\\Delta x}_{0} x dx = \\frac{\\Delta x^2}{2} & & I_N[x] &= w_1 \\frac{\\Delta x}{2} + w_2\\Delta x\\\\\n",
    "    &\\text{if}~f = x^2:  &I[f] =  \\int^{\\Delta x}_{0} x^2 dx = \\frac{\\Delta x^3}{3} & & I_N[x^2] &= \\frac{\\Delta x^2}{4} w_1 + w_2\\Delta x^2\\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "We then have the system of equations:\n",
    "$$\\begin{aligned}\n",
    "    w_0 &+& w_1 &+& w_2 &=\\Delta x \\\\\n",
    "        &\\quad& \\frac{\\Delta x}{2} w_1 &+& \\Delta x w_2  &= \\frac{\\Delta x^2}{2} \\\\\n",
    "        &\\quad& \\frac{\\Delta x^2}{4} w_1 &+& \\Delta x^2 w_2 &=\\frac{\\Delta x^3}{6} \\\\\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "    1 & 1 & 1 \\\\\n",
    "    0 & \\Delta x / 2 & \\Delta x \\\\\n",
    "    0 & \\Delta x^2 / 4 & \\Delta x^2 \\\\\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "    w_0 \\\\ w_1 \\\\ w_2\n",
    "\\end{bmatrix} = \\begin{bmatrix} \n",
    "    \\Delta x \\\\ \\Delta x^2 / 2 \\\\ \\Delta x^3 / 3\n",
    "\\end{bmatrix} \\Rightarrow \\begin{bmatrix}\n",
    "    1 & 1 & 1 \\\\\n",
    "    0 & 1 / 2 & 1 \\\\\n",
    "    0 & 1 / 4 & 1 \\\\\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "    w_0 \\\\ w_1 \\\\ w_2\n",
    "\\end{bmatrix} = \\begin{bmatrix} \n",
    "    \\Delta x \\\\ \\Delta x / 2 \\\\ \\Delta x / 3\n",
    "\\end{bmatrix} \\Rightarrow \\begin{bmatrix}\n",
    "    1 & 1 & 1 \\\\\n",
    "    0 & 1 / 2 & 1 \\\\\n",
    "    0 & 0 & -1 \\\\\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "    w_0 \\\\ w_1 \\\\ w_2\n",
    "\\end{bmatrix} = \\begin{bmatrix} \n",
    "    \\Delta x \\\\ \\Delta x / 2 \\\\ -\\Delta x / 6\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "Leading to \n",
    "\n",
    "$$ w_2 = \\frac{\\Delta x}{6} \\quad w_1 = \\frac{2}{3} \\Delta x \\quad w_0 = \\frac{\\Delta x}{6}$$\n",
    "\n",
    "Another way to write Simpson's rule is to use intervals of three points (similar to one of the ways we did this last time).  The formulation here effectively has a $\\Delta x$ half of what the intervals show but is easier to program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Note that this calculates the cummulative integral from 0.0\n",
    "\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 10\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.subplots_adjust(hspace=.5)\n",
    "axes = fig.add_subplot(2, 1, 1)\n",
    "\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "axes.plot(x, f(x), 'b')\n",
    "\n",
    "for i in range(num_partitions):\n",
    "    axes.plot([x_hat[i], x_hat[i]], [0.0, f(x_hat[i])], 'k--')\n",
    "    axes.plot([x_hat[i + 1], x_hat[i + 1]], [0.0, f(x_hat[i + 1])], 'k--')\n",
    "    coeff = numpy.polyfit((x_hat[i], x_hat[i] + delta_x / 2.0, x_hat[i + 1]), \n",
    "                          (f(x_hat[i]), f(x_hat[i] + delta_x / 2.0), f(x_hat[i+1])), 2)\n",
    "    x_star = numpy.linspace(x_hat[i], x_hat[i+1], 10)\n",
    "    axes.plot(x_star, numpy.polyval(coeff, x_star), 'k--')\n",
    "    \n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Partition and $f(x)$\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "I_hat = numpy.zeros(x_hat.shape)\n",
    "I_hat[0] = delta_x * (1.0 / 6.0 * (f(x_hat[0]) + f(x_hat[1])) + 2.0 / 3.0 * f(x_hat[0] + delta_x / 2.0))\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + delta_x * (1.0 / 6.0 * (f(x_hat[i]) + f(x_hat[i+1])) + 2.0 / 3.0 * f(x_hat[i] + delta_x / 2.0))\n",
    "    \n",
    "axes = fig.add_subplot(2, 1, 2)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_hat + delta_x, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Error Analysis\n",
    "\n",
    "From before we have a particular quadrature scheme $I_N$ which we can also write as\n",
    "$$\n",
    "    I_N[f] = \\sum^{N-1}_{i=0} w_i f(x_i).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Define the error $E[f]$ such that\n",
    "\n",
    "$$I[f] = I_N[f] + E_N[f]$$\n",
    "\n",
    "The degree of $I_N[f]$ is the integer $N$ such that $E_N[p_i] = 0 \\quad \\forall i \\leq n$ and $\\exists p_{n+1}$ such that $E[p_{n+1}] \\neq 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Quadrature Accuracy\n",
    "\n",
    "We can also use our polynomial analysis from before to analyze the errors made using both of the aforementioned methods.  From Lagrange's theorem we have the remainder term as before which we can use to look at the error\n",
    "\n",
    "$$R_N(x) = (x - x_0)(x - x_1) \\cdots (x- x_N) \\frac{f^{(N+1)}(c)}{(N+1)!}$$\n",
    "\n",
    "and integrate it to find the form and magnitude of the error on a single interval.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To find the total error we must sum the error over all the intervals:\n",
    "\n",
    "$$I[f] = \\sum_{i=0}^N \\int^{x_{i+1}}_{x_i} P_N(x) dx + \\sum_{i=0}^N \\int^{x_{i+1}}_{x_i} R_N(x) dx = I_N[f] + E_N[f]$$\n",
    "\n",
    "as we defined before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Trapezoidal error\n",
    "\n",
    "With $N=1$ we have\n",
    "\n",
    "$$R_1(x) = (x - x_i) (x - x_{i+1}) \\frac{f''(c)}{2}$$\n",
    "\n",
    "Integrating this leads to\n",
    "\n",
    "$$\\int^{x_{i+1}}_{x_i} (x - x_i) (x - x_{i+1}) \\frac{f''(c)}{2} dx = \\frac{\\Delta x^3}{12} f''(c)$$\n",
    "\n",
    "giving us a form for the error.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we sum up across all the intervals the total error is\n",
    "\n",
    "$$E_N[f] = -\\frac{\\Delta x^3}{12} \\sum_{i=0}^{N} f''(c_i)$$\n",
    "\n",
    "or more illustrative\n",
    "\n",
    "$$E_N[f] = -\\frac{1}{2} \\Delta x^2 (b - a) \\left [ \\frac{1}{N} \\sum^{N-1}_{i=0} f''(c_i) \\right ]$$\n",
    "\n",
    "where the expression in the brackets is the mean value of the second derivative over the interval $[a,b]$.  This also shows that the trapezoidal rule converges quadratically as $\\Delta x \\rightarrow 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Simpson's Rule Error\n",
    "\n",
    "Similarly here we have $N = 2$ and \n",
    "\n",
    "$$R_2(x) = (x - x_i) \\left(x - x_i - \\frac{\\Delta x}{2} \\right) (x - x_{i+1}) \\frac{f'''(c)}{3!}$$\n",
    "\n",
    "Integrating and summing the error contributions we find\n",
    "\n",
    "$$E_N[f] = -\\frac{1}{180} (b - a) \\Delta x^4 f^{(4)}(c)$$\n",
    "\n",
    "Interestingly we have gained two orders of accuracy by increasing the polynomial order by only 1!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Example 1:\n",
    "\n",
    "If $f(x) = \\sin \\pi x$ look at the relative accuracy of midpoint, trapezoidal and simpson's rules for a single interval $x\\in[0,1]$.\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\text{Exact:}  &I[f] = \\int^1_0 \\sin \\pi x = \\left . \\frac{-\\cos \\pi x}{\\pi} \\right |^1_0 = \\frac{2}{\\pi} \\approx 0.636619772 \\\\\n",
    "    \\text{Midpoint:}  &I_1[f] = \\Delta x f(1/2) = \\sin (\\pi / 2) = 1 \\\\\n",
    "    \\text{Trapezoid:}  &I_1[f] = \\frac{\\Delta x}{2} (\\sin(0) + \\sin(\\pi)) = 0 \\\\\n",
    "    \\text{Simpson's:}  &I_1[f] = \\frac{\\Delta x}{6} \\sin(0) + \\frac{2 \\Delta x}{3} \\sin(\\pi / 2) + \\frac{\\Delta x}{6} \\sin(\\pi) = \\frac{2 \\Delta x}{3} = \\frac{2}{3}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compute the error as a function of delta_x for each method\n",
    "f = lambda x: numpy.sin(numpy.pi * x)\n",
    "\n",
    "num_partitions = range(50, 1000, 50)\n",
    "delta_x = numpy.empty(len(num_partitions))\n",
    "error_mid = numpy.empty(len(num_partitions))\n",
    "error_trap = numpy.empty(len(num_partitions))\n",
    "error_simpson = numpy.empty(len(num_partitions))\n",
    "\n",
    "for (j, N) in enumerate(num_partitions):\n",
    "    x_hat = numpy.linspace(0.0, 1.0, N + 1)\n",
    "    delta_x[j] = x_hat[1] - x_hat[0]\n",
    "\n",
    "    # Compute Midpoint\n",
    "    x_star = 0.5 * (x_hat[1:] + x_hat[:-1])\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N):\n",
    "        I_hat += f(x_star[i]) * delta_x[j]\n",
    "    error_mid[j] = numpy.abs(I_hat - 2.0 / numpy.pi)\n",
    "    \n",
    "    # Compute trapezoid\n",
    "    I_hat = 0.0\n",
    "    for i in range(1, N):\n",
    "        I_hat += (f(x_hat[i + 1]) + f(x_hat[i])) * delta_x[j] / 2.0\n",
    "    error_trap[j] = numpy.abs(I_hat - 2.0 / numpy.pi)\n",
    "    \n",
    "    # Compute simpson's    \n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N):\n",
    "        I_hat += delta_x[j] * (1.0 / 6.0 * (f(x_hat[i]) + f(x_hat[i+1])) + 2.0 / 3.0 * f(x_hat[i] + delta_x[j] / 2.0))\n",
    "    error_simpson[j] = numpy.abs(I_hat - 2.0 / numpy.pi)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, error_mid, 'ro', label=\"Midpoint\")\n",
    "axes.loglog(delta_x, error_trap, 'bo', label=\"Trapezoid\")\n",
    "axes.loglog(delta_x, error_simpson, 'go', label=\"Simpson's\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error_trap[0], 2.0) * delta_x**2.0, 'b--', label=\"2nd Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error_simpson[0], 4.0) * delta_x**4.0, 'g--', label=\"4th Order\")\n",
    "axes.legend(loc=4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Recursive Improvement of Accuracy\n",
    "\n",
    "Say we ran the trapezoidal rule with step size $2 \\Delta x$, we then will have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\t\t\t\\int^{x_2}_{x_0} f(x) dx  &= \\frac{2 \\Delta x}{2} (f_0 + f_2) =  h (f_0 + f_2) \\Rightarrow \\\\\n",
    "\t\t\t\\int^b_a f(x)dx &\\approx I_{2\\Delta x}[f] = \\sum^{N/2-1}_{j=0} \\Delta x (f_{2j} + f_{2j+2}) \\\\\n",
    "\t\t\t&= \\Delta x (f_{0} + f_{2})  + \\Delta x (f_{2} + f_{4})  + \\cdots + \\Delta x (f_{N-2} + f_{N}) \\\\\n",
    "\t\t\t&= \\Delta x\\left ( f_0 + f_N +  2 \\sum^{N/2-1}_{j=1} f_{2j} \\right )\n",
    "  \\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now compare the two rules for $\\Delta x$ and $2 \\Delta x$:\n",
    "\n",
    "$$I_{\\Delta x}[f] = \\frac{\\Delta x}{2} \\left (f_0 + f_N + 2 \\sum^{N-1}_{j=1} f_j \\right)~~~~~~~~~ I_{2 \\Delta x}[f] = \\Delta x \\left ( f_0 + f_N +  2 \\sum^{N/2-1}_{j=1} f_{2j} \\right )$$\n",
    "\n",
    "$$I_{\\Delta x}[f] = \\frac{1}{2} I_{2\\Delta x} + \\Delta x(f_1 + f_3 + \\cdots + f_{N-1})$$\n",
    "\n",
    "Here we see we can actually reuse the work we did to calculate $Q_{2 \\Delta x}[f]$ to refine the integral."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Arbitrary Intervals (Affine Transforms)\n",
    "\n",
    "Mapping $\\xi \\in [-1,1] \\rightarrow x \\in [a,b]$ can be done through an *affine transform* or *affine map* which is a linear transformation.\n",
    "\n",
    "$$x = \\underbrace{\\frac{b - a}{2}}_{\\text{scaling}} \\xi + \\underbrace{\\frac{a+b}{2}}_{\\text{translation}} ~~~~~ \\text{or} ~~~~~ \\xi = \\left( x - \\frac{a + b}{2}\\right) \\frac{2}{b-a}$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    I[f] &= \\int^b_a f(x) dx = \\int^1_{-1} f(x(\\xi)) \\frac{dx}{d\\xi} d\\xi = \\frac{b - a}{2} \\int^1_{-1} f(x(\\xi)) d\\xi\\\\\n",
    "    I_N[f] &= \\sum_i w_i f(x(\\xi_i)) \\left . \\frac{dx}{d\\xi}\\right|_{\\xi_i}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Newton-Cotes Rules\n",
    "\n",
    "We can rewrite our previous quadrature rules so that they are given on the interval $\\xi \\in [-1, 1]$ instead of $x \\in [x_i, x_{i+1}]$.  Recall that a general quadrature rule can be written as\n",
    "$$\n",
    "    \\sum^N_{i=1} w_i f(\\xi_i)\n",
    "$$\n",
    "where $w_i$ are the weights and $\\xi_i$ are the points specified to evaluate the function at.  For Newton-Cotes rules we know that the points $\\xi_i$ are uniformly distributed on $[-1, 1]$ but we still need to define the weights.  For trapezoid rule we can do this to find that\n",
    "$$\n",
    "    \\int^1_{-1} f(x) dx \\approx f(-1) + f(1)\n",
    "$$\n",
    "so that $\\xi_0 = -1$, $\\xi_1 = 1$, and $w_0 = w_1 = 1$.  Note that if we map this using our affine transform we would get back the original trapezoid rule:\n",
    "$$\n",
    "    I_N[f] = \\sum_i w_i f(x(\\xi_i)) \\left . \\frac{dx}{d\\xi}\\right|_{\\xi_i} = (f(-1) + f(1)) \\frac{\\Delta x}{2}\n",
    "$$\n",
    "Similarly for Simpson's rule we have\n",
    "$$\n",
    "    \\xi = [-1, 0, 1] \\quad \\text{and} \\quad w = \\left[\\frac{1}{3}, \\frac{4}{3}, \\frac{1}{3} \\right].\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtoAAAEWCAYAAABYLDBhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XecVNX5x/HPs4DAilgQG7IzdgU1ItjFYNTECraocTRR0YkVsaEGCxrX3sW2aiw/RgSjWDG2iBE7WBLsojsLWABF2oICe35/nLsyu+xsY2fulO/79drX7t5pz9w597nP3HvOueacQ0RERERE2lZJ2AGIiIiIiBQiFdoiIiIiIhmgQltEREREJANUaIuIiIiIZIAKbRERERGRDFChLSIiIiKSASq0M8TMRpjZqCy9VqWZ7Z2N18okM3vOzP4SdhwNMbMJZnZi2HGsLDOLmpkzs/ZhxxIGMxtgZh81cvsoMxuRode+wsweyMRzy8pRvm455evMy2S+NrMyM1tgZu3a+rkLSVu0pbwqtFuSoAplQ2srwY7EmdlOYceSjnNuP+fcg239vJkuLluykzaz48xsYibiyLRgx7og+FliZr+k/H9X2PE1h3NugnOud2sea2Ynmtmy4P3OM7P3zWy/to6xUChft57ytfJ1WzCzDc3sMTObbWZzzWyKmR0H4Jyrcs51cc4tCznMZgnaRWUrHvdAyr7qRzN70cy2zECIaeVVoZ0r8u1ooJkZ8Gfgx+B3pl4nr9ZLIcrk0Ylgx9rFOdcFSADX1v7vnDu5gVgKsT28Frz/NYB7gbFmtlrIMUkj8q0dKl8XjywcTf4/YBoQAboBxwLfZ/g1c9G1Qd7uAcwA7svmi+dtoV37TdPMrjezOWb2de3RJTMrB/oDI4NvMSOD5VsG32Z+NLPPzOyIlOfrZmZPB0eq3g1O805Mud2Z2Wlm9gXwRbDsFjObFjxmspn1b2bsa5rZM2Y2K4j9GTPbMOX2CWb2dzN73czmm9kLZrZ2yu3HmlnSzH4ws+HNeMn+wPrAEOAoM1ul3np83cxGBt94PzWzverFcpWZvRO8zyfNbK3gttojD4PNrAr4d7B8oJl9ZGY/BY/fKli+SbDutw/+3yBYBwNSXuvEenHdFDzPV2a2a7B8mpnNtJTTlmZ2QHCEcV5w+4iU9/+f4PdPQXvYJXjMCWb2SfAZPG9mkZTn2ydYF3OD9mPNWM+1j3VmdrKZfRHEfrt5WwF3AbsEcfwU3L9j0I6rzOx7M7vLzDqnPN8wM/vWzL4xf1TVmdmmwW0PmNmdZjbezBYCezaxLjLGzPY2fxTzb2b2HXBPsF2NT2nrT5tZj5THTDSzcjObFKzrcWa2Zsrtu5nZW8F6/MDM9giW97flR9MXmNliM/syuK2Tmd0arLMZZnZjbZuvjTHl+fsGzzvfzEYDHZvzXp1zNfidWBeg9rOo89zBsum17buB9dXgeytEpnydNOVr5evs5+sdgAeccwudc0udc+87554LYqlz5iD4PK8wszeC9/u0+e0sYcu3s2i99TYk+Kxnm9l1ZlYS3Lapmb0afB6zzWxMyuN2DZ5rbvB715TbGt2WUpnZ+ebz+3zz+WGvhu6Xyjm3CBgLbJfyPHXOcNRfLw28btp22NgL580PUAnsHfx9HLAEOAloB5wCfANYcPsE4MSUx66K/2Z3PNAe6APMBnoFtz8S/JQCvYL7Tkx5vANeBNYCOgfLjsF/S2wPnAN8B3QKbhsBjErzProBhwWvtRrwKPBEyu0TgKnA5kDn4P+rg9t6AQuAPfBFwY3A0tr1kub17sM3rg7AD8BhKbcdFzz+rOD2I4G5wFopscwAtg7W4WO17wuIBuvloeC2zkHMC4F9gucbBnwJrBI85iTg4+C9Pw9cX+99n1gvruODz/cKoAq4PXjfvwfmA12C+w8AtsF/edwW/6394Hpxtk95rUFBXFsFn99FwBvBbWsHz3148B7OCmI5Mc36rfNZB6/1DP6oZxkwC9g35X1NrPf4m4Cn8G1rNeBp4Krgtn3x7ap3sM5GBc+/aXD7A8HntVvw3ju1dF20clt8ALii3rK9g/V0JbBK0B66A4cEf3cFHgf+mfKYifhtrRe+DT2B3zEA9MS31z8E72Vf/Dbbrd7rrhI8z9+D/68E3gheex3gbeDSlBgrg787AtPxBU0H4Ch8ThmR5j2fCEwI/m4PnAn8DKxd/7lTHjMdGBD8fUVL31s+/6B8DcrXytch5mvgJeB1fG4rq3dbnecOPs8vgU2A1YPP/XN8XmuPbzf311tvrwTroSy4b217GA0MT3mPuwfL1wLm4I+stwf+FPzfraltqV7sW+C3+Q1S3ssmadbBAwT7Kny7/z/gw0baQ0PrpfZ9pW2HjX4OYSXhVjaaSuom7i9TbisNVs569VdO8P+R+NO+qc93N3ApPjEsAbZIue0KVkzcv2sivjnAbxr68Jp43HbAnJT/JwAXpfx/KvCv4O9LgEdSblsV+IU0iTtYL/NYvuHeDTyZcvtxpOzwgmXvAMemxHJ1ym29gtdrl9IgN065/WJgbMr/JfjEPyBl2VPA/4D/Ah3rve/UxP1Fym3bBK+1bsqyH4Dt0rzvm4GbGtpwgmXPAYPrxVmNP8X2Z+CtlNsMXzC1JHHvnvL/WOCClPc1sd5zLyQlSQC7AF8Hf/+DIIkH/2/Kion7oSbaV6PropXb4gM0XGgvJthJp3lcP2BWyv8TU58Hv6NZHKyX4aQk9uD2l4FYvWUVwJMsL9qSwO9Tbj+AIFdQt9D+HT5Z12/7I9LEfiJ+B/4TPl9UU7cIakmh3az3ls8/KF+D8rXydYj5GlgTuBr4CFgGfADs0NBzB5/n8JTH3gA8l/L/QcAH9dbbvvXa/cvB3w/h8/KG9eI5Fnin3rI3geOa2pbqPWZTYCY+53ZoYh08gN+n/ATUAF8D2zbSHhpaL7XtPG07bCyGvO06Eviu9g/nXHXwZ5c0940AOwWnhn4KTgPFgPXwR77a43e6taY18Bx1lpnZucEphLnB862O/3bdKDMrNbO7zZ9OnIc/VbaG1e2v9V3K39Up72uD1DiccwvxCSydQ/DFwfjg/wSwn5l1T7nPDBe0mkAyeJ1a0+rd1oG67zP19g2C+9TGVxPc3iPlPvfgj7jc5pz7uZHYU/uSLQqer/6yLgBmtpOZvWL+1OZc4GQa/ywiwC0pbeFHfBLtwYrr2NFwe2hMus+vvu74nevklFj+FSynfixp4qjfLlu6LmofF7PlXTGea+r+aXzvnPsl5Tm7mNm9wWnWefjT1fVjqd++OuKPfESAP9XbZncmpW2a2Wn4o0PHpLThOm0w+Du1/ZFyv+kNtP3GTHTOrRHENx7YvYn7p9PkeytAytfK18rXWczXzrk5zrkLnB8Avi6+0H7CzNJ1ran/eTX4+aV5L6ntcBj+83nHfLekE4Ll9XNz7eNS21uTn4Vz7ktgKL5Inmlmj5hZY7nz+iBvR4P3sUUj921MY+0wrXwvtBvj6v0/DXjVObdGyk8X59wp+FNFS4ENU+7fs7HnNN+/bxhwBLBm8CHOpXl9w87Bf9A7Oee64k8r0szHfpsam5mV4k9tpvMXfEOtMt9v9lF84j065T496m14ZfijJrV61rttCf40bq3Udf0NvjHWxmfB42cE/3fBf2O/DxhhQf/BNvAw/shLT+fc6vi+dbXvqX5bAN8e/lqvPXR2zr3BiuvYaLg9tEb9WGbjN/zeKXGs7vzADYJYmt0uA42ti/SBOZdwywc3tnY2jfqxnAdsBOwYtPXfNfCY+u3rZ3wCm4Y/6pv6Ga3qnLsOwHxf0UuAgc65+SnPUacNBs85o4HXrb9ua+/bpOD1TgEGm9m2weKF+J0wQXztSb9tNvreipDytad8rXydkXztnJsNXI8vdtvqc6zf1r4JXus759xJzrkNgL8Cd5jvp14/N9c+rqH83Cjn3MPOud2D53PANc14TBW+y98ttrxffZ28jf8yn05j7TCtQi60vwc2Tvn/GWBz8wNTOgQ/O5jZVs5Pb/M4PpGUmp/6panR3qvhk/0soL2ZXYLvg9ocq+E31p+CxHVpC97XP4EDzWx384NkLifN52h+0NlewIH4053bAb/BN8jU97cOMCRYJ3/E9z8an3L7MWbWK9hJXI7vY5tuSqCxwAFmtpeZdcDvpH7G95kFuAWY5Jw7EXgWn1TawmrAj865xWa2I3V3TLPwp4xS28NdwIVm1hvAzFYP3jtBXL3N7NCgWBpC4xtfS3wPbBh8drVHkO4BbjKzdYJYepjZH4L7jwWON7OtgvV/cTNeo7F1kW2r4Y9KzDGzbvjCuL4/mx/4tipwGf5UtsP3pTvE/ECnduYHOe5pflBWBN8P8Bjn3NR6zzcauMTM1g6OBF6M7ytZ30SgxMxON7P25gfbbd/cN+acm4U/VVz7mXwKrGZmfwja/qX4Iqkhad9bc1+/wChfK1/XUr5uI2Z2jZltHeS31fAHB750zjV2VqUlzjM/WLgnvoAdE7zuH235gOE5+EK4Bt9ONzezo4OYjsR3b3qmhe9rCzP7nZl1xHcLWRQ8f5Occy/iC/54sOgDYA/z84qvDlzYyMMba4dpFXKhfQtwuPmRobcGR6B+jx8U8A3+9MQ1LJ9l4HT8qcTv8DvB0fiEk87z+FNGn+NPfSym+aerbsZ39J8NvBU8T7M45z4CTsN/C/4W34inp7n7sfg+VS8E3zC/c859B9wKbGtmWwf3exvYLIinHDi83ob4f/h+Tt/hBzYMaSS+z/CDjm4Lnu8g4CDn3C9mNgg/WOSU4O5nA9ubWay5778RpwKXm9l8fDE3NiWm6uB9vW7+lM/Ozrlx+M//EfOng6cA+wX3nw38Ed+37Qf8unm9DWIE33XiI+A7M6s9ynQ+foDFW0EsLxGc2nJ+hPit+EEnX+LbCzTeNtOuixDciN+ufsDvvBs6xfl/+EL4W3xf0qEAzrlK/Kn0i/E73yp8IVCCH7y1DjDOlp8+/TB4vsuAD/Gf6X/x7fuq+i8anAY/BD/ga07w9xMtfH83AQPNrJdzbg5wBvAg/gjNj9Q9DZr62o29t2KkfK18XRuT8nXbKQXG4fsnf4U/+juwDZ//SWAyvlh9luXT5u0AvG1mC/BH6890zn0VtNMD8bnuB/xZpgODz7AlOuI/79n4dr4OjRfI9V0HDDOzjkHhPQa/r5hMI0V/Y+2wMbWDh6QeM7sGP1DnL2HHkknmJ68/MTgF09DtE/ADBe7NZlzSMPNTTk3BD0paGnY8K8v8lGz3OuceCDsWyV/K17/ePgHl65xRaPk6lZk5YDPn+0tLI4r16MkKglPX25q3IzAY/01QJFRmdoj5uVvXxH+bfrrQkrZISyhfS65Svpb6slJom1lP86NqPzY/AvXMBu4zwPxo8A+Cn4b6cmbSavh+fwvxpxFuwJ8WEQnbX/FTGU3FT9F0SuN3F1k5eZCzla8lVylfSx1Z6TpiZusD6zvn3gs65E/GzxP6ccp9BgDnOucOzHhAIiKSlnK2iEjbyMoRbefct86594K/5wOf0MS8gyIiEg7lbBGRttHgtdwzycyi+Mvpvt3AzbsEMwd8gz9S8lEDj48TTMuy6qqr9t1yyy0zF6yISAZNnjx5tnOue9P3DI9ytoiI15qcndVZR8xPfv8qUO6ce7zebV2BGufcAjPbH7jFObdZY8/Xr18/N2nSpMwFLCKSQWY22TnXL+w40lHOFhFZrjU5O2uzjpifDP8xIFE/YQM45+Y55xYEf48HOphZk5chFRGRtqecLSKy8rI164jhJzL/xDl3Y5r7rBfcj2C6phL8hOYiIpJFytkiIm0jW320d8Nf9ep/ZvZBsOxv+Gvc45y7CzgcOMXMluIvp3mU09V0RETCoJwtItIGslJoO+cmAtbEfUYCI7MRj4iIpKecLSLSNnRlSBERERGRDFChLSIiIiKSASq0RUREREQyQIW2iIiIiEgGqNAWEREREckAFdoiIiIiIhmgQltEREREJANUaIuIiIiIZIAKbRERERGRDFChLSIiIiKSASq0RUREREQyQIW2iIiIiEgGqNAWEREREckAFdoiIiIiIhmgQltEREREJANUaIuIiIiIZIAKbRERERGRDFChLSIiIiKSASq0RUREREQyQIW2iIiIiEgGqNAWEREREckAFdoiIiIiIhmgQltEREREJANUaIuIiIiIZIAKbRERERGRDFChLSIiIiKSASq0RUREREQyQIW2iIiIiEgGqNAWEREREckAFdoiIiIiIhmgQltEREREJANUaIuIiIiIZIAKbRERERGRDMhKoW1mPc3sFTP72Mw+MrMzG7iPmdmtZvalmf3XzLbPRmwiIlKXcrbkm0QiQTQapaSkhGg0SiKRCDskEQDaZ+l1lgLnOOfeM7PVgMlm9qJz7uOU++wHbBb87ATcGfwWEZHsUs6WvJFIJIjH41RXVwOQTCaJx+MAxGKxMEMTyc4Rbefct86594K/5wOfAD3q3W0Q8JDz3gLWMLP1sxGfFLmaGvjsM/jXv+Dee0kcdhjR1VenxIzo6quTOOIIeOgh+Pe/YcaMsKMVyTjlbMl5M2b4nPzQQww//fRfi+xa1dXVDD/7bJ/ba2pCClIke0e0f2VmUaAP8Ha9m3oA01L+nx4s+7be4+NAHKCsrCxTYUohcw6mTIEnnoCJE+Htt2HuXAAS+MZVm7KT8+YRf/RRePRRfj0usuGGsPPO8NvfwiGHQI/69YdI4VDOlpwwYwaMGwevvgpvvQXTp/96U1Wah1TNnAlbbglrrAE77QS77QYHHwxbbw1m2Ylbil5WB0OaWRfgMWCoc25ea57DOVfhnOvnnOvXvXv3tg1QClsyCRddBFtsAdtuC5deCt99B0ceCffdB6+/zvAePaiu97BqYPgGG8DLL8PNN0P//jBpEpxxhi+6d90VbrsN5rWqSYvkLOVsyZYG+1jPm+dz6667+lx7xhk+9/bv73Pxyy/D559T1rNng89Ztu66PrcfcQR8+63P+dtu6/cBF10EVVXpX1ukrTjnsvIDdACeB85Oc/vdwJ9S/v8MWL+x5+zbt68TadJrrzl32GHOlZT4n733du6uu5z77rsV7mpmDljhx8xWfN5PPnHuiiuc+81vnAPnunRx7owznPviiyy8KSkEwCSXpRzc0h/lbMmWUaNGudLS0jo5t7R9ezeqY0efW3/zG59rP/mk+Y8vLXWjRo2qe8dvv3Xuzjv9PiDYH4zaYQdX2rFj048Vca3L2dlK2AY8BNzcyH0OAJ4L7rsz8E5Tz6ukLemMGjXKRdZbzxm4CLhRpaXOnX++c8lko4+LRCINFtqRSKTxF3znHeeOOca5Dh3cKDMX6dLFmZmLRCJK2JJWrhbaytmSTWnz7qqrOvfuu816jlGjRrlIJNL8vJtMOjdsmIuUlLQu50tRak3ONv+4zDKz3YHXgP8BtaMS/gaUATjn7jIzA0YC++LP1h/vnJvU2PP269fPTZrU6F2kCCVuu434WWdRvWzZr8tKO3em4p57mhyBXn/0OkBpaSkVFRXNGr2euP124kOHUr10aYtfW4qPmU12zvULO476lLMlm0pKSmioFjEzajI8kDHM15b805qcnZVCO1OUtKWOmhq4+26ip51GsoF2HYlEqKysbPJpEokEw4cPp6qqirKyMsrLy5tdJEejUZLJ5Iqvvc46VH7/fbOeQ4pHrhbamaKcLSt46SWi++1HMuXgRK3m5uyVkTZnr7UWlbNmQYmu6yfLtSZnqwVJYZg61c8CcuqpVKX58lhVlW5sel2xWIzKykpqamqorKxs0ZHodK9RNXMmDB4M8+c3+7lERArW/Pk+J+6zD+XdulHasWOdm0tLSykvL894GOXl5ZSWltZ97ZISyn/80e9Tpk7NeAxS2FRoS/4bMwb69PFT9t1/f9opxLIxtVja1+7aFR54APr2hfffz3gcIiI56733YPvtfU684AJiX39NxX33EYlEMDMikUizu+utrFgsRkVFRd3XfughYvff7/cp228PY8dmPA4pXCq0JX8tWgR//SscdZSfF/WDD+C44yi/8soVj1CEeXSktJTyO+6AV16B6mo/B/ftt/v5vEVEioVzMHIk7LKLz98TJsBVV0Hnzit1JnFlNfjaxx3n9ym9evkpYE8+2ccs0kIqtCU/ffstDBgAFRVw/vn+IgaRCJDmCEWYR0dqX3uPPXzi3ntvOP10f9r0l18yHpOISOh+/hlOOMHPhb3PPj4X9u8fdlSNi0TgP/+BYcPg7rthzz39tRdEWkCDISX/fPABHHQQ/PgjJBL+Sl/5pKYGLrsMLr/c72gefxzWXjvsqCQEGgwpRWH2bDj0UHjtNX/RmEsuyb9Bhk88AbEYdOsGTz8Nv/lN2BFJCDQYUgrfs8/C7rv7v19/Pf+KbPA7mMsug4cfhnfe8ZcG/uKLsKMSEWl7n3/uc9w778Do0TBiRP4V2eD3NRMn+u4vu+0G48eHHZHkiTxs7VK0Hn7YJ7stt/RJe7vtwo5o5fzpT77Ly/z5/svDhx+GHZGISNup7R4yf77PdUcdFXZEK6dPH7/v2XJLGDTIf3EQaYIKbckPd94Jxxzjk/Yrr8D664cdUdvYaSd/OrVjRz+V1Ouvhx2RiEirJBIJotEoJSUlRNdbj8Suu/rc9tprPtcVgvXXh3//2x8cicVInHDC8vccjZJIJMKOUHKMCm3JWb8mbTOip55Kok8ff7putdXCDq1tbbGFPyW57rok9tyT6LrrKmmLSF6pvapuMpnEOUfy+++JL15M4txzfY4rJF27wvjxJLbbjvj99y9/z8kk8XhceVvqUKEtOalO0gaSQPzTT0k89ljYoWVGWRmJoUOJL11KcuZMJW0RySvDhw+nurq6zrJq5xh+440hRZRhnTsz/IcfqK63uLq6muHDh4cSkuQmzToiOSntZXGzcEnesBTjey52mnVECkVJSQkN1RNmRk1NTQgRZV4xvudip1lHpGBUNVBwQvMvo56P0l6+vYDfs4gUhrJ11214eRauyBuWMK9CLPlDhbbknieeIF2aKuQEljZpF8rATxEpTJ9/TvmCBZSa1VmcrSvyhqXBKwED5YccEk5AkpNUaEtuee01OOooyjfZhNLOnevcVLRJ+5df/JUwRURyzbffwh/+QKxzZyquvz6UK/KGZYUrAffsScUmmxC76y4/wF0EFdqSS6ZMgYEDIRol9tZbVNxzT3En7UiEissvJ7ZoEey3H8ydG3aIIiLLzZ3rc9OsWTB+PLGzz6ayspKamhoqKysLOl/XisViy99zVRWxt96CsjJ/9eIpU8IOT3KABkNKbvjuO9hxR1i6FN54A6LRsCPKHc8/DwceCL/7nb8yZvv2YUckbUSDISVvLV0KBxzg55R+5hn4wx/Cjih3VFbCrrtChw7+Ajdp+q9L/tFgSMlPixb5Kz7+8INP2Cqy6/rDH/wFe154Ac46K+xoRERg6FCfk+66S0V2fdEoPP20P9J/8MGweHHYEUmIVGhLuJyDwYPh7bdh1CjYfvuwI8pNJ54I55wDI0fCHXeEHY2IFLPbb/c/557r87esqG9fv0976y044QS/r5OipEJbwnXFFTB6NFx1FWikduOuucZ3IRkyBF58MexoRKQYvfginHmmH09z9dVhR5PbDj0UrrzS7+OuuCLsaCQkKrQlPM88A5dcAsceC+efH3Y0ua9dO3j4YdhqKzjqKN8PUEQkW77+Go48Enr18kdr27ULO6Lcd8EFfh93ySV+nydFR4W2hOOLL+CYY3xXkbvvhnrzr0oaq60G48bBsmX+aMmiRWFHJCLFoLra5xznfA5abbWwI8oPZn4f16eP3+d9+WXYEUmWqdCW7FuwwCfsdu3gsceg3nzZ0oRNN4VEAt5/H04+WX3/RCSznPO55sMPfe7ZZJOwI8ovnTvD44/7fd4hh8DChWFHJFmkQluyyzk46ST4+GN45BHNMNJaBxwAl14KDz3kZyQREcmUO++E//s/n3P23z/saPJTNOr7an/0kR/crgMkRUOFtmTX3Xf7Avvvf4d99gk7mvx2ySV+p3fWWfDee2FHIyKF6L33fI7Zf3+4+OKwo8lvv/+9HxT5yCNQURF2NJIlKrQlez74wM+9uu++foCIrJySEnjwQejeHY44AubNCzsiESkk8+b53LLOOv7sWYlKhpV2wQV+3vEzz/RdcaTgaauR7Jg/3yfsbt2UsNvS2mv7oyOVlRCP63SkiLSN2m5+lZU+x3TrFnZEhaGkxO8D11rL7xPnzw87IskwVTuSebUDaaZO9X3UuncPO6LCsvvuvivOmDE6HSkiKyWRSBCNRikpKSE6diyJww6D3XYLO6zCss46MHo0ic8/J7reen5dR6MkEomwI5MMUKEtGZVIJIh2707Jww8TXW01EtOmhR1SYTr/fBLbbEP05JOVtEWkVRKJBPF4nGQyiQOSQPyZZ5RLMiAxfTrxDh1IVlfjnCOZTBKPx7WuC5AKbcmYRCJB/KSTSP7wg0/ac+cqkWRIYvRo4l9+SRKUtEWkVYYPH051dXWdZdXV1QwfPjykiArX8OHDqV6ypM4yrevCZC6P+3T269fPTZo0KewwJI1oJEKyqmqF5ZFIhEpd1bBNRaNRksnkCsu1rnObmU12zvULO45sUc7ObSUlJTRUE5gZNTU1IURUuLSu81NrcraOaEvGVDVQZDe2XFpP61pEVlZZmvEzZWVlWY6k8KVbp1rXhUeFtmTGO++QLl0okbS9tEm7R48sRyIieWnOHMp/+YVSszqLS0tLKS8vDymowlVeXk5paWmdZaVA+QknhBOQZIwKbWl7ixbBn/9M+VprUVrv8upK2pmRNmlHIuEEJCL5ZcgQYgsWUHH55UQiEcyMSCRCRUUFsVgs7OgKTiwWo6KiYvm63nBDKtZck9jo0X4fKgUjK4W2mf3DzGaa2ZQ0tw8ws7lm9kHwc0k24pIM+dvf4LPPiI0ZQ8U99yhpZ8EKSTsSoeLQQ4m9/jqMGxd2eJKHlLeLyOOPw6hRcNFFxC66iMrKSmpqaqisrFS+zqBYLLZ8XU+bRmzsWPj0U9CAyIKSlcGQZrYHsAB4yDm3dQPnFD9dAAAgAElEQVS3DwDOdc4d2JLn1cCaHDRhAuy5J5x2GowcGXY0xW3JEth5Z5g2DaZM8XO3Sk7J5cGQmcjbytk5aOZM6N0bIhF4803o0CHsiIrbaafBnXfCK6/Ab38bdjRST84OhnTO/Qf4MRuvJSGaPx+OPx423RSuuSbsaKRDB38Fsnnz/AWD8niGIck+5e0iUHsxsfnzfa5QkR2+a6+FTTaB446DBQvCjkbaQC710d7FzD40s+fMrHe6O5lZ3MwmmdmkWbNmZTM+acqFF0IyCQ8+CKuuGnY0Av5I1eWX++4jjz4adjRSeJrM28rZOWzsWJ8b/v536NUr7GgE/L7zgQf8vvTCC8OORtpA1ubRNrMo8EyaU5BdgRrn3AIz2x+4xTm3WVPPqdOQOeTVV2HAABg6FG66KexoJNXSpbDrrvD11/Dxx5BmCi/JvlzuOgJtn7eVs3PIrFm+uN54Y3jjDWjXLuyIJNXQoXDLLX7fusceYUcjgZztOtIU59w859yC4O/xQAczWzvksKS5qqth8GB/uksziuSe9u3h/vth7lwYMiTsaKRAKG/nuTPO8N3K/vEPFdm5qLzcfwkaPNjvYyVv5UShbWbrmfnJO81sR3xcP4QblTTbxRfD1Klw771Qb4o5yRG9e8Mll8Ajj8ATT4QdjRQA5e08Nm4cjBnjc0LvtD01JUyrrur3qV9+6T8nyVvZmt5vNPAmsIWZTTezwWZ2spmdHNzlcGCKmX0I3Aoc5fL52vDF5O23fVeRU07xXUckd51/Pmy3nf+s5swJOxrJccrbBWrOHDj1VOjTB4YNCzsaacyee/rBqjfe6Pe1kpey1kc7E9TfL2S//AJ9+/rE/fHH0LVr2BFJU957D3bcEU44ASoqwo6m6OV6H+22ppydA046yXcle/ddX2xLbps3z/elX2stmDxZM8OELG/7aEueuu46Pz/zHXeoyM4X228PZ50F99zjB9mISPF49VXfHeHss1Vk54uuXf0+9n//8/tcyTs6oi2t8/nnsO22MHCgnyJK8kd1NWy9tT8y8uGH0KlT2BEVLR3RlqxZvBh+8xs/C9H//qfxNPnmj3+Ep5+G//4XNt887GiKlo5oS3bU1EA8Dp07w623hh2NtFRpKdx9t/+ypFliRIrDFVf4bf6uu1Rk56Nbb/UHRf76V118LM+o0JaWe+ABfwryuutgvfXCjkZaY5994Nhj4eqr4aOPwo5GRDLpo4/81Xr//Ge/7Uv+WX99v8+dMMHvgyVvqNCWlpk9G847D3bf3Q+ok/x1ww2+/9/JJ/uzFCJSeGpq/Dbetavf5iV/DR4Mu+3m98GzZ4cdjTSTCm1pmfPO86Og77oLStR88lr37nDttTBxoo6QiBSq++/32/h118Haup5QXisp8d3+5s7V1Ix5RJWSNN+rr/qC7LzzdJGDQnH88f7sxHnn+Usyi0jhmDXLF2T9+/ttXfJf795w7rn+C9R//hN2NNIMKrSleX7+2Z9+3GgjuOiisKORtlJS4s9OzJvni20RyXuJRIJoNErJOusQ/fFHEgceCP4inlIILr6YxNprE91rL0pKSohGoyQSibCjkjRUaEuTEokE0fXWo+TTT4kuXEhi3LiwQ5K21Ls3if33J/rgg0raInkukUgQj8dJJpM4IAnEL7tM23QBSYwbR3z+fJJLl+KcI5lMEo/H9RnnKM2jLY1KJBLETzqJ6kWLfl1WWlpKRUUFsVgsxMikregzDo/m0Za2Fo1GSSaTKyyPRCJUVlZmPyBpc/qMw9OanK1CWxqlDbrw6TMOjwptaWslJSU0tF83M2o0u1BB0GccHl2wRtpcVVVVi5ZL/tFnLFI4yjbYoOHlZWVZjkQyJd1nqc84N7W40DazVc2sXSaCkRxTXU1Zmin8tEEXDiXtwqacXVzKe/Sg/nUfS0tLKddVYAtGeXk5pfWu7llqRvmll4YUkTSmyULbzErM7Ggze9bMZgKfAt+a2cdmdp2ZbZr5MCUUV19N+bJllHbsWGexknZhaTBpA+VHHRVOQLJSlLOL2AsvEHvnHSoOP5xIJIKZEYlENN6iwMRiMSoqKpZ/xuuuS4VzxNTVLyc12UfbzF4FXgKeBKY452qC5WsBewJHA+Occ6MyHOsK1N8vg6ZO9fN1HnYYif33Z/jw4VRVVVFWVkZ5ebmSdoFJJBLLP+OePSlfvJjYmmvCf/8Lq6wSdngFKxN9tJWzi9Qvv8A22/grQU6ZAvUOkEiBO/poePxx+Phj2HjjsKMpWBkZDGlmHZxzS1b2PpmgpJ1BBx0EEybAZ59Bmj5/UsDGj4cDDvBXjtT82hmToUJbObsYXXstnH++33b32y/saCTbZsyALbeEPfeEp54KO5qClZHBkLXJ2MxuMWt4xvswErZk0DPP+J9LL1WRXaz23x8GDoTLLvMJXPKGcnYRmj4dLr8cBg1SkV2sevSASy6Bp5+GZ58NOxpJ0ZLBkPOBp8xsVQAz+4OZvZ6ZsCQ0ixfDmWf6b8ZDhoQdjYTppptg6VId0c5fytnF4rzzYNkyv81K8ardd595pt+XS05odqHtnLsIGA1MCJL12cAFmQpMQnL99fDVVzBypPrmFruNN4YLLoDRo303IskrytlFYsIEeOQRv61utFHY0UiYVlkFbrvNj7G6/vqwo5FAsy9YY2Z7ARcBBqwPDHTOfZbB2Jqk/n5trKrKfxvef3/45z/DjkZywaJFsNVW0LUrvPcetG8fdkQFJZMXrFHOLgJLl0KfPrBggR8E17lz2BFJLjj8cHjuOfj0U+jZM+xoCkqmL1gzHLjYOTcAOBwYY2a/a8mLSY4791z/+4Ybwo1DckfnznDjjfC//8Fdd4UdjbSMcnahu/NOP8PIjTeqyJblrr/ezz5Tu0+XULWk68jvnHMTg7//B+wHXJGpwCTLXnkFHn3Un36MRMKORnLJIYfAXnvBxRfDrFlhRyPNpJxd4GbN8oPf9t4bDj447Ggkl0Sjfl8+dqy6/eWA5lywJt2o9W+BvRq7j+SJpUv9wMdoVAPfZEVmcOutMH8+DB8edjTSBOXsIvG3v/kuI7fe6rdRkVTDhvmDZkOG+H28hKY5R7T/bWZnmFmd6zGb2SrALmb2IPCXjEQn2aHTj9KUXr3gjDPg3nth8uSwo5HGKWcXukmT4L77/Da51VZhRyO5KLXb3513hh1NUWvOBWsqgA+AwfgBNT8BnYB2wAvAHc659zMcZ4M0sKYNzJoFm28O/frBCy/oyIikN3cubLaZ/5k4UW2lDWTogjXK2YXMOdhtNz+zxOefw+qrhx2R5CrnYJ99/MGRL76AtdcOO6K8l6nBkDs45+7Aj1wvw5963N45F3HOnRRWwpY2cvHFvkvALbeocJLGrb46XHUVvPEGPPxw2NFIesrZhSyRgDff9NuiimxpjJnft8+f7/f1EormFNovm9mbwLrAn4ENgEUZjUqy4/33oaICTj/ddw0Qacrxx0Pfvr7/34IFYUcjDVPOLlQLFvjLrPfrB8cdF3Y0kg9694bTTvP7+g8+CDuaotScS7CfCxwDLAM2Ai4GppjZR2Y2JsPxSaY45wdJdOsGI0aEHY3ki5ISP/jqm2/8ETXJOcrZBezKK/22d+utflsUaY4RI2DNNf0+v5nXTpG206yrTzjnpprZ3s65z2uXmVkXYOuMRSaZNWaM72dbUQFrrBF2NJJPdt0VjjnGz9U6eLC/gqTkFOXsAjR1qr/GwbHHwi67hB2N5JM11/Rf0v76Vz/l35FHhh1RUWn2lSFzkQbWtNLChf4KkN27w7vvQrt2YUck+eabb/wg2n32gXHjwo4mb2XyypC5SDl7JRx8MLz0kh8AucEGYUcj+WbZMthhB5g9218xsrQ07IjyUqavDCmF4tprYfp0f/pRRba0xgYb+Dm1n3jC7/xFJHNefBGefBIuukhFtrROu3Z+YOS0ab4GkKxRoV1skkm/kR11FOy+e9jRSD476yzfbWToUF0QQSRTli7121jttibSWv37+24j11wDVVVhR1M0VGgXm/PO81P+6ButrKxOnXyf0Y8+grvuCjsakcJ0553w8cf+4iOdOoUdjeS7a6/1NYCuAp01WSm0zewfZjbTzKakud3M7FYz+9LM/mtm22cjrmKRSCSIRqOUlJQQffRREvvtBz17hh2WFIJBg0j06kX0zDN9+4pGSSQSYUclbUB5O1yJRIJoz56UDBlCtFMnEppOU9pCWRmcfz6JsWOJrree8nYWZOuI9gPAvo3cvh+wWfATB3S90DaSSCSIx+Mkk0mccySB+HPPaaOSNpF4+GHiX31FsqbGt69kkng8rvZVGB5AeTsUv+bt6dNxQHLxYm1X0mYSZWXEzUh+/73ydhZkbdYRM4sCzzjnVpheyszuBiY450YH/38GDHDOfdvYc2oEe9Oi0SjJZHKF5ZFIhMrKyuwHJAVF7Wvl5PqsI22dt5Wzm0fblWSS2lfr5fOsIz2AaSn/Tw+WrcDM4mY2ycwmzZo1KyvB5bOqNAMe0i0XaQm1r6LWrLytnN1y2q4kk9S+sitXCu1mc85VOOf6Oef6de/ePexwcl5ZWVmLlou0hNqXNEU5u+XKunVreLm2K2kDytvZlSuF9gwgdXTehsEyWUnlp51G/WnpS0tLKS8vDyUeKSzl5eWU1rvwQakZ5SNGhBOQZJPydiYsWkQ5fjtKpbwtbaXBvN2pk9pXhuRKof0U8OdgFPvOwNym+mdLMzhH7OWXqejcmciGG2JmRCIRKioqiMViYUcnBSAWi1FRUUEkEvHta511qHCO2AzVW0VAeTsTbryR2OzZVFx44fLtSnlb2tAKeduMis02I3b00WGHVpCyMhjSzEYDA4C1ge+BS4EOAM65u8zMgJH4Ee7VwPHOuSZHzGhgTROefhoGDoSbbtKFDiR7DjsMnn8ePvsMejQ41EICuTwYMhN5Wzm7CdOnwxZbwL77wmOPhR2NFIubb/YXIHv6aTjwwLCjyWmtydlZm3UkE5S0G/Hzz9C7N6yyCnz4IXToEHZEUiy+/hq22goOPxxGjQo7mpyWy4V2JihnNyEW8wX2J5/ARhuFHY0UiyVLYNtt/VVIp0yBjh3Djihn5fOsI9LWbr4Zpk71v1VkSzZttBGcey4kEvDGG2FHI5IfXn8dHn7YX7FPRbZkU4cOvlb48ku45Zawoyk4OqJdiL79FjbfHH73O3jyybCjkWK0cKE/Bb7eevDOO1Ci7/QN0RFtAaCmBnbYAb7/3ne5WnXVsCOSYjRwILzyCnzxhc/dsgId0Rbvwgvhl1/ghhvCjkSK1aqrwjXXwOTJcP/9YUcjktvuvx/eew+uvVZFtoTnhht8t9MLLgg7koKiQrvQvP02PPigH9iw6aZhRyPF7OijYddd4W9/g7lzw45GJDf99JM/OLLbbvCnP4UdjRSzzTbztcODD/paQtqECu1CUlMDZ5wB668Pw4eHHY0UOzO47TaYNQsuvzzsaERy0+WXw+zZflupN3e2SNZddJHvNjJkiK8pZKWp0C4kDz4I777rTz+utlrY0YjA9tvDiSfCrbf6mRREZLmPP/YF9kknQZ8+YUcj4muHa6/1Y2seeijsaAqCCu1CMXeu71e1yy5+iiiRXFFe7vudDh0KeTz4WqRNOQdnngldusAVV4QdjchysZivJS64AObNCzuavKdCu1D8/e/+FL1OP0qu6d7dnx5/4QXNgiNS68kn4aWX/LbRvXvY0YgsV1Liz0LOnKluf21AhXYh+OQTP/fl4MHQt2/Y0Yis6JRT/AWUzj4bFi0KOxqRcC1a5Aed9e7ttw2RXNOvH5xwgq8tPv007GjymgrtfOecH7TQpQtceWXY0Yg0rEMHf7bl66/huuvCjkYkXNdeC5WVMHIktG8fdjQiDbvySt/t74wz1O1vJajQznePP+5PP/797zr9KLltzz3hiCPgqqt8kSFSjCor4eqr4cgjYcCAsKMRSW+ddXxt8dJLMG5c2NHkLRXa+ay62p+K33ZbOPnksKMRadr11/v+f+ecE3YkIuE4+2y/DVx/fdiRiDTtlFNgm218u62uDjuavKRCO59dfTVUVflT8jr9KPmgZ08/x/vjj8OLL4YdjUh2vfCCPzJ40UWw4YZhRyPStPbtfY2RTPqr/UqLmcvjfjf9+vVzkyZNCjuMcEyd6gfSHHooPPxw2NGINN/PP8PWW0O7dvDf/8Iqq4QdUWjMbLJzrl/YcWRLUefsX37xZx+XLYMpU6Bjx7AjEmm+o4/2B0g++gg22STsaELTmpytI9r5qHYAZIcOOv0o+adjRz911Gefwc03hx2NSHbcdJNv87fdpiJb8s911/maY+jQsCPJOyq089HTT8P48XDZZbDBBmFHI9Jy++0Hgwb5OVqnTQs7GpHMmjbNt/WDD4Z99w07GpGW69EDRoyAZ57xNYg0mwrtfLNokb+aWO/efsodkXx1883+NLoGRkqhO/tsfybyppvCjkSk9YYMgV69/G9dD6HZVGjniUQiQTQapaS0lGhlJYmDD/ancUTyVTQKw4eTePRRouuuS0lJCdFolEQiEXZkIivt15xdUkL0n/8kccABvs2L5KsOHeD220lUVhJdf33l7GZSoZ0HEokE8XicZDKJA5JA/Kab1Lgl7yU23JC4GcmZM3HOkUwmicfjatuS1+rkbOd8zh4/Xu1a8l5ixgzi7dqRnDtXObuZNOtIHohGoySTyRWWRyIRKnXhD8ljxd62NetIYSr2di2Fq9jbtmYdKVBVVVUtWi6SL9S2pRCpXUuhUttuORXaeaAszYUNysrKshyJSNtK14bVtiWflfXs2fBytWvJc8rZLadCOw+U9+pFab1lpaWllJeXhxKPSFspLy+ntLRu6y5dZRW1bclr5QcfrJwtBanBnN2+vdp2I1Ro57r33iP24otU7LUXkUgEMyMSiVBRUUEsFgs7OpGVEovFqKioWN62O3SgorSU2EEHhR2aSOvMm0fssceoKCsjUlamnC0FZYWc3aULFUuXEuvVK+zQcpYGQ+aypUth551hxgz45BNYY42wIxLJrLffhl128XPE33JL2NFknAZDFqAhQ2DkSHjzTdhpp7CjEcmsn36CLbeEnj3hrbegXbuwI8ooDYYsNCNHwuTJvuBQkS3FYKed4NRT/WWq33037GhEWuadd3zePu00FdlSHNZYw9cokyb5ti8r0BHtXDVtGmy1FeyxBzz7LJiFHZFIdsyb59v+Ouv4Yrt9+7Ajyhgd0S4gS5bADjvA7Nnw8cfQtWvYEYlkh3Ow//4wcaJv+2kGAxcCHdEuFM7B6af733fcoSJbikvXrv6I9gcfFEX3ESkQt9wCH37o266KbCkmZr5WWbbMd/uTOlRo56LHHoOnnoLLLtMle6U4HXIIDBwIl1wCX38ddjQijfvqK7j0Uhg0yLddkWKz0Ua+ZnnySV/DyK9UaOeaH3/0R7O33x6GDg07GpFwmMHtt/uBNfG4P7sjkouc8220XTv1UZXidtZZ0KePH6MwZ07Y0eQMFdq55pxzfB+/++4r6L6pIk3acEO49lp46SV44IGwoxFp2P33w8svw3XX+TYrUqzat/e1y+zZvpYRQIV2bnnxRV9QnH8+bLdd2NGIhC8e9wOCzz4bvvsu7GhE6vr2W19Q7LEHnHRS2NGIhK9PHxg2zH8BfemlsKPJCVkrtM1sXzP7zMy+NLMLGrj9ODObZWYfBD8nZiu2nLBggS8qttgCLr447GhEckNJCdxzDyxa5E9HStYoZzfD6af7tnnPPb6tiogfW7P55v7L58KFYUcTuqxkBjNrB9wO7Af0Av5kZg1dRmiMc2674OfebMSWMy68ECorfcLu1CnsaERyx+abw4gR8Pjj8OijYUdTFJSzm2HsWN8mR4zwbVREvE6d4N57fU1z4YVhRxO6bH0F3xH40jn3lXPuF+ARYFCWXjv3TZjgB9GceSb07x92NCK559xzoV8/fzGbmTPDjqYYKGc35vvvfVvcYQffNkWkrv79/VVSb7sNXn017GhCla1CuwcwLeX/6cGy+g4zs/+a2T/NrHBnPE+1YAEcfzxsuilceWXY0Yjkpvbt4cEH/cVsTj1Vs5BknnJ2Os75Njh/vh9To0HrIg278krYZBNf4yxYEHY0ocmlTmVPA1Hn3LbAi8CDDd3JzOJmNsnMJs2aNSurAWbE+edDMukHDpSWhh2NSO7q1cvP0/rYY/60vYStOHP2mDG+y8jll/s2KSINW3VVX9tUVsIFKwzzKBrZKrRnAKlHOzYMlv3KOfeDc+7n4N97gb4NPZFzrsI5188516979+4ZCTZr/v1vfzWloUNh993DjkYk9517Luy4ox8YqVlIMkk5uyHffecHQO60k6YvE2mO/v19t9jbb/c1TxHKVqH9LrCZmW1kZqsARwFPpd7BzNZP+Xcg8EmWYgvHTz/BX/7iZxm54oqwoxHJD+3b+9P1CxfCiSeqC0nmKGfX5xwMHuzbnrqMiDRfebkfMHzccb72KTJZKbSdc0uB04Hn8cl4rHPuIzO73MwGBncbYmYfmdmHwBDguGzEFprTT/dHR0aNUpcRkZbYait/IZtnn/Wz9EibU85uQEUFjB/vL0yz5ZZhRyOSP0pLfa3zzTdwxhlhR5N1Weuj7Zwb75zb3Dm3iXOuPFh2iXPuqeDvC51zvZ1zv3HO7emc+zRbsWVLIpEgGo1SUlJCNJEgMWiQn0lBRFrmtNNgn31InH460R49/DYVjZJIJMKOrGAoZ9fL2aecQmLrrf1ASBFpmR12gEsuITFqFNHu3YsqZ+fSYMiClkgkiMfjJJNJnHMkgfhzzxVFIxNpcyUlJAYOJL5kCclvvvHbVDJJPB7XNiVtYoWc7RzxqVNJjB4ddmgieSmx0UbES0pIzp5dVDnbXB73cezXr5+bNGlS2GE0SzQaJZlMrrA8EolQWVmZ/YBE8lwhbFNmNtk5VzSntZSzRYpXIWxTrcnZOqKdJVVVVS1aLiKN0zYlmaT2JdK2inWbUqGdJWXrrNPw8rKyLEciUhjSbTvapqQtlPVo6Po8al8irVWsOVuFdjb89BPlS5dSalZncWlpKeXl5SEFJZLfysvLKa03Y09pu3bapmTlOUd59+7Unw9KOVuk9RrM2WaUDx8eUkTZoUI704K5V2M//UTFpZcSiUQwMyKRCBUVFcRisbAjFMlLsViMioqK5dvUmmtSsWwZsZkzww5N8t3NNxN7/30qjjlGOVukjayQs9dbjwog9vzzBX1NBA2GzLSbboKzz4brr9eVxEQyyTk49FB45hl49VXYddewI2qSBkPmoNdfhwED4MAD/aXW652JFJE2dP31cN55vlYaOjTsaJqkwZC55vXXYdgwOOQQX2yLSOaYwf33Q1kZHHEEzJoVdkSSb2bOhCOPhEjEtyUV2SKZdc45cPDBvth+442wo8kIFdqZMnOm39krYYtkzxprwD//CbNnQywGy5aFHZHki2XLfJuZPdu3oTXWCDsikcJXe4AkEvE1UwF2/VOhnQlLlvijIj/+6BP26quHHZFI8ejTB0aOhBdfhIsvDjsayRcXXwwvvQS33w7bbRd2NCLFI/UAyVFH+RqqgKjQzoSzz4YJE6CiQglbJAyDB8NJJ8FVV8GYMWFHI7luzBjfVk46CU44IexoRIrPdtvBPffAK68U3Hg2Fdpt7b77/NG0c86BY48NOxqR4mTmt8PddoPjj4cPPgg7IslV77/v28huu/k2o25+IuE49lh/oPK22+Af/wg7mjajQrstvfEGnHIK/P73cPXVYUcjUtxWWQUeewy6dYNBgwqy75+spJkz/UCstdf2bWWVVcKOSKS4XXMN7LOPr6XefDPsaNqECu22MnWqT9hlZfDII9C+fdgRici668ITT/iCatAgWLQo7IgkVyxaBAMH+rYxbpxvKyISrvbtfQ3Vs6fP2V99FXZEK02Fdlv44QfYf38/an38eFhzzbAjEpFafftCIgFvv+1PTdbUhB2RhK2mBo45Bt55Bx5+2LcREckNa60Fzz7ra6r99/cTS+QxFdora/FifyS7shKefBI23zzsiESkvkMPhRtu8N0Dhg0LOxoJ23nn+YvR3Hijv86BiOSWLbbwZyO//trXWIsXhx1Rq6nQXhnLlsFf/gITJ8JDD8Huu4cdkYikM3QonH66L7hvvTXsaCQst9ziC+wzzoAzzww7GhFJp39/ePBBeO01X2vl6XURVGi3lnO+s/7YsXDddX7ebBHJXWZw883+COaZZ/ovx1LwEokE0WiUkpISomuvTWLoUN8GbrpJM4yI5LqjjvI11tixJPbZh2gk4rflaJREIhF2dM2iQruFEonE8g/6nntIDBwI554bdlgi0hzt2sHo0bD33iT+8hei66yTd0lbmi+RSBCPx0kmkzjnSP7wA/GSEhKDBvm2ICK579xzSQwcSPyVV0hWVfltOZkkHo/nRd5Wod0CvybtqiockATiL72UFx+0iAQ6diTxpz8RLykhOWtW3iVtab7hw4dTXV1dZ1l1TQ3DL700pIhEpDWGf/AB1fWWVVdXM3z48FDiaQlzzoUdQ6v169fPTZo0KWuvF41GSSaTKyyPRCJUVlZmLQ4RWTm5si2b2WTnXL+svWDIsp2zS0pKaGgfZ2bUaPYZkbyRK9tya3K2jmg3l3NUNbBjBqiqqspyMCKyMtJts9qWC0vZ2ms3vLysLMuRiMjKSLfNlvXsmeVIWk6FdnM4BxdeSLrUrKQtkl/SJu1u3bIciWTMY49RPns2pSV1d3OlpaWUl5eHFJSItEZ5eTmlpaV1lpUC5dtu62u0HKZCuylLlsDxx8M111C+114rftBK2iJ5p8GkXVJC+ezZcPvtIUUlbWbkSPjjH4ntsgsVd99NJBLBzIhEIlRUVBCLxcKOUERaIBaLUVFRsXxbLiujYq+9iD3zjK/RliwJO8T0nHN5+8fQAx4AAAvdSURBVNO3b1+XUfPmObfPPs6Bc5df7lxNjRs1apSLRCLOzFwkEnGjRo3KbAwikhErbMv/+IdzAwf67f38851btizjMQCTXA7k0mz9ZDxnL1vm3LBh/jMcNMi5hQsz+3oiEp6aGucuu8xv77//va/ZMqw1OVuDIdOprPRXI5oyBe65x39jEpHCtnSpv5DJXXf5ufHvuw9WXTVjL6fBkG1o4UIYPBjGjPHXOLjtNk3hJ1IM/vEPiMdhm21g3DiIRjP2UhoM2VZefBH69vXFdu1pCREpfO3bwx13wDXX+ItR7borTJ0adlTSlKlTYZdd4NFH/Wd3++0qskWKxQkn+Frt66+hXz9fw+UQFdqpamrgqqtg331hgw1g0iT/t4gUDzMYNgz+9S+YNs0n7mefDTsqSeeZZ/xnNGMGPPec/+x0xUeR4rLvvvDuu7Deev7vq6/2NV0OUKFdK5mEvfaCv/0N/vhHePNN2HTTsKMSkbD8/vcwebI/DXnggb47wsKFYUcltRYs8J/JQQf5z2jSJP+ZiUhx2mwzeOstX8NdeKGv6dJMy5xNRVdoJxIJotFo3csuP/ggbLutT9T33ecv0dylS9ihikjYNtrIf+k+5xy4+27Ybjt4442G84hkzArre8QI6NPHfybnnOM/o402CjtMEQlbly6+hrvvPl/Tbbutr/GcCy9vt3T0ZC79tHQE+6hRo1xpaakDfv0pLSlxo8C5/v2d++qrFj2fiBSRV15xrqzMjQJX2q5d3TxSWtqqGYjQrCNNajBvgxvVrZtzEya0+PlEpEh89ZWv7cCN2nJLV9qp00rn7dbk7KKadSTtZZfXWovKmTM1eEZEGjdvHtGePUnOm7fCTa25fLtmHWla2rzdsyeVupKniDRm2TK4916ip55KsoE+2y3N25p1pAlpL7s8Z46KbBFpWteuVM2f3+BNunx7BixYQFWaPpZV06dnORgRyTvt2sFf/0pVmoPK2cjbWSu0zWxfM/vMzL40swsauL2jmY0Jbn/bzKJt8sJLl/oZAw4/nLI0K1qXUBeR5kp7+XaA00/3AyibOFNY21cQ6NvG4bWZ0HK2c34dnnYa9OxJuuysvC0izZU2b3fqBOPH+1qxESuTs7NSaJtZO+B2YD+gF/AnM+tV726DgTnOuU2Bm4BrWv2CP/4I//ynn8C8rMzPGPDqq5Tvuy+lnTvXuasuoS4iLdHg5ds7daJ8553h3nv9VHPbbOOnmXvpJVi8uM59E4kE8Xi8we4QuSLrOXvxYr+uhg3z665fP38Riv32o3zEiBXXt/K2iLRAg3m7fXvK27WDAw7wtWI87mvHH3+sc7+Vztkt7dTdmh9gF+D5lP8vBC6sd5/ngV2Cv9sDs8H3IU/303fjjZ27917nbrzRXzL5oIOc22QT58ycA+e6dnXu0EOde/xx537+2TnXwGWXdQl1EWmhtHlkzhzn7rzTuQEDnOvQweehVVZxbpttnDvySOcuvdRF1lyzzoAclwODFOv/ZCxnb7KJcw884Nyttzp36aV+nWyzjV9H4NfZgAF+Hc6Z0/T6FhFppgbzyM8/+xrx0EN9zQi+htxkE19Tnn/+SufsrAyGNLPDgX2dcycG/x8L7OScOz3lPlOC+0wP/p8a3Gd2veeKA3GAvtD312E1HTrAFltAr16w9dZ+/sQdd/RXehMRybYFC+A//4FXX4WPP/Y/X39NiXOkZl3nXM5dXSUrOdvMT8nXq5f/+e1vYY89NLWqiIRj6VJ45x14+WWYMsXn7M8+o2TJkpXK2XlXhTrnKoAKgH69ezv+9S/o2tUnZw1oFJFc0aUL7L+//6m1ZAllG29MsogG8q2Qs5991ufs1VbzB0hERHJB+/aw667+p9ayZZRttBHJadNa/bTZGgw5A+iZ8v+GwbIG72Nm7YHVgR8afdbOnaFnT1h9dRXZIpL7OnSg/OqrV+grmIMyl7OjUVhrLRXZIpL72rWj/KqrVipnZ6vQfhfYzMw2MrNVgKOAp+rd5yngL8HfhwP/dtno1yIikkWxWIyKigoikUjYoTRGOVtEhJXP2VkptJ1zS4HT8YNnPgHGOuc+MrPLzWxgcLf7gG5m9iVwNrDCdFIiIoUgFovVXiRhcsihNEg5W0RkuZXJ2Vnro+2cGw+Mr7fskpS/FwN/zFY8IiKSnnK2iMjKK6orQ4qIiIiIZIsKbRERERGRDFChLSIiIiKSAVm5YE2mmNl84LOw41hJa+OvqJav8j1+yP/3oPjD19r3EHHOdW/rYHKVcnZOyPf4If/fQ77HD/n/HrKWs/PugjX1fOac6xd2ECvDzCbl83vI9/gh/9+D4g9fIbyHLFHODlm+xw/5/x7yPX7I//eQzfjVdUREREREJANUaIuIiIiIZEC+F9oVYQfQBvL9PeR7/JD/70Hxh68Q3kM2FMJ6yvf3kO/xQ/6/h3yPn/9v725CPZvjOI6/P2ZoGE+FNBllFrKxQCgPSUREWFIoGxtEFsJGStnJzmaGyMMkDyXJQ1FIGON5DJKUKxolMUqevhb3LAYz5dzp3N/5zbxfdbv/81/cPrd7+9zvPb/fOYf+v4dly9/1xZCSJEnSXPV+RluSJEmaJQdtSZIkaQLdDtpJLkjyWZIvktzaOs9YSe5Psi3Jx62zLEWSo5O8kuSTJFuS3Ng60xhJViV5O8kHQ/47W2daiiQrkryX5NnWWZYiyVdJPkryfpJ3WucZK8mhSZ5I8mmSrUlOa51pruzstuzs+ei5t3vvbFj+3u5yj3aSFcDnwHnAArAJuKKqPmkabIQkZwHbgYeq6vjWecZKsgZYU1XvJjkI2Axc1svPIEmA1VW1Pcm+wOvAjVX1ZuNooyS5GTgZOLiqLm6dZ6wkXwEnV1WXDz5I8iDwWlWtT7IfcEBV/dg619zY2e3Z2fPRc2/33tmw/L3d6xntU4EvqurLqvoN2Ahc2jjTKFX1KvBD6xxLVVXfVtW7w+ufga3AUW1T/X+1aPtwuO/w0dV/nUnWAhcB61tn2RslOQQ4C9gAUFW/OWTvkp3dmJ09D/Z2Wy16u9dB+yjg6x2OF+ioMPY0SY4BTgTeaptknGH57n1gG/BSVXWVH7gXuAX4q3WQ3VDAi0k2J7m2dZiR1gHfAw8My8Drk6xuHWqm7OwZsbOb6r23e+5saNDbvQ7amokkBwJPAjdV1U+t84xRVX9W1QnAWuDUJN0sBye5GNhWVZtbZ9lNZ1bVScCFwHXD8nwvVgInAfdV1YnAL0B3e4+1d7Gz29lDervnzoYGvd3roP0NcPQOx2uH97SMhn1yTwKPVNVTrfMs1bBs9ApwQessI5wBXDLsl9sInJPk4baRxquqb4bP24CnWdxi0IsFYGGHs2pPsFjg+i87ewbs7Oa67+3OOxsa9Havg/Ym4Ngk64aN7JcDzzTOtFcZLkzZAGytqnta5xkryRFJDh1e78/iRVqftk31/1XVbVW1tqqOYfH3/+WqurJxrFGSrB4uymJYujsf6OaODlX1HfB1kuOGt84FuriwrAE7uzE7u73ee7v3zoY2vb1yyi8+lar6I8n1wAvACuD+qtrSONYoSR4DzgYOT7IA3FFVG9qmGuUM4Crgo2HPHMDtVfVcw0xjrAEeHO6GsA/weFV1d6ulzh0JPL3495+VwKNV9XzbSKPdADwyDI9fAtc0zjNLdvYs2NnaXXtCZ8My93aXt/eTJEmS5q7XrSOSJEnSrDloS5IkSRNw0JYkSZIm4KAtSZIkTcBBW5IkSZqAg7YkSZI0AQdtSZIkaQIO2tK/JDklyYdJVg1PwtqS5PjWuSRJ/2Vna858YI20E0nuAlYB+wMLVXV340iSpF2wszVXDtrSTgyPZt0E/AqcXlV/No4kSdoFO1tz5dYRaecOAw4EDmLxLIkkab7sbM2SZ7SlnUjyDLARWAesqarrG0eSJO2Cna25Wtk6gDQ3Sa4Gfq+qR5OsAN5Ick5Vvdw6myTpn+xszZlntCVJkqQJuEdbkiRJmoCDtiRJkjQBB21JkiRpAg7akiRJ0gQctCVJkqQJOGhLkiRJE3DQliRJkibwNzQVVD5SDyIGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Note that this calculates the cummulative integral from 0.0\n",
    "\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 20\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "# Trapezoid\n",
    "xi_map = lambda a, b, x: (b - a) / 2.0 * x + (a + b) / 2.0\n",
    "\n",
    "# Trapezoid/Simpson's\n",
    "xi_methods = [[-1, 1], [-1, 0, 1]]\n",
    "w_methods = [[1, 1], [1/3, 4/3, 1/3]]\n",
    "titles = [\"Trapezoid Rule\", \"Simpson's Rule\"]\n",
    "fig, axes = plt.subplots(1, 2)\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "for k in range(2):\n",
    "    xi = xi_methods[k]\n",
    "    w = w_methods[k]\n",
    "    I_hat = numpy.zeros(x_hat.shape)\n",
    "    for j in range(len(xi)):\n",
    "        I_hat[0] += f(xi_map(x_hat[0], x_hat[1], xi[j])) * w[j]\n",
    "    I_hat[0] *= delta_x / 2.0\n",
    "    for i in range(1, num_partitions):\n",
    "        for j in range(len(xi)):\n",
    "            I_hat[i] += f(xi_map(x_hat[i], x_hat[i+1], xi[j])) * w[j]\n",
    "        I_hat[i] *= delta_x / 2.0\n",
    "        I_hat[i] += I_hat[i - 1]\n",
    "\n",
    "\n",
    "    axes[k].plot(x, I(x), 'r')\n",
    "    axes[k].plot(x_hat + delta_x, I_hat, 'ko')\n",
    "    axes[k].set_xlabel(\"x\")\n",
    "    axes[k].set_ylabel(\"$f(x)$\")\n",
    "    axes[k].set_title(\"Integral and Approximated Integral - %s\" % titles[k])\n",
    "    axes[k].set_xlim((0.0, 2.0 * numpy.pi))\n",
    "    axes[k].set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimal Quadrature Methods\n",
    "\n",
    "Can we determine $I_{\\Delta x}[f]$ to maximize degree for a given number of function evaluations (points)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Generalized Gaussian Quadrature\n",
    "\n",
    "Given $g(x) \\in P_N(x)$ with roots $\\{x_i\\}^N_{i=1}$ we have\n",
    "$$\n",
    "    \\int^1_{-1} w(x) x^i g(x) dx = 0 \\quad \\forall i < N,\n",
    "$$\n",
    "i.e. $g(x)$ is orthogonal to the $x^i$ with respect to the weight function $w(x)$.\n",
    "\n",
    "Recall something similar:\n",
    "$$\n",
    "    <x, y> = \\sum^N_{i=1} x_i y_i = ||x|| \\cdot ||y|| \\cos \\theta.\n",
    "$$\n",
    "If $<x, y> = 0$ then the vectors $x$ and $y$ are orthogonal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given the above $g(x)$ there then exists $\\{w_i\\}$ such that\n",
    "\n",
    "$$\\int^1_{-1} w(x) P_j(x) dx = \\sum^N_{i=1} w_i P_j(x_i) \\quad \\forall j \\leq 2 N - 1$$\n",
    "\n",
    "In other words, given a polynomial basis function and weight and orthogonality to all polynomials of order $i < N$ we can exactly integrate polynomials of order $2 N - 1$.  Choosing the correct weighting function and basis leads to a number of useful quadrature approaches:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Gauss-Legendre\n",
    "\n",
    "General Gauss-Legendre quadrature uses $w(x) = 1$ and $g(x) = \\ell_N(x)$ which can be shown to have weights\n",
    "\n",
    "$$w_i = \\frac{2}{(1-x_i^2)(P'_n(x_i))^2}$$\n",
    "\n",
    "and $x_i$ is the $i$th root of $\\ell_N$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The first few rules yield\n",
    "\n",
    "<table width=\"80%\">\n",
    "    <tr align=\"center\"><th>$$N$$</th> <th align=\"center\">$$x_i$$</th> <th align=\"center\"> $$w_i$$ </th></tr>\n",
    "    <tr align=\"center\"><td>$$1$$</td>           <td> $$0$$ </td> <td> $$2$$ </td> </tr>\n",
    "    <tr align=\"center\"><td>$$2$$</td>           <td> $$\\pm \\sqrt{\\frac{1}{3}}$$ </td> <td> $$1$$ </td> </tr>\n",
    "    <tr align=\"center\"><td rowspan=2>$$3$$</td> <td> $$0$$ </td> <td> $$8/9$$ </td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\sqrt{\\frac{3}{5}}$$ </td> <td> $$5/9$$</td> </tr>\n",
    "    <tr align=\"center\"><td rowspan=2>$$4$$</td> <td> $$\\pm \\sqrt{\\frac{3}{7} - \\frac{2}{7} \\sqrt{\\frac{6}{5}}}$$</td> <td> $$\\frac{18 + \\sqrt{30}}{36}$$ </td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\sqrt{\\frac{3}{7} + \\frac{2}{7} \\sqrt{\\frac{6}{5}}}$$</td> <td>$$\\frac{18 - \\sqrt{30}}{36}$$ </td> </tr>\n",
    "    <tr align=\"center\"><td rowspan=3>$$5$$</td> <td> $$0$$ </td> <td> $$\\frac{128}{225}$$ </td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\frac{1}{3} \\sqrt{5 - 2 \\sqrt{\\frac{10}{7}}}$$</td> <td> $$\\frac{322 + 13\\sqrt{70}}{900}$$</td> </tr>\n",
    "    <tr align=\"center\">                     <td> $$\\pm \\frac{1}{3} \\sqrt{5 + 2 \\sqrt{\\frac{10}{7}}}$$</td> <td> $$\\frac{322 - 13\\sqrt{70}}{900}$$</td> </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Example 2:  2-Point Gauss-Legendre Quadrature\n",
    "\n",
    "Let $N=2$ on $x \\in [-1,1]$\n",
    "\n",
    "$$I_2[f] = w_0 f(x_0) + w_1 f(x_1)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using undetermined coefficients again we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    &\\text{if}~f = 1:  &I[f] =  \\int^{1}_{-1} 1 dx = 2 & & I_2[1] &= w_0 + w_1\\\\\n",
    "    &\\text{if}~f = x:  &I[f] =  \\int^{1}_{-1} x dx = 0 & & I_2[x] &= w_0 x_0 + w_1 x_1\\\\\n",
    "    &\\text{if}~f = x^2:  &I[f] =  \\int^{1}_{-1} x^2 dx = \\frac{2}{3} & & I_2[x^2] &= w_0 x_0^2 + w_1 x_1^2\\\\\n",
    "    &\\text{if}~f = x^3:  &I[f] =  \\int^{1}_{-1} x^3 dx = 0 & & I_2[x^3] &= w_0 x_0^3 + w_1 x_1^3\\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    &w_0 + w_1 = 2\\\\\n",
    "    &w_0 x_0 + w_1 x_1 = 0\\\\\n",
    "    &w_0 x_0^2 + w_1 x_1^2 = \\frac{2}{3}\\\\\n",
    "    &w_0 x_0^3 + w_1 x_1^3 = 0\\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "Note that we need to solve for 4 unknowns $x_0$, $x_1$, $w_0$, and $w_1$.  Solving these equations leads to\n",
    "\n",
    "$$x_0 = -\\sqrt{\\frac{1}{3}}, x_1 = \\sqrt{\\frac{1}{3}}  \\quad \\text{and} \\quad w_0 = w_1 = 1 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compute Gauss-Legendre based quadrature and affine transforms\n",
    "f = lambda x: numpy.sin(x)\n",
    "I = lambda x: 1.0 - numpy.cos(x)\n",
    "x = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "num_partitions = 10\n",
    "x_hat = numpy.linspace(0.0, 2.0 * numpy.pi, num_partitions + 1)\n",
    "delta_x = x_hat[1] - x_hat[0]\n",
    "\n",
    "xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "xi_0 = -numpy.sqrt(1.0 / 3.0)\n",
    "xi_1 =  numpy.sqrt(1.0 / 3.0)\n",
    "I_hat = numpy.zeros(x_hat.shape)\n",
    "I_hat[0] = (f(xi_map(x_hat[0], x_hat[1], xi_0)) + f(xi_map(x_hat[0], x_hat[1], xi_1))) * delta_x / 2.0\n",
    "for i in range(1, num_partitions):\n",
    "    I_hat[i] = I_hat[i - 1] + (f(xi_map(x_hat[i], x_hat[i+1], xi_0)) + f(xi_map(x_hat[i], x_hat[i+1], xi_1))) * delta_x / 2.0\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, I(x), 'r')\n",
    "# Offset due to indexing above\n",
    "axes.plot(x_hat + delta_x, I_hat, 'ko')\n",
    "\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$f(x)$\")\n",
    "axes.set_title(\"Integral and Approximated Integral\")\n",
    "axes.set_xlim((0.0, 2.0 * numpy.pi))\n",
    "axes.set_ylim((-0.1, 2.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Other Quadrature Families\n",
    "\n",
    " - Clenshaw-Curtis (Gauss-Chebyshev): If $w(x) = \\frac{1}{\\sqrt{1 - x^2}}$ and $g(x)$ are Chebyshev polynomials then we know the roots of the polynomials to be $x_i = \\cos\\left(\\frac{2i-1}{2N} \\pi \\right)$ (the Chebyshev nodes) and we can derive that $w_i = \\frac{\\pi}{N}$.\n",
    " - Gauss-Hermite:  If $w(x) = e^{-x^2}$ and $g(x)$ are Hermite polynomials $H_i(x)$ then\n",
    "   $$w_i = \\frac{2^{N-1} N! \\sqrt{\\pi}}{N^2 (H_{N-1}(x_i))^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Example 3:\n",
    "\n",
    "If $f(x) = e^x$ look at the relative accuracy of midpoint, trapezoidal, Simpson's and 2-point Gauss-Legendre quadrature for a single interval $x \\in [-1,1]$.\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\text{Exact:} &I[f] = \\int^1_{-1} e^x = \\left . e^x \\right |^1_{-1} = e - \\frac{1}{e} \\approx 2.350402387 \\\\\n",
    "    \\text{Midpoint:} &I_2[f] = 2 e^0 = 2 \\\\\n",
    "    \\text{Trapezoid:} &I_2[f] = \\frac{2}{2} (e^{-1} + e^1) = e + \\frac{1}{e} = 3.08616127 \\\\\n",
    "    \\text{Simpson's:} &I_2[f] = \\frac{2}{6} e^{-1} + \\frac{4}{3} e^0 + \\frac{2}{6} e^1 = \\frac{4}{3} + \\frac{1}{3} (e^{-1} + e^1) \\approx 2.362053757 \\\\\n",
    "    \\text{Gauss-Legendre:} &I_2[f] = e^{-\\sqrt{\\frac{1}{3}}} + e^{\\sqrt{\\frac{1}{3}}} \\approx 2.342696088\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compute the error as a function of delta_x for each method\n",
    "f = lambda x: numpy.sin(numpy.pi * x)\n",
    "I = 2.0 / numpy.pi\n",
    "\n",
    "# num_partitions = range(50, 1000, 50)\n",
    "num_partitions = range(5, 50, 5)\n",
    "delta_x = numpy.empty(len(num_partitions))\n",
    "error_trap = numpy.empty(len(num_partitions))\n",
    "error_simpson = numpy.empty(len(num_partitions))\n",
    "error_2 = numpy.empty(len(num_partitions))\n",
    "error_3 = numpy.empty(len(num_partitions))\n",
    "error_4 = numpy.empty(len(num_partitions))\n",
    "\n",
    "for (j, N) in enumerate(num_partitions):\n",
    "    x_hat = numpy.linspace(0.0, 1.0, N)\n",
    "    delta_x[j] = x_hat[1] - x_hat[0]\n",
    "    \n",
    "    # Compute trapezoid\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        I_hat += (f(x_hat[i + 1]) + f(x_hat[i])) * delta_x[j] / 2.0\n",
    "    error_trap[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute simpson's    \n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        I_hat += delta_x[j] * (1.0 / 6.0 * (f(x_hat[i]) + f(x_hat[i+1])) + 2.0 / 3.0 * f(x_hat[i] + delta_x[j] / 2.0))\n",
    "    error_simpson[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute Gauss-Legendre 2-point\n",
    "    xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "    xi = [-numpy.sqrt(1.0 / 3.0), numpy.sqrt(1.0 / 3.0)]\n",
    "    w = [1.0, 1.0]\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        for k in range(len(xi)):\n",
    "            I_hat += f(xi_map(x_hat[i], x_hat[i+1], xi[k])) * w[k]\n",
    "    I_hat *= delta_x[j] / 2.0\n",
    "    error_2[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute Gauss-Legendre 3-point\n",
    "    xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "    xi = [-numpy.sqrt(3.0 / 5.0), 0.0, numpy.sqrt(3.0 / 5.0)]\n",
    "    w = [5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0]\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        for k in range(len(xi)):\n",
    "            I_hat += f(xi_map(x_hat[i], x_hat[i+1], xi[k])) * w[k]\n",
    "    I_hat *= delta_x[j] / 2.0\n",
    "    error_3[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "    # Compute Gauss-Legendre 4-point\n",
    "    xi_map = lambda a,b,xi : (b - a) / 2.0 * xi + (a + b) / 2.0\n",
    "    xi = [-numpy.sqrt(3.0 / 7.0 - 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0)), \n",
    "           numpy.sqrt(3.0 / 7.0 - 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0)),\n",
    "          -numpy.sqrt(3.0 / 7.0 + 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0)),\n",
    "           numpy.sqrt(3.0 / 7.0 + 2.0 / 7.0 * numpy.sqrt(6.0 / 5.0))]\n",
    "    w = [(18.0 + numpy.sqrt(30.0)) / 36.0, (18.0 + numpy.sqrt(30.0)) / 36.0,\n",
    "         (18.0 - numpy.sqrt(30.0)) / 36.0, (18.0 - numpy.sqrt(30.0)) / 36.0]\n",
    "    I_hat = 0.0\n",
    "    for i in range(0, N - 1):\n",
    "        for k in range(len(xi)):\n",
    "            I_hat += f(xi_map(x_hat[i], x_hat[i+1], xi[k])) * w[k]\n",
    "    I_hat *= delta_x[j] / 2.0\n",
    "    error_4[j] = numpy.abs(I_hat - I)\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "# axes.plot(delta_x, error)\n",
    "axes.loglog(delta_x, error_trap, 'o', label=\"Trapezoid\")\n",
    "axes.loglog(delta_x, error_simpson, 'o', label=\"Simpson's\")\n",
    "axes.loglog(delta_x, error_2, 'o', label=\"G-L 2-point\")\n",
    "axes.loglog(delta_x, error_3, 'o', label=\"G-L 3-point\")\n",
    "axes.loglog(delta_x, error_4, 'o', label=\"G-L 4-point\")\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error_trap[0], 2.0) * delta_x**2.0, 'r--', label=\"2nd Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[0], error_simpson[0], 4.0) * delta_x**4.0, 'g--', label=\"4th Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[1], error_3[1], 5) * delta_x**5, 'b--', label=\"5th Order\")\n",
    "axes.loglog(delta_x, order_C(delta_x[1], error_4[1], 7.0) * delta_x**7.0, 'g--', label=\"7th Order\")\n",
    "\n",
    "axes.legend(loc=2)\n",
    "axes.set_xlim((5e-3, delta_x[0]))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## SciPy Integration Routines\n",
    "\n",
    "SciPy has a number of integration routines that we have derived here including general purpose integrators that can control error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.integrate as integrate\n",
    "# integrate?"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.1"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}