{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.2"
  },
  "name": "",
  "signature": "sha256:447faf4dda91c9e17acd7d2066d62f7a06350bdbbd962dfeef81b6ace71c4e9a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "\\newcommand{\\dt}{\\Delta t}\n",
      "\\newcommand{\\udt}[1]{u^{({#1})}(T)}\n",
      "\\newcommand{\\Edt}[1]{E^{({#1})}}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the first in a series of posts on testing scientific software. For this to make sense, you'll need to have skimmed [the motivation and background](http://ianhawke.github.io/blog/close-enough.html)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We're starting from Euler's method, which is first order. That means that the difference between the exact solution $u(T)$ and the numerical approximation $\\udt{\\dt}$ should be proportional to $\\dt$: that is, \n",
      "\n",
      "$$\n",
      "\\begin{equation}\n",
      "  u(T) - \\udt{\\dt} = c_1 \\dt + {\\cal O}(\\dt^2) \\simeq c_1 \\dt^{s_e}\n",
      "\\end{equation}\n",
      "$$\n",
      "\n",
      "with $s_e=1$. By performing some simulations, we can measure (to some accuracy) the slope $s_m$. The problem is that the measured slope $s_m$ will be \"contaminated\" by the higher order terms, proportional to $\\dt^2$ and higher powers of $\\dt$. For sufficiently small $\\dt$ the measured value $s_m$ will be indistinguishable from the expected value $s_e=1$. However, we rarely get to use a sufficiently small $\\dt$. So how close does $s_m$ need to be to $s_e$ to be close enough?\n",
      "\n",
      "One argument that can be justified: $0.585 \\lesssim s_m \\lesssim 1.585$ is close enough when $s_e = 1$."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Error bars and Richardson extrapolation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The argument starts by assuming that what we care about is not the algorithm, but the solution $u(T)$. From the solution of *any* numerical algorithm we cannot compute $u(T)$ exactly, but only estimate it with some error. Given a set of numerical calculations - $\\udt{\\dt}, \\udt{2\\dt}$ for example - and a model of how the algorithm behaves, then both the exact solution *and* the error bars can be estimated. \n",
      "\n",
      "**Claim**: if two models estimate exact solutions that lie within each other's error bars, they are indistinguishable. Hence, if the model corresponding to the *measured* slope is consistent with that of the *expected* slope, the algorithm implemented is close enough."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To make this concrete, we need a model of the algorithm. The simplest models are\n",
      "\n",
      "$$\n",
      "\\begin{align}\n",
      "  u(T) - \\udt{\\dt} & = c_e \\dt^{s_e}, \\\\\n",
      "  u(T) - \\udt{\\dt} & = c_m \\dt^{s_m}.\n",
      "\\end{align}\n",
      "$$\n",
      "\n",
      "Using [Richardson extrapolation](http://en.wikipedia.org/wiki/Richardson_extrapolation) we can estimate both the solution $u(T)$ and also the error $\\Edt{\\dt} = u(T) - \\udt{\\dt}$. For Richardson extrapolation we perform two numerical calculations $\\udt{\\dt}$ and $\\udt{2\\dt}$ and combine them as\n",
      "\n",
      "$$\n",
      "\\begin{align}\n",
      "  \\udt{e} & = \\frac{2^{s_e} \\udt{\\dt} - \\udt{2\\dt}}{2^{s_e} - 1}, & \\Edt{e, \\dt} & = \\frac{\\udt{\\dt} - \\udt{2\\dt}}{2^{s_e} - 1}, \\\\\n",
      "  \\udt{m} & = \\frac{2^{s_m} \\udt{\\dt} - \\udt{2\\dt}}{2^{s_m} - 1}, & \\Edt{m, \\dt} & = \\frac{\\udt{\\dt} - \\udt{2\\dt}}{2^{s_m} - 1}.\n",
      "\\end{align}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The models are consistent if\n",
      "\n",
      "$$\n",
      "\\begin{equation}\n",
      "  \\udt{e} \\in \\udt{m} \\pm \\Edt{m, \\dt} \\quad \\text{and} \\quad \\udt{m} \\in \\udt{e} \\pm \\Edt{e, \\dt}.\n",
      "\\end{equation}\n",
      "$$\n",
      "\n",
      "That is, if both models predict the same exact solution, within their own error bounds, then we can't tell the difference between them, and they're close enough."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Checking the MOOC solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can use this criteria to check the result from [the full phugoid model notebook](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_03_PhugoidFullModel.ipynb).\n",
      "\n",
      "Let's take the code from there directly."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "from math import sin, cos, log, ceil\n",
      "import numpy\n",
      "from matplotlib import pyplot\n",
      "%matplotlib inline\n",
      "from matplotlib import rcParams\n",
      "rcParams['font.family'] = 'serif'\n",
      "rcParams['font.size'] = 16"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# model parameters:\n",
      "g = 9.8      # gravity in m s^{-2}\n",
      "v_t = 30.0   # trim velocity in m s^{-1}   \n",
      "C_D = 1/40.  # drag coefficient --- or D/L if C_L=1\n",
      "C_L = 1.0    # for convenience, use C_L = 1\n",
      "\n",
      "### set initial conditions ###\n",
      "v0 = v_t     # start at the trim velocity (or add a delta)\n",
      "theta0 = 0.0 # initial angle of trajectory\n",
      "x0 = 0.0     # horizotal position is arbitrary\n",
      "y0 = 1000.0  # initial altitude"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def f(u):\n",
      "    \"\"\"Returns the right-hand side of the phugoid system of equations.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    u : array of float\n",
      "        array containing the solution at time n.\n",
      "        \n",
      "    Returns\n",
      "    -------\n",
      "    dudt : array of float\n",
      "        array containing the RHS given u.\n",
      "    \"\"\"\n",
      "    \n",
      "    v = u[0]\n",
      "    theta = u[1]\n",
      "    x = u[2]\n",
      "    y = u[3]\n",
      "    return numpy.array([-g*sin(theta) - C_D/C_L*g/v_t**2*v**2,\n",
      "                      -g*cos(theta)/v + g/v_t**2*v,\n",
      "                      v*cos(theta),\n",
      "                      v*sin(theta)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def euler_step(u, f, dt):\n",
      "    \"\"\"Returns the solution at the next time-step using Euler's method.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    u : array of float\n",
      "        solution at the previous time-step.\n",
      "    f : function\n",
      "        function to compute the right hand-side of the system of equation.\n",
      "    dt : float\n",
      "        time-increment.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    u_n_plus_1 : array of float\n",
      "        approximate solution at the next time step.\n",
      "    \"\"\"\n",
      "    \n",
      "    return u + dt * f(u)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we're going to compute the result at a specific time and do the convergence test. We'll use an end time of `T=100` as in the original notebook, and convergence test on `v` as an example."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T = 100.0\n",
      "dt_values = numpy.array([0.001*2**i for i in range(3)])\n",
      "v_values = numpy.zeros_like(dt_values)\n",
      "for i, dt in enumerate(dt_values):\n",
      "    N = int(T/dt)+1\n",
      "    t = numpy.linspace(0.0, T, N)\n",
      "    u = numpy.empty((N, 4))\n",
      "    u[0] = numpy.array([v0, theta0, x0, y0])\n",
      "    for n in range(N-1):\n",
      "        u[n+1] = euler_step(u[n], f, dt)\n",
      "    v_values[i] = u[-1,0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We now have three calculated values of `v` using parameter values `dt=0.001, 0.002, 0.004`. That is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i, dt in enumerate(dt_values):\n",
      "    print(\"Value of v({}) = {:.7g} when dt={}.\".format(T, v_values[i], dt))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Value of v(100.0) = 29.86798 when dt=0.001.\n",
        "Value of v(100.0) = 29.86667 when dt=0.002.\n",
        "Value of v(100.0) = 29.864 when dt=0.004.\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As in the original notebook, we can measure the order of convergence using\n",
      "\n",
      "$$\n",
      "\\begin{equation}\n",
      "  s_m = \\log_2 \\left| \\frac{\\udt{4\\dt} - \\udt{2\\dt}}{\\udt{2\\dt} - \\udt{\\dt}} \\right|.\n",
      "\\end{equation}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s_m = log(abs((v_values[2]-v_values[1])/(v_values[1]-v_values[0]))) / log(2.0)\n",
      "print(\"Measured convergence rate is {:.4g}.\".format(s_m))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Measured convergence rate is 1.023.\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As in the original notebook, this is close to one, but not exactly one. So we now use the two most accurate data points with our two models, with *expected* convergence rate $s_e = 1$ and *measured* convergence rate $s_m \\simeq 1.023$. We compute our exact solutions and error intervals."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "v_exact_e = (2.0*v_values[0] - v_values[1])\n",
      "v_exact_m = (2.0**s_m*v_values[0] - v_values[1]) / (2.0**s_m - 1.0)\n",
      "error_e = abs(v_exact_e - v_values[0])\n",
      "error_m = abs(v_exact_m - v_values[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Then we check if they are consistent:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"Does the expected solution {:.7g} lie within the interval \"\n",
      "      \"[{:.7g}, {:.7g}]\\nfrom the measured model? {}\".format(\n",
      "        v_exact_e, v_exact_m-error_m, v_exact_m+error_m, \n",
      "        v_exact_m-error_m <= v_exact_e <= v_exact_m+error_m))\n",
      "print(\"Does the measured solution {:.7g} lie within the interval \"\n",
      "      \"[{:.7g}, {:.7g}]\\nfrom the exact model? {}\".format(\n",
      "        v_exact_m, v_exact_e-error_e, v_exact_e+error_e, \n",
      "        v_exact_e-error_e <= v_exact_m <= v_exact_e+error_e))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Does the expected solution 29.8693 lie within the interval [29.86798, 29.87053]\n",
        "from the measured model? True\n",
        "Does the measured solution 29.86925 lie within the interval [29.86798, 29.87061]\n",
        "from the exact model? True\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So the two models are indistinguishable: $s_m \\simeq 1.023$ is *close enough*."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Finding the limits"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So, how different can $s_m$ be from one and still be close enough? First, let's do this experimentally, by plotting the intervals as we vary $s_m$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s_m_values = numpy.linspace(0.4, 2.0)\n",
      "v_exact_m = numpy.zeros_like(s_m_values)\n",
      "interval_lower = numpy.zeros_like(s_m_values)\n",
      "interval_upper = numpy.zeros_like(s_m_values)\n",
      "for i, s_m in enumerate(s_m_values):\n",
      "    v_exact_m[i] = (2.0**s_m*v_values[0] - v_values[1]) / (2.0**s_m - 1.0)\n",
      "    error_m = abs(v_exact_m[i] - v_values[0])\n",
      "    interval_lower[i] = v_exact_m[i]-error_m\n",
      "    interval_upper[i] = v_exact_m[i]+error_m\n",
      "pyplot.figure(figsize=(8,6))\n",
      "pyplot.plot(s_m_values, v_exact_e*numpy.ones_like(s_m_values), 'g-', linewidth=2, \n",
      "            label='Exact model')\n",
      "pyplot.plot(s_m_values, v_exact_m, 'r-', linewidth=2, label='Measured model')\n",
      "pyplot.fill_between(s_m_values, v_exact_e-error_e, v_exact_e+error_e, facecolor='green', \n",
      "                    alpha=0.5, edgecolor='None', label='Exact model interval')\n",
      "pyplot.fill_between(s_m_values, interval_lower, interval_upper, facecolor='red', \n",
      "                    alpha=0.5, edgecolor='None', label='Measured model interval')\n",
      "pyplot.legend()\n",
      "pyplot.xlabel(r\"$s_m$\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAGYCAYAAAAnchq9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl81NW9//HXmWxk35NJCAnIIqCAVepuC1qo3tZaXO51\nqaKtrdpbq7/WSq29F5dqe71e27q0akXAaq1arbsVF0Bt3RdAQXbClhVCEpKQ9fP745uM2TNAYCbJ\n+/l4zCPkzPmeOTMhmfec7znn68wMERERkfZ8oe6AiIiIhB8FBBEREelCAUFERES6UEAQERGRLhQQ\nREREpAsFBBEREelCAWEIcs79m3Ou2Dk3P9R9ERGR8KSA0I5zLtY59wPn3FvOuc+ccyudc284507f\nizYudc4tc86tcs5tcs7d65xL6abejNa2P229LXXOfa2HNs9qrfuBc269c+5959x39vH53QPcAGQB\n2gRDRES6NeQCQuub9uwe7j4D+D1wnZkdZmYTgaeBZ5xz3w2i7euBu4EfmtkEYDIwBXjROedrV+9o\n4EXgTTM73MwOB/4JvOCc+3KnNv8fcDNwoZlNBQ4F1gAn79UT95wFVADH78OxIiIyhAy5gID3qbmn\nT84GPGFmbwYKzO4AtgBX99aocy4euB540sz+2XpsFfBfwLHAOe2qnw1EALe1K7sNiALObNfmSODX\nwJVmVtjaZhNwDXBP70+zW38xs1+2tiEiItKjoRgQevM4cEk35cVAl9MEnRwGDAM+61S+ovXrt9qV\nNbZ+jWpX1vbv9m/eF7Z+/3r7Bs2syMw+bF/mnPtK66mRtc65Dc65h51z/k7HtfTxHERERAAFhA7M\n09y+zDkXARwCLOnj8JrWr51f07Y35UPbld0HFAE3OecinXNReKcRtgH3tqt3PLAV+Dfn3OLWORH/\n6ny6wzl3AvAq8LSZjQXGAbHAa8656D76LSIi0sVQDQhuL+qeifdme1Mf9dYBu4EjO5Uf0fo1qa3A\nzDYDXwOmAjuBcuBLwAwz29bu2BGttxuAC1rnRNwO3Oec+0W7ev8DbDOz21vbbwJ+AUwAzgvqWYqI\niLQzqAOCc26qc+7j9jcgF++Te/vy+3s4Pgf4LXCpma3r7bHMrB7vjfx059yZzpML3ArUAnXt2j0R\n+BfwHJAKpAGPAe8456a3a3YYXji53sy2tz7OU8AzwC+cc8Occ3HAccA7nbq0Bu/0xLQ+XiYREZEu\n3FC73LNzbiMw18we6qNeCrAYmGdmd+9F+98Hvgck451GuA3vU/8mM/tma533gBQzG9fp2I+A2NYV\nELQGmslAjpmVtqv3S7wRjalAKbAZbxSi/egDQCaw1Mwu6KafLcACM+tzdYaIiAw9kaHuQDhqDQeL\ngPl7Ew4AzOxPwJ/ateUDCoC/tKs2CS98dLYWOMc5F2VmjcDneMskO4/0tM2T8OGdomgBnjWzS/em\nryIiIj0Z1KcY9oVzLhkvHDxsZne2K3+2U71o51x6p7Kpracl2puGF8T+3K6sBC80dFYAVLSGA4C2\nx5zcqd7heKctPjOzWrzTFUd0qoNz7qfOuX/v5nFERER6pYDQTrtwUAvscM59p+2Gt5dBe88BW5xz\n+e3KfgT8b+vKh7Z9DO4C5pjZlnb17gQmOOd+0O6xLwCObq3f5jHgfWCucy6ptd6JeBse3WJmbfMa\nrgUOd85d1q69acBPgfd6e8q93CciIkOYTjF0NBv4Mt6GSV/pdF/nyRpFeOf969qVvQpcBax1zlXg\nBY3/ap1Y+EVDZnc457YDVznnftxaXAN818wWtKvX4pw7FW+VwgrnXB1QD/ynmc1rV69tcuOtzrnr\nWvtVAvybmW1qq9c6obFtnwYDzm4NEs+a2VV9vTgiIjJ0DLlJiiIiItI3nWIQERGRLhQQREREpItB\nPQfBOafzJyIiMqSYWb9MQB/0IwhmNmhvc+fODXkf9Pz03PT8Bt9Nz2/g3vrToA8IIiIisvcUEERE\nRKQLBYQBbNq0aaHuwgE1mJ/fYH5uoOc30On5CQzyfRCcczaYn5+IiEh7zjlMkxRFRETkQFFAEBER\nkS4UEERERKQLBQQRERHpYlDvpCgiwXNOV/8WCXcHc+K9AoKIBGjVj0j4OtghXqcYREREpAsFBBER\nEelCAUFERES6UEAQERGRLhQQREREpAsFBBEZ8iZPnkxWVhY+n4+YmBhycnK6vcXExLBw4cJQd3dQ\n8vv9pKSk4PP59vk1fuihh8jJySE6OhqfT29v+0uvoIgMecuXL+f9998H4IQTTqCoqKjb2/HHHx+W\n+0XccMMN+Hw+3njjjVB3ZZ8VFxfz+9//Htj35XwXXXRRWP+cBhoFBBERtAdEONDPILxooyQRkSC9\n9tprYT10rTdY6U/h+z9dRCRMXHzxxdx4440dwsEZZ5xBYmIiPp+PxMREZs+eDcCPfvQjUlNTiYqK\nIjc3lxUrVgCwZcsWfvGLX3DUUUeRl5dHSkoKhx12GL/5zW9oamrq9nE//vhjzjrrrMAciLFjx/Kt\nb32LBQsW0NzcTHNzM36/n//7v/8D4MwzzwzU/fWvf93rc5o8eTLp6en4fD5uvPFG5s+fz5QpU0hI\nSOCII47gmWeeAeDee+9l8uTJJCcnc/TRR/Ovf/2r2/bWrl3LRRddRF5eHn6/n4KCAq644gqKioq6\n1G1sbOSGG25g5MiRpKenM3HiRG6//XZaWlp67O/q1as5//zzyc7OJj09nZEjR/LDH/6QkpKSXp+n\n7AczG7Q37+mJSDCG+u/Lxo0bzTln06dP73Lf7Nmz7cYbb+xSvnr1aouPj7fRo0dbdXV1oPy0006z\nO+64o0Pd+fPnW1xcnL3wwguBsmeeecbi4uLsvPPO69L2Sy+9ZDExMXb++edbZWWlmZlt3rzZpk2b\nZs45KywsDNS94YYbzDlnS5cu3avnvGTJEnPO2ZQpU+x3v/udNTc3W0VFhX35y1+2iIgI+8lPfmK/\n/e1vrbm52Xbu3GlHHXWUpaamWlVVVYd2PvjgA0tKSrKZM2daSUmJmXmv59SpUy0nJ8c2btzYof65\n555rPp/PHnjgATMzq6urs1tuucUmTpxozjlbuHBht+2feuqpVlpaamZmK1assEMPPdTy8/OtrKys\nQ/2vfvWr5vP59uq1GAiC+R1trdM/76H91VA43ob6HzyRvTHUf1/aAkJ0dLT5/f4Ot9jY2G4DgpnZ\nvffea845u/DCC83M7Pe//73NnDmzS71nn33W5s6d26X8mmuuMeecrV69OlBWU1NjWVlZlp2dbXv2\n7OlQv7Cw0Hw+X4eAMHfu3H0KCIsXLzbnnB155JEdyufNm9dt+QMPPGDOOXvssccCZS0tLTZp0iSL\njIy0zZs3d6j/9ttvm3POTj311EDZa6+9Zs45+/d///cu/TnppJO6BISWlhY7/PDDLT4+vksQeOml\nl8w5Z1dddVWHcgWE/nkP1RwEEdln7saDP1Pc5h7Y8+wnnHACr7/+eoeySy65pO1DRxeXXXYZL7zw\nAg8//DCjRo3iwQcf5L333utS7/TTT+f000/vUj5hwgTAW0kxbtw4AF555RXKyso4//zziYmJ6VA/\nPz+f1157jaysrH16ft057rjjOnyfl5fXa/mWLVsCZcuWLePTTz/ly1/+MiNGjOhQ/9hjjyUrK4tX\nXnmFnTt3kpaWxlNPPQXAqaee2qUfM2fO5K233upQtmzZMj777DO++tWvkpGR0eG+o48+GoDnn3+e\n3/3ud0E/XwlO2AQE55wPuAa4GfiBmWmxsYiEjd6Wzc2bN49JkyZx8803M2/ePHJycrrUaWpqYuHC\nhfz5z3+msLCQuro6nHPU1dUBUFtbG6i7evVqgC5vuG2mTZu2H8+kq85vvNHR0d2Wt4WVmpqaQNma\nNWsAGD58eLdtDx8+nNLSUtasWcOxxx7LunXrAMjNze1St7vXre21ePfdd7u9Pz4+noqKiu6fmOyX\nPgOCcy4L+C1wVGvRCuBqM9sWxLFRwH8DZwNNQBVwrZn9s1O9fOAhIB6IAnr8iOCcywRuAo4FIoAk\n4F28UFHZV59EpP8c6E/z4WL+/Pm93p+Zmcnhhx/O66+/zkMPPcTFF1/cJVB897vf5eGHH+aWW27h\nyiuvJCEhAYCFCxdyySWXdNtufX19/zyBPvS0MuNArtgI9rm1vY6nnnpqYPRBDo5ef/rOuWjgFbwg\nMbH1VgMsds7FB9H+XcA5wIlmNgl4EFjknJvSqd7VwB/xRhB6608G8Daw3sy+ZGaTgVnAt4DkIPoj\nItLv7rnnHhoaGrjyyitZunQpt912W4f7Kysrefjhhxk/fjzXXXddIBxA90sT2047bNvW/eew7du3\nd/gUH0rjx48HYOvWrd3ev23bNnw+X+D0ydixYwPlnW3fvr3H9tuf1mhv7dq13Z7Skf3XVzycDUwC\n5phZi5m1AHOAQ4ArejvQOXco8H3gN2a2A8DM5gEbgVs6Vb/GzB4D+jqheQtQZGa3txWY2cfAN4Cy\nPo4VEelRMDvvfe973+Pmm2/uULZq1SpuvfVWHnnkEf73f/+XSZMm8d///d989NFHgTpRUVFERER0\n+xgbN27sUjZjxgyys7N5+eWXO5x6APjss8/Iy8tj1apVgbL4eO/zWmNjY6BP11zT6+etfjN58mQm\nTZrERx99xKZNmzrc9/bbb1NaWsrXv/510tLSAG8pJsBzzz3Xpa1XX321x/Y//vhj1q9f3+E+M+M/\n/uM/+Otf/9pPz0ba6ysgnAUUmtmmtgIzKwFWtt7Xm1l4b/iLO5UvBmY65+Latdnz4tdWzrlY4ALg\n+c73mdnrZlbX7YE7dvTVtIhI4JN8T5MRAZqbmzus1W9oaOCCCy7gjjvuYMSIEURHR/PII4/g8/m4\n4IILAvML4uLimDVrFqtWreK2224LvJEvXryYO++8s8vjDhs2jAULFlBfX8/ll19OVVUVAOvXr+fi\niy/m7LPPZurUqYH6U6Z4g7KffPIJAI888ghLly7d6+e+r+Xz588nISGBSy+9NLAvwaZNm7jyyivJ\nycnhD3/4Q6Du9OnTOffcc3n55ZeZN28eLS0t7Nmzh1//+tesXbu21/YvueQSCgsLAW9U5oc//CEV\nFRVce+21Qfdd9kJvSxyA7cDr3ZQ/C+zu49i/AI2A61T+E6AFmNrNMdNa77uom/uOab3vcuBO4BNg\nDbAAyO+hD2avvNLnshARGdrLHAsKCiwmJsZ8Pp9FRETYsGHDur1FRkYGljted911lp6ebhEREZaT\nk2Mff/yxmZllZ2cH2kpLS7Orr77azMx2795tc+bMsdGjR1t8fLyNHDnSzj33XLv++uvNOWfJycl2\n+OGHd+jXRx99ZGeeeaZlZ2eb3++38ePH2y233GL19fVdnsNPf/pTy83NNb/fb8ccc4x9+OGHvT7n\nU045xdLS0szn81lCQoKNGDHCzMzOOeecoMpzcnI67P2wdu1au/DCCy03N9eys7MtPz/fLr/8ctu+\nfXuXx25sbLQbbrjBCgoKLC0tzcaMGWO//OUv7f777w+8Fjk5OdbQ0NCl/ZycHPP7/TZ69Gi74oor\nbOvWrYE6CxcuNL/fb9HR0ebz+czv99uVV17Z6+swkATzO0o/LnN01kvKcs41AC+Z2Rmdyh8Gzgdi\nzazbmSbOuUXAMWaW3Kn8UuB+4DQze7nTfdOA14GLzeyhTvedDTyONwfip8CfgFTgCeBQ4AgzK+90\njNntt8P/+38QxtujioQD55w+dYmEsWB+R1vr9Mv6477eNcPpr8Ww1q8rzOz+1rC0E7gKyAV+2O1R\n1dXQuqxGREREgtPXMsdyILGb8iSgpqfRg3bHxjvnnHWMPEmtX/d2ckB169dPOpV/hreEcirduGHJ\nEvj0Uzj8cKZNm9bv64dFRERCZcmSJSxZsuSAtN1XQFgOjO+mfBTefgi9WQacC4wANnc6thFvouPe\naJuy292oR0sP5dwwbRpERMBPfgLxwazMFBERGRg6f/C98cYb+63tvk4xPAUUOOcK2gqcc9l4oeHJ\n9hWdc9mu4xqev+Odopjeqc3pwCIzq2UvmNkavEmJkzvdNQaIBt7v8eDmZli+fG8eTkREZEjrKyAs\nwBsp+B/nXETrdsi/ATbgbWwEgHPuBLwVD/e0lbW+od8PXOecS2+tdwneCML1fTxuTxMsrgWOcc6d\n09peNHAbsA24u9cWP/64j4cUERGRNr0GBDNrBGYAzXinBFYCCcDJnUYAqoEKvJDQ3pV4qwz+6Zxb\nAVwKzDSzDh/nnXMnO+c2Ao/ijTrc7pzb6Jw7r1N/nsXbf+Hnzrk1wGqgFjjeWjdj6lFpKfSwK5mI\niIh01Osyx4HOOWc2d+4XBVOnwje/GboOiYQxLXMUCW/htsxxcPn0U2jdwUxERER6NrQCwp490G7/\nchEREene0AoIAO0uoCIiIiLdG3oBobAQKipC3QsREZGwNvQCgpmWPIqIiPRh6AUEgGXLvKAgIiIi\n3RqaAaGyEtavD3UvRCQMNDc34/f7ycrKwufz4fP5ePHFF/s87plnnsHn8xEZGUlOTg5nn332Qejt\nwHf77beTk5NDREQEo0aNCnV39tqWLVvIyckhMTERn8/H0qVL96mdW265Jexfh6EZEECnGUQEgIiI\nCIqLi3n//S92a7/pppv6PK6tTn5+PkVFRfztb387YH0cTK655hqKiorIz8+n4+78A8OIESMoKiri\nmmuuAdjn53D99deH/eswdAPC6tVQVxfqXohImGjbgOZLX/oS7733Hi+//HKPdZ9//vmD1a1Ba6Bv\nytVf/Q/n12HoBoSmJl3ASUS6+OUvfwn0flW8m2++OVBPZLAaugEBdJpBRLqYNWsWhx12GO+88w6v\nvPJKl/tfeukl9uzZw6xZs3ptZ/v27Vx22WXk5eWRnp7O8OHDueiii1jfaf5TfX09d9xxB9OnT6eg\noID09HQKCgq4/PLLKSsr69JuRUUF11xzDWPHjiU3N5eCggK+/vWv84c//IGGhgYAvva1r5Geno7P\n5+sQdNrO//t8PqZP/+JCu+vWretwXv3VV1/lxhtvZMKECcTHx+Pz+Vi4cGGg/qJFizjllFNIT08n\nLS2NSZMm8etf/5qmpqYu/f3www+ZMWMGSUlJ5Ofn881vfpPPP/+819eu8+uTk5NDSkoKPp+PBQsW\ncNtttzFhwgQSEhI4/vjjeeutt2hpaeFXv/oV48ePJyUlhZNPPplVPWyM9+GHHzJr1ixycnLw+/2M\nHTuWOXPmUFlZ2aVudXU1V111Fbm5uWRlZTFlyhQWLFjQa5/ff/99vvWtb5GRkUF6ejrjxo3j5z//\nObt37w76eYcFMxu0N8Bs7tzeb9u2mYiYeX8Ohq6NGzeac87MzB599FFzztmJJ57Ypd6xxx5rjz32\nmJmZOeds1KhRXeps2rTJcnJybOrUqbZhw4ZA2XHHHWcpKSm2evXqLo/7q1/9yhobG83M7JNPPrHx\n48fb2LFjraampkPbX/va12zy5Mm2ZcsWMzPbs2ePXXfddeacs8LCwkC9JUuWmHPObrzxxi79c87Z\n9OnTu5TfcMMN5pyzo48+2hYsWGBNTU22a9cuGzNmjC1cuNDMzO69917z+Xx23XXX2Z49e6ylpcWe\neeYZS0hIsG984xsd2vvwww8tLi7Ojj32WNu6dauZmX366ac2bdo0S09P7/a168mCBQvMOWdHHXWU\nPf7442ZmtmXLFhs5cqQlJCTYlVdeGfi5bN682QoKCmzs2LHW3NzcoZ3nn3/eoqOj7aKLLrKqqioz\nM1u2bJmNGjXKJkyYYBUVFYG6jY2NduKJJ1psbKw999xzZmZWXV1tV1xxhU2cONGcc7Z06dJu2589\ne3ag/TfffNP8fr8deeSRVldX16F+QUFB0K9DML+jrXX65z20vxoKx1tQAeFvf+vzBRcZCvYpIHgL\nhg/u7QBpHxBaWlrs0EMPNeecvfbaa4E6L7/8so0fPz7wfU8B4Rvf+Ib5fD777LPPOpSvWrXKnHN2\nxhlnBMq2b9/e5Y3VzHujcc7ZfffdFyjbtWuXOefspz/9aZf6U6dOte3btwe+X7x48V4HhLlz55pz\nzr773e92KF+0aJGtXr3atm7datHR0Xb00Ud3OXbOnDnmnLO///3vgbKTTjrJfD6frVy5skPd1157\nrcfXrifz588355ydeeaZHcr/67/+q9vyX/7yl+acs3fffTdQVlNTY1lZWZaSkmK1tbUd6reFwssv\nvzxQNm/ePHPO2bXXXtuhbnNzs40cObJLQKipqbHMzEwbMWKE1dfXdzjmj3/8oznn7Le//W2H8nAO\nCEP7FAPAypVQVRXqXohIGHHOcd111wEd5yLcdNNN/OIXv+j12J07d/Liiy+Sn5/PxIkTO9w3fvx4\nEhMTWbRoEc3NzQDk5OR0O+lxwoQJACxvN1cqLi6O5ORk5s+fz8MPP8yePXsC973//vvk5OTs5TPt\n3owZM7p8P27cOJ544gkaGxs57bTTuhxz9NFHA19M4CwrK+Ott94iLy8v8FzaTJs2jYiIiH3q23HH\nHdfh+7y8vF7Lt2zZEihbtGgRZWVlzJw5k9jY2A71204ZPfbYY4Gyp556CoBTTz21Q12fz8fJJ5/c\npW+LFi2ivLycU045hejo6A73dX59BoLIUHcg5Jqb4d13odMvhIgEwcJ3Bvb++s53vsNNN93Em2++\nyZIlS2hubqaoqIgLLrig1+PWrl0LeHMQenrDjoyMpKysDL/fD8Crr77KPffcw8qVK9m1axc+ny8Q\nIGprawPHRUVF8fjjj3PxxRdz0UUXccUVVzBjxgzOO+88Zs2aRWRk//xJb+tXZ6tXrwbgd7/7Hffd\nd1+H+5qbm0lISKC8vBwgMNciNze3Szs+n6/Hx+hLRkZGh+/b3og7l8fExABQU1MTKFuzZg0Aw4cP\n79JuTEwMqamp7Nq1i+LiYvx+P+vWrcM51+1z6K6s7fX529/+xj/+8Y8O95kZCQkJ7Nixo8/nGC4U\nEAA+/BC++lXolPhEZOiKiIjg5z//OZdddhk33ngjZsYvfvELfL7eB17b1rQfdthhfBTExeHmzZvH\n97//fc4++2xef/31wJtXYWFhtxvozJgxg8LCQl566SWefPJJnnrqKZ5++mkmTZrE66+/Tnp6eq+P\n19LS0mefenqObc/tpptu4sc//nGf7YA3ybA/9dS3vn4u+yPY59D2+nz/+9/njjvuOGD9OVh0igG8\ny0BrRYOIdHLxxReTl5fH0qVL2bRpE7Nnz+7zmHHjxuGcY+vWrd3ev2XLFpYsWRL4/p577sE5x913\n393hk631MDrT3NxMZGQkp59+OgsWLKCoqIhLLrmEFStWcNdddwXqRUVFAdDY2Njh+JKSkj6fQ0/G\njx8PwObNm7u9/7333guMoIwePRrwRlI6a2lp2a9+7Ku2/nf3s6mrq6OiooLk5OTA6MaYMWMwM7Zt\n29alfnfPq6/XZ/ny5R1OGYU7BYQ277wzqIdLRWTvRUVFMWfOHJxzzJkzJ6gh/JSUFL7xjW9QXl7O\n4sWLu9z/4x//mDvvvDPwfUxMTLdhYOPGjV3KNm3aRHJycof68fHxXH311QDs2rUrUN52Dr6wsLBD\nG/u6NTDAOeecQ0xMDE899VSXkYjS0lK+8pWvBJYWZmZmctJJJ1FaWsp7773XpQ/dLYk80GbOnElW\nVhYvv/xylyWHTz/9NADnnXdeoOyss84C4LnnnutQt7m5udufbVv7ixYtoqrT3Lb6+nq+/vWvdwiH\n4U4BoU1FBezF2lwRGRr+8z//k8bGRq644oqgj7n77rvJzc3lRz/6EZ9++ikAe/bs4eabb+aNN97g\nlltuCdS98MILAfjRj35EReul6NetW8dPf/pToOtIQm1tLT/72c8C59arq6u56667iIqK6vDmVlBQ\nwFFHHcXzzz/PJ598AsDKlSsDW0L3NELR2305OTnceeedFBYWcuWVV1JdXQ3Ahg0bOOusszjllFM4\n/fTTA/V/97vfERcXx1VXXRX4FL5y5Up+/vOfk5CQ0Gsf9rZvwZQPGzaMBx98kPr6ei677LLAm/jy\n5cu5/vrrmTBhArfeemug/uzZsznxxBNZsGBBYHLh7t27ueqqqwITRLtrv6Ghge9+97uBfSxKSko4\n77zzGD58OJdeemnQfQ+5/loOEY43glnm2P72wAO9rx8RGcQYwvsgTJo0yTIzM83n85nf77fTTjut\n1/pf//rXze/3m8/ns8jISPP7/faDH/ygQ52ioiK74oorbMSIEZadnW2jRo2yCy64wD7//PMu7d1/\n//02ZcoUS0hIsLy8PDv55JPt/vvvN+ecxcbGWk5Ojq1Zs8YaGhrsT3/6k5166qk2cuRIy8nJsfz8\nfPv2t79t77zzTpd2t2zZYrNmzbKMjAwbPny4nX/++VZaWmrOOYuJiTG/32+PPvqoVVVVWXZ2tiUk\nJJjP57O0tDTz+/321ltvdfv8X331VZs5c6alp6fb8OHD7bDDDrNbb721yxp/M28vhBkzZlhiYqIN\nHz7cTjnlFHv33Xdt5MiRgdfu1ltv7fX1HjdunCUnJ5vP57Pk5GQ7/vjjzczbkyKY8vZLU83MPvjg\nA5s1a5ZlZ2dbdna2HXLIIXbttddaZWVll8eurq62q666ynJyciwzM9MmTJhgd9xxR2BJaFpamk2Z\nMqXb9jMzMy0nJ8fGjRtnc+bM6bDHwq9+9Svz+/0WGRkZeB1uu+22Xl+HYH5H6cdljs7CNbn0A+ec\n2dy5e3fQpZdC69CcyFDinAvfTzIiEtTvaGudfrn6k04xdPb226HugYiISMgpIHS2ahW0m+gjIiIy\nFCkgdNbS4m2cJCIiMoQpIHTno4+gnzf3EBERGUgUELpTX++FBBERkSFKAaEn77zjnW4QEREZghQQ\nelJZ6V3pUUREZAhSQOiNljyKiMgQpYDQm23boIeLboiIiAxmCgh90SiCiIgMQQoIffn8c9i5M9S9\nEBEROaj6vnbpUGcGb74JZ5wR6p6IHHDO9csW7iIyCCggBGPZMjjxREhPD3VPRA4YXahJRNrTKYZg\ntLTAkiVtEggoAAAgAElEQVSh7oWIiMhBo4AQrE8/hdLSUPdCRETkoFBACJaZRhFERGTIUEDYG6tW\nQVFRqHshIiJywCkg7A0zWLw41L0QERE54BQQ9taaNbB1a6h7ISIickApIOyL118PdQ9EREQOKAWE\nfbFhAxQWhroXIiIiB4wCwr7SKIKIiAxiCgj7qrAQ1q8PdS9EREQOCAWE/aFRBBERGaQUEPbHtm2w\nenWoeyEiItLvFBD21+LF3v4IIiIig4gCwv4qLvZ2WBQRERlEFBD6g0YRRERkkFFA6A9lZbB8eah7\nISIi0m8UEPrLq69CQ0OoeyEiItIvFBD6S3U1vPFGqHshIiLSLxQQ+tPbb8POnaHuhYiIyH4Lm4Dg\nnPM55651ztU752YHUf8U51yLc27+wehfUJqb4R//CHUvRERE9lufAcE5l+Wce8Q593nr7Qnn3PBg\nGnfORTnnbnbOrXLOrXDO/dM5d0I39fKB14FzgCig1yUBzjkf8H+t34bX8oE1a2Dt2lD3QkREZL/0\nGhCcc9HAK0AkMLH1VgMsds7FB9H+XXhv+iea2STgQWCRc25Kp3pXA38Ergmy3xcDW4Kse/D94x/e\naIKIiMgA1dcIwmxgEjDHzFrMrAWYAxwCXNHbgc65Q4HvA78xsx0AZjYP2Ajc0qn6NWb2GOD66nBr\nMJkL/KyvuiGzY4c3H0FERGSA6isgnAUUmtmmtgIzKwFWtt7Xm1l4b/iLO5UvBmY65+LatdkSbIeB\na4EXzOzzvTjm4HvjDW9lg4iIyADUV0CYjPeJv7NNeCMLfR3bDGzuVL6RL05Z7BXnXC7wA7wRhPDW\n0ACvvBLqXoiIiOyTvgJCBtDdx+AqIM45F9PHsbVmXfYgrmr9mh5cFzu4BbjLzMr24diDb/ly2BK+\nUyVERER60ldACJsVAs65I4DpfLF6YWB48UVdp0FERAacyD7uLwcSuylPAmrMrL6PY+Odc67TKEJS\n69cdwXcT8ILB3D4es4sbliwJ/HvayJFMGzlyLx92PxUVwYcfwtSpB/dxRURk0FuyZAlL2r3P9SfX\n9QxAuzudewkYb2ajOpWvAKrN7Phejp0D/BoYaWab25XfhTePINXMajsdMw1vP4SLzeyhduVJwGqg\nhI6jGlOACrx5DoVm9u1O7ZnNDYPpCnFxcOWVEBsb6p6IiMgg5pzDzPpcERiMvk4xPAUUOOcK2j14\nNjAeeLJTp7Kdc+079Xe8N/PpndqcDizqHA56Y2ZVZpZjZkeY2Zfabq13P9P6/bd7bSSUamu9S0KL\niIgMEH0FhAXACuB/nHMRrTsY/gbYgLexEQCtuyNuB+5pKzOzNcD9wHXOufTWepcAo4Dr+3jcvUk/\n/ZKUDrgPPoCSklD3QkREJCi9BgQzawRm4C1XXNl6SwBO7jQCUI031L+9UxNXAk8A/2w9LXEpMNPM\nlrev5Jw72Tm3EXgUb9ThdufcRufced31yzm3uLW+AWe31v1xUM84VFpa4PnnNWFRREQGhF7nIAx0\nYTMHob2ZM+H4HqduiIiI7LODOQdB+tvixd5WzCIiImFMAeFga2yEp5/WqQYREQlrCgihsGULvPNO\nqHshIiLSIwWEUHn9dZ1qEBGRsKWAECqNjfDMMzrVICIiYUkBIZQ2b4Z33w11L0RERLpQQAi1116D\nnTtD3QsREZEOFBBCTasaREQkDCkghAOdahARkTCjgBAudKpBRETCiAJCuNCqBhERCSMKCOGksFCn\nGkREJCwoIISbV1+F4uJQ90JERIY4BYRw09QETzwB9fWh7omIiAxhCgjhaMcOeO65UPdCRESGMAWE\ncPXpp/Dhh6HuhYiIDFEKCOHspZegpCTUvRARkSFIASGcNTXB449DQ0OoeyIiIkOMAkK403wEEREJ\nAQWEgWDFCvjoo1D3QkREhhAFhIFC8xFEROQgUkAYKBobvf0RNB9BREQOAgWEgaS8HF54IdS9EBGR\nIUABYaBZtgw+/jjUvRARkUFOAWEgeuEF2LYt1L0QEZFBTAFhIGpqgkcfhaqqUPdEREQGKQWEgWr3\nbvjLXzRpUUREDggFhIGsuBieegrMQt0TEREZZBQQBrrPP4fXXgt1L0REZJBRQBgM3nrLW90gIiLS\nTxQQBovnnoMtW0LdCxERGSQUEAaLpib4619h165Q90RERAYBBYTBpKbGW/5YXx/qnoiIyACngDDY\nlJTAk09qZYOIiOwXBYTBaM0aeOWVUPdCREQGMAWEwepf/4K33w51L0REZIBSQBjMFi2CTz4JdS9E\nRGQAUkAYzMzg2We9zZRERET2ggLCYNfSAn/7G2zcGOqeiIjIAKKAMBS0Xf1Rl4gWEZEgKSAMFQ0N\n8MgjUFYW6p6IiMgAoIAwlNTWwp//rN0WRUSkTwoIQ01VFTz0EOzeHeqeiIhIGFNAGIp27oSHH4Y9\ne0LdExERCVMKCENVcTH85S/e3AQREZFOFBCGss2bvZEEXdxJREQ6UUAY6jZv9iYu6nSDiIi0o4Ag\nsHUrLFwIdXWh7omIiIQJBQTxFBXBggVQUxPqnoiISBhQQJAvlJR4IUFLIEVEhjwFBOmorAzmz/f2\nSxARkSFLAUG62rHDCwnacVFEZMhSQJDuVVR4IaGiItQ9ERGREAirgOCc8znnrnXO1TvnZoe6P0Ne\nZaUXEnSBJxGRISeogOCcy3LOPeKc+7z19oRzbniQx0Y55252zq1yzq1wzv3TOXdCN/XygdeBc4Ao\nwLqpE+uc+4Fz7i3n3GfOuZXOuTecc6cH0xfZB1VV8OCDUFgY6p6IiMhB1GdAcM5FA68AkcDE1lsN\nsNg5Fx/EY9yF96Z/oplNAh4EFjnnpnSqdzXwR+CaXto6A/g9cJ2ZHWZmE4GngWecc98Noi+yL+rq\nvM2UPvss1D0REZGDJJgRhNnAJGCOmbWYWQswBzgEuKK3A51zhwLfB35jZjsAzGwesBG4pVP1a8zs\nMcD10qQBT5jZm4ECszuALXgBQw6Upib429/g7bdD3RMRETkIggkIZwGFZraprcDMSoCVrff1Zhbe\nG/7iTuWLgZnOubh2bbYE0ZfHgUu6KS8GUro9oiWYZiUoZvDyy/CPf3j/FhGRQSuYgDAZ7xN/Z5vw\nRhb6OrYZ2NypfCNfnLIImnma25c55yLwRjOWdHvQm292Wyz74Z13vNGEpqZQ90RERA6QYAJCBlDd\nTXkVEOeci+nj2FqzLh8323bhSQ/i8ftyJhAL3NTtvUuWwPr1/fAw0sFnn3nzEnT9BhGRQSmYgBC2\nY8nOuRzgt8ClZraux4pPPaWdAQ+EwkJvhYM2VBIRGXQig6hTDiR2U54E1JhZfR/HxjvnXKdRhKTW\nrzuC62ZXzrkU4EW8CZB/7aneDamp3mY/DzzAtDPOYNro0fv6kNKdsjKYNw/OOw9yc0PdGxGRIWXJ\nkiUsWbLkgLTtuo7+d6rg3EvAeDMb1al8BVBtZsf3cuwc4NfASDPb3K78LuAHQKqZ1XY6ZhrefggX\nm9lDPbSbAiwCHjazO3t5fLOf/Qzuu88bQTjmGDj11F6fr+yjqCg4/XSYPDnUPRERGbKcc5hZb6sB\ngxbMKYangALnXEG7DmQD44EnO3Us2znXvmN/xztFMb1Tm9OBRZ3DQTCcc8l0Ew6cc892e0BcHJxz\nDvh88O67sHLl3j6kBKOx0TuV8/LLWjkiIjIIBBMQFgArgP9xzkU453zAb4ANeBsbAdC6O+J24J62\nMjNbA9wPXOecS2+tdwkwCri+j8ftkoDahYNaYIdz7jttN+DYHlvKy4OZM71/P/OMdzEiOTDefhse\neUSTF0VEBrg+A4KZNQIz8JYrrmy9JQAndxoBqAYq8EJCe1cCTwD/bD0tcSkw08yWt6/knDvZObcR\neBRv1OF259xG59x57arNBr4MnAT8GXio3a33FRFHHw0TJ0JDAzz+uPeJVw6M9evh/vuhpCTUPRER\nkX3U5xyEgcw5ZzZ37hcF9fXwpz95IwhTpsAZZ4Drl1M10p3oaPj2t71gJiIiB9zBnoMweMTEePMR\nIiNh2TL4+ONQ92hwaxutee017bwoIjLADK2AAJCdDd/8pvfvF1+EoqLQ9mcoePNNePRR2LMn1D0R\nEZEgDb2AAN7phSOPhOZmeOIJqN3rxRSyt9as8Zabbu88RUVERMLR0AwIAKedBjk53iZKjz6qSYsH\nQ0WFt6nS22/rlIOISJgbugEhMtLb/S85GbZuhSef1Pr9g6G52dsr4dFHNXIjIhLGhm5AAEhMhAsu\ngGHDYPVqeOEFfbI9WNasgXvv9a7nICIiYWdoBwSAzExvJCEyEj76CN54I9Q9GjqqqmDhQli6VMFM\nRCTMKCAA5OfDWWd5eyIsWeIFBTk4Wlpg8WJ46CGo7u6q4iIiEgoKCG3Gj4d/+zfv388/7w2By8Gz\ncaN3ymHt2lD3REREUEDoaOpUOOkkb7j7iSe8yYty8NTUeNdxePZZb9dLEREJGQWEzqZPhyOOgKYm\nb6a9Lux08H30Efzxj96ogoiIhIQCQmfOeTstjhnjLcN7+GHYvTvUvRp6du3y5iW89JL2qBARCQEF\nhO5ERHjXbMjN9d6o/vIXbRMcCmbw7rve3IQtW0LdGxGRIUUBoSfR0XD++ZCa6l2v4c9/hrq6UPdq\naNqxA+bPh1df9TZaEhGRA04BoTfx8XDRRZCS4l1D4KGHtPtfqLS0wFtveddz0AW2REQOOAWEvqSk\nwMUXQ1oaFBd7G/toTkLolJbCn/7kbdfc0BDq3oiIDFoKCMFITvZCQkaG9wa1cKE29Qmllhbvgk9/\n+IP2qxAROUAUEIKVmOiFhKwsKC+HBQugsjLUvRra2iaQPv64ApuISD9TQNgb8fEwezb4/bBzpxcS\ndu0Kda9k5Uq4+2547z1d00FEpJ8oIOytuDhv4uLw4V44mD/fCwsSWvX18OKL8MAD3lwRERHZLwoI\n+yI2Fi68EEaM8K5IOH++d9pBQm/bNrj/fli0SNs1i4jsBwWEfRUTA9/5Dowc6a1qWLBAy+/CRUsL\n/OtfcNdd8PHHOu0gIrIPFBD2R9tmSocc4l1oaP58WL061L2SNrt3wzPPeMsitROjiMheUUDYX1FR\ncN55MHmyd82Av/7VW4KnT63hY/t2mDcPnnzSOyUkIiJ9cjaI38icc2Zz5x6cBzODN9+ExYu97486\nCk47zbuug4SPqCg48UQ4/njv3yIig4hzDjNz/dKWAkI/+/RTePpp75oBo0fD2WfDsGEHtw/St+Rk\nmDEDDj881D0REek3CghBCklAAO9891//6l23ITPTm6eQknLw+yF9Gz4cvvY1GDUq1D0REdlvCghB\nCllAAKiogEcfhbIyb4Olc8+FvLzQ9EX6NmaMFxT8/lD3RERknykgBCmkAQFgzx544gnYsMGbi/Dt\nb2tIO5w55/18Tj7Zu8y3iMgAo4AQpJAHBPDmIrz0Enz4off9iSfC9Ong0wKSsBUR4U0y/epXvdEf\nEZEBQgEhSGEREMBb4fDOO/DKK96/CwrgrLO8C0BJ+IqOhuOO81Y8xMSEujciIn1SQAhS2ASENps2\neWvxd+/2rulw5pneSgcJb7GxcOyxcMwxWpEiImFNASFIYRcQwAsHf/+7Ny8B4KSTYNo0nXIYCIYN\n84LCsccqKIhIWFJACFJYBgTwrhXw5puwdKlOOQxECgoiEqYUEIIUtgGhzcaN8NRT3qhCfDzMmqVT\nDgPJsGHeaYfjjlNQEJGwoIAQpLAPCOCFg6ee8sICwFe+4s2e1ymHgWPYMJg61RtRSEgIdW9EZAhT\nQAjSgAgI8MUphyVLvO/z8rw9E9LTQ9ot2UuRkTBlCpxwAqSlhbo3IjIEKSAEacAEhDYbN3oTGKur\nvTebk0/2hrA1mjCwOAcTJnh7XuTmhro3IjKEKCAEacAFBIC6Onj5ZVi2zPt+xAg44wyNJgxUo0Z5\nQUFzS0TkIFBACNKADAht1qyB557z5ihERsIpp3ijCa5ffu5ysPn93oZLhx2mS4CLyAGjgBCkAR0Q\nwBtN+Mc/YPly7/v8fG80Qee3B66EBG9C49SpmtAoIv1OASFIzjlbPPuroe7GfkvfUs64t9cSU9dA\nc4SPDUeOYtuE4RpNGMBaInyUjsxk68Q8dqdr/wsR6R83Tr+x3wJCZH80IgfWjhEZvJ+VzJj31uHf\nUMrY99eTtamUtceM1ZvLAOVrbsG/vgT/+hIqs5LZOjGP8vwMzKfQJyLhQQFhgGiKieLzkyZQVpDJ\nuHfWklxWzVHPf8T2Q3PZ+KWRNMVEhbqLso+SSytJLq1kT3wM28cPp2iMn8bY6FB3S0SGOAWEAWZH\nfgbv+VMYuayQvFVbGb56O5mbyth45CiKxvp12mEAG1ZTzyEfbmDkxxspz89g+6G57MpJDXW3RGSI\nUkAYgJqjI1n/5dEUjfEz9r21pBZXcujba8hZW8TaY8ZSnaHTDgOZr8XI2lRG1qYyapPj2D4uh+Ix\nfo0SichBpUmKA515byaj319PTF0DBhSNy2HDl0bRNExvKINF26TG7YfmUpWVHOruiEiY0iRF+YJz\nlI7KYkdeGgXLCslbuY3cNUVkFpax8YiRFI3LwbQT44DXflJjTUo8xWP8FI/O1lwFETlgNIIwyMTt\nqmHse+tILdoFQG1iLBuPHElZQabmJwwyLT7Hzrx0isf42ZGXrhUQIqIRBOlZbUo8y2ZMJmNzOYd8\ntJG4qjoOW7qKqvStbDhqlCa9DSK+FiNjczkZm8tpiI2m5JBsisf4qUmND3XXRGQQUEAYjJyjvCCT\nHSMy8K8rYuQnhSTtqOaIRcvZmZvKhiNHaf+EQSa6roERn21hxGdbqMpIpHiMn7KRWTRqHoqI7CMF\nhEHMfI6icbmUHJJN3qpt5K/YTNr2CtK2V1AyKouNXxrJnsTYUHdT+llSeTVJ5dWMeW8dFblplIzO\npnxEOi2RugaEiARPAWEIaImMYPOkfLaPzaFgxWaGf76N7I2lZBaWUTQ2h82TRlAfPyzU3ZR+5msx\n0rfuIH3rDpqiIijPz6DkkGx25aRqvoKI9ClsJik653zANcDNwA/MbGE/tDnkJikGI2b3HkZ9sons\n9SU4vMluxaOz2TwpXyMKQ0BDbDSlo7IoOSRbe2aIDDIHdZKicy4L+C1wVGvRCuBqM9sWxLFRwH8D\nZwNNQBVwrZn9s1O9fOAhIB6IAnpMLc65C4CfAQ7wAb83swf66ot8oT5hGJ+fOJ7Nh42gYHkhWZvK\nyF1bTM66YkoO8YJCbXJcqLspB0h0XQN5K7eSt3IrdQnDKBuVRenITM1LEZEOeg0Izrlo4BXgc2Bi\na/GDwGLn3JfMrKaP9u8CpgEnmNkO59z3gEXOuePNbFm7elcDfwSKgcW99OdcYB5wopl94JybBLzt\nnPOZ2f199EU6qU2NZ9VXJ7LpiFryV2zGv8FbZ5+9voSykZkUTs6nJlWXJB7MYnfvIX/FZvJXbKYu\nMZbSkZmUjcpid5p+7iJDXa+nGJxz3wfuAw4xs02tZdnANuDnZnZ7L8ceCqwEvmdmC9qVfwpsMrNv\ntivzmVmLc24a8DpwsZk91Kk9H1AIvGZmF7crvxs4F8g1s4ZOx+gUw14YVl1H/qdb8K8rxtfi/b8o\nG5HO5skFGooeYmqTYikb6Y0s1CgsiAwYB/MUw1lAYVs4ADCzEufcytb7egwIwCy80wCdRwQWA5c5\n5+LMrLa1zZYg+no0MLyH9n4ITAdeDqId6cGexFjWHDeOwskFjPh0Czlri8jcsoPMLTvYld16SeK8\ndNAEt0EvrqqOguWFFCwvpC4xlvL8DMrzM6jMStKGWyJDRF8BYTLe6YXONgEnB3FsM7C5U/nG1sed\nCHzQdxc7tNd2fOf2ACahgNAv6uNjWHfMGDZPzifvsy3krikipaSSlJJK6hKGsXXicIrH+GmO0iKY\noSC2ui6wx0JDbDTlI9IpL8ikwp+CRWgbb5HBqq+/8BlAdTflVUCccy7GzOp7ObbWup7DqGr9mh58\nNwPt0U1/9rU96UNDbDQbpo6mcEoB/rXF5K3aRuzuPYx9bz2jPt5E0dgctk4YTn2ClkgOFdF1DeSu\nKSJ3TRFNURHsyEunPD+DncPTaI5WYBQZTPr6jQ6PNZASUs1RkWybmMe28cPJ2LKDvJVbSSmtZMTK\nreSt2kpZfiZbJw6nKlPDz0NJZGMz2RtLyd5YSovPUZmdwo4R6ezIS6cuSctlRQa6vgJCOdDd7LQk\noKaX0YO2Y+Odc67TKEJS69cdwXcz0B7d9KfX9hZ8sinw7yP8KRzhT9nLh5UAn6O8IIPyggwSy6vJ\nW7mVzE1lZBV6t92p8Wwfl0vJIVn6NDnE+FqM1KIKUosqGPPeOmqT4yhvDQtVWcnamEnkANn0ySY2\ntXuf6099/RVfDozvpnwU3n4IvVmGt7pgBB3nIYwCGvFWOOyNtmWRI4E3OrXX1tcuLj5i5F4+jASj\nOiORVV+ZwIajDiH3823krC0moaKGce+uZfSH6ykZlUXRuFytfhii4iprya+sJf/TLTTGRLIzN42d\neenszE3VJapF+tHII0Yyst373NKFS/ut7b4CwlPAfc65AjMrhMAyx/HAz9tXbC0vbTda8HfgVrzV\nBe13RZwOLGpbwbAX3gO2th7ffgnkdLzRgyV72Z70g/r4GDYedQibjhhJ5uZyctZsJ7W4kty1xeSu\nLaY6PYHt43IpHZVFc5SuBTAURdU3BU5FmIPdaYnsHJ7GzuFpVGUmaXRBJEz1tQ9CFN5Kg1XABXhz\nEuYBxwNfanuTd86dgPep/j4z+2G74/+I9wbetlHSJcDdwHFm1uUTf7t9EC7pbqtl59x/AAvwNkr6\nsHWjpH8BPzGzP3VTX/sghEBcZS05a7bjX19CVH0TAE1REZSMyqJ4jN8bVdBcBQGaoiOpyEkNBIb6\n+JhQd0lkQDto+yCYWaNzbgbeVssr8QLCCuDkTiMA1UAFsL1TE1cCc4F/Ouca8VYczOwcDpxzJ+MF\nj2Gtj3G7c+4G4Bdm9mi7/jzmnIsAHnTeG4wP+H/aajm81CbHsf7LY9h45CFkbiojZ812UkqrGL6m\niOFriqhNiqV4dDYlo7N1kaghLrKhiczCMjILywCoSY6jIjeVitw0dvlTNOokEkJhc7GmA0EjCOEj\nvqIG/7pisjeUEL2nEfCS4C5/CsVjsinPz9SbgXTQ4nNUZyR5gSEnVacjRILQnyMICghyULkWI3X7\nTvzrS8jYXB7Y0rk50kdZQSbFo7PZlZ2i3Rqli6aoCCqzU6jITWWXP4XdqfE6VSXSyUG9mqNIfzKf\n82az56V7w8ubSvGvLyG5tAr/eu9iUQ3DoigryKR0VCaVWcl6ExDA23chfesO0rd6K5obYyLZlZ3C\nLr93q1FgEOlXCggSMk3RkRSNy6VoXC6xVXVkbyghe0MJsdV7GL56O8NXb6c+LprSAu8Kg1Wa3Cjt\nRNU3kbm5nMzN3hYpDcOiqGwfGFLi9P9FZD/oFIOEFzMSdu4ma2MpWZvKGFbzxV5ce+JjKB2ZSWnb\n5Yj1x1960RgTRWVWkhcaspPZnZ6oOQwy6GkOQpAUEAY4M5LKq8nc5IWFmNovrua9Jz7Gu8LgiHQq\ns1P0h1/61BwZQVVmEruyk6nMTqYqM4mWSE2MlcFFcxBkaHCOqswkqjKTWD91NMmlVWRuKiWzsJxh\nNfXkrdpG3qptNMZEehcNGpHBztxUWrQaQroR0dQc2A4avFUSNakJVGYlUZWVTGVWsvZhEGlHIwgy\n8JiRWF5N5uZyMjaXE1dVF7irOcJHRW4q5SMy2JGXpm19Za/siY/xwkKmFxp2pyVodEoGFI0gyNDm\nHNWZSVRnJrHhqEOIq6wlozUsJJVXk7FlBxlbvJnuVemJ7MxLY8fwNO3gKH0aVlPPsI2lZG0sBbzl\nt9XpiVRnJAVGszTKIEOFRhBkUImurfcCwuZyUop3BfZZAG+Wu3fRoDR25qbSFBMVwp7KQFUfF01V\nhhdQqzKTqE5P1CZfEjY0giDSg4a4GLYfmsv2Q3PxNTaTUrKL9K07Sd+6g2E19fg3lODfUII5qMpM\nYmduGhW5qVRrhrsEKaa2ocPySnOO2uQ4qjMSqcpIpDo9kZq0BFoifCHuqcj+0QiCDA1mxFXWkrZt\nJ+lbd5JcWtlhdKEpKoIKf0pgW9+6pFidjpB91jYBsro1MFRlJFKbEq8QKgecRhBE9pZz1KbEU5sS\nz9bDRhDR0ERq0a7ArPa4qjoyt+wgs3Xuwp74GCpyUr3A4E/RZEfZK74WI3FHNYk7qgNlLRE+dqfG\nszs9keq0BHanJ1KTGq+RBglbCggyJDVHR1JekEF5QQYAMbv3kFpUQdp2LzAMq6knZ10xOeuKAe8q\ng94OfcnsylZgkL3na24hqbyapPJ2ocHnBde2wLA7LYHdqfE0R+tPs4Se/heKAPUJwygem0Px2JzA\nbo6pRRWkbq8gubSK+Mpa4itrGb7au6J5TXIcu7KT2eVPodKfQoMCg+wDX4v3fy1h525oDaPmYE9C\nrDfakJYQuNUn6NLocnApIIh05pz3aS49kS2H5+OaW0jcUU1K8S5SiitJLqv8IjCsKQKgNimWyqxk\nb2vfrGTNYZB95gxiq+uIra4LTIQE7+JUu1MTqGkdZahJTaAmJU67QcoBo4Ag0geL8FGVlUxVVjKb\nJ/NFYCipJKV4F8mllcRV1RFXVRc4JdEwLCqw2U5llrcUznSuWfZDVH0TqcW7SC3eFSgzB3WJsdSk\ntoUGLzjUJQ5TQJX9poAgspc6BIZJ+biWFhJ27ia5pIrkskqSSyqJ3tPYYdJjc4SP6vQEb918hjbc\nkf7hjEA4zSwsC5Q3R0ZQmxznBYaU1uCQHKfTFLJXBv0yR24IdS9kyDE4pAJO2AwnboYTtsBhZV2r\nbUuEd/Lg3eHwbh58kAu1msogIvvjBnQ1x2AoIEi4SKuFo7fBsVvhmK1wzDZI3dOxTpODz7Lg/Vwv\nLOtdn/QAABZVSURBVHyQCyuyoUHjfCISrBsUEIKijZIkbJkRW1VHUlmVt/StrIqEit24Tr+OLT5H\nTUp8YMOd6vQEalLiNZ9BDoqWCB91ibHUJsd1uWl76fCkjZJEBjrnqEuOoy45jpIxfgB8Tc0k7Nzt\nbbBT7n2Nq6wlceduEnfuBrwVE22hYXd6ArtTE7zQkJqgP9jS73zNLcTvqiF+V02X++rjoqlN8sJC\nXXIctUmx1CXFsSdhmHaMHCQUEETCREtkRGDyY5uIxiYSdrSGhtavcVV17UKDx4C6pNgO6+Z3p8Z7\n+zNoNrscADG1DcTUNnRYVQFegN2TGBsIDLVJsV4YToylPk7/HwcSBQSRMNYcFUll62ZMbSIamkio\n2N26wU4NCTt3E7+rJjCbPWvTFzMiA2vnU+O/WAaXEq+183LA+Fq8657EVdYCOzrc1xzpnbKoS4yl\nLqnd16Q4hYcwpIAgMsA0///27j1IsrMu4/j3N33v2d1skpXFJGQTtAAjbKCkEAilJJEECizlIoKi\nEQRTKCkoAQEhRSBoEcFCsLgZLikuFpAKFw0aIlZSWAENSNjN1QDZBJQQctns7sz0vX/+8b592T49\ntz49092zz6fq1Onz9ukz77vd8+4z73nP6XyWQ7t3cmh3LzRYq838oSW2PbjQDQ/zBxeHXztPGG3o\nXgK3s8zicfNUdpQ0t0E2VKbZZtvBRbYdTJ6yaGXnqM4XqewoUe2EiO1FKttL4bSFPpubTgFBZAvw\nzFz31EKv0Cks1Zk/uMC2g4vMx6V8aKnv2vnenfraZiE47JxnaWe5Gx4qO0r4nDpn2ViZZrt7h9JB\nblArF0JYiIGhsr0YgsS2or4bZYMoIIhsVWbU5gvU5gs8dMqJveJWm/KhJbYdXKT88CLzDy8x//Ai\nxYVqr4O+p3eYthnV7UUWOzPYd5a7k9P0pUKyGcyhuFijuFiDgTkPEG4M1R8aqtuKRy1NfU5Hon81\nkWOMZ+ZYPCHc07/fXLNF+dBSnLUe1uWHlyguVLsjDvz46HPKtVK+Fxh2hPPJS3EmO5rJLpsk02wt\ne7UFQDOfPSowVLYVqW0rUp0vhABRyG1yjWeDAoKIAOEqis6XVPWba7YoHa4wHyee9S+FSp1Cpc7x\n9yZnsncmoB01m317ibomo8kmy9abvW/NHKKVzXTDQq0TJOYL1ObDul4uHJOXbiogiMiK2tnM0BEH\n2k5xsRrDQoXS4TC3oXS4QnGptuz55FZmrjeDfXux73EpfD+FwoNsskyzteznFcDNqJfyITT0hYfa\nfCE8Lhe25DwIBQQRGU283r26vcRDpww81WhRPhJDQwwPpSNVSocr5GuNZWeyt+csnktOTkSrbi/S\nyqnLks1n7hSWahSWanD/4aH7tDNz1Mphzk//uhrnAdXLBerF3EwFYP22icjYtXOZ5FUVUbbepHgk\njDSU+tblwxXy1UZvvsMQ9WIuMQGt/5xyW5fCyYTMtdrh83xk+GcXQgCulwvUSvlugKj3h4pynnq5\nMDWfYwUEEdlUzXx26FwHCCMPpYVqCBALVYpHqpQWKmF9JASIfLXBjgeODD12rZTvnT/eVqQ6X6S6\nLZ5L3lbQCIRM1FzbKS5UKS5UYcg3vHY0ClnqpRAYuiEihodaKU89Lht9bwj9tojI1GjnMuEGTsfP\nJ5+M93XodLCDS2Gx1p00edwynW8jnw3njbvnkQvUyr1zyfVyXvd8kInL1Zrkas1lr8qAcG+IRiEX\nwkJfcBgnBQQRmQ1993U4tPu45PPtcJ64uFCluFiluFCjENfFxRAgcvUmuXpz6PwHCJ1uvZgPow4x\nOITzyfk4BKwQIdPBnO6IGst8ntNSQBCRrWHOqMW5CIeGPe9OrtqguNgJDiE0FJdqYfRhsUY+jkAU\nKnVg+GkMB+qlfDcsdM4dd7dLYd3MZ2dqQprIIAUEETk2mNEo5WmU8hzZlZz/AOEuk4VKnUIccSgs\n1uK3FoYZ7IkQ8eDQwwDhcs7+wFAr56mXCjFc9M4jK0jItFJAEBGJPDPXvTJiOdb2EBI6gWEpzn1Y\nqpFf6q2zzVa4tPNIdcWf2crMdcNCcslRL27epDSRfgoIIiLr4HO9uRD83PL7ZRrNbmAoLNXJV+rk\nKzFEVOJ2J0gsVCktrBwkIEyyrJfyNIqd4JCjUewFiUYpTz0+187pK70lHQUEEZEN0MplqRyXpXJc\necX9Mo1WLzxUGvHxkKVa706yHD7JYuDnZ+dCaCjkaHQCRDFHvRhCRfdxIUejmKOdVaCQoykgiIhM\nUCuXoZIL31uxIndytQb5SoNctR7CRLVOrlKPs9nr5GJZvlIn02yveWQCBgJFMS6FvnUhR6OY7ZY1\n87lj8vsJjiUKCCIis8As/uWfB4bcJ6KfO5lmqxsYcvFyuFw1hIlcpR7CRrVBrtogV2usO1BAOOXR\nCQ7NTogohLLh29lwl0BNypwJCggiIluNGa1cllYuS3W1kQkIgaLRCkGiFoJErtokV6uHm/bEENFd\n1xpka83eKY/hV4QO1Z6zGBiyNPIhNHTCQzMfH+cHHsfnFCw2lwKCiMixzoxWPksrn6XKGgIFQNtD\nQKg1yHbDQ992fNwNE/HxXNv77jWxdk64TXd3KWRjkMjRzGdoFnI089leuBhYdAXI+ikgiIjI+s1Z\nd67Cul7WbPUCQ73ZDRPZzna973EtbOdqDbKNVm/EYgStzNxAaMjQyg0EiVym+zg8l6GZ62xnjrkR\nDAUEERHZNO1shno2Q32+sK7XWdtDeKg3ydYbvfDQLesEikbcbvXt3yTTapMZYeSiXzMXAkNrIDh0\nwkUnVLRyR283c2G/Vi5LKzs7czAUEEREZOr5USMWazwN0n2xM9dsd8NCrt4k02j2gkR8nOkEikYv\nXGQa8flGq7uwlKIdhCtXeqFh8HGmGya6j7PZgecytLJhn428kkQBQUREtjYz2rkM9dz6Ry662k62\n2QohohGCRabRCxiZblncjkGju91ZN9u9oDEGrcxcDAthGScFBBERkdXMWXd+Qi3FYawdLkHtBYhO\neBh43AzhItPoBY1MLO8EjkyzFU6dtNpQa4ytqR0KCCIiIpvExxQ0wsGcuVb7qIDB1d8dRzUBBQQR\nEZHZZEY7m6GdzdBY57SMtdCFoSIiIpKggCAiIiIJCggiIiKSoIAgIiIiCasGBDN7hJl91szuiMuV\nZnbyWg5uZjkzu9TMbjezm83sBjM7a5l9X2dmt5rZPjP7bzP7rSH75M3s4r79bjezy83skWupj4iI\niKzNigHBzPLAvxGudjgjLovAdWa2yveNAvD3wO8Az3D3JwCfAK41szMHfs6bgbcCz3P3M4E3AVea\n2bMHjvce4M3A78X9zgKeCnzVbEbuXSkiIjIDVhtBuAB4AvAmd2+7e5vwn/ejgVev9EIzeyzwKuDd\n7v4ggLt/HDgA/FXffjuBi4EPuvuBuN/XgWuB9w4c9qXAte6+L+73EPBx4EnAY1ZtrYiIiKzJagHh\nhcA97n53p8Dd7wNui8+t5PmAAdcNlF8HnGdm5bj9bMKNtYftd0YMGh0NYPCrwzrbo33Fl4iIiCSs\nFhD2Ev7iH3Q3YWRhtde2gB8NlB+gd8qis1+nfHA/Bn7OO4FzzewcADM7DbgQ+JS7/3CV+oiIiMga\nrXYnxV3AkSHlh4GymRXcfbm7Re4Cltzdh7wW4MS+/Rjycwb3w90/amYt4Cozq8TXfhh4/SrtEBER\nkXVYbQRh8D/3iTKz9wB/DZzv7icBpxBGIP6j75SFiIiIpLTaCMIDwPYh5TuAxRVGDzqvnTczGxhF\n2BHXD/btR/w5B5fbz8x+mTBScLG73wjg7j8zs9cC3wMuAi4brMTHbv1x9/Hek07gzJNPWKHKIiIi\ns2Pf/z3E/p88tCHHXi0g7AceN6T8dODmVV67D3gJ8CiOnodwOmGy4W19+wGcNmS/Th0AHh/X3x/4\nOZ3tvQzxmW9raoKIiGxN5w9sf3aMV/yvdorhi8AeM9vTKTCz3YTQcFX/jma2e+BeBF8inKI4e+CY\nZxMuVVyK29cAS8vsd6u73xm374vrPQP77Rl4XkRERFJaLSBcQRgpuMzMMmY2B7wbuIswORCAeHfE\nnwAf7JTF/9j/AXiLmZ0Y93s5YWTgrX37HQIuBf7MzE6P+/0GcB7whr66fIN4KsHMfiHuVyacVlgC\nLl9n20VERGQZK55icPeGmT0LeB/hlIATAsM5fSMAEK5AOEgICf0uAt4O3GBmDcKVCee5+/7+ndz9\nMjOrAlebWZNweeSL3P1rffu0zexc4C3AV+LVDHngVuBp7n77OtsuIiIiy7DkVYhbR3J+pIiIyNZl\nZrj7WCYi6NscRUREJEEBQURERBIUEERERCRBAUFEREQSFBBEREQkQQFBREREEhQQREREJEEBQURE\nRBIUEERERCRBAUFEREQSFBBEREQkQQFBREREEhQQREREJEEBQURERBIUEERERCRBAUFEREQSFBBE\nREQkQQFBREREEhQQREREJEEBQURERBIUEERERCRBAUFEREQSFBBEREQkQQFBREREEhQQZtj1118/\n6SpsqK3cvq3cNlD7Zp3aJ6CAMNO2+od8K7dvK7cN1L5Zp/YJKCCIiIjIEAoIIiIikmDuPuk6bBgz\n27qNExERGcLdbRzH2dIBQUREREajUwwiIiKSoIAgIiIiCcdsQDCzd5lZ28wumHRdZOsws583s2vM\nrD3pumwEtW9Nx5javmUrv39buW2TMnMBwcweYWafNbM74nKlmZ28zmOcAvw54HGZGmnbZ2ZnmtlX\nzOy7ZnZ7PMZlG1nn9UjTvtgBfCy2a5+Z3WpmbzOz3EbXey3M7AXADcBprPNzZWY5M7s0tu1mM7vB\nzM7aiHqOatT2xfftHfE9uzm28Soze/xG1XUUad6/vmNMc9+Sqn3T3Lek/N2b6n4FwMyeaGaXm9lt\nZrY/1vH9ZrZrDa8dvW9x95lZgDywD/g8IdzMAVcAdwLz6zjOp4B/BtrAH066XeNqH/B04H7g7L6y\nPwXumnTb0rYv7nsTsB84PpY9EVgC/m7SbYv16XRQVwDtdb72I8AdwIlx+4+BReDMSbcrbfti2/4H\nODluF4AvxPY9ftLtGsf713eMqexb0rZvBvqWUT+bU9+vxDrdAVwJlOL2ScDt8fequMprR+5bJt7w\ndf4jvSr+4p3WV7YbaAJvWOMxfgX4AXDetP0Sp2kfYPEDc+lAeRY4f9JtG0P7zoivfe1A+ZeBn066\nbZ33IK7X20k9FmgBfzRQfgtw9aTbNYb2fRh4xUDZo+P7+YFJtytt+/peP7V9S8r3bxb6llHbNvX9\nSqzPbcCjB8peEev+ghVel6pvmbVTDC8E7nH3uzsF7n4f4R/vhWs8xt8CfwnUx1679NK07xmED8PV\n/YXu3nT3r425nqNK075mXA8O++WAxrgqmIbH37wRPJ/QCV83UH4dcJ6ZlVNVbExStO81wCcHyu6N\n652j12i8UrSvY5r7ljTtm/q+JUXbpr5fifa6+10DZWv5HUrVt8xaQNgLHBhSfjfwhNVebGa/DRTc\n/Qtjrte4pGnf0+O6HM/v3hLPN/2VmRXHWckURm6fu98J/CNwoZntATCzc4BzgLePt5qbbi8h5f9o\noPwA4a+0Mza9RmPk7q0hHfhj4vr6Ta7OhpiBviWNWehbRjIr/Yq7N4cUP4Yw3+IbK7w0Vd+SXUcd\np8Eu4MiQ8sOED2/B3WvDXhgnnLybMCwzrUZuH/CouP4U8GJ3/1acBPavwJOB88de2/VL0z6AC4D3\nAd83s/uBMvAad//E+Ku6qXYBS0P+Ez0c1yducn02w58Qhjk/PemKpDUjfUsas9C3pDFz/YqZZQhz\nCT7m7j9YYddUfcusjSCkGQJ8NXCLu39zXJXZAGna10nyn3b3bwG4+y3AZcCzzOzX0lZuDEZuX/xL\n5XpCh7TH3U8mpPx3mNn7x1M92Qxmdi7wYsJ/NtM0jDuqWehb0piFvmUkM9yvXAzUgNdt5A+ZtYDw\nALB9SPkOYHGF0YOdwJvjknh6fNVLbaT2RZ2/zL83UN7ZfnLKuo1Dmva9AjgLeKO73wvg7jcBfwNc\nZGa/Ou7KbqIHgHkzG/ws7ojrBze5PhvGzM4kTCT7TXe/Y8LVSW2G+pY0ZqFvGdXM9Stm9nLgRcBz\n3L2yyu6p+pZZCwj7gdOHlJ8O3LzC655KmIxypZndZGY3AZfH594Zy9423qqOZNT2QZhlDMn3tLVM\n+SSkaV9njsL3B8o723tT1GvS9hHen0cNlJ9OmCh126bXaAOY2V7gS8Dvuvt/Tro+YzIrfUsas9C3\njGqm+hUz+wPCfTbOcfcH1vCSVH3LrL2xXwT2dCaTAJjZbuBxwFX9O5rZ7k5qcvdr3P1Ud39SZwFe\nGXe9OJa9a5PasJKR2hf9C+EXdvAD3bkZzbfHX911S9O+++J6D0fbM/D8tFj2dMqQtn0p7n/2wK5n\nA9e6+9IG1C+t9bSvEw6+DLysMxQfb1DzkY2t5sjW1L4Z6lsGref9m4W+pd962jYz/YqZvQz4C+Bc\nd/9ZLHuemb2qb5/x9i2beS1n2oVw6ck+4HNAhhBwPkm4WUS5b7+zCB/oD61wrGcSriG9YNLtGlf7\nCJdZ/Qz4pbh9EiEJXzPptqVtH+EmKIeArwHbYtmphOvO7wDyk25fX12vYJlrsVd47z7M0TczeTnh\nZiZ7J92etO0j/JV2f2zjy/qW1wHXTbo943j/BvaZur4lbfumvW8ZtW2z0q8Avw9UgNcP/A59FHj7\nKu/dyH3LTF3F4O4NM3sWYcbpbYRkdDNhuKU/CR0BDgI/GTyGmT0S+BZh4o0D7zWzSwjDnjdubAtW\nNob2vZFwzumfzKxJuErlKqbkcp007XP3u83sKcAlwLfNrE4IHF8l3MBl4teem9mHgOcQZg67mR0g\ntPGx3puMt9x7dxHhfbrBzBqEWcbnufv+Tan8GqRo3yXACcCFcel3/QZWeV1Svn9T3bdA6vZNdd8y\nattmoV+JPkC4E+17BsodeEd8vMCY+5bO3adEREREumZtDoKIiIhsAgUEERERSVBAEBERkQQFBBER\nEUlQQBAREZEEBQQRERFJUEAQERGRBAUEERERSVBAEBERkQQFBBEREUlQQBAREZGEmfqyJhGZLma2\nC3gTvS+IudPdvzrBKonImGgEQUTS+AjwOXd/H7AEvHTC9RGRMdG3OYrIyMzsJuBu4Argv4Aldz88\nyTqJyHhoBEFE0rgQOB74HCEgFCZbHREZFwUEERmJmR3v7je6+zOBk4D7gWdMtlYiMi4KCCKybmZW\nBv7XzJ4Tiw4B9wE3TK5WIjJOmoMgIiMxs0uAe4B54FTgC+7+HTN7OnAu8CDwU+CJwAFCiDgVuNfd\nPz+RSovImukyRxEZibtfssxTZeAB4Dh3/5CZOfBcd3+lmT0FeAmggCAy5XSKQUTGyt2/DjwTuDIW\nnQV8Jj7+deCbE6iWiKyTAoKIbIRfdPcfxsdPoxcKngv8u5k9bTLVEpG1UkAQkbEys1OBG+PjInCf\nu9fj03cBzwO+M6HqicgaaZKiiIiIJGgEQURERBIUEERERCRBAUFEREQSFBBEREQkQQFBREREEhQQ\nREREJEEBQURERBIUEERERCRBAUFEREQS/h+6sCyZNa1IxAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106227110>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We see a few things. Firstly the lower limit is identical for both models for all $s_m$. This shouldn't be surprising: one of the edges of the interval will always be $\\udt{\\dt}$, from the definitions of the error and the interval.\n",
      "\n",
      "Second, we see that for small enough $s_m$ the predicted solution for the measured model falls outside the interval for the exact model: this occurs when $s_m \\simeq 0.6$. For large enough $s_m$, the predicted solution for the exact model falls outside the interval for the measured model: this occurs when $s_m \\simeq 1.6$. \n",
      "\n",
      "This gives a very wide range of values for $s_m$ that are \"close enough\". Also, it's specific to this data. However, we can do a general calculation."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The general interval"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have a definition for our two models, defining two intervals and the solutions that should fall within them. We can do a symbolic calculation to check when this happens. First, we write out the inequalities that need to be satisfied:\n",
      "\n",
      "$$\n",
      "\\begin{equation}\n",
      "\\left| \\frac{2^{s_e} y^{(h)} - y^{(2h)}}{2^{s_e} - 1} - \\frac{2^{s_m} y^{(h)} - y^{(2h)}}{2^{s_m} - 1} \\right| \\le \\left| \\frac{y^{(2h)} - y^{(h)}}{2^{s_e} - 1} \\right|\n",
      "\\end{equation}\n",
      "$$\n",
      "\n",
      "and\n",
      "\n",
      "$$\n",
      "\\begin{equation}\n",
      "\\left| \\frac{2^{s_e} y^{(h)} - y^{(2h)}}{2^{s_e} - 1} - \\frac{2^{s_m} y^{(h)} - y^{(2h)}}{2^{s_m} - 1} \\right| \\le \\left| \\frac{y^{(2h)} - y^{(h)}}{2^{s_m} - 1} \\right| .\n",
      "\\end{equation}\n",
      "$$\n",
      "\n",
      "Noting as above that only two need solving (as one edge of the inequality will always agree), we can solve:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sympy\n",
      "sympy.init_printing()\n",
      "yh, y2h, se, sm = sympy.symbols('y^h, y^{2h}, s_e, s_m')\n",
      "Eq1 = sympy.Eq((2**sm*yh-y2h)/(2**sm-1)-(2**se*yh-y2h)/(2**se-1) , (y2h-yh)/(2**sm-1))\n",
      "sympy.solve(Eq1, sm)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\left [ \\frac{\\log{\\left (2^{s_{e} + 1} - 1 \\right )}}{\\log{\\left (2 \\right )}}\\right ]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAJkAAAA/BAMAAAABN3d6AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARM1UEDKrie+7Inbd\nmWYDwv/QAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADw0lEQVRYCe2YTWgTQRTH//lsssmGguC1tSge\nPDRQDwqWBoXqQTAIFj0UI3gQFY14yEGC8RMPHgKi0JMLnpTW1lPRIuTgB1SxEaTiodKDKCoWqy1V\nUdZ5szubTXZ38tEc+w4zb9/Hb9/Mzm4eQZeud6MdclvXk+jq6+tsBwzr+gYYrSnUfll0xJMWci14\n+KQbbcI0etMOOdIyZNnhMDPDetPoTdvoSCuSxZXWMU8uwJOmiuo7ppJGKLxpeGKECNrQZjNFTOMp\nYOzdTeB1iFOYXUJ7bKQJWvCqwJhzHghmcL4TI1tTpsmdFrjM3ONGiKDhlpkipktAOIVIAUN6mmzh\nnp6LPVuc+zY4u8S88XmKqexbDS3EYiIlhJcDWngCnxYo1r22GNFU7vOkJQqAf4nR4mk1nSiqFm0b\naXbhtFiBm+wr3T11hh2d/NeHzKN2c3d8KbHwAP6zH+iK7j82c4I7KgOn+ajAqpWGphGdYOZXMWaP\nGqV3ZSkoqgVpMkyk2YXTQsvcZKstnkVoOVwCf3f8Ge5+xkflfWXfuME+cJqyyE022uEksKIWDFq8\nTG6VjyJXEwqw8yXJJjIYtJ/cZ6MNdAK/leUgOxvsCWRoNA8lqZ4iqe2fMjmZpkROC5Rxz5MiHMa+\n/eGXttoiZfgWA6xAknCRDcPAI34lGzit9pmOIHQF/qxvwwGNcukAKc9zsyW6kAqn0fFkImob0o9j\nLHcE6NV/lZmdbubXdb1EQTJRj/045XwXzIxoMnjnBekXZAiHzz/BTaI24d/OlM908VZYGpr3pXhY\nLW1vGrhBnoPc3ejw0QispQVzU3kGZA812SiJ4p4awbU0CxHsttT6SiJrxHjS0F8fYkXcNTVvmjjH\nVopEeVOXJkn2dHnX5pkicQgaO/arFXYXQZPcsAnXGs11s+4b1i+uTjJK9s3R6ChFKHMjaYQ1L5yE\n5mh0/BoG0bEI71dYQnM0Ot+Ao8A14FwLtTlo08B14HQae9pAU0rsVyNFNH+nB67OSu2NTqzMGTOp\nmhbARpbTqhudecrzrQCJbtJcRE6ranTi/NsezwCBkguJTHJaVaMT4bs1x5JC31uiVTU6nMY3r0Ua\n1WY1Onyl/VA0+Eot1VbV6ETn2RrLiGmtPYWaRqejDIzm8qz1VLPN1+ZodNj6enX9L2t30s3TRIbV\n6LA3y5BRodTOshNixlqNDnvrDWnprTdzrUYnqhmW1r5IJs1qdNjXkougmn7b1MBKK9Gr+ZJXKA1r\nTdVWl7pGq7tFrgG0b239V2pX+/6V6tP+A/XNQ3RzkyOMAAAAAElFTkSuQmCC\n",
       "prompt_number": 11,
       "text": [
        "\u23a1   \u239b s\u2091 + 1    \u239e\u23a4\n",
        "\u23a2log\u239d2       - 1\u23a0\u23a5\n",
        "\u23a2\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u23a5\n",
        "\u23a3     log(2)     \u23a6"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Eq2 = sympy.Eq((2**sm*yh-y2h)/(2**sm-1)-(2**se*yh-y2h)/(2**se-1) , -(y2h-yh)/(2**se-1))\n",
      "sympy.solve(Eq2, sm)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\left [ \\frac{\\log{\\left (\\frac{2^{s_{e}}}{2} + \\frac{1}{2} \\right )}}{\\log{\\left (2 \\right )}}\\right ]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAIkAAAA/BAMAAAAmmfaSAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARM1UEDKrie+7Inbd\nmWYDwv/QAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD40lEQVRYCe2YTUhUURSAz+j8z5tpQKpdTla4\nKHBgNgmZD4JsITEISS3EKVpZkK2EYsSKokWQUAQG0YO2mi4CCze2NdBZBOFCcFcmiYJiRfE659z7\n3ryfueNMugm6i3vPPefc751773vnDAPNppmB3bSXppmF5lwuXQVyaLWKkU1NuU6kqL2KL3QYOF/Z\n3uNQR5WUUBoiRjgDfccgNbzsWCLEvus4JqRaTbkEoBmhnxGjC1ag5KPAKVQFCkKvphxFSjq+kVp7\nA62ndeHt7IkCw0KjpCSnyKFxk/rnoNPgbkzRdFZalN5Wtw9Msl1j1sTasseKU6YkxFYtSvCBx6/I\nc9F7TGLKlPg6TywKPPO43qN5SjzKYWrk4EjBFJhnm4oS4gM5CPsdBBK9lMWqlNQImkNHck/Zq9yV\nKSdZ+ZZ7ZyxdMzcBDhRX36ElmcEuappi3+zKnU2ZWLhGimadegclNAexqcAmzNMbGRslq7/ZFGnq\nz5PgoGiDENoKzwK92tBQoN7fvJTJNPk4KP1ZgO3kiKBoJTSajraO80hLy4m7LS0FFPeRCcdJXOSi\ndKYBfsS3gnzHu4rld3x6Ok/0MoWyQ7nZO4osfiRtv0G9Y0fREgQ2GjEgauFRMYrsIGUcbMotuEJa\nzx2NQeg+NAwGDl8wyJoYoR4bZQchcW9T5kQYX1hrxdJrDsDE0GWANvN7CS1446JRdpAiDTblPfD1\neN5d6RnLBl99IPmOVFjZQU5tCkCnjrrjrLdikU7QjsIKTT5ZGtCmbBEgaNgTelhQhOyldOcBnpDj\nRdtbkR0itPFUgd28lODQTBFBeElZtlfKDsLwjQaNfcs3LdfYQzAjRX92YENgMICEz8LJG4tcikOH\nECtkBzb05Np1iI8KJzVFvn8VsgOvfEyfUVLfiSLstfXqWGpb747FkQPqFRG0t7HUE7/f95+M5bXY\nx1f/dlCj2JGv+ONbGl8ay0PYqIRRUHzFv8GAsxDZAPvzcsEUFF/xxy/4KsBDgNuu5XJSK2UO4BHA\njTyc2wUlPosZWSdKQ7oCpkoszuKfKPHaBR2SQnKj1BR38V+mZYFtTH0ZkjxNTXEVfy1L67QCVoRZ\nkjxNTXEV/yifxhIuDq17CDRVU1zFnyl8OHVSKBa7+POOOiBuQGC2rliizuIfW8a9lCBh1He6nuJP\nJWx8qIg/s5KDtcfiK/64jzbT/IU/bPK1UyxPu/jjFyDauCU4R9UdSR+7+HM9JWVdX6OkdGP8XPxj\nhtDUlxkkxS7+mKW4WTRpl8MOOyo7/03GLK+uSao5lqq0/5TKx0Pnsif/Vpyp/m9F5ac7tU25nPEH\nmzNHZ8UJAgUAAAAASUVORK5CYII=\n",
       "prompt_number": 12,
       "text": [
        "\u23a1   \u239b s\u2091    \u239e\u23a4\n",
        "\u23a2   \u239c2     1\u239f\u23a5\n",
        "\u23a2log\u239c\u2500\u2500\u2500 + \u2500\u239f\u23a5\n",
        "\u23a2   \u239d 2    2\u23a0\u23a5\n",
        "\u23a2\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u23a5\n",
        "\u23a3   log(2)   \u23a6"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So, for Euler's method where $s_e=1$ we find the interval must be:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"Acceptable interval is [{:.6g}, {:.6g}].\".format(log(1.0+0.5)/log(2.0), \n",
      "                                                        log(2.0**2-1.0)/log(2.0)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Acceptable interval is [0.584963, 1.58496].\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This matches well with the results above, and is completely independent of the data, and indeed the method: the only thing that matters is the expected convergence rate $s_e=1$. As we vary the expected convergence rate, the acceptable *measured* convergence rate also varies, and in fact the acceptable interval around $s_e$ gets larger:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s_e = numpy.arange(1,10)\n",
      "upper_limit = numpy.log(2.0**(s_e+1)-1.0)/numpy.log(2.0)\n",
      "lower_limit = numpy.log(2.0**(s_e-1)+0.5)/numpy.log(2.0)\n",
      "pyplot.figure(figsize=(8,6))\n",
      "pyplot.fill_between(s_e, lower_limit-s_e, upper_limit-s_e, facecolor='green', alpha=0.5)\n",
      "pyplot.xlabel(r'$s_e$')\n",
      "pyplot.ylabel(r\"Acceptable $s_m-s_e$\")\n",
      "pyplot.ylim(-1.5, 1.5);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAhIAAAGQCAYAAAD7rYnxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcnXV99//X58yWyb6zk8QNREAoFmxBDVZt0aq3S7Vq\nF/uz3txdbO3Pta0VrNXKrdW72mptq7e1/VWptlZFVECiQqSoQFgEkhQzSSBAAtkmmSWzfH5/XGcm\nk8mEzDlzZs6cmdfz8Thc51zbfE4S5rzP9/p+v1dkJpIkSdUo1bsASZLUuAwSkiSpagYJSZJUNYOE\nJEmqmkFCkiRVzSAhSZKq1vBBIiJOiohvRcRgvWuRJGm2aa53ARMREa8EPgIcAiqaECMiOoA9Y2x6\nW2beOPHqJEma+Ro6SABvA54PXAk8rcJjMzPPr3lFkiTNIo0eJC7JzIyIetchSdKs1NB9JNL5vSVJ\nqquGDhITFRFXRcT6iNgYEd+OiJfWuyZJkhrJbA4SO4HbMvNi4BnAV4GvRsTv1bcsSZIaR8yEqwMR\n8TngNzJzQsEoIq4BngOszMzeWtQmSdJM1uidLWvth8CLgbOAO4ZWRkTjpy1JkiqUmccdzTArL21E\nxJyImDfGpoHysmn0hsycsY8rrrii7jX4/nx/s+29+f4a/zHT3994zaQgccx3HREnxJFjRH8V+Ksx\ndr0A6AHurXFtkiTNSDMpSIzZ/BIRFwM7gL8dtel1EfGsEfu9Fng58L8zs2vSqpQkaQZp6D4SEfFJ\n4DJgOZARsYWiZeKMzOwr79ZJMRX2jhGHXgucCnwyIlqAxcBu4PLM/Mepqn+6WLt2bb1LmFS+v8Y1\nk98b+P4a3Ux/f+M1I0ZtTLaISP+cJEmzSUSQdraUJEmTySAhSZKqZpCQJElVM0hIkqSqGSQkSVLV\nDBKSJKlqBglJklQ1g4QkSaqaQUKSJFXNICFJkqpmkJAkSVUzSEiSpKoZJCRJUtUMEpIkqWoGCUmS\nVDWDhCRJqppBQpIkVc0gIUmSqmaQkCRJVTNISJKkqhkkJElS1QwSkiSpagYJSZJUNYOEJEmqmkFC\nkiRVzSAhSZKqZpCQJElVM0hIkqSqGSQkSVLVDBKSJKlqBglJklQ1g4QkSaqaQUKSJFXNICFJkqpm\nkJAkSVUzSEiSpKoZJCRJUtUMEpIkqWoGCUmSVDWDhCRJqppBQpIkVc0gIUmSqtbwQSIiToqIb0XE\nYL1rkSRptmnoIBERrwTWA6uBrPDYloh4f0TcFxF3R8T6iLh4MuqUJGmmauggAbwNeD7wX0BUeOwn\ngF8BLsnMc4DPAtdFxDNrW6IkSTNXoweJSzKzo9KDIuIM4M3AhzLzcYDM/AywBfhATSuUJGkGa+gg\nkZkVXc4Y4RUULRjrRq1fB7woIuZOqDBJkmaJhg4SE3AuMABsG7V+C9AMnDXlFUmS1IBma5BYDnSN\n0aKxv7xcNsX1SJLUkJrrXYAkqXFl5vBjcHDwiOXQY/R+T/S62m312Hfo+cg/i2Otq8W+1Rx3RM0c\n/Z6OtW9LSwvjNVuDxGPAvIiIUa0SC8vLx0cfcOWVVw4/X7t2LWvXrp3M+iRNU4ODgwwMDBzxGGvd\nEz1G7j/yQ3fo+cDAAAOD5f3GWGbm4dcDAwzm4e0j9x3MweF1g3l43+F1I7aN+Xys1wwyOHA4MAxm\neQqfgIggI4kIImJ4/RGG9iuP2B9+HhDFf47Yl6AY3D/q+dD5R56HKD4ERz4fz74ZefhnjKg5SYKj\n6zviODjy2PLxQx8tQ+9vuIbMI+sZek9j/FkN/fzh5yPOMdbPG3ne4T//0X/8Y/y9DK17ePPDPLL5\nEQAG9g6MefxYZmuQuBP4VeA0juwnsQboA+4dfcDIICFpcgx9iB46dGj4Ue0Hdn9/P339ffQP9HOo\n/9Dw65GP/oHyPuVtAwMD9A/009/fP/ZyoL/4JV6CKBUfKlGK4nmJYv3QB0wUr0cuM5IolbeXDn+Q\nDH+gxOEPgYjyeSn/jKEP5/IH4+jnR+zbXMG+Y2wvRYmmaBpz32Mdp8Z3+jmnDz/fft127v6vu8d1\n3EwKEsccwRERJwA7R7Q+fAX4IHAp8E8jdr0UuC4zuyatSmmGGBwcPOIDf6xHX18fhw4dorunm66e\nLrp6uornvcWyu7ebnp4eeg710Huol57enuJDtQmiKYimOPyhPcaH8vAHdkCW8oj1UQpKTaViWSod\nfl3+cCw1l4jWUdtG7NvU1ERzqXn49chtkg6bSUFizP+7y7NVfh/4NPC7AJm5KSL+HvjjiLgmMx+P\niN+iaJF4/VQVLE2FzBz+QH+iD/tDhw7R29s7/IE//MHf2013Tzc9vT1093Zz6NAhenp76BvoG/6w\nj+aAJo58lCCbksHSINEcNDU3FY+WJprammiaf/j1nOY5zGuZR6mpRKlptvYBlxpTQweJiPgkcBnF\nKIyMiC0ULRNnZGZfebdOYA+wY9ThbwGuANZHRB/FiI0XZeZdU1K8VIGBgQEOHjzIgQMHhh+dnZ08\nvu9xDnYdHPMbfk9vsTx06FDxrb4poLm8HPGBn6UcXmZTFh/0Iz/05zTRtLx43dLSwpzmOSxpWXJE\nU7ek2Suqn9Np9ji6T6Y0cZlJT0/PEeHgwIED7N23l8f2Psbj+x5n997d7Nm/h86uTkqtJaItoBWy\nJRlsGaR5TjPNrc2HP/THWJaaS37gS6rI9uu289m//CyZedxfHg3dIiFNR/39/UeFg87OzsPhYN9u\ndu/dzb7OffTTT6mtVISD1iRbkmgLWttbaZ3fSuuKVha0L2Bp21KvzUualgwS0jhkJl1dXUcFhD37\n9gwHhD379rBn/x66uruItqDUVhoOB9maNM9pLgLCSa20PqmVE9tPpKm5qd5vTZImxCChWa2vr4/O\nzs4jwsH+zv1FONhbtB7s2b+HfZ37GCwNDrce0AqDLcXr1vZWWhe10npiK4vbF7O8bbmXEiTNGgYJ\nzVg9PT08/PDD7N27lwMHDrB7324e2/vY8KWFvfv30t3XTamtGAY4FA5ohZb2liIgnNZKa3srp7Sf\n4mgCSRqDQUIzQnd3Nw8//DA7duxgU8cmNm/dzCO7HyHmB7QXAaFpThOt7a20LWuj9dRWlrUvo6ml\nydYDSZoAg4QazlBoeOihh9i8dTMbOzayc89OYkGQC5K2RW0sOGsBpy863Q6KkjTJDBKa1rq6uoZb\nGjZ2bGTz1s3s3LuT0oISg/MGmbNkDvPPns/pCw0NklQPBglNG11dXezYsaMIDVs2snnbZnbt20XM\nL1oa5iyew4JzFhShwcsRkjQtGCRUFwcPHhy+PLGpYxMbOzbyeOfjxeWJ+cmcJXNYcK6hQZKmO4OE\nJt3BgwfZsWMHD+14qGhp2LqZxzsfp7SwuDzRvrSd+efNNzRIUgMySKimDhw4wMMPP8yDDz3Ipo5N\nbOrYxO4Du4dbGtqXtrPg/AWcvsDQIEkzgUFCVTtw4MBRLQ27D+4e7gjZvrSdBT9jaJCkmcwgoXHp\n7Ow8qqVhT9eewy0NS9pZcMECTp9vaJCk2cQgoSNk5hGhYWPHRjZ1bGJf1z5Ki4qWhrlL5zL/WfMN\nDZIkg8Rs19PTQ0dHRxEaykMu93XvKzpCzi9Cw4ILF7Bo3iJDgyTpKAaJWWrPnj18f/33+db6b9E7\ntxcWUFyeMDRIkipgkJhFMpPt27dz/feuZ/3d64mTg5WXrKRtXlu9S5MkNSiDxCwwMDDAfffdxzXf\nuYb7H7mf1lWtnPLCU2hqaap3aZKkBmeQmMF6enr48W0/5mvf+Ro7B3cy/0nzWXXuKi9bSJJqxiAx\nA43s/9C9sJtlZy1j9fLV9S5LkjQDGSRmiNH9HzgJTrjkBE6Yd0K9S5MkzWAGiQY3ODjIvffeO9z/\noeX0Fvs/SJKmjEGiQR3R/2FgJ/OfPJ9V56wiSvZ/kCRNHYNEg9mzZw83/eAmvnnzN+lZ2MPSs5ba\n/0GSVDcGiQaQmTz44INc/73rufmum4f7Pzj/gySp3gwS09jg4GAx/8ON13Dfjvvs/yBJmnYMEtPQ\nmP0fnmH/B0nS9GOQmEZG9n/oXuD8D5Kk6c8gMQ0Mzf9w8103w4mw8uKVnDDf+R8kSdOfQaJOjur/\ncJr9HyRJjccgMcV6enq47fbb+OoNX7X/gySp4RkkpsjevXuL/g83fZOuBV0sffpSVq9YXe+yJEma\nEIPEJBvu/3BnMf/DyotXsnL+ynqXJUlSTRgkJsHI/g/3PnQvbavaOPmFJ9Pc6h+3JGlm8ZOthnp7\ne4fnf3i0/1HmPWkeq5+x2v4PkqQZyyBRAyP7P3Qv7GbJmUvs/yBJmhUMEhPw4IMPcsP3buD7G75v\n/wdJ0qxU0yARES8FlgBfAE4Clmfm7bX8GfU21P/h2huv5ScP/YTWVa32f5AkzVq1/vQL4GrgpZn5\nHxFxETAjgsRQ/4evf+frPNL/CPOeNM/5HyRJs16tg8SzgRuAg+XXXTU+f9380ZV/VMz/cKbzP0iS\nNKTWQeIaYANwf0Q8GTgF+EaNf0ZdLLx4of0fJEkapVTLk2XmD4BLgeuBQ8AHa3n+epozf069S5Ak\nadqpqkUiIl4CvA74XaANuARYn5k7M/Mh4BO1K1GSJE1X1V7aeBpwEzAXuBXYBTwaEVdk5o9rVdzx\nRMRK4GPABeVVdwNvLYeZ4x3bAewZY9PbMvPGmhUpSdIMVnUficz8dES8HFgJnAfsAz4ETEmQiIhW\nikso9wNnlVd/FlgXEedn5sFjHlzIzDx/MmuUJGmmq7aPxMkRcQHwRuDWzNyTmYPAtppVdny/CZwD\nvCszB8s//13Ak4DfmcI6JEmataoNElcBfwysBt4JEBGXcXjY51R4FbA1MzuGVmTmo8C95W2SJGmS\nVRUkMvOxzHx1Zp6fmT+MiCXAl4FVtS3vCZ0LbBljfQdFS8VxRcRVEbE+IjZGxLfLM3NKkqRxqsnw\nz8zcA6zIzCtrcb5xWg50jrF+PzA3ItqOc/xO4LbMvBh4BvBV4KsR8Xu1LVOSpJmrZvNIZOZUz2KZ\nEzo486LM/Lfy8/7M/CRwLfDBcYQQSZLEBGe2jIjWzDxUq2Iq9BiwYIz1C4GDmdlbxTl/CLyYYhTI\nHSM3fPdz3x1+vvq81aw+b3UVp5ckaXrq2NBBx4YOAPY9sG/cx010iuzPAa+f4DmqdRdw5hjr11DM\nJ3FMETEHaBpjiOhAedk0+pi1b1xbRYmSJDWGkV+St1+3nQ03bxjXcRO9tNEyweMn4j+AVREx3MEz\nIk6gCBf/PnLHiDghIkbepvNXgb8a45wXAD0UIz8kSdJx1PReG1PscxQtD1dFRFNElCgmxPop8Kmh\nnSLiYmAH8Lejjn9dRDxrxH6vBV4O/O869PeQJKkh1frun1MmM/si4oUUU2TfS9H58m7g+aOCQCfF\nVNg7Rqy7FjgV+GREtACLgd3A5Zn5j1NRvyRJM0HDBgmAzNwJvOE4+9xFMVR09HF/UX5IkqQqNfKl\nDUmSVGcGCUmSVDWDhCRJqppBQpIkVW2iQeJbNalCkiQ1pAkFicz8TK0KkSRJjcdLG5IkqWoGCUmS\nVDWDhCRJqppBQpIkVW3SgkRE/Fz5PhaSJGmGqmmQiIjXRsRHI+LVwP0Ut+uWJEkzVK1bJE4Cvgis\nAb4CXFrj80uSpGmk1nf/3JaZPwR+CHy4xueWJEnTTK1bJBZFxD9GxNqIaK3xuSVJ0jRT6yCxBvgG\n8MvAuoj4jxqfX5IkTSO1vrSxITO/QtE/goiYU+PzS5KkaaTWLRInRsS7I+LJAJnZU+PzS5KkaaTW\nQWIZsA/4YETcFhGfqvH5JUnSNFLrSxvfAeZm5qcAIuL0Gp9fkiRNI1W1SETESyLiXyJiYUSsiIhX\nRMTKzPxBZt4wtF9mbqtdqZIkabqptkXiacBNwFzgVmAX8GhEXJGZP65VcZIkaXqruo9EZn4auAhY\nCbwQeCnwmhrVJUmSGkC1LRInR8QFwBuBWzNzD0BEeClDkqRZpNogcRXwd8Bq4HKAiLgMOFibsiRJ\nUiOoKkhk5mPAq4deR8QS4Mt4fw1JkmaVmgz/zMw9EbEiM7tqcT5JktQYajYhlSFCkqTZp9YzW0qS\npFnEICFJkqpmkJAkSVUzSEiSpKoZJCRJUtWqGv4ZEc8D1gLtmfnuiFgL3JaZnTWsTZIkTXMVtUhE\nxIKIuAFYB1wB/EZ502XAXRGxpsb1SZKkaazSSxsfAtopgsPpwE6AzHwX8Mfl7ZIkaZao9NLGZcAz\nhy5hREQObcjML0bEO2pZnCRJmt4qbZHoO04/iMUTKUaSJDWWSoPEwYh49VgbIuIlwOMTL0mSJDWK\nSi9tvB/4ckSsB9YDKyLiPcB5wEuBV9W4PkmSNI1VFCQy8ysR8XqK24VfUl7958A24PWZeU2N65Mk\nSdNYxfNIZObVEfFvwBnAcuAxYGNm5hMfKUmSZpqqJqQqh4b7R6+PiE9m5u9OuCpJktQQnjBIRMRv\nAuNtaQjgJROuqAIRsRL4GHBBedXdwFsz86FxHNsCvBd4NdAP7AfemZnrJ6lcSZJmnOO1SPzfCs83\nZZc3IqIVuJ6iZeSs8urPAusi4vzMPHicU3yCYprvizPz8Yh4E3BdRPx8Zt45WXVLkjSTHC9I3Ae8\nmKK1YTy+MbFyKvKbwDnAyzNzECAi3gU8BPwO8JFjHRgRZwBvBt6UmY8DZOZnIuKPgA8AvzzJtUuS\nNCMcL0h8PDO3jvdkEfHxCdZTiVcBWzOzY2hFZj4aEfeWtx0zSACvoAhH60atXwdcHhFzM7OrxvVK\nkjTjPGGQyMxPj7U+Ip4KXAicDOwAfpiZm4+1/yQ5lzE6fAIdwPPHcewAxbDVkbZQ/JmcBfx4gvVJ\nkjTjVTRqIyLmAp8GXseRs2IORsS/ApdnZncN63siy4GxpuveD8yNiLbM7H2CY7vGGLK6v7xcVqMa\nJUma0Sod/vnXwPOAdwO3A3uApRSjJn4f+DhF34OpMKXzVnz3c98dfr76vNWsPm/1VP54SZImVceG\nDjo2dACw74F94z6u0iDxCuD8zNw+av13IuKLFOFiqoLEY8CCMdYvBA4+QWvE0LHzIiJGtUosLC+P\numfI2jeurbZOSZKmvZFfkrdft50NN28Y13GV3rRr2xghAoDM3AY8WOH5JuIuYM0Y69dQzCfxRO6k\neO+njXFsH3DvhKuTJGkWqDRI3B4R5421obx+rM6Pk+U/gFURsWpEDScAZwL/Pqq2EyJi5BDWr1Bc\nGrl01DkvBa5zxIYkSeNT6cyWd1Dc/fN64CcUnRMXAmcDL6O4mddU+RxFv4yrIuIN5To/BPwU+NTQ\nThFxMfB9ik6ivwuQmZsi4u+BP46Ia8oTUv0WRYvE66fwPUiS1NCqndny8mOs/wjFlNWTLjP7IuKF\n5Z93L0WQuBt4/qgWhU6KTqE7Rp3iLcAVwPqI6KMIRS/KzLsmvXhJkmaIRp7ZkszcCbzhOPvcRTHc\nc/T6fuDPyg9JklSFWs9s+YkJ1iNJkhrIE3a2rGKmys9MoBZJktRgKh21cTy31vh8kiRpGqt0Qioi\n4gLgTRQjHNpGbgKeWqO6JElSA6ioRSIiLqMYSnk+cAlFeChR3LzreUztPBKSJKnOKm2RuAJ4QWbe\nEhF3ZObwhE4R8SvAz9W0OkmSNK1V2kdiTmbeMtaxmfklipYKSZI0S1QaJAZGPO+LiJOHXkTEYorp\nqSVJ0ixRaZDYERHvj4g5wPeAb0fEWyPireXX3uxKkqRZpNI+Eh8FXgOsAD5I0Sfio+VtPwVeV7vS\nJEnSdFdRkMjMdcC6odcRcQnwFOAEIDLTFglJkmaRSod/3jzydWYOZuam8svPR8SnxjhMkiTNUJVe\n2pg31srMvCkingzcOfGSJElSozhukIiIVcAqismn5kfEc8faDTgNWFDb8iRJ0nQ2nhaJ3wLeO+L1\nd4+xXwLvn2hBkiSpcYwnSHyOw+HhH4DfpmiBGKkP6MjMh2pWmSRJmvaOGyQyswPoAIiI/5OZ35vk\nmiRJUoOodPjn3wJExFOBCylu1rUD+GFmbq59eZIkaTqrKEhExFzg0xQTT40cOjoYEf8KXJ6Z3TWs\nT5IkTWOVDv/8a4rbhb8buB3YAywFLgB+H/g48OZaFihJkqavSoPEK4DzM3P7qPXfiYgvUoQLg4Qk\nSbNEpTft2j5GiAAgM7cB2yZekiRJahSVBol7I+K8sTaU198+at0nqi1MkiRNf5Ve2rgNuDYivgr8\nBNgPLATOBn4B+FhE/EZ53wBeBrylRrVKkqRpptIg8ZHy8vJjbP+bUa+zwvNLkqQGUmmQuA94MUfP\nbHks36jw/JIkqYFUGiQ+nplbx7tzRHy8wvNLkqQGUlFny8z89GTuL0mSGkulLRIARMTzgLVAe2a+\nOyLWArdlZmcNa5MkSdNcRS0SEbEgIm4A1gFXAEMjNC4D7oqINTWuT5IkTWOVziPxIaCdIjicDuwE\nyMx3AX9c3i5JkmaJSi9tXAY8c+gSRkQMD+/MzC9GxDtqWZwkSZreKm2R6DtOP4jFEylGkiQ1lkqD\nxMGIePVYGyLiJcDjEy9JkiQ1ikovbbwf+HJErAfWAysi4j3AecBLgVfVuD5JkjSNVRQkMvMrEfF6\n4MPAJeXVf05x18/XZ+Y1Na5PkiRNYxXPI5GZV0fEvwFnAMuBx4CNmel9NSRJmmWqmpCqHBruH3od\nEQsp7gQqSZJmkUonpPq1iNgTEdtGbfp2RPxzRMypYW2SJGmaq3TUxq8DnwGePmr9S4CDzOAJqQ7s\nPlDvEiRJmnYqvbRxEvBLo/tDZObuiPh94LaaVTbNDGwYoKOtg5Vnr2Tuorn1LkeSpGmh0haJtmN1\nqszMforps2ekq95zFW+85I10/bCLrbdupfdgb71LkiSp7ioNEvsj4mVjbYiIXwb2Tryk8YuIt0bE\nTyLizoi4LSJePs7jroyIrRFxx6jH/znWMa2trax93lo+/J4P85rzXsOem/aw9batHOo+VLs3JElS\ng6n00sb7gK9ExDrgR8BuYAnws8ClwCtrW96xRcS7gbcBF2bmloh4AXBtRLwsM791nMMT+LPM/Hyl\nP7e9vZ3LfvEyLvn5S7j+xuu5dt215CnJSU8/iebWqgbBSJLUsCqdkOqaiHgd8FfAC0Zs2ga8bqom\npIqIxcCfAR/OzC3l2m6IiOuAjwDHCxIAMZEaFixYwCtf/koufe6lfPP6b3LDDTfQtLqJE592Ik3N\nTRM5tSRJDaPSSxtk5peB1cBZwHPLyzWZ+e+1Le0J/RJFf4x1o9avA86KiDOmqpAlS5bw+te8ng+9\n/UNcNP8itl+3nYc3PszgwOBUlSBJUt1UHCSgmJAqM+/PzJsz835gQY3rOp5zy8sto9YPvT5nHOf4\npYi4ISLuKveveF9EVN1ZdOXKlfz2b/w2f/EHf8EzeAbbbtjGzi07ccJPSdJM1qgTUi0vL0ff0nxo\nds1lxzm+CzgAvDIzzwXeBPwa8J2ImFBHh1NPPZU/uPwPeO+b38uqzlV03NDB4w8+bqCQJM1I02JC\nqoh4QUQMjuNxYzXnHy0zP5yZb87M/eXXG4B3Ac8GXlOLn/GkJz2Jd77lnbz7197Nsh3L6FjXwd5H\npnRQiyRJk266TEi1HjhzHPt1lZePlZcLgD0jti8sLx+vooYflpcXAf86euOVV145/Hzt2rWsXbv2\nuCeMCJ7+9Kfz3jPey1133cUXr/kiWzZuYfnZy1mwbKqvBkmSdGwdGzro2NABwL4H9o37uEqDxBNO\nSFVtH4PM7AY2VXDIneXlaooRI0PWlJd3PdHBEbEiM3eNWj1QXo455GJkkKhUqVTivPPO45xzzuFH\nP/4R/3btv9HR6iyZkqTpY/V5q1l93moAtl+3nQ03bxjXcY06IdW3KFonLh21/lLgJ5k5HEoiYm5E\nLBq139aIGP3eLygvb69ppSM0NTXx7IuezYf+9EPDs2R23NpBz4GeyfqRkiRNqoackCoz90XE+4G3\nRcTnR0xI9SLgpaN2vwNYEhGryi0fAHOA90XEFZk5GBGrKPp33M8YlzVqbWiWzIsuvIjv3fQ9/vPG\n/6R3RS8nnXUSre2tk/3jJUmqmYackKpcy1UR0QNcExH9FJcmXp2Z3x616w7gENA/Yt0bgNcDGyKi\nCZgLfJNitsspax5ob2/nl170S8OzZH5j3TfIk5OTznKWTElSY4hqhiVGRABnUAzDfKw8lwQR8YnM\nfEttS6y/iDhW15Ca2rNnTzFL5o9voGlVEyee4SyZkqSpt/267Xz2Lz9LZh53FuiqgsRRJyn6G7wE\n+IfMPHHCJ5xmpipIDNm5cydf//bXuemem2h7UhsnPOUESk1VzR0mSVLFKgkSE/p0iohnRsRHKS4f\nfBVYMZHzqbBy5Ure9Otv4gN/8AHOjrPZev1WHv3po05qJUmadiq+EB8RKyn6GPwmxVTVfcBNwDeA\n361pdbPcKaecwlv+51vYsmULX/r6l/jJ9T9h4dMXsvTUpRRXlyRJqq9xBYmIaAVeRhEefrF83Pco\nJoZ6ytAMkeVOj6qxNWvW8I63vIONGzdy9dev5qebf8qSs5aw+MTF9S5NkjTLHbePRER8EngtxTDP\nPcA/AX+XmZsiYktmrnnCE8wAU91H4okMDg5y9913c/U1V/NQ/0POkilJqrlK+kiMp0Xicoqhle8A\n/iYzeydaoKpXKpV45jOfydlnn82Pb/sxV197NVuat3DCOSc4S6YkacqNJ0icTnFnzDcA50fE32Xm\nzZNblo6nqamJiy68iPPPO59b/usWvvTtL7Fz4U5OfMaJzJk/VTdhlSTNdhUN/4yIn6XoJ3EO8GXg\nHZl5+ojtTxs5PfVMMZ0ubRxLT08P3/3+d50lU5I0YZM+j0S58+UvU9xWvEQxK+Q1wNcy82cqPuE0\n1whBYsiBAwe4Yd0NXHPzNQyeNMjJZ53sLJmSpIpM6YRUEbGC4rLHbwBnZ+aM+xrcSEFiyN69e/nm\n9d/k+h9dT+n0EiedeZKzZEqSxmXKZ7YcPlnEnZn5zJqdcJpoxCAxZNeuXXztW1/jprtvovXJrZz4\nlBOdJVN9elyXAAAV8ElEQVSS9ISmbGbLMYy+rbfqbMWKFcUsmX/4AZ7Z9Ey23bDNWTIlSTVT0xaJ\nmaqRWyRG27JlC1++5svc8/A9LDhzActOW+YsmZKkI9R6HgnNIGvWrOHtv/92Nm3axNVfv5oHNj/A\n/KfMZ+nJS2lqsQ+FJKkyBolZKCI444wzeM9T38M999zDjbfcyN3X3U0uTJpXNrPstGW0zW2rd5mS\npAZgkJjFSqUS5557Lueeey49PT3893//N7fffTu33nwr3c3dxPJgyalLmLdknpc/JEljso/EOMyk\nPhLjMTAwwLZt27j73rv5wR0/4JHOR4jlwYKTF7D4xMWO+pCkGc4+EpqQpqYm1qxZw5o1a3jpi1/K\nrl27uO/++7hlwy1s2rCJXJLMOWEOS09dSktbS73LlSTVkUFCTygiWLlyJStXruR5z30eBw4cYPPm\nzfzozh9x24230dfeR2lFiaWnLvWmYZI0CxkkVJH58+dz/vnnc/7559PX18eWLVu48547ueX2W9ja\nt5VYHiw6ZRELVywkSvarkKSZzj4S4zDb+khUIzPZsWMH9953L+vvWM/WXVthKcw9aa5DSyWpwdhH\nQlMuIjjllFM45ZRTeOELXsjevXvZuHEjt955K3dfdzcDCweKoaWnLqNtnkNLJWmmMEhoUixevJiL\nLrqIiy66iJ6eHh544IFiaOn6W3m06VFYDktOXcL8pfMdWipJDcxLG+PgpY3aGRwcZNu2bdxz7z2s\nv2M9j3Q+AkthwSkOLZWk6cJLG5q2SqUSq1evZvXq1bzkspfw2GOPDQ8t3XjHRnJJ0nZCG8tOXUbL\nHIeWStJ0Z5BQ3UQEK1asYMWKFTz3Oc/l4MGDbNq0iR/f/WNuW3cbh9oPEcuDZactc2ipJE1TXtoY\nBy9tTL3+/n62bNnCXT+5i/Ub1rOndw8sh0WnLGLRykUOLZWkSVTJpQ2DxDgYJOorM3n44Ye59757\n+cEdP2DLo1tgGcw9cS5LTl5Cc6sNa5JUS/aR0IwSEZx88smcfPLJvOAXXsC+ffuKoaUbbuXO6+9k\ncMEgTSubWHbqMubMn1PvciVpVjFIqOEsWrSICy+8kAsvvJDe3l4eeOAB7rjnDm75wS08Wnp0+K6l\n85c5tFSSJpuXNsbBSxuNYXBwkO3bt3PPvffwgw0/YMfeHbAM2pa2MW/xPOYumusMm5I0Dl7a0KxU\nKpVYtWoVq1atGh5aunHjRjZv3cyWLVt46NGHGGwZJOYFg+2DtC1qY+6iucxdNNd+FpJUJX97asZa\nvnw5y5cv5+KLLwaKFos9e/awa9cudu7cScdDHWzZWgSM/qZ+mAc5N2ld0Mq8JfMMGJI0Dv6W1KxR\nKpVYtmwZy5Yt48wzz+S5PBcoAsbevXuHA8bWHVvZsn0LD/7oQfpKfcMtGC0LWopLJIvn0tLmZFmS\nBAYJiVKpxNKlS1m6dClnnHEGz+E5QDHsdN++fUe0YHTs6GD77dvpzV5K80sMzi0CxtxFc5m3eJ6z\ncUqadQwS0jFEBIsXL2bx4sU89alP5WKKSySZyf79+9m1axe7du1i245tbHloC9s2bKN7sJuYF2R7\n0rygefgSScucFkeQSJqRDBJShSKCRYsWsWjRIp7ylKfwc/wcUASMzs7OowPGndvo6u8qWjDaB2ma\n3zR8iaS1vdWAIamhGSSkGokIFi5cyMKFC3nyk5/Ms3k2UASMgwcPHg4YD21jy44tbPvJNg4cOkBp\nXomcm5Tml4YvkbTONWBIagwGCWmSRQTz589n/vz5rFmzhgu5cHjbyIDx4MMPFi0YG7fx8MGHKc0/\nHDCG5sFom9dmwJA0rRgkpDqaN28e8+bNY/Xq1fwsPzu8vqurazhgPPTIQ/z0wZ+ybdM2Hj3waBEw\n2hPmMdx60TqnlZY5LZSaSnV8N5JmI4OENA3NnTt3eHKtZ/Gs4fXd3d089thj7Nq1ix2P7OCnD/2U\n3Y/uZt/+fTx64FGylJTaSkRrkC1JtiSDLYO0zGkZDhut7cWyubXZ1g1JE2aQkBpIe3s7p512Gqed\ndho/w88csS0z6enp4eDBg0c8Dhw4wO59u9m9fze7d+1mb+de9nbupauni2gNojWgleHgEa0xZvBo\nanZ6cUlHa+ggEREl4O3A+4H/mZn/VOeSpLqJCNrb22lvb2f58uXH3X9gYOCo0HHw4EH27d9XhI79\nu9m7qwgdj3c+Tn/2E21BtJSDR3PR2tE8p/no4NHWQpRs7ZBmg4YNEhFxOvB5YB7QAlR0V62IeBHw\nF8Cc8vH/BFzl3bk0WzQ1NQ2PMjmezOTQoUNjtnbs2b+H3ft3s2ffHvZu38u+zn10dnVCE5TaSsX/\nXSNaPEa2cgyFj6aWJi+zSA2qYYME8FbgU8AjwLpKDoyIS4CvA7+SmV+LiFOBW4GFwJ/UulCp0UUE\nbW1ttLW1sXTp0uPuPzg4SHd391HBY3/n/sPB45E97O3cy2P7H6Onr4dSa+nwZZZya0eptURTSxNN\nzeVH+XmpuXTE66bmJltApDpp5CDx9swcjIi1VRx7FXBLZn4NIDMfjIiPAR+IiE9k5sO1LFSabUql\n0vCIlPHo6+ujq6truJVj+DJL5z66erro7umm+1A33Qe66T3US3dvN9293fT09hSPQz0QFK0gTSVo\nKp5Hc5ClLF6XYLA0SJaSaI7DYeQ4AaXUXLK1RHoCDRskMnOwmuMi4iTg54D3jdq0jqIR9mXApydW\nnaRKtLS0DM8WWo3MZGBggEOHDnHo0CH6+vqOej5yXe+h3sMBpbebrp4ueg4WgaS7pwgo3YeKZV9/\nXzEaphxQoimOCihZSrKpeERTjCugNLUUrSiGFDW6hg0SE3BOebll1Poto7ZLahARQXNzM83Nzcyd\nO7em587MMcPIsZ739PYUIaW3ezio9HQeDig9h3ro7O2kt7eX/oF+CIrLMqXifQw9p0TRyhKHnw9t\ny0iIw8uh7RlJkkXwKZWIUlBqKi/Lr0c+H8+2iKP3k0aajUFiqDt756j1+8vLZVNYi6RpLiJobW2l\ntbW15uceHBxkcHCQgYGBIx5jrat026G+Q/T199E/0E9/fz99A3309fUVy/4++vv66e/pP7x9xL79\nA/0MDAwMPx9eDvSTmYfDThR/PkPBZ2QgGgo5EUFSHJPk8OuhkDRy3fC28rxqI58f9XroZ5eXwyEM\njnxeDkOj9x3zuBEtROM5brgGDi9HrhteHxy5Dxx1/BPtO91NiyARES8ArhvHrt/NzOdPdj2SNBVK\npRKlUonm5mnxq3hcjhV+xgozQ/tm5lHLsdZVsm1wcJCBwcM/bzCL1wODA+RgFtty8HBNWSyP2NZf\nLIfON/L50LmH1w+UlyP2yUwGcxCyCDk5WG4RGnperjcpL8fYDhx9nvJxR4ST4j/DISw4vA0OB7ah\n5wTFWMY4ctvo10cElqEWsISmQ+OfN2a6/OtdD5w5jv26avCzHisvF4xaPzQG7vGxDrryyiuHn69d\nu5a1a9fWoBRJaiyNGH4a0VCYONbz472uZt+bb76Z9evXA5W1iESjT5tQHrVxI/DGzPz8OPY/CXgI\nuDIz/3zE+guAHwH/KzP/ftQxTi8hSZpVIoLMPG6imPF3+ImIuREx3BW8PLTzFuDSUbteCvRRzC8h\nSZLGYSYFiWOlpjuAzRHRPmLdO4Gfj4iXApQnpHor8BHnkJAkafwa9tJGRDwf+AzFFNcrgd3AAeBP\nMvMLI/ZbRzFS42cys2/E+qEpstuAVuCfMvNDx/hZXtqQJM0q47200bBBYioZJCRJs419JCRJ0qQz\nSEiSpKoZJCRJUtUMEpIkqWoGCUmSVDWDhCRJqppBQpIkVc0gIUmSqmaQkCRJVTNISJKkqhkkJElS\n1QwSkiSpagYJSZJUNYOEJEmqmkFCkiRVzSAhSZKqZpCQJElVM0hIkqSqGSQkSVLVDBKSJKlqBglJ\nklQ1g4QkSaqaQUKSJFXNICFJkqpmkJAkSVUzSEiSpKoZJCRJUtUMEpIkqWoGCUmSVDWDhCRJqppB\nQpIkVc0gIUmSqmaQkCRJVTNISJKkqhkkJElS1QwSkiSpagYJSZJUNYOEJEmqmkFCkiRVzSAhSZKq\nZpCQJElVM0hIkqSqNXSQiIhSRLwzInoj4jfrXY8kSbNNc70LqFZEnA58HpgHtABZwbFXAr8F7B61\n6XuZ+dZa1ShJ0kzXsEECeCvwKeARYF2FxybwZ5n5+ZpXJUnSLNLIQeLtmTkYEWurPD5qWYwkSbNR\nw/aRyMzBetcwU3z3u9+tdwmTyvfXuGbyewPfX6Ob6e9vvBo2SNTAL0XEDRFxV0TcFhHvi4j2ehdV\nDzP9fwbfX+Oaye8NfH+Nbqa/v/Fq5EsbE9EFHAAuz8z9EXEe8O/ACyPiuZnZX9/yJElqDNOiRSIi\nXhARg+N43FiLn5eZH87MN2fm/vLrDcC7gGcDr6nFz5AkaTaIzHGPmpy8IopLCqeNY9euzHxw1LFr\ngRuBN05kFEZ5OGkH8InM/MNR2+r/hyRJ0hTLzOMOTJgWlzYysxvYNFU/LyJWZOauUasHysum0fuP\n5w9SkqTZaFpc2phMETE3IhaNWr01Ika/9wvKy9unoCxJkmaEmRQkjtVqcAewedSIjDnA+4bCRESs\nAj4E3A/866RWKUnSDNKwQSIinh8RW4AvUMxU+ZGI2BIRrxu16w7gUWDkSIw3AOcBGyLiJ8B3y4/n\nZGbPiJ9xUkR8KyKcs0LTVkTcVO6MfHq9a5E0+0yLzpbTUUS8EvgIcAh4amYe1XeiUZWHu/4ecDFF\nwGoCbgDen5mP1bO2WoiIJwO/A6wtr1pAESY/lJnX1quuyRARrwK+RBGm12TmtjqXNGERsRq4B9g8\nxua1mblvSguaBOW/tz8E5gJLKO7789eZ+S91LWyCIuJzFL9XDozatBQ4AViUmb1TXVctRcSzgPcD\nqyh+f/YAf5WZV9e1sBqJiOcC7wNOpbiP1U8oZpK+71jHNGyLxBR4G/B84L+YedNpfxFYDFyQmecC\nLwReBKyPiDl1raw2LgNeC7wmM58FnAncDHyt/D/JjBARrRSX5K5l5v0b/VFmnj/GYyaEiD+i+CD6\n9fK/zzMoOps/v66F1UYCbxr99wZ8E/jKDAgRq4HvADuBs8u/Pz8DfCEi/kcdS6uJiHgexfv7RmY+\nNTNXU4xm/H5EnHqs4wwSx3ZJZnbUu4hJMgi8qzxahszcAXwYeCrw4noWViMPAldk5k8Bsmh2u4ri\n3/vL6llYjf0ecCvwo3oXovEpfxD9JfCWzNwKUJ4A7+3A39avspr5NEVfs2ERMQ/4VeDv61JRbb2Y\nooXzo0O3acjMTwP7Kd5jo3s/8HBmfmTEundT3GX7T4910LQY/jkd5cy+5nPuGLN3PlxeLp7qYmot\nM/9zjNVDI3dGD/ttSBGxlOLD59nAm+pcjsbv1ymaw4+YXC8zH+bw/4MNKzP/a4zVrwEezcxK79I8\nHQ393mwZWhERQXF5uK8uFdXWsyhab4dlZmdEbKX4EvY7Yx1ki8QsdIwpwJ9G0Sz5/SkuZ9JFxCkU\n3/ZuY2Z86wN4L/DPmbm93oVMkhMi4p8j4taI2BgR/19EnF3vomrg5ylazF4cEesi4t6I+EFE/D/1\nLmwSvRn4h3oXUSNfoGhxeU9EzCuP/PsTihDxV3WtrDYOMHYuGAROjIj5Yx1kkBAR0UTxrfYfM/O/\n611PrUTEkyPiv4HtFH0IXpGZozuBNZyIeCrwK8AH6l3LJBmg+Ob30cy8iOJbUh9wa7mjWyM7rfy4\nEnhDZp5F0an70xHxJ/UsbDJExFnAzwCfq3MpNZGZncAvAO3AYxSduH8DeEn5VguN7g7g7IgYvlpR\nnodpTfnlwrEOMkgI4M+AXuCt9S6kljLzgcx8CsVljc3AnRFxcZ3LqoWrgL8s/1KbcTJze2aem5l3\nlF93Av8LOAh8sK7FTdwcig+hPy33TSIz/wP4KvAnM/AOxL8N/OdMGA0GEBFnUPRJ2gIsycwVFK2D\nN4wx9UAjeg+wDPhgRLSU/z1+jMOXdLrHOsggMctFxG8BrwYuG+p8OdNkZmdm/hHFt4dP1rueiYiI\n5wDPAP5urM1TXM6UKc/vcg9Fn5BG1klxCXH0t9cNFENBnz7lFU2S8qiiX2NmdLIc8n6Kb+V/ODTn\nUHnY57eBT0VEyxMdPN1l5o+AFwDnUgz7/B5wF/AvwMHM3DPWcXa2nMUi4teB/xd4/kz5xgAQEXNG\nTiw2wj3AqyKiJTMbtWPUCyg6dv2o6OMFwInl5bURcQj448z8Vj2Kq4WIWAh0j/F3NEDjf/m5H3gm\nR7+PoXv9NPr7G+l/AHszsyZ3bZ4mzgEeHGMY62aK97uasec/aRiZ+T2KADEsIq6lmAphTDPpH+1k\nmnEjOCLi14B3Ar+QmTvL6345It5c38pq4lsRMdY319XAvgYOEWTmFZn5lFFj9IdaJy4rr2vYEFH2\nceBVI1eUv92eQ+PfC+dr5eW5o9afDXRRfAucKX6bmdPJcsijwMnlfmUjraLokNjQo8Ii4vSIOHfU\nukXAJTzB36VBYnxmVJNxRLyB4h/F54AXRcSvlYPFS4GT61lbjSTFvVSWQjE8KyL+gKLT3sfrWtnk\niFHLRpfAOyLiRBjuDPxhimu376tnYTVwNcU19ivKLS9ExCUUwekDM+XyYvn+Rc8B/m+9a6mxT1Bc\n2vjzoRURcSnwCuBfMnNvvQqrkecDX4yIJTA8B8g/ANdn5peOdZBTZB9DRHySYobE5RTXLrdR/II7\no5G/0QJExOMU80WM/uBJ4H2Z+edHH9U4IuLnKb4N/SxFJ6E5FD2sP5mZX6hnbbVUnknvYxR/lwsp\n7ivTl5lPqmthE1Qe5nk5xQcRFP8P3kvxQfu9Yx7YIMq/pK8CfpGi81ov8PHM/ExdC6uhiHgfcGZm\nvrbetdRaRLyIYpKmEyl+vwxSfCn7m2MMrW8YEXEBxWy5TwX2UFxy+xLwkcwcOOZxBglJklQtL21I\nkqSqGSQkSVLVDBKSJKlqBglJklQ1g4QkSaqaQUKSJFXNICFJkqpmkJAkSVUzSEiSpKoZJCRJUtUM\nEpIkqWrN9S5AkoZExHLgXRQ3IAPYlJnfqGNJko7DFglJ08nfAV/MzI8BXcDr6lyPpOPw7p+Spo2I\nuAPooLgt861AV2bur2dNkp6YLRKSppPLgSXAFymCRFt9y5F0PAYJSdNCRCzJzB9m5lrgZGAXcEl9\nq5J0PAYJSXUXEXOBByPisvKqfcCjwPr6VSVpPOwjIWlaiIgrga3APOB04N8y88flbT8P/DLwX0Bm\n5tfrVaekIxkkJE1rEXEKsA44DxgAXpuZn69vVZKGOI+EpOnulcD9wHMoWiu+UN9yJI1kkJA03XUD\n38rMbwNExOnAtvqWJGmIlzYkTWsRMQ/4E+AHwFzggcy8vb5VSRpikJAkSVVz+KckSaqaQUKSJFXN\nICFJkqpmkJAkSVUzSEiSpKoZJCRJUtUMEpIkqWoGCUmSVDWDhCRJqtr/D4IPdtnkgjypAAAAAElF\nTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106b32d10>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If all you care about is the answer that the algorithm gives you, then the measured convergence rate doesn't matter very much. Use Richardson extrapolation to estimate the exact solution, given multiple numerical solutions, and check that the estimated error is in a range that you're comfortable with."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "More reading"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See also [this page from NASA on grid convergence](http://www.grc.nasa.gov/WWW/wind/valid/tutorial/spatconv.html). Also [this paper by Liu](http://ocw.mit.edu/courses/mathematics/18-304-undergraduate-seminar-in-discrete-mathematics-spring-2006/projects/xtrpltn_liu_xpnd.pdf), which is rather mathematical, on Richardson extrapolation and its extensions."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}