lessons/05_relax/05_03_Iterate.This.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, G.F. Forsyth.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Relax and hold steady"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ready for more relaxing? This is the third lesson of **Module 5** of the course, exploring solutions to elliptic PDEs.\n",
    "In [Lesson 1](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/05_relax/05_01_2D.Laplace.Equation.ipynb) and [Lesson 2](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/05_relax/05_02_2D.Poisson.Equation.ipynb) of this module we used the Jacobi method (a relaxation scheme) to iteratively find solutions to Laplace and Poisson equations.\n",
    "\n",
    "And it worked, so why are we still talking about it?  Because the Jacobi method is slow, very slow to converge.  It might not have seemed that way in the first two notebooks because we were using small grids, but we did need more than 3,000 iterations to reach the exit criterion while solving the Poisson equation on a $41\\times 41$ grid. \n",
    "\n",
    "You can confirm this below: using `nx,ny=` $128$ on the Laplace problem of Lesson 1, the Jacobi method requires nearly *20,000* iterations before we reach $10^{-8}$ for the L2-norm of the difference between two iterates. That's a *lot* of iterations!\n",
    "\n",
    "Now, consider this application:  an incompressible Navier-Stokes solver has to ensure that the velocity field is divergence-free at every timestep.  One of the most common ways to ensure this is to solve a Poisson equation for the pressure field. In fact, the pressure Poisson equation is responsible for the majority of the computational expense of an incompressible Navier-Stokes solver. Imagine having to do 20,000 Jacobi iterations for *every* time step in a fluid-flow problem with many thousands or perhaps millions of grid points! \n",
    "\n",
    "The Jacobi method is the slowest of all relaxation schemes, so let's learn how to improve on it. In this lesson, we'll study the Gauss-Seidel method—twice as fast as Jacobi, in theory—and the successive over-relaxation (SOR) method. We also have some neat Python tricks lined up for you to get to the solution even faster. Let's go!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the same example problem as in [Lesson 1](./05_01_2D.Laplace.Equation.ipynb): Laplace's equation with boundary conditions\n",
    "\n",
    "\\begin{equation}\n",
    "  \\begin{gathered}\n",
    "p=0 \\text{ at } x=0\\\\\n",
    "\\frac{\\partial p}{\\partial x} = 0 \\text{ at } x = L\\\\\n",
    "p = 0 \\text{ at }y = 0 \\\\\n",
    "p = \\sin \\left(  \\frac{\\frac{3}{2}\\pi x}{L} \\right) \\text{ at } y = H\n",
    "  \\end{gathered}\n",
    "\\end{equation}\n",
    "\n",
    "We import our favorite Python libraries, and also some custom functions that we wrote in [Lesson 1](./05_01_2D.Laplace.Equation.ipynb), which we have saved in a 'helper' Python file for re-use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "from matplotlib import pyplot, cm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "%matplotlib inline\n",
    "from matplotlib import rcParams\n",
    "rcParams['font.family'] = 'serif'\n",
    "rcParams['font.size'] = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from laplace_helper import p_analytical, plot_3D, L2_rel_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have the analytical solution in the array `p_analytical`, and we have the functions `plot_3D` and `L2_rel_error` in our namespace.  If you can't remember how they work, just use `help()` and take advantage of the docstrings.  It's a good habit to always write docstrings in your functions, and now you see why!\n",
    "\n",
    "In this notebook, we are going to use larger grids than before, to better illustrate the speed increases we achieve with different iterative methods. Let's create a $128\\times128$ grid and initialize. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nx = 128\n",
    "ny = 128\n",
    "\n",
    "L = 5\n",
    "H = 5\n",
    "\n",
    "x = numpy.linspace(0,L,nx)\n",
    "y = numpy.linspace(0,H,ny)\n",
    "\n",
    "dx = L/(nx-1)\n",
    "dy = H/(ny-1)\n",
    "\n",
    "p0 = numpy.zeros((ny, nx))\n",
    "\n",
    "p0[-1,:] = numpy.sin(1.5*numpy.pi*x/x[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We said above that the Jacobi method takes nearly 20,000 iterations before it satisfies our exit criterion of $10^{-8}$ (L2-norm difference between two consecutive iterations). You'll just have to confirm that now. Have a seat!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def laplace2d(p, l2_target):\n",
    "    '''Solves the Laplace equation using the Jacobi method \n",
    "    with a 5-point stencil\n",
    "    \n",
    "    Parameters:\n",
    "    ----------\n",
    "    p: 2D array of float\n",
    "        Initial potential distribution\n",
    "    l2_target: float\n",
    "        Stopping criterion\n",
    "        \n",
    "    Returns:\n",
    "    -------\n",
    "    p: 2D array of float\n",
    "        Potential distribution after relaxation\n",
    "    '''\n",
    "    \n",
    "    l2norm = 1\n",
    "    pn = numpy.empty_like(p)\n",
    "    iterations = 0\n",
    "\n",
    "    while l2norm > l2_target:\n",
    "        pn = p.copy()\n",
    "        p[1:-1,1:-1] = .25 * (pn[1:-1,2:] + pn[1:-1,:-2] +\\\n",
    "                              pn[2:,1:-1] + pn[:-2,1:-1])\n",
    "        \n",
    "        ##Neumann B.C. along x = L\n",
    "        p[1:-1,-1] = .25 * (2*pn[1:-1,-2] + pn[2:,-1] + pn[:-2, -1])\n",
    "        \n",
    "        l2norm = numpy.sqrt(numpy.sum((p - pn)**2)/numpy.sum(pn**2))\n",
    "        iterations += 1\n",
    "     \n",
    "    return p, iterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Jacobi method took 19993 iterations at tolerance 1e-08\n"
     ]
    }
   ],
   "source": [
    "l2_target = 1e-8\n",
    "p, iterations = laplace2d(p0.copy(), l2_target)\n",
    "\n",
    "print (\"Jacobi method took {} iterations at tolerance {}\".\\\n",
    "        format(iterations, l2_target))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Would we lie to you?  19,993 iterations before we reach the exit criterion of $10^{-8}$.  Yikes!\n",
    "\n",
    "We can also time how long the Jacobi method takes using the `%%timeit` cell-magic.  Go make some tea, because this can take a while—the `%%timeit` magic runs the function a few times and then averages their runtimes to give a more accurate result.  \n",
    "\n",
    "- - -\n",
    "##### Notes\n",
    "\n",
    "1. When using `%%timeit`, the return values of a function (`p` and `iterations` in this case) *won't* be saved.\n",
    "\n",
    "2. We document our timings below, but your timings can vary quite a lot, depending on your hardware. In fact, you may not even see the same trends (some recent hardware can play some fancy tricks with optimizations that you have no control over).\n",
    "- - -\n",
    "\n",
    "With those caveats, let's give it a shot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 5.92 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "laplace2d(p0.copy(), l2_target)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The printed result above (and others to come later) is from a mid-2007 Mac Pro, powered by two 3-GHz quad-core Intel Xeon X5364 (Clovertown). We tried also on more modern machines, and get conflicting results—like the Gauss-Seidel method being slightly slower than Jacobi, even though it required fewer iterations. Don't get too hung up on this: the hardware optimizations applied by more modern CPUs are varied and make a big difference sometimes.\n",
    "\n",
    "Meanwhile, let's check the overall accuracy of the numerical calculation by comparing it to the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pan = p_analytical(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.1735513352884566e-05"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L2_rel_error(p,pan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a pretty small error.  Let's assume it is good enough and focus on speeding up the process.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss-Seidel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will recall from [Lesson 1](./2D_Laplace_Equation.ipynb) that a single Jacobi iteration is written as:\n",
    "\n",
    "\\begin{equation}\n",
    "p^{k+1}_{i,j} = \\frac{1}{4} \\left(p^{k}_{i,j-1} + p^k_{i,j+1} + p^{k}_{i-1,j} + p^k_{i+1,j} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "The Gauss-Seidel method is a simple tweak to this idea: use updated values of the solution as soon as they are available, instead of waiting for the values in the whole grid to be updated. \n",
    "\n",
    "If you imagine that we progress through the grid points in the order shown by the arrow in Figure 1, then you can see that the updated values $p^{k+1}_{i-1,j}$ and $p^{k+1}_{i,j-1}$ can be used to calculate $p^{k+1}_{i,j}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/solvepath.svg\" width=350>\n",
    "\n",
    "#### Figure 1. Assumed order of updates on a grid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The iteration formula for Gauss-Seidel is thus:\n",
    "\n",
    "\\begin{equation}\n",
    "p^{k+1}_{i,j} = \\frac{1}{4} \\left(p^{k+1}_{i,j-1} + p^k_{i,j+1} + p^{k+1}_{i-1,j} + p^k_{i+1,j} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "There's now a problem for the Python implementation.  You can no longer use NumPy's array operations to evaluate the solution updates.  Since Gauss-Seidel requires using values immediately after they're updated, we have to abandon our beloved array operations and return to nested `for` loops. Ugh.  \n",
    "\n",
    "We don't like it, but if it saves us a bunch of time, then we can manage. But does it? \n",
    "\n",
    "Here's a function to compute the Gauss-Seidel updates using a double loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def laplace2d_gauss_seidel(p, nx, ny, l2_target):\n",
    "    \n",
    "    iterations = 0\n",
    "    iter_diff = l2_target+1 #init iter_diff to be larger than l2_target\n",
    "    \n",
    "    while iter_diff > l2_target:\n",
    "        pn = p.copy()\n",
    "        iter_diff = 0.0\n",
    "        for j in range(1,ny-1):\n",
    "            for i in range(1,nx-1):\n",
    "                p[j,i] = .25 * (p[j,i-1] + p[j,i+1] + p[j-1,i] + p[j+1,i])\n",
    "                iter_diff += (p[j,i] - pn[j,i])**2\n",
    "        \n",
    "        #Neumann 2nd-order BC\n",
    "        for j in range(1,ny-1):\n",
    "            p[j,-1] = .25 * (2*p[j,-2] + p[j+1,-1] + p[j-1, -1])\n",
    "            \n",
    "        iter_diff = numpy.sqrt(iter_diff/numpy.sum(pn**2))\n",
    "        iterations += 1        \n",
    "        \n",
    "    return p, iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would then run this with the following function call:\n",
    "\n",
    "```Python\n",
    "p, iterations = laplace2d_gauss_seidel(p,1e-8)\n",
    "```\n",
    "\n",
    "<br>\n",
    "But **don't do it**.  We did it so that you don't have to!  \n",
    "\n",
    "The solution of our test problem with the Gauss-Seidel method required several thousand fewer iterations than the Jacobi method, but it took nearly *10 minutes* to run on our machine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### What happened?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you think back to the far off days when you first learned about array operations, you might recall that we discovered that NumPy array operations could drastically improve code performance compared with nested `for` loops.  NumPy operations are written in C and pre-compiled, so they are *much* faster than vanilla Python.  \n",
    "\n",
    "But the Jacobi method is not algorithmically optimal, giving slow convergence.  We want to take advantage of the faster-converging iterative methods, yet unpacking the array operations into nested loops destroys performance. *What can we do?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use Numba!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Numba](http://numba.pydata.org) is an open-source optimizing compiler for Python. It works by reading Python functions that you give it, and generating a compiled version for you—also called Just-In-Time (JIT) compilation.  You can then use the function at performance levels that are close to what you can get with compiled languages (like C, C++ and fortran).\n",
    "\n",
    "It can massively speed up performance, especially when dealing with loops.  Plus, it's pretty easy to use. Like we overheard at a conference: [*Numba is a Big Deal.*](http://twitter.com/lorenaabarba/status/625383941453656065)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Caveat\n",
    "\n",
    "We encourage everyone following the course to use the [Anaconda Python](https://www.continuum.io/downloads) distribution because it's well put-together and simple to use.  If you *haven't* been using Anaconda, that's fine, but let us **strongly** suggest that you take the plunge now.  Numba is great and easy to use, but it is **not** easy to install without help.  Those of you using Anaconda can install it by running <br><br>\n",
    "\n",
    "`conda install numba`<br><br>\n",
    "\n",
    "If you *really* don't want to use Anaconda, you will have to [compile all of Numba's dependencies](https://pypi.python.org/pypi/numba).\n",
    "\n",
    "\n",
    "- - -"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intro to Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's dive in!  Numba is great and easy to use.  We're going to first walk you through a simple example to give you a taste of Numba's abilities.  \n",
    "\n",
    "After installing Numba (see above), we can use it by adding a line to `import numba` and another to `import autojit` (more on this in a bit)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numba\n",
    "from numba import jit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You tell Numba which functions you want to accelerate by using a [Python decorator](http://www.learnpython.org/en/Decorators), a special type of command that tells the Python interpreter to modify a callable object (like a function). For example, let's write a quick function to calculate the $n^{\\text{th}}$ number in the Fibonacci sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fib_it(n):\n",
    "    a = 1\n",
    "    b = 1\n",
    "    for i in range(n-2):\n",
    "        a, b = b, a+b\n",
    "        \n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several faster ways to program the Fibonacci sequence, but that's not a concern right now (but if you're curious, [check them out](http://mathworld.wolfram.com/BinetsFibonacciNumberFormula.html)).  Let's use `%%timeit` and see how long this simple function takes to find the 500,000-th Fibonacci number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 4.22 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "fib_it(500000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try Numba!  Just add the `@jit` decorator above the function name and let's see what happens!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def fib_it(n):\n",
    "    a = 1\n",
    "    b = 1\n",
    "    for i in range(n-2):\n",
    "        a, b = b, a+b\n",
    "        \n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 138.50 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "1000 loops, best of 3: 499 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "fib_it(500000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Holy cow!* In our machine, that's more than 8,000 times faster!\n",
    "\n",
    "That warning from `%%timeit` is due to the compilation overhead for Numba.  The very first time that it executes the function, it has to compile it, then it caches that code for reuse without extra compiling.  That's the 'Just-In-Time' bit.  You'll see it disappear if we run `%%timeit` again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 499 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "fib_it(500000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would agree if you think that this is a rather artificial example, but the speed-up is very impressive indeed. Just adding the one-word decorator!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Running in `nopython` mode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba is very clever, but it can't optimize everything.  When it can't, rather than failing to run, it will fall back to the regular Python, resulting in poor performance again.  This can be confusing and frustrating, since you might not know ahead of time which bits of code will speed up and which bits won't.  \n",
    "\n",
    "To avoid this particular annoyance, you can tell Numba to use `nopython` mode. In this case, your code will simply fail if the \"jitted\" function can't be optimized. It's simply an option to give you \"fast or nothing.\"\n",
    "\n",
    "Use `nopython` mode by adding the following line above the function that you want to JIT-compile:\n",
    "\n",
    "```Python\n",
    "@jit(nopython=True)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- - -\n",
    "\n",
    "##### Numba version check\n",
    "\n",
    "In these examples, we are using the latest (as of publication) version of Numba: 0.22.1.  Make sure to upgrade or some of the code examples below may not run.  \n",
    "\n",
    "\n",
    "- - -"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.22.1\n"
     ]
    }
   ],
   "source": [
    "print(numba.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Back to Jacobi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We want to compare the performance of different iterative methods under the same conditions. Because the Gauss-Seidel method forces us to unpack the array operations into nested loops (which are very slow in Python), we use Numba to get the code to perform well. Thus, we need to write a new Jacobi method using for-loops and Numba (instead of NumPy), so we can make meaningful comparisons. \n",
    "\n",
    "Let's write a \"jitted\" Jacobi with loops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def laplace2d_jacobi(p, pn, l2_target):\n",
    "    '''Solves the Laplace equation using the Jacobi method \n",
    "    with a 5-point stencil\n",
    "    \n",
    "    Parameters:\n",
    "    ----------\n",
    "    p: 2D array of float\n",
    "        Initial potential distribution\n",
    "    pn: 2D array of float\n",
    "        Allocated array for previous potential distribution\n",
    "    l2_target: float\n",
    "        Stopping criterion\n",
    "        \n",
    "    Returns:\n",
    "    -------\n",
    "    p: 2D array of float\n",
    "        Potential distribution after relaxation\n",
    "    '''\n",
    "        \n",
    "    iterations = 0\n",
    "    iter_diff = l2_target+1 #init iter_diff to be larger than l2_target\n",
    "    denominator = 0.0\n",
    "    ny, nx = p.shape\n",
    "    l2_diff = numpy.zeros(20000)\n",
    "    \n",
    "    while iter_diff > l2_target:\n",
    "        for j in range(ny):\n",
    "            for i in range(nx):\n",
    "                pn[j,i] = p[j,i]\n",
    "                \n",
    "        iter_diff = 0.0\n",
    "        denominator = 0.0\n",
    "        \n",
    "        for j in range(1,ny-1):\n",
    "            for i in range(1,nx-1):\n",
    "                p[j,i] = .25 * (pn[j,i-1] + pn[j,i+1] + pn[j-1,i] + pn[j+1,i])\n",
    "                \n",
    "        \n",
    "        #Neumann 2nd-order BC\n",
    "        for j in range(1,ny-1):\n",
    "            p[j,-1] = .25 * (2*pn[j,-2] + pn[j+1,-1] + pn[j-1, -1])\n",
    "            \n",
    "            \n",
    "        for j in range(ny):\n",
    "            for i in range(nx):\n",
    "                iter_diff += (p[j,i] - pn[j,i])**2\n",
    "                denominator += (pn[j,i]*pn[j,i])\n",
    "        \n",
    "        \n",
    "        iter_diff /= denominator\n",
    "        iter_diff = iter_diff**0.5\n",
    "        l2_diff[iterations] = iter_diff\n",
    "        iterations += 1    \n",
    "        \n",
    "    return p, iterations, l2_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numba Jacobi method took 19993 iterations at tolerance 1e-08\n"
     ]
    }
   ],
   "source": [
    "p, iterations, l2_diffJ = laplace2d_jacobi(p0.copy(), p0.copy(), 1e-8)\n",
    "\n",
    "print(\"Numba Jacobi method took {} iterations at tolerance {}\".format(iterations, l2_target))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 2.41 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "laplace2d_jacobi(p0.copy(), p0.copy(), 1e-8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our old machine, that's faster than the NumPy version of Jacobi, but on some newer machines it might not be. Don't obsess over this: there is much hardware black magic that we cannot control.\n",
    "\n",
    "Remember that NumPy is a highly optimized library.  The fact that we can get competitive execution times with this JIT-compiled code is kind of amazing.  Plus(!) now we get to try out those techniques that aren't possible with NumPy array operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Note\n",
    "\n",
    "We're also saving the history of the L2-norm of the difference between consecutive iterations.  We'll take a look at that once we have a few more methods to compare.\n",
    "- - -\n",
    "\n",
    "\n",
    "##### Another Note\n",
    "\n",
    "Why did we use \n",
    "\n",
    "```Python\n",
    "l2_diff = numpy.zeros(20000)\n",
    "```\n",
    "\n",
    "Where did the `20000` come from? \n",
    "\n",
    "We cheated a little bit.  Numba doesn't handle _mutable_ objects well in `nopython` mode, which means we can't use a *list* and append each iteration's value of the L2-norm.  So we need to define an array big enough to hold all of them and we know from the first run that Jacobi converges in fewer than 20,000 iterations.\n",
    "\n",
    "\n",
    "- - -\n",
    "\n",
    "##### Challenge task\n",
    "\n",
    "It is possible to get a good estimate of the number of iterations needed by the Jacobi method to reduce the initial error by a factor $10^{-m}$, for given $m$. The formula depends on the largest eigenvalue of the coefficient matrix, which is known for the discrete Poisson problem on a square domain. See Parviz Moin, *\"Fundamentals of Engineering Numerical Analysis\"* (2nd ed., pp.141–143).\n",
    "\n",
    "* Find the estimated number of iterations to reduce the initial error by $10^{-8}$ when using the grids listed below, in the section on grid convergence, with $11$, $21$, $41$ and $81$ grid points on each coordinate axis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Back to Gauss-Seidel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you recall, the reason we got into this Numba sidetrack was to try out Gauss-Seidel and compare the performance with Jacobi.  Recall from above that the formula for Gauss-Seidel is as follows:\n",
    "\n",
    "\\begin{equation}\n",
    "p^{k+1}_{i,j} = \\frac{1}{4} \\left(p^{k+1}_{i,j-1} + p^k_{i,j+1} + p^{k+1}_{i-1,j} + p^k_{i+1,j} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "We only need to slightly tweak the Jacobi function to get one for Gauss-Seidel.  Instead of updating `p` in terms of `pn`, we just update `p` using `p`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def laplace2d_gauss_seidel(p, pn, l2_target):\n",
    "    '''Solves the Laplace equation using Gauss-Seidel method \n",
    "    with a 5-point stencil\n",
    "    \n",
    "    Parameters:\n",
    "    ----------\n",
    "    p: 2D array of float\n",
    "        Initial potential distribution\n",
    "    pn: 2D array of float\n",
    "        Allocated array for previous potential distribution\n",
    "    l2_target: float\n",
    "        Stopping criterion\n",
    "        \n",
    "    Returns:\n",
    "    -------\n",
    "    p: 2D array of float\n",
    "        Potential distribution after relaxation\n",
    "    '''\n",
    "    \n",
    "    iterations = 0\n",
    "    iter_diff = l2_target + 1 #initialize iter_diff to be larger than l2_target\n",
    "    denominator = 0.0\n",
    "    ny, nx = p.shape\n",
    "    l2_diff = numpy.zeros(20000)\n",
    "    \n",
    "    while iter_diff > l2_target:\n",
    "        for j in range(ny):\n",
    "            for i in range(nx):\n",
    "                pn[j,i] = p[j,i]\n",
    "                \n",
    "        iter_diff = 0.0\n",
    "        denominator = 0.0\n",
    "        \n",
    "        for j in range(1,ny-1):\n",
    "            for i in range(1,nx-1):\n",
    "                p[j,i] = .25 * (p[j,i-1] + p[j,i+1] + p[j-1,i] + p[j+1,i])\n",
    "                \n",
    "        \n",
    "        #Neumann 2nd-order BC\n",
    "        for j in range(1,ny-1):\n",
    "            p[j,-1] = .25 * (2*p[j,-2] + p[j+1,-1] + p[j-1, -1])\n",
    "            \n",
    "            \n",
    "        for j in range(ny):\n",
    "            for i in range(nx):\n",
    "                iter_diff += (p[j,i] - pn[j,i])**2\n",
    "                denominator += (pn[j,i]*pn[j,i])\n",
    "        \n",
    "        \n",
    "        iter_diff /= denominator\n",
    "        iter_diff = iter_diff**0.5\n",
    "        l2_diff[iterations] = iter_diff\n",
    "        iterations += 1    \n",
    "        \n",
    "    return p, iterations, l2_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numba Gauss-Seidel method took 13939 iterations at tolerance 1e-08\n"
     ]
    }
   ],
   "source": [
    "p, iterations, l2_diffGS = laplace2d_gauss_seidel(p0.copy(), p0.copy(), 1e-8)\n",
    "\n",
    "print(\"Numba Gauss-Seidel method took {} iterations at tolerance {}\".format(iterations, l2_target))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool!  Using the most recently updated values of the solution in the Gauss-Seidel method saved 6,000 iterations!  Now we can see how much faster than Jacobi this is, because both methods are implemented the same way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 2.09 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "laplace2d_gauss_seidel(p0.copy(), p0.copy(), 1e-8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get some speed-up over the Numba version of Jacobi, but not a lot. And you may see quite different results—on some of the machines we tried, we could still not beat the NumPy version of Jacobi. This can be confusing, and hard to explain without getting into the nitty grity of hardware optimizations.\n",
    "\n",
    "Don't lose hope!  We have another trick up our sleeve!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Successive Over-Relaxation (SOR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Successive over-relaxation is able to improve on the Gauss-Seidel method by using in the update a linear combination of the previous and the current solution, as follows:\n",
    "\n",
    "\\begin{equation}\n",
    "p^{k+1}_{i,j} = (1 - \\omega)p^k_{i,j} + \\frac{\\omega}{4} \\left(p^{k+1}_{i,j-1} + p^k_{i,j+1} + p^{k+1}_{i-1,j} + p^k_{i+1,j} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "The relaxation parameter $\\omega$ will determine how much faster SOR will be than Gauss-Seidel. SOR iterations are only stable for $0 < \\omega < 2$.  Note that for $\\omega = 1$, SOR reduces to the Gauss-Seidel method. \n",
    "\n",
    "If $\\omega < 1$, that is technically an \"under-relaxation\" and it will be slower than Gauss-Seidel.  \n",
    "\n",
    "If $\\omega > 1$, that's the over-relaxation and it should converge faster than Gauss-Seidel.  \n",
    "\n",
    "Let's write a function for SOR iterations of the Laplace equation, using Numba to get high performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def laplace2d_SOR(p, pn, l2_target, omega):\n",
    "    '''Solves the Laplace equation using SOR with a 5-point stencil\n",
    "    \n",
    "    Parameters:\n",
    "    ----------\n",
    "    p: 2D array of float\n",
    "        Initial potential distribution\n",
    "    pn: 2D array of float\n",
    "        Allocated array for previous potential distribution\n",
    "    l2_target: float\n",
    "        Stopping criterion\n",
    "    omega: float\n",
    "        Relaxation parameter\n",
    "        \n",
    "    Returns:\n",
    "    -------\n",
    "    p: 2D array of float\n",
    "        Potential distribution after relaxation\n",
    "    '''    \n",
    "    \n",
    "    iterations = 0\n",
    "    iter_diff = l2_target + 1 #initialize iter_diff to be larger than l2_target\n",
    "    denominator = 0.0\n",
    "    ny, nx = p.shape\n",
    "    l2_diff = numpy.zeros(20000)\n",
    "    \n",
    "    while iter_diff > l2_target:\n",
    "        for j in range(ny):\n",
    "            for i in range(nx):\n",
    "                pn[j,i] = p[j,i]\n",
    "                \n",
    "        iter_diff = 0.0\n",
    "        denominator = 0.0\n",
    "        \n",
    "        for j in range(1,ny-1):\n",
    "            for i in range(1,nx-1):\n",
    "                p[j,i] = (1-omega)*p[j,i] + omega*.25 * (p[j,i-1] + p[j,i+1] + p[j-1,i] + p[j+1,i])\n",
    "        \n",
    "        #Neumann 2nd-order BC\n",
    "        for j in range(1,ny-1):\n",
    "            p[j,-1] = .25 * (2*p[j,-2] + p[j+1,-1] + p[j-1, -1])\n",
    "            \n",
    "        for j in range(ny):\n",
    "            for i in range(nx):\n",
    "                iter_diff += (p[j,i] - pn[j,i])**2\n",
    "                denominator += (pn[j,i]*pn[j,i])\n",
    "        \n",
    "        \n",
    "        iter_diff /= denominator\n",
    "        iter_diff = iter_diff**0.5\n",
    "        l2_diff[iterations] = iter_diff\n",
    "        iterations += 1\n",
    "        \n",
    "    \n",
    "    return p, iterations, l2_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That wasn't too bad at all.  Let's try this out first with $\\omega = 1$ and check that it matches the Gauss-Seidel results from above.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numba SOR method took 13939 iterations at tolerance 1e-08 with omega = 1\n"
     ]
    }
   ],
   "source": [
    "l2_target = 1e-8\n",
    "omega = 1\n",
    "p, iterations, l2_diffSOR = laplace2d_SOR(p0.copy(), p0.copy(), l2_target, omega)\n",
    "\n",
    "print(\"Numba SOR method took {} iterations\\\n",
    " at tolerance {} with omega = {}\".format(iterations, l2_target, omega))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the exact same number of iterations as Gauss-Seidel.  That's a good sign that things are working as expected.  \n",
    "\n",
    "Now let's try to over-relax the solution and see what happens.  To start, let's try $\\omega = 1.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numba SOR method took 7108 iterations at tolerance 1e-08 with omega = 1.5\n"
     ]
    }
   ],
   "source": [
    "l2_target = 1e-8\n",
    "omega = 1.5\n",
    "p, iterations, l2_diffSOR = laplace2d_SOR(p0.copy(), p0.copy(), l2_target, omega)\n",
    "\n",
    "print(\"Numba SOR method took {} iterations\\\n",
    " at tolerance {} with omega = {}\".format(iterations, l2_target, omega))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow!  That really did the trick!  We dropped from 13939 iterations down to 7108.  Now we're really cooking!  Let's try `%%timeit` on SOR.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 1.18 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "laplace2d_SOR(p0.copy(), p0.copy(), l2_target, omega)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Things continue to speed up. But we can do even better!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tuned SOR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, we picked $\\omega=1.5$ arbitrarily, but we would like to over-relax the solution as much as possible without introducing instability, as that will result in the fewest number of iterations.  \n",
    "\n",
    "For square domains, it turns out that the ideal factor $\\omega$ can be computed as a function of the number of nodes in one direction, e.g., `nx`.  \n",
    "\n",
    "\\begin{equation}\n",
    "\\omega \\approx \\frac{2}{1+\\frac{\\pi}{nx}}\n",
    "\\end{equation}\n",
    "\n",
    "This is not some arbitrary formula, but its derivation lies outside the scope of this course.  (If you're curious and have some serious math chops, you can check out Reference 3 for more information).  For now, let's try it out and see how it works.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numba SOR method took 1110 iterations at tolerance 1e-08 with omega = 1.9521\n"
     ]
    }
   ],
   "source": [
    "l2_target = 1e-8\n",
    "omega = 2./(1 + numpy.pi/nx)\n",
    "p, iterations, l2_diffSORopt = laplace2d_SOR(p0.copy(), p0.copy(), l2_target, omega)\n",
    "\n",
    "print(\"Numba SOR method took {} iterations\\\n",
    " at tolerance {} with omega = {:.4f}\".format(iterations, l2_target, omega))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow!  That's *very* fast.  Also, $\\omega$ is very close to the upper limit of 2.  SOR tends to work fastest when $\\omega$ approaches 2, but don't be tempted to push it.  Set $\\omega = 2$ and the walls will come crumbling down.  \n",
    "\n",
    "Let's see what `%%timeit` has for us now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 184 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "laplace2d_SOR(p0.copy(), p0.copy(), l2_target, omega)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regardless of the hardware in which we tried this, the tuned SOR gave *big* speed-ups, compared to the Jacobi method (whether implemented with NumPy or Numba). Now you know why we told you at the end of [Lesson 1](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/05_relax/05_01_2D.Laplace.Equation.ipynb) that the Jacobi method is the *worst* iterative solver and almost never used.\n",
    "\n",
    "Just to convince ourselves that everything is OK, let's check the error after the 1,110 iterations of tuned SOR:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.7927433550683514e-05"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L2_rel_error(p,pan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking very good, indeed.  \n",
    "\n",
    "We didn't explain it in any detail, but notice the very interesting implication of Equation $(5)$: the ideal relaxation factor is a function of the grid size.  \n",
    "Also keep in mind that the formula only works for square domains with uniform grids. If your problem has an irregular geometry, you will need to find a good value of $\\omega$ by numerical experiments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decay of the difference between iterates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [Poisson Equation notebook](./05_02_2D.Poisson.Equation.ipynb), we noticed how the norm of the difference between consecutive iterations first dropped quite fast, then settled for a more moderate decay rate. With Gauss-Seidel, SOR and tuned SOR, we reduced the number of iterations required to reach the stopping criterion. Let's see how that reflects on the time history of the difference between consecutive solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAICCAYAAAAd2s0iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xdc1dX/wPHXuQgoKqigiKaAMxkyTFNbSNqwRLPxpRyp\nqWVZjkoblmhlS7O+lWlpmvnVTBtqmTYUtVwpOMn1yz0DceBAxvv3B3oT4cJFwMt4Px+P+9B7Pudz\nzvtzQXh7PuecjxERlFJKKaUcxeLoAJRSSilVvmkyopRSSimH0mREKaWUUg6lyYhSSimlHEqTEaWU\nUko5lCYjSimllHKocpOMGGPqGmO+NcZ8YoxZYIxp4eiYlFJKKQWmvOwzYoyZB8wWkZnGmObAd0Aj\nKS8fgFJKKVVCleiREWOMjzFmkTEms5DtVAfuBX4AEJFNgBNwS+GjVEoppVRhlNhkxBjTFVgJNABs\njl4YY2oaY2YYY7YZY/4yxswxxtS9opofkCoipy4rOwL4F3XcSimllCqYEpuMAM8D7YE/bFUwxjgD\nvwLOQDMgADgDLDXGuF2LIJVSSilVOCU5GblJRP4vnzq9gCBgmFwEDCdrNGXAZfX2Ai7GGPfLymoD\ne4ouXKWUUkpdjRKbjIiIPfNEugL7RGTvZecdBRKA+y8rOw78SNa8EYwxIUAGsLwoY1ZKKaVUwVVw\ndACF1BzYnkv5biDyirKngA+NMbcA9YBoXUmjlFJKOV5pT0a8gHW5lJ8C3IwxriKSCiAiB4D7rmVw\nSimllMpfaU9GioUxRkdMlFJKlSsiYhzVd2lPRhKBqrmUuwNnL42KXA29g1M6xcTEEBMT4+gw1FXS\nr1/ppl+/0ssYh+UhQAmewGqnTWTtIXIlf2DztQ1FKaWUUlejtCcj3wK+xpj6lwqMMd5k7Tky12FR\nKaWUUspupSEZyWvsaBpZIyBvG2OcjDEW4C3gb2DiNYhNlTARERGODkEVgn79Sjf9+qmrVWIflGeM\neQfoQNYy3OrAxouHWolI+mX1agLjgZZAJrAFGCwiBwvRt4wcOZKIiAj9x6WUUqrMio2NJTY2llGj\nRjl0AmuJTUYcyRijW5AopZQqN4wxDk1GSsNtGqWUUkqVYaV9aa9SShUpPz8/9u7dm39FpUoRX19f\n9uzZ4+gwbNLbNLnQ2zRKlV8Xh6sdHYZSRSq/72u9TaOUUkqpck2TERtiYmKIjY11dBhKKaVUsYmN\njS0Ru+bqbZpc6G0apcovvU2jyiK9TaOUUkoplQdNRpRSSinlUJqMKKWUKnZhYWF4enrSoEEDR4ei\nSiBNRpRSqgxbtWoVYWFh+Pj4YLFYCAoK4vXXX7/mccTHxxMVFWVX3SFDhtCyZctijkiVJJqMKKVU\nGdamTRvi4+N54oknMMbw008/MWLECEeHlafatWvj5+fn6DDUNaQ7sCqlVDlSGlYKDR8+3NEhqGtM\nR0aUUqocSkpK4umnnyYsLIzw8HBCQ0N5+eWXSU1NzVF36tSpNG/enMDAQEJCQujSpQs//vhjtjpT\npkwhODiYZs2a0bBhQ4YMGcLZs2dz7funn37illtuwc/Pj+DgYBYvXmw9NnDgQHx9fbFYLOzbt69o\nL1qVWJqM2KCbnimlyrJdu3YRGxvLsmXLiIuLY8WKFaxYsSLHqMS4ceN45pln+Oyzz9i6dSvr1q2j\nUqVKvPLKK9Y677zzDkOGDGHKlCn89ddfrF+/nuXLl9OpU6cc/SYmJjJ//nyWL1/Onj17eOCBB4iK\nimLbtm0AfPTRR4wePRpjHLblRblSUjY9Q0T0dcUr62NRSpVH9v77B4r1VdRiYmLEYrHI3r17RUTk\n7NmzcuDAgWx1Jk6cKFWqVLG+P3XqlFSpUkX69++frd727duldevWIiJy8uRJqVy5sjz++OPZ6syf\nP1+MMfLdd99Zy3r16iUWi0WOHj1qLUtNTZUaNWpIjx49rGXTpk3LFqsqvPy+py4ed9jvXZ0zopRS\n5VClSpVYunQpU6ZMISkpiQoVKpCUlMTZs2c5evQo3t7erFy5kjNnznDDDTdkO7dJkyasWrUKyFqt\nc/bs2Rx1WrVqBcAvv/xCly5drOXVq1enVq1a1vcuLi4EBQVZ21Plk96mUUqpq1Dc/1MsbhMnTuTR\nRx9l8ODBbNq0ibi4OEaPHg1gnTeSmJiIMYYaNWrYbOdSnerVq2crv3ROYmJitnJ3d/ccbVSvXp1D\nhw4V6npU6abJiFJKlUNffvklwcHBdO7c2WYdLy8vRITk5OR86xw/fjxb+aX3Xl5e2cpPnTqVo43j\nx49Tp06dgoSvyhhNRpRSqhxKTU3NMUn08OHD2d63bduWypUrs27dumzlW7dupX379kDWPiZubm78\n+eef2eqsXbsWYwx33HFHtvLk5GSOHDlifX/hwgW2bNlC27ZtC31NqvTSZEQppcqhe+65h82bN/Pr\nr78CcOTIESZPnpytTtWqVRk1ahSzZs1i7dq1AJw7d46XXnrJmjy4u7szcuRIvvrqK1avXg1kJRyj\nRo2iXbt2OUZeXFxcGDFiBJmZmQCMGTOGs2fP8tJLL1nrXIvbVKpk0QmsSilVDr300kucPn2a3r17\nU7t2berUqUNUVBQffPABHTt2ZPTo0XTt2pWhQ4dSo0YNHnvsMTIyMnBxceH+++/Ptovr888/j6en\nJ/369SMtLY0LFy7QpUuXbNvOh4WFsX//furUqUP79u1p3bo1R48excPDgwULFtC0aVMga5+RH374\nAYCOHTsydOhQ+vTpc20/HHXNGc1AczLGiH4uSpVPxpgy+T/zV199lTfeeIODBw9Su3ZtR4ejrrH8\nvq8vHnfY5i56m0YppcqghIQE/vjjD+v7pKQkXFxcqFmzpgOjUip3mowopVQZ9Oeff1pvk6SkpPDr\nr7/SqVMnnJycHByZUjnpnBGllCqDgoKCOHnyJAEBAUDWqpfx48c7OCqlcqdzRnKhc0aUKr/K6pwR\nVb7pnBGllFJKqTxoMqKUUkoph9JkxIaYmBhiY2MdHYZSSilVbGJjY4mJiXF0GDpnJDc6Z0Sp8kvn\njKiySOeMlFLn0887OgSllFKqXNBkxIaks0mODkEppZQqFzQZsSE9M93RISillFLlgiYjNmgyopRS\n6kpbtmwhLCwMV1fXAj/Ab9KkSQQGBmKxWJg+fXoxRVg6lZtkxBhTwRgzzBiTYoypn199TUaUUmVN\nQkICjz32GKGhobRo0YLmzZvTqlUrhg4dypIlSxwdXoGJCJMnT6Z169a0aNGC0NBQgoODiY6OZurU\nqQVub+3atXh6erJu3TqbdYKCgoiPj6dOnToFbv/xxx9n4cKFBT6vPCg3yQjwOLAMqGRPZU1GlFJl\nyf/+9z9uvvlmOnToQFxcHOvXr2fTpk1MnDiRRYsW0aFDB44fP+7oMAvkxRdfJCYmhunTp7N+/Xo2\nbNjAkiVL+Oeffxg5cmSB26tcuTJ+fn5Urly5GKJVeSlxyYgxxscYs8gYk1mU7YrIxyKyBrBr6VJl\nF/1mVEqVDXFxcfTp04f333+f6OhoLJZ/f/SHh4czZ84cB0Z39aZMmcJ//vMfmjRpYi2rWbMm48aN\nu6r2AgMDWb9+Pc2aNSuqEJWdSlQyYozpCqwEGgA2F0QbY2oaY2YYY7YZY/4yxswxxtQtylj8qvkV\nZXNKKeUwr7/+OlWrVqVbt265Hg8MDGTixInWEYFPP/2Um2++mVatWtG8eXM6derE9u3brfV///13\nwsLCsFgsjB49GoCTJ0/anEsxbdo0wsPDadGiBSEhIfTo0YONGzdaj69atYqIiAhatGhBWFgYd999\nN99//32+15WWlsa+fftylIeGhrJ69epsZefPn+e5556jQYMGNGvWjNDQUGbMmGE9/vPPP+e4pksW\nLlxIUFAQ/v7+3HrrrSxatCjXeJKTk+nXrx9+fn5cf/31tG7d2mZddQURKTEvYBXQEJgKZNio4wxs\nBGaTNcphgGnADsDtsnpPAzsvlkdeVp4J1M8nDlFKlU9l7d9/RkaGVKlSRTp06GD3Oc2aNZNffvnF\n+v69996T6667TlJSUrLVM8bIqFGjspX5+flJ7969re+XL18uFStWlD179oiIyNmzZyUiIsJ63unT\np6V69eoyc+ZM6znDhg2Tdu3a5Rtnly5dxBgj/fr1ky1btuRZt2PHjtK4cWM5cuSIiIj88ccfUrFi\nRfnyyy/zvKZNmzaJs7OztSwzM1P69esnVatWzXadqampEhYWJm3atLF+TnPmzJEKFSpIbGystd6e\nPXvEGCNffPFFvtdXlPL7vr543HG//x3ZeY5gwHLxz7ySkX5ABuB7WZk3kA48a0cfmowopWwqyL9/\nskZwc7yKqn5ROHbsmBhjpFu3bnafs3379mzvz507J8YY+frrr7OV25OMjB07VqpVqyapqanWst9/\n/10WL14sIiLr1q0TY4ysWrXKevzIkSPy8ccf5xvnoUOHJCIiQiwWixhjpHHjxjJkyBDZuHFjtnq/\n/PKLGGPk888/z1YeHR0tDRo0yPOaHn74YalevbpcuHDBWnbw4EExxmS7zs8++0wsFossWbIkW3ut\nW7fOllhpMpL7q0TdphERe+aJdAX2icjey847CiQA9xdXbEopVdbExsYSFhZG8+bN8fHx4b333gMg\nPT2dbt26ERISQnh4OG3atMEYw99//13gPm6++WZSUlJo2bIln332GYmJidx0003ccccdADRt2hRv\nb286d+7M66+/zs6dO/H29ubJJ5/Mt20fHx+WLl3KmjVrGDJkCMYYPvjgA0JDQxk0aJC13m+//YYx\nhrZt22Y7PzAwkD179uR6q+eS1atXExgYiLOzs7WsTp06VKtWLVu9S320adMmRx8rV64kIyMj3+sp\nz0pUMmKn5sDuXMp3A8G2TjLG3GKM+ZCs/4mMMsZEFVN8Sqlywtb/8oqqflHw9PSkSpUqHDlyJMex\niIgI4uPjWbBgAUePHiUlJYVDhw5x6623ArBmzRri4uKIj49HREhNTS1w/zfeeCPLli3Dz8+PgQMH\n4uPjwwMPPMChQ4cAqFKlCmvXrqVr166MGzeOpk2bcuONN7Jy5UoADh8+TFhYGOHh4YSFhXHvvffm\n6OOGG25g3LhxbN++nd9//52QkBA++ugjfvvtNwASExMRER566CHCw8OtrxkzZlC7dm2SkmzvuH3k\nyJEciQeAh4dHtveJiYnWhOfyPpYtW4anpyfJyckF/uzKkwqODuAqeAG5LQI/BbgZY1xFJMe/GBFZ\nAawgay5Jvo6kHKF2ldqFClQppRzNYrHQoUMHlixZQlpaWrb/4efmxx9/JDk5mWHDhlGxYsV8274y\nmTpz5kyOem3btmXevHn8888/TJkyhddee43o6GiWL18OQL169fjkk094//33+eabb3jppZfo2LEj\ne/bswcfHh/j4+Fz7nzlzJg8//DDG/LtIsk2bNkyYMIGbbrqJ+Ph4br/9dry8vDDGsHDhQurWLdha\nBx8fn1wTiRMnTmR77+XlhcViIS4urkDtqyylcWTkmhj64lBiYmKIiYkhNjbW0eEopdRVGzFiBOfO\nnWPy5Mm5Hr88ocht9CO3URWAWrVqZftFnZSUlGOUYdasWfzwww9A1rLbF154gb59+7Jp0yYga0fT\nN998EwBXV1ceeeQRxo8fz+nTp9mzZ0+e1/Xyyy+zc+fOHOWXkpOaNWsC0KFDBwA2bNiQrd7BgweJ\njo4mPd32vlJt2rQhISGBtLQ0a9mBAwc4efJktnodOnQgPT2drVu3ZiuPj4/niSeeyPM6HCE2Ntb6\nOy4mJsbR4ZSsCayXDVnmNYH1ILAkl/J5wOki6l++/+t7UUqVP5TRCewzZswQd3d3mTJliqSlpVnL\nd+zYIY8++qhYLBYZP368bNu2TVxcXKR3796SkZEhIiIDBw7MdbJqdHS0BAQEWFePjBgxQjw8PLJN\n7IyJiZFbbrlFTp06JSJZq2luvvlm6dq1q4iIxMbGipeXl+zcuVNEslarDB06VOrUqSPnz5/P85r8\n/PzkrrvukqNHj1rLDh06JO3atZN69erJyZMnreVRUVESHh5uXU2TkpIiDz74oAwaNChbm1deZ0JC\ngri6ukpMTIyIiKSnp8ujjz4qFStWzHadFy5ckJYtW8qdd94pp0+fFhGRpKQkufnmm+X999+31tu9\ne7dOYM3t964jO7cZVN7JyE/A37mUbwJWFlH/Mnfr3Dy/cEqpsqmsJiMiIlu3bpVevXpJYGCghIWF\nyfXXXy8BAQHy2GOPZVsFMm/ePAkODhZ/f3+JjIyUCRMmiMViER8fH3n00Uet9fbv3y8dOnSQevXq\nSWRkpCxevFj8/f3F09NTWrZsKSIiGzdulB49elj7DAwMlCeeeEKSk5NFRCQxMVGee+45ad68uYSF\nhUlQUJB06tQp36W6l+Ls3bu3NG/eXMLDw6Vp06bSsGFD6dOnj+zduzdb3QsXLsiLL74o/v7+1vqv\nvfaaZGZmiojI4sWLJTQ01HqdnTt3tp67ePFiCQ4OFl9fX7nxxhtl9uzZ1usMCwuz1jt58qQMGDBA\nfH19JTQ0VFq2bCkTJ060Hp84caIEBASIxWIRX19fGTBgQAG+eoVT0pMRkxVDyWKMmQr0FBGnXI71\nAyYC/iKy72KZN3AAGC4i7xVB//LV5q/4T9B/CtuUUqqUMcZQEn8uKlUY+X1fXzxu1w7lxaGkzhnJ\n6wOZBmwG3jbGOBljLMBbwN9kJSlFQp9No5RSSl0bJSoZMca8Y4yJB+69+D7u4su66kdE0oAOZG18\nlgBsBaqQtcvq2aKK5ZfPf9GJq0oppcq0SxNZHa1E3qZxNGOM6OeiVPmkt2lUWaS3aZRSSiml8qDJ\niFJKKaUcSpMRpZRSSjmUJiNKKaWUcihNRmx46vmndDWNUkqpMk1X05RgxhgZuXQkMRExjg5FKXWN\n6WoaVRbpappS6kLGBUeHoJRSSpULmozYoMmIUkopdW1oMmJDWkZa/pWUUkopVWgV8q9SPunIiFKq\nrIiLi2PMmDH8/fffAKSlpVGzZk1uueUWnnrqKWrVqmWtu3HjRt599102bNiAi4sLaWlpBAcHM2zY\nMEJDQ631/vnnH+644w727dtHcnIyoaGhiAipqamcP3+eli1b8vbbb+Pn53etL1eVQjoyYkOwd7Cj\nQ1BKqULbuHEjN910E5GRkcTFxREXF8emTZvo1KkTr732Glu2bLHWnTVrFu3ataNLly5s2bKFuLg4\nNm/eTFRUFBEREXz55ZfWujVr1iQ+Pp6oqCiMMcTFxREfH09CQgJr165l27Zt3H333aSmpjrislUp\no8mIDf1b9Hd0CEopVWgzZsygcuXKPPnkk9YyYwxDhgwhJCTEWrZx40Z69+7N+++/zwMPPJCtjejo\naMaNG0ffvn3ZsGFDvn16eXnRu3dvduzYwZo1a4ruYlSZpcmIDemZ6Y4OQSmlCi0tLY0zZ86QnJyc\n49ivv/7KzTffDMBrr71G5cqV6datW67t9OzZEzc3N9544w27+wU4fvz4VUauyhNNRmwYPWq0bnqm\nlCr12rVrR2pqKu3atWPevHnWJAHA09MTFxcXRISff/6ZFi1a4OTklGs7zs7OhIeH8/PPP+fb565d\nu/jss89wdXWlZcuWRXYtquiVlE3PdAKrDcNeHka1itUcHYZSqoQyo4p3fygZWTQbr3Xu3JmRI0fy\n5ptvct999+Hu7k779u15+OGHiYqKwtnZmcTERFJSUrJNZM2Nt7c3KSkpJCYm4uXl9W+sIoSHh5OZ\nmcnBgwdJSkqiSZMmfP3119StW7dIrkMVj4iICCIiIhg1apRD49CRERtOnznt6BCUUqpIjBw5kgMH\nDjB+/HhuuOEG5s+fz4MPPkiLFi04cOBAgdszxuR4HxcXx4YNG/jrr79o1aoVTz31FPfee29RXYIq\n43RkxIble5fTzTP3e6dKKVVUIxfXipeXF4MGDWLQoEEcO3aM119/nY8++ogXX3yRL774gsqVK3Ps\n2LE82zh69ChVq1bF09Mzz37GjBlD+/btCQgI4Pbbby/qS1FlkI6M2DB2zVhHh6CUUoW2fv16/vrr\nr2xltWrV4r///S9NmjQhPj4ei8XCnXfeyfr168nIyMi1nbS0NOLi4ujYsWO+fUZGRhIeHs7o0aOL\n5BpU2afJiA2pabo2XilV+v3www9Mnz4912PGGOvcj1deeYWzZ88yY8aMXOtOmzaN8+fP8/LLL9vV\n76BBg/j9999ZtWrV1QWuyhVNRmxIzdBkRClVNkycOJEVK1ZY32dkZDB27Fh27NjBwIEDAQgJCWHq\n1KkMHTqUOXPmWJ/wKiJ89dVXvPDCC0yfPp3AwMBsbdt6Emx0dDS1atXizTffLKarUmWJzhmxQbeD\nV0qVBY888gjp6ekMHz6cc+fOkZGRwalTp2jYsCHff/89nTp1staNjo4mICCAN998k1GjRuHi4sKF\nCxcICQlh2bJlBAUFWete2g5+//79AISHh3PXXXcxZswYIGsp8BNPPMHo0aMJDw9n6NChdO/e/dpe\nvCo1jK2stjwzxkidd+tw8LmDjg5FKXWNGWNs/m9fqdIqv+/ri8eLd716HvQ2jQ3h3uGODkEppZQq\nFzQZsWFA8ABHh6CUUkqVC5qM2KCraZRSSqlrQ5MRGzQZUUoppa4NTUZsuJCmq2mUUkqpa0GTERu2\n7djm6BCUUkqpckH3GbFh7ty53N7wdm6P1OcqKKWUKptiY2OJjY11dBi6z0hujDFCDJx56Qxuzm6O\nDkcpdQ3pPiOqLNJ9Rkox3YVVKaWUKn6ajORBkxGllFKq+GkykgdNRpRSSqniVy6SEWNMdWPMFGPM\nBxdfC4wxjfM7T5MRpZQqf7Zs2UJYWBiurq706dPH0eGUC+UiGQHqA+dEZJCIDAIWA5/ndYJ/JX8q\nVqh4TYJTSqnitnHjRrp3705QUBDh4eEEBwfzyCOPsGHDhkK1O2rUKJYvX56jfMiQIbRs2bJQbV9u\n7dq1eHp6sm7duiJr05agoCDi4+OpU6eOXfVFhMmTJ9O6dWtatGhBaGgowcHBREdHM3Xq1FzrT5s2\njdtuu43Q0FDCw8Np0aIFI0eO5NSpU9nqTpkyhbCwMCwWC56enoSHhxMWFkajRo24/vrreeutt8jI\nyCiS63YoESlRL8AHWARkFmMfHYH/y+O4dP+4uyilyp+sH4tly8yZM6V69eoyZ86cbOWzZs0SDw8P\nmT59+lW3bYyRUaNG5Sh/66235IEHHrjqdq+0ZcsWCQ8Pl4SEhCJrMz9+fn7Su3fvfOsNHz5c6tat\nK9u3b7eWHTt2TCIjI6VevXrZ6qanp0uXLl2kTZs2smvXLmv5iRMnpFevXtK4cWPZt29fjj6MMdKn\nT59sZQsXLhQnJycZOXJkvjHm93198bjjfvc7svMcwUBXYDewA8jIo15NYAawDfgLmAPULUA/nwBD\n8zguXd/vmucXTilVNpW1ZGTDhg3i6uoqX3zxRa7HJ0+eLC4uLhIfH39V7dtKRsoCe5MRLy8vGTp0\naI7y+Pj4HMnIsGHDxNPTU5KTk3Nt65ZbbpFWrVrlKDfG5BpLSEiI+Pv75xtjSU9GStptmueB9sAf\ntioYY5yBXwFnoBkQAJwBlhpj3C6r97QxZqcxZocxJvKy8nuAqiLyXl6BnEs7V6gLUUqpkuC1116j\ncuXKdOvWLdfjPXv2xM3NjTfeeAOASZMmERgYiMVi4b333qNbt26EhYXh5eVFv379OHcu62fj2rVr\nCQsLwxjDxIkTCQsLIzw8nG3btjFw4EB8fX2xWCzs27cvR7uffPIJ/fv3Jzg4mIYNG/LNN9+QkZHB\n0KFDCQkJoVGjRnz33XfWGH/++WfrrYrRo0dby/38/AgPD7e+wsLCcHJyombNmtY6IsKYMWNo0qQJ\nAQEBBAQE8N57OX/8L1y4kKCgIPz9/bn11ltZtGiR3Z9xWlqa9TovFxoayurVq63vk5KS+O9//0u3\nbt2oVq1arm09/fTT/Pnnn/z444929Z2ens7x48ftjrXEcmQmdOULsFz8cyo2RkaAfkAG4HtZmTeQ\nDjybT/sdgYlc3Owtj3oSMSYizyxSKVU2UYZGRjIzM6Vq1arSoUOHPOtFRkaKu7u79f2ePXvEGCM+\nPj6yYcMGERHZv3+/1KtXT3r06JHtXGOMjB49Okeb06ZNE4vFInv37s3RbmhoqBw4cEBERF566SVx\ncXGRkSNHyv79+0Uk67ZHlSpV5OTJkzn6unwU5soRgSlTpojFYsk2CjRgwADx8vKy3kLZtm2b1KxZ\nU1577TVrnU2bNomzs7O17czMTOnXr59UrVrVrpGRLl26iDFG+vXrJ1u2bLFZ7+uvvxaLxSL/+9//\nbNbZt2+fGGPkqaeeynHtV8Yybdo0McZIp06d8o0xv+9rdGTkXyKSaUe1rsA+Edl72XlHgQTgflsn\nGWMeBO4QkScufmHfz6uTHf+3w86olVLllRllcn0VVf3CSkxMJCUlhVq1auVZz9vbm5SUFBITE7OV\nd+nShZCQEACuu+46nnnmGWbNmsWuXbuy1cv6XWa/9u3bU7duXQDuv/9+0tLSSElJ4brrrgPgwQcf\n5OzZs/z55595tvP2229b//73338zZMgQunbtSs+ePQHYtWsXkyZN4sknn6RJkyYANG3alMcee4x3\n3nnHOsrz5ptvUqVKFV588UUgazfSmJgYUlJS7LqeCRMmcNtttzFlyhSCg4Np0qQJQ4cOZdOmTdnq\n7dmzByDPr4e3t3e2upebP38+4eHhBAQEUKlSJR5//HE6derEpEmT7IqzJCtRyYidmpM1r+RKu4Hg\n3E4wxgQDM4GHjDGHjTGHgb55dZJSMYX9J/cXNlallCo1jMmeGAUEBGR736JFCzIyMlizZk2h+mnY\nsKH17zVq1ACgUaNG1jJPT09EhCNHjuTZzoMPPghAZmYmPXv2pEqVKtl+Mf/2228AtG3bNtt5gYGB\npKSkWJOd1atXExgYiLOzs7VOnTp1bN5KuZKPjw9Lly5lzZo1DBkyBGMMH3zwAaGhoQwaNMiuNq50\n5dcCICplFNp2AAAgAElEQVQqiri4OBISEpg7dy7NmjXjtddew8fH56r6KElK44PyvIDc1nadAtyM\nMa4iknr5ARHZTNYcE7udqnmKP/b/QbRH9NVHqpQq02RkwUYEClq/sDw9PalcuTLHjh3Ls97Ro0ep\nWrUqnp6e2crd3d2zva9evToAhw4dKlRcbm7/PvPr0i/d3MrsXbL61ltvsWrVKhYuXGhNbiBrZEhE\nGDx4MJUqVbKWp6am4uPjw4kTJwA4cuQIgYGBOdr18PAowFXBDTfcwA033MC4ceNYtWoVTz75JB99\n9BFRUVHcfvvt+Pr6IiJ5fj2OHj0KgL+/f5593XPPPcyePZsuXbqwY8cOKlQojb/O/1UaR0aumfPp\n5x0dglJKXTWLxcKdd97J+vXrbf5iT0tLIy4ujo4dO+Y4duWeF5cmStq7/8a1EBcXx6hRoxgwYAB3\n3nlntmNeXl4YY5g8eTJxcXHW19atWzl48CBRUVFA1shGcnJyjrYvJSv5mTlzZo5bVW3atGHChAmI\nCPHx8QBERkbi6uqabVLrlVatWoUxhnvvvTfffocNG8aePXv48ssv7YqzJCuNyUgiUDWXcnfg7JWj\nIldtKcydMJeYmJgS8XhlpZS6Gq+88gpnz55lxowZuR6fNm0a58+f5+WXX85xbOvWrdner1u3Dicn\nJ1q1amUtq1ChgvUX8bZt29i4cWMRRp+38+fP0717dxo0aMDYsWOt5RMnTuTkyZO0b98eIMfGbqmp\nqTz44IMkJSUBWYlDQkICaWlp1joHDhzg5MmTdsXx8ssvs3Pnzhzll0Z4vLy8rH8OGjSImTNn5pr8\nAHz44Ye0bduWO+64I99+g4KCaNeuHe+++65dcV4uNjaWmJgY68vhHDl71taLvFfT/AT8nUv5JmBl\nEfUvxCDvr3o/lznHSqmyjDK0muaSWbNmSY0aNeTrr7+WzMxMEclaMXJ5+eUurXpp1KiRdTXN3r17\npX79+tKzZ89sdZs0aSJ9+/YVEZFu3bpZV6lMnTo1x2qa3bt3izEm22oXe8tEcq6mefrpp8XFxUXW\nrl2brV5ERIS136efflp8fX1lx44dIiKSlpYmgwYNks6dO1vrJyQkiKurq8TExIhI1sZkjz76qFSs\nWNGu1TR+fn5y1113ydGjR61lhw4dknbt2km9evWyrQpKT0+Xrl27Sps2bawxiYgkJydLr169JCAg\nQA4fPpyjD1v7jMybN0+MMTJ37tw8Y8zv+xrd9KzAycilpb31LyvzBtLIYyOzAvYvxCBv//52nl88\npVTZUxaTERGRjRs3SnR0tAQGBkpYWJgEBgbKI488Ips3b85R91Iy8vHHH8tjjz0mYWFh4uXlJf36\n9ZOzZ89mqztv3jxp1KiRhISESPv27SUxMVGeeuop8fX1FYvFIoGBgTJlyhSZPXu2BAQEiMViEV9f\nXxkzZowsXrw437Jhw4bJ4sWLJTQ0VCwWi/j4+Ejnzp3l8OHDYrFYpGLFiuLv7299+fn5SaVKlbIl\nQe+++640adLEeu3PPvtsjutYvHixBAcHi6+vr9x4440ye/Zs8ff3F09PTwkLC8vzs503b5707t1b\nmjdvLuHh4dK0aVNp2LCh9OnTJ1scl/v888/lpptukubNm0toaKiEhYXJqFGjJCUlJVu9yZMnW6/9\nUiyrVq2yHs/MzJRGjRqJl5eXhIWF5VgOfUlJT0ZMVgwlizFmGtBDRJxyOeYM/EnWzqvdAQGmAG2B\nMBE5WwT9i9MIJ+Y8Mof7mt1X2OaUUqWIMYaS+HPxWtq7dy/+/v5MmzbNukxWlW75fV9fPF5868zz\nUaKm3xpj3gE6APUuvo+7eKiViKQDiEiaMaYDMJ6svUUygS1AZFEkIlZLofrt1bP2eFVKKaXKoNjY\n2BIxL7JEjow4mjFGnF5yIv2NdEeHopS6xnRkJGvDrQYNGujISBlS0kdGSuNqmmsio0JGuf+BpJQq\nfyZNmsQ999yDMYZXX32Vp59+2tEhqXJAR0ZyYYwRYuDMi2dwc3HLt75SquzQkRFVFunISCl2MPGg\no0NQSimlyjxNRvLw41b7HuGslFJKqaunyUgepm2d5ugQlFJKqTKvRC3tLVGWwj9N/nF0FEoppVSx\n0aW9JdilCay1LbU5/MphR4ejlLqGdAKrKotK+gRWHRnJgz61V6nyx9fX1/qAM6XKCl9fX0eHkCdN\nRvJgnPUHklLlzZ49exwdglLljk5gzUNVqeroEJRSSqkyT5ORPHhe8HR0CEoppVSZp8lIHk6knnB0\nCEoppVSZp8lIHk5nnHZ0CEoppVSZpxNYbVkKp2qecnQUSimlVLHRfUZKsEv7jHAS5D39fJRSSpVt\njt5nRG/T5KUyJCUnOToKpZRSqkzTZMQWASpA7PpYR0eilFJKlWmajNhy8e7Mis0rHBuHUkopVcZp\nMmJLZtYfWw5scWwcSimlVBmnyYgNTplOABw+qw/KU0oppYqTJiM2uLm4AeBc3dnBkSillFJlmyYj\nNtSvVB+A5LPJDo5EKaWUKts0GbHhgcYPAHAqRTc+U0oppYqTJiM2hDcKB+CUnCIzM9PB0SillFJl\nlyYjNjT1aQpAZpVMDhw44OBolFJKqbJLkxEbfKv5Zv3FA1auWunYYJRSSqkyTJMRGypWqEiFsxXA\nCaZ9O83R4SillFJlliYjNuxM2omrxRWA7ce2OzgapZRSquyq4OgASqpnhj/DmfQz4A9HUo84Ohyl\nlFKqyMXGxhIbG+voMDAi4ugYShxjjExeP5m+C/pmFSyDC4sv4OysG6AppZQqe4wxiIhxVP96m8aG\nyi6V/31THZYtW+a4YJRSSqkyTJMRGyo7Z09G/vzzT8cFo5RSSpVhmozY4Obs9u8bT0hMTHRcMEop\npVQZpsmIDXXd6/JQwENUMBXADZavW+7okJRSSqkyqVxMYDXGjAeqAieAUGCqiPwvj/py6XNpMbEF\ncUfjcJruxKnNp3Bzc7N1mlJKKVUq6QTWa+OCiPQVkeeAl4DJxhi7rj2odhAAGdUzWL16dTGGqJRS\nSpVPJSoZMcb4GGMWGWOK9Ml0IjL8srfXA1tExK4+mnk1y/pLTVixYkVRhqWUUkopSlAyYozpCqwE\nGgA27x0ZY2oaY2YYY7YZY/4yxswxxtS1o/3mxpivgaeBB+2Ny5qMeMGcOXPsPU0ppZRSdioxyQjw\nPNAe+MNWBWOMM/Ar4Aw0AwKAM8BSY4zbZfWeNsbsNMbsMMZEAojIJhF5CBgIrDDGVLcnqGY1/x0Z\n2bp1K6dPn76aa1NKKaWUDSUpGblJRP4vnzq9gCBgmFwEDCdrNGXApUoi8qGINBaRJmQlKpUvO7YG\nOAfcll9AP+z4geV7l+NscQYPwBXmzp1b4AtTSimllG0lJhmxcw5HV2CfiOy97LyjQAJwv41z6gOT\nL70xxngC3kB+iQ/P/PQM/Rb0o7Fn46wCb5g+fbodYSqllFLKXqXtQXnNgdweobsbiLRxznEAY8zn\nQDJZE1gHicjm/Dpzd3UHoHGNxiT8kwA+sGbNGkQEYxy2AkoppZQqU0pbMuIFrMul/BTgZoxxFZHU\nyw+IyGng4avpzKOiBwC+Hr5ZBbXh3JpzrFmzhtatW19Nk0oppZS6QmlLRq6ZmJgYDm8+DElwrsY5\nACr6V+Q855k1a5YmI0oppUqt2NhYYmNjHR2GVYnbgdUYMxXoKSJOuRw7CGwXkcgryucBkSJStYhi\nEBGh27fdmLl5Jp91+oz+C/pjMRYyRmfgU8uH/fv34+SUI0SllFKq1NEdWAtmE+CXS7k/kO8ckIK6\ntf6t9AzpSRPPJjTxbEKGZFAntA6HDx/WDdCUUkqpIlLakpFvAV9jTP1LBcYYb7L2HCnyNbeP3/A4\nX3T5glt9byXMJwyAkLtCAF1Vo5RSShWVkpiM5DVMNI2sEZC3jTFOF58v8xbwNzCxOINqVacVAFWu\nrwLAV199xYkTJ4qzS6WUUqpcKDHJiDHmHWNMPHDvxfdxF1/WSbYikgZ0ADLI2ltkK1CFrPkiZ4sy\nnpiYmGyTe1pflzVhNeFUApGRkZw7d44ZM2YUZZdKKaXUNRUbG0tMTIyjwyh5E1hLgksTWC+Xmp6K\n+1vupGWkcevKW1n28zKaNGnCtm3bdM8RpZRSpZpOYC0lXCu4Eu4TjiBY6md9bDt27ChRS6OUUkqp\n0kiTkTycTz/P/zb9j8lxWbvJt7muDQD12tSz1nnrrbccEptSSilVVmgykoeMzAy6f9edZ356Bvh3\n3shRl6P4+PgA8PPPP7NhwwaHxaiUUkqVdpqM5MHN2Y2KFSpyLv0cZy6csY6MrDm4hp6P9rTWe/vt\ntx0VolJKKVXqaTJiQ0xMDMuWLcPLzQuAxLOJ1POoRz33epw4f4JbH7yVSpUqYYxh9uzZbN261cER\nK6WUUgVTUlbTaDJiQ0xMDBEREdZk5J+z/wAQ4RcBwM4LOzl27BgDBgxARBgxYoSjQlVKKaWuSkRE\nhCYjpUFNt5pA1sgIQDu/dgDE7o2lSpUqjBgxAjc3N77//ntWr17tsDiVUkqp0kqTkXzc0/geBtww\nAJ8qWRNWL42MLNuzjEzJxMfHh8GDBwMwfPhwdN8WpZRSqmB007Nc5Lbp2eX83vdj78m9xPWPI8wn\njJMnT9KwYUOSkpKYNWsW0dHR1zBapZRSqnB007NSqJ1/1q2apXuWAuDh4WFdUTN06FBOnz7tsNiU\nUkqp0kaTkasQ4RsBQOyeWAD27NnDnDlzcHd35/Dhw4waNcpxwSmllFKljN6myUV+t2n2ndyH7/u+\nuLu6kzQsiTOnz1C/fn1OnToFgJOTE2vWrKFFixbXKmSllFLqqultmlKovkd9mng24VTqKdYcWIOH\nhwcDBw7MOla/PhkZGTz66KOkpqY6OFKllFKq5NNkxIaYmBhiY2M5nXqaT9d/yqR1k7Idv6vhXQAs\n2rUIgOeeew4PDw/27dtH3bp12bp1a4lYu62UUkrZUlI2PdPbNLm4/DbN0ZSj1B5Xm5puNTn2/DFr\nnZ92/kTHmR25oc4N/NnvTwDGjBnDyy+/TLNmzdi2bRvGGH7//XfatGnjkOtQSiml7KG3aUq4GpVq\nAJB0LolMybSW3+Z3G65Orqw7tI5jZ7KSlEGDBtGkSRMefvhhhg4dSmZmJtHR0Rw/ftwhsSullFKl\ngSYj+XB2cqZaxWpkSibJ55Kt5W7ObtzmdxsAv/zfLwBUrlyZhIQEXnnlFcaMGUOrVq3Yt28fvXr1\n0s3QlFJKKRs0GbHDlVvCX2KdN/J/i6xlTk5OALi4uDB79myqVavGggULeO+9965RtEoppVTposmI\nHa58WN4ldzXKSkYW71qc7RbOJX5+fkybNg3I2ip+yZIlxRuoUkopVQppMmKHhwIf4rk2z+Fd2Ttb\n+fVe1+Pr4cs/Z/9h/aH1uZ7buXNnhg0bRkZGBg888AC7du26FiErpZRSpYYmI3YY3How797xLo09\nG2crN8ZwT+N7AJi/fX6O80SEOXPmcN9993HvvfeSnJxMp06dOHHixDWJWymllCoNNBkppM7XdwZg\n3vZ5OY59+umnPPTQQ/Tv358vvviCoKAgtm3bxn/+8x/S0tKudahKKaVUiaTJSCFF+EXg7urO5mOb\n2Z28O9uxHj160KBBAzZv3szEiRNZsGABNWvW5Oeff6ZPnz5kZuacZ6KUUkqVN5qMFJKLkwt3N7ob\nyDk64ubmxqeffgrA6NGjSU1N5ccff6Ry5crMmDGDZ599Vpf8KqWUKvc0GSkCnZvavlVz++2307t3\nb1JTU+nfvz8tWrTgu+++w9nZmffff5+33377WoerlFJKlSi6HXwurnxq78nzJ5m2YRqCMLj14Bz1\nT5w/Qc13ayIiHHv+mHXX1kuOHz9Os2bNqFmzJkuWLKFWrVp8/fXXREdHIyJ88MEHPPPMM8V+XUop\npVRudDv4UuB8+nkGLx7MmBVjcj1erWI1IvwiyJAMftzxY47jNWrUYPHixaxbt45atWoB8NBDD/HJ\nJ58AWdvI//e//y2+C1BKKaVKME1G7ODp5onBkHg2kYzMjFzrXLpVM39HziW+AKGhoVSsWDFb2eOP\nP86ECRMATUiUUkqVX5qM2BATE0NsbCwAFSwVqFGpBoKQdC4p1/pRTaOArKf5nks7Z3c/AwYMyJaQ\njB8/vnCBK6WUUnaKjY0lJibG0WHonJHcXDlnBCDg4wD+SvyLzQM2E1QrKNfzWn7WknWH1vHdf76j\ny/VdCtTnJ598wpNPPgnAiBEjGD16NMY47PadUkqpckTnjJQStSpnzfU4duaYzToPNHsAgLkJc/Nt\n7/Tp0wwePJjk5KwnAQ8YMICpU6fi5OTE66+/zoABA8jIyP2WkFJKKVWWaDJip+7NuzPytpHUc69n\ns84DAVnJyPzt8zmffj7P9vr3788HH3xA7969rXuN9OrVi2+//RZXV1cmTZrEww8/TGpqatFdhFJK\nKVUC6W2aXOR2m8Ze4ZPCiT8Sz/zo+XRq2slmvd27dxMeHs6JEycYN24cQ4cOtR5bvnw5nTp14tSp\nU0RERPDNN99Qo0YNm20ppZRShVFmbtMYY14vqrZKswcDHgRgTsKcPOv5+/szdepUAIYPH87q1aut\nx2699VaWLVtG7dq1iY2NpU2bNuzcubP4glZKKaUcqMDJiDHGyxhT/8oX0K8Y4itSxpjnjDHF+kCY\ny2/VpKbnfYulS5cuDBkyhPT0dB566CGSkv5dqRMaGsratWsJCQlhx44dtG7dmmXLlhVn6EoppZRD\n2JWMGGM8jTEzjDFngaPA7lxeXsUWZREwxgQCEUCx3pdq7NmYEO8QTqae5Ne/f823/ltvvUWrVq24\n6667qFKlSrZj9erVY8WKFdx7770cP36cDh06WEdTlFJKqbLCrjkjxphvgduAH4ADwIUrqwBDRcS9\nUMEY4wNMBe4QkaK8hVQB+Bp4FdgoIk751L/qOSMAry9/nVeWvsKjIY8yrcu0fOunpKTkSEQul5GR\nwfPPP2/dg+Spp57ivffew8XF5apjVEoppS5x9JwRe5ORZKCliOzKo856EWlx1YEY0xUYB6QBDW0l\nDMaYmsB44AayRjm2AINF5GAebb8O/ALsAf6+mmTkVOopPlzzIZmSySu3vZLntWxP3M71H19PtYrV\nOPrcUVyciiZpmDx5Mk899RQXLlygbdu2zJkzhzp16hRJ20oppcovRycj9o4+HMgrEQEoTCJy0fNA\ne+APWxWMMc7Ar4Az0AwIAM4AS40xbpfVe9oYs9MYs8MY8wpQSUSWkTWCc1UyJZMRS0cwdtXYfOs2\n9WpKUK0gTpw/wW9//3a1XebQt29fVqxYwXXXXcfKlStp0aIFK1asKLL2lVJKKUewNxl50xjTN68K\nxph1hYzlJhH5v3zq9AKCgGFyETAcaAAMuFRJRD4UkcYi0gSoDFQ3xkwAXr8Y6wRjzP0FCc7D1QNn\nizOnUk/lu4cI/Luq5uuErwvSjVVycjLHjuXcYK1Vq1asX7+edu3aceTIESIjIxk/fjyFua2klFJK\nOZLd+4wYY3qS9Qt/PZAIXLkq5bnCzhm52M9UoGdut1KMMT8B14uI/xXlm4AUEWmbT9u+wO785qPY\nmjNS9726HDp9iH2D91HPw/bmZwDbErfR7ONmeLh6cPS5o7hWcM2z/uV27NjBPffcg6enJ0uXLqVS\npUo56qSnp/Piiy8ydmzWSM29997L1KlT8fIq0fOIlVJKlUCl4jaNMaYjMAm4EXiSrImgMVe8KhdD\nfFdqTtbKnSvtBoLzOtEYcxswChBjzH+NMbk/YCYP9mwJf8n1XtcTWjuUk6knWbRrUYH68fDwIC0t\njTVr1tCjRw8yM3OuRq5QoQLvvvsu3377LdWqVeOHH34gJCTE+nA/pZRSqrSw9zbNO2StpGkHNAX8\nr3g1AHJ/nG3R8gJO51J+CnAzxtgcfhCRZSLSS0ScROQZEdlS0M4vJSNHzxy1q350YDQAs7bMKlA/\n3t7e/Pjjj3h4ePDNN9/w4osv2qx73333sWHDBtq2bcuhQ4eIjIxk5MiRpKenF6hPpZRSylHsXU1z\nHKgpIjaf3GaMGSYi7xQ6oLxv06QCi0Sk8xXlXwKPAG4iUuiHuRhjZOTIkdb3ERERREREMHPzTA6c\nOkDXZl1pVKNRvu3sPbEXvw/8qFShEseeP0YVF9vLd3Pz66+/cvfdd5Oens6kSZPo37+/zbrp6enE\nxMQwZswYRISbb76Z6dOn4+/vb/McpZRS5VNsbGy2kfRRo0aViqW9sUAnEcltVOJSnZtExOZKGLsD\nyjsZOQhsF5HIK8rnAZEiUrWw/V9sr1D7jFzups9vYuX+lczsOpOHgx8u8Pmff/45jz32GMOHD+et\nt97Kt/6SJUvo1q0bR44coUqVKrz//vv06dMHYxz2PaaUUqqEKxVzRoCBwEfGmGZ51Mn7YSxFYxPg\nl0u5P7D5GvRfYFd7q+aSPn36sGbNGrsSEYDIyEg2b97M/fffT0pKCn379iUqKoojR45cVf9KKaVU\ncbM3GZkP3A1sMcacNsbsMcb8ffkLqFl8YVp9C/hefBYOAMYYb7L2HJl7DfovsIcCH8JiLCzatYjj\n545fVRutWrUqUH0vLy/mzJnDjBkz8PDw4IcffiAoKIi5c0vkR6SUUqqcszcZqQtsBZYD68havbL3\nstc+wOZ8kgLKa5hoGlkjIG8bY5yMMRbgLeBvYGIR9V+kvKt4E+kfSVpmGt/99d0169cYQ7du3diy\nZQsdOnQgKSmJBx98kG7dupGYmHjN4lBKKaXyY28yclxE2uXxigCSCxOIMeYdY0w8cO/F93EXXxUu\n1RGRNKADWYlPAlkJUhWy5oucLUz/V4qJiSmyZbKFvVWTmwMHDnDq1Kl861133XUsXryYjz/+GDc3\nN2bOnElAQACzZs3SjdKUUqqci42NJSYmxtFh2D2BdT/wMzBLRHJ9FK0xprmIbCri+BzC1gTWMxfO\nMGbFGFIzUhl7R/7bwl+SfC4Z77HeZEgGB4cepHaV2oWK76+//uKOO+6gSZMmLFy4EFdX+zZU27Vr\nF/3792fp0qVA1kZpEyZMoF69vDdwU0opVbaVlgmsdcl6HsxhWxXKSiKSFyeLE2N+H8MHaz4o0KhC\n9UrVubvx3WRKJnO2Fn6er6urK+np6SxZsoTu3buTkWHfHbJGjRrx22+/8dlnn1nnkgQGBvLJJ5/k\nurGaUkopdS3Ym4wcE5GeIrK1WKMp4SpWqIi7qzvpmemcOH+iQOcW5a2aBg0asGjRItzd3Zk7dy4D\nBw60OzkyxtC3b18SEhLo0qULp0+f5sknn+S2225j27ZthY5NKaWUKih7k5HNxhifvCoYY74ognhK\nvIJsCX+5qKZRuDm7serAKvac2FPoOEJCQpg/fz6urq5MnDiRF154oUDn16lTh2+//ZY5c+bg7e3N\n77//TvPmzXn55Zc5e7ZIp98opZRSebI3GRkEfGKMaZFHnQ5FEE+J513ZG7B/S/hLKrtUJqppFACz\nt8wuklhuu+025syZQ4UKFahRo0aBzzfG8MADD5CQkEDfvn1JS0tjzJgxBAQEsGDBgiKJUSmllMqP\nvcnID0BbYK0x5owD9xm5ZmytprnakREonlU1nTp14q+//mL48OFX3UaNGjX47LPPWLlyJaGhoezd\nu5eoqCiioqLYs2dPkcWqlFKqZCltq2lSgZV5VQFai0jFogrMkfLaDn7B9gXsO7mPOxvdadfzaS6X\nmp5K7XG1OXH+BFuf3EpAzYCiCLdIpaenM2HCBEaMGMHp06epVKkSI0aM4Nlnn7V71Y5SSqnSxdGr\naexNRg6LSH5zRvKtU1oU5bNprtRvfj8mx0/mpZtf4o3b3yiWPorC4cOHefbZZ5k1K2sUp2nTpowf\nP567777bwZEppZQqao5ORuy9TdPdjjp3FiaQ8qJ786yP8n+b/0emFN9y2u3btxdq3oePjw8zZ87k\nt99+o2nTpmzfvp2OHTtyzz33sH379iKMVCmlVHlnVzIiIr/ZUafM7zNSFG7xvYV67vXYe3Ivf+wr\n9EOOc3XkyBFuvfVW7r//fn766adCtRUZGcmmTZsYO3Ys7u7uLFy4kKCgIJ599llOnCjY8mallFIq\nN/aOjABgjOlpjPnBGLPt4muBMcaeURN1kcVYeCT4EQBmbJpRLH14e3vTrVs30tLSuO+++1i8eHGh\n2nNxceHZZ59l586d9OvXj4yMDN577z0aN27Mp59+avema0oppVRu7J0z4gJ8D9yVy2EBFgOdLz47\nptQrzjkjAFuObSH4k2CqV6zO4WcP41qh6CeGiggDBw5kwoQJVKxYkQULFtC+ffsiaTs+Pp5Bgwax\nYsUKIGvPkw8++IDbbrutSNpXSil1bZWWOSMvAOHA80BToOrFV1Ng+MVjBdt1q5RKTU9l8KLB9JnX\n56rbCKoVRIh3CMnnk/lpV+Fuo9hijOHDDz/k8ccf5/z580RFRXHw4MEiaTssLIxly5Yxe/Zs6tev\nz8aNG4mIiOD+++9n586dRdKHUkqp8sPekZHtQLSIxNs43oKsh+g1KeL4HCKvkZFMycT1dVfSM9M5\n//L5qx7VGLtyLM//8jz3N7ufuQ/NLUy4ecrMzGTAgAEEBwczcODAIm//3LlzjB07ljfffJNz585R\noUIFnnjiCV599VVq1ixTW88opVSZ5eiREV3amwtjjIwcOZKIiAgiIiJyHK8zrg6HUw6zf8h+rnO/\n7qr6OHjqIPXG18PZyZmjzx2lWsVqhYzaNhHBmOL9Hjt06BCvvvoqU6dOJTMzE3d3d1544QUGDx5M\npUqVirVvpZRSVyc2NpbY2FhGjRpVKpKR/UALEcl121FjTG1gnYhc3W/mEia/OSOhE0PZeHQj6/uv\nJ9wn/Kr7uX367SzZvYTJnSbzWPhjV91OSbJ582aGDx9uXcVz3XXX8frrr9O9e3ecnJwcHJ1SSqnc\nOJm4DsIAACAASURBVHpkxN45I4uAb4wxYVceMMaEA3OAhUUZWEnmXeXi82lSCvZ8mit1D85aiDRj\nc/GsqslPWlrRzzcODg5m4cKF/Prrr4SGhnLgwAF69epFixYt+OWXX4q8P6WUUqWfvcnICKA+sM4Y\nc8gYs/7i6zDwJ3Ad8EpxBVnSFOb5NJfr2qwrFStUJHZPLPtO7iuK0Oy2adMmmjZtysqVee3yf/Vu\nv/121q9fz/Tp06lXrx4bN27kjjvu4K677mLDhg3F0qdSSqnSyd5Nz44CNwBTgUpA2MVXReBzoNXF\nOuVC79DeTO08lZvr31yodjwqetCpSScAZm0uuofn2WPSpEns3r2bu+66q9gSEovFQo8ePdi+fTtv\nv/02Hh4eLF68mLCwMB5++GFdeaOUUgqwc85IthOyZkJeWibxT7FuyOEgxb3PyOXmb59P5686E1wr\nmE0Drt0mtunp6fTo0YOvvvqKKlWq/D979x2f8/X+cfx1MiQIYjSxRRDU6ECNojGC2qNmlbTf2lUU\nNdJW0krMqlHa2ruK2nvGpq3dIoQkRlE7BCHJ+f0R8lMNgjs5951cz8cjj29z587n826/HnLlfM65\nLlavXk3lyi9XXD3L5cuXCQoKYsKECURHR2Nvb8///vc/vvrqK/LkyZOs9xZCCPFkpveMPHcx8sQL\nKVVca33UIhczLCWLkXux98j1bS6u3rnKwc4HKe1eOkXuC/EFSfv27Zk7dy4ZM2Zk5cqVKdK47MyZ\nMwQEBCScvHF2dqZ79+7069eP7NmzJ/v9hRBC/JvpYuS52sE/wzPn14j/SmefjhavtgCSrz38kzg4\nODBz5kzatWtHVFQUoaGhKXLffPnyMXnyZP766y/ee+897t69y4gRI/D09GTw4MHcunUrRXIIIYSw\nDkleGVFKNQIaA7kBx0TeUklr7WzBbMak5MoIwI7TO6g8rTJ5MuXhdK/T2ClL1ojPFhsby9atW6lW\nrVqK3vehvXv3MnDgQNatWweAm5sbX3zxBR07dsTJyfKt8oUQQvyb6ZWRpPYZ6QmMAmKAS8C9RN6W\nV2udWJFic1K6GNFa4znWk/Dr4Wxqt4lqBc0UBaZt3ryZAQMGsGfPHgAKFChAQECA9CgRQohkZroY\nSeqv4J8AAwEXrXUerXXBxz+Ay8kX07rE6Tg6Lu9I01+aYomiRSn1/z1HUvhRjTWpVq0au3btYsmS\nJZQoUYKIiAh8fX0pUaIE8+bNIy4uznREIYQQySCpKyOXtNZPHTSilPLWWgdbKphJSVkZyTI0C5HR\nkVz9/CpZ02d96Xseu3yM4uOLk9kpMxd6XyC9o/kW6nv37uXcuXM0bNgwxe8dGxvLnDlzCAgI4NSp\nUwCUKFGCgIAAmjRpgp1dyj7KEkKI1MxWVkb+VEo965hDmvrpYKnGZw8Vy1GMsrnLEhkdyZJjSyxy\nzZdx9uxZfHx8aNasGQsXJt8gvyext7enXbt2HDt2jEmTJpE/f/6EDa9lypRh2bJlFlmVEkIIYV5S\nC4juwASlVKmnvCdNPV+wdDEC4PuaLwAzDs6w2DVfVJ48eejYsSMxMTG0atWKX375xUgOR0dHPv74\nY44fP86ECRPIkycPBw4coFGjRpQvX541a9ZIUSKEEDYuqcXIMqAacEApdUspFa6UOvXoB//fCC1V\n8Pf3Jzg4+Ilfd8/4YD5NlOUaz7Yq2Yp09ulYf2o95yLPWey6L0IpxZAhQ/Dz8yM2NpY2bdowZ84c\nY3mcnJzo0qULoaGhjB49Gnd3d37//XfeffddKleuzMaNG6UoEUKI5xQcHIy/v7/pGEneMxINPK1n\nuAIqpKWjvZ1XdOanvT8xvu54upbrarF7N1/QnIVHFjKkxhD6V+5vseu+KK01X3/9Nf7+/jg4OHD8\n+HEKFixoOha3b99mwoQJDBs2jMuX4/dOV61alW+++YaqVasaTieEELbF9J6RpBYj57XWuV72PbYi\nKcXInrN7OH3jNGVzl6VgVsv9cF55fCX1f65P0exFOdrtKPHd980LCgoid+7c+Pr6mo7yLzdv3mTc\nuHGMHDmSa9euAVCzZk2+/vprKlasaDidEELYBlspRmporZ/aYVUpVVprnXLDVZJRSvcZeVRMXAx5\nR+XlYtRFdv1vFxXyVjCSw9bcuHGD0aNHM2rUKCIjIwF49913GTRoEOXLlzecTgghrJvpYiSpU3v/\nVYgopf7zEzK1FCKmOdg50LZ0fM+R6Qemmw1jQ7JkycKgQYMIDw/Hz88vYfhfhQoVqFOnDrt27TId\nUQghxBO80KA8pdTfWuvcyZDHKphcGQE4fPEwpX8sTRanLJzvfd4qeo48ydmzZ8mbN6/pGP9x+fJl\nRo0axbhx4xJm3fj4+DBo0CDefvttw+mEEMK62MTKSCKsYyNDEimlpiml/n7kY57pTE9Tyr0UZXKV\n4Ub0DZaFLDMd54l+++03Xn31Vb788kurO8mSI0cOgoKCElZKMmXKxPr166lcuTI+Pj5s377ddEQh\nhBAPvGgxYl0/eZ5Na61zP/LRynSgZ/F93ReA6QenG83xNGFhYdy+fZvBgwfTt29fqytIALJnz87g\nwYMJDw/niy++IHPmzGzYsIEqVapQo0YNtm7dajqiEEKkeS/6mCZZTs4opXIB04BaWmuLdXRVSk0D\n/gHsiS/ARmqt/37K+5P0mMZ3iS/nbp5jResVODlYdrrsldtXyPVtLmJ1LKd7niZP5jwWvb6l/Prr\nr7Ru3Zr79+/TrVs3xo4da9Wt2q9du8bo0aMZM2YMN27cAMDb25tBgwbh7e1tNpwQQhhiq49pmlo0\nBaCUakp8LxNPnrLyopR6RSk1Wyl1TCl1VCm1QCn1rJ/US4HvtNZ9gNXAdqXUS/dEWXdyHRtObeDS\n7Usve6n/yJ4hOw2LNiROx1n18LxmzZqxaNEi0qVLx/jx4+na1XI9V5JD1qxZCQgIIDw8HH9/f7Jk\nyUJwcDDVqlXD29ubzZs3W+UKjxBCpGYvVIxorZPjaEJfoCaw40lvUEo5AhsAR6A48CoQBWxWSmV4\n5H3dlVInlFLHlVLVtdZLtNYXHmRfDzgAlV82cHK0hH9U+9faA/GPaqz5B2T9+vVZsWIFLi4uVKtW\nzXScJHF1dU04fRMQEICrqytbtmyhevXqvPPOO2zatMmq/5sLIURq8sLr6Uqp1kqp75VSnyil0j94\nrbBSqvODVY7n9bbW+uQz3uMLlAQ+1w8A/YhfTeny8E1a63Fa6yJaay+t9SalVJHHrhMNvPQRleQu\nRuoUroN7RneOXT7GzjNPa4Brno+PD2FhYbRs2dJ0lOfi6urKV199RXh4ON988w1Zs2Zl27Zt1KhR\ng6pVq7JhwwYpSoQQIpm9UDGilBoEDAXcgTbEz6wpoLUOBZYDC573mlrruCS8rSlwWmsd8cj3XQSO\nAM2e8n0JQ1WUUp5AdmDP82Z8nLvLg/k0tyw3n+ZRjvaOCRtZJ+2blCz3sKQcOXKYjvDCsmTJwhdf\nfEF4eDiBgYFky5aN7du34+PjQ+XKlVm3bp0UJUIIkUxedGWkOFBUa91ca10JaAmMVUrlBWJIvqO/\npYGwRF4PA542UfiQUmqOUmoU8B3QXGv90ssZbhmSd2UE4OM3PwZg/l/zuX73erLdR8TLnDkzAwcO\nJDw8nKCgILJnz87OnTupXbs2FStWZMWKFVKUCCGEhb1oMbJHa3334Sda6wNAC6ATz9iA+pJyADcT\neT0SyKCUSvRIi9b6Y631+1rrz7TWjZ7V2j6pfF/3ZUXrFbQsmXyPJgpnK0yNgjW4E3OHOYfMTc19\nUVu3bqV+/foJjcdsRaZMmRgwYABhYWEMHTqUHDlysGfPHho0aECZMmVYtGgRcXFJWcwTQgjxLC96\ntLcp4Ar4A3W11n8+8rWOwASttcMLBYo/httOa22fyNeigTVa60aPvT6L+MdFGbTW0S9y38eupwcN\nGpTwube3t9Fjn/P/mk/LhS0p5VaKg50PWs3wvGeJjY2lRIkShISEULFiRVavXk2WLFlMx3ohUVFR\n/PTTT4wYMYILFy4AULJkSfz8/GjevDn29v/54yqEEFYrODiY4ODghM8DAgKsf1Beot8Yv/eiJLBa\na33/sa+9rbV+4qmYZ1z3acXIOSBEa139sdeXAtW11ple5J6J3MdoO/jH3Yu9R95Rebl0+xK7/7eb\n8nltZ/BbaGgoNWrU4PTp05QtW5a1a9eSLVs207Fe2N27d5kyZQrDhg3jzJkzAHh5eTFw4EDef/99\nHBxeqAYXQgijbLXPCEAh4veHxD7+hRctRJLgEOCRyOsFgcPJdE/j0tmnS9jIOnHvRLNhnlPhwoXZ\nunUrnp6e/PHHH1SrVo1//km+PTbJzdnZmW7duhEaGsrEiRMpWLAgx48fx9fXFy8vLyZNmsS9e/dM\nxxRCCJvywsXIg34ddsBypdRwpVRxy8V6okVAAaVU/ocvKKXcid9QuzAF7m/Mw42s8/6ax427Nwyn\neT4FChRg69ateHl5ceTIEQ4cOGA60ktLly4dHTp0ICQkhOnTp+Pl5UVYWBgdO3akcOHCjB8/nrt3\n7z77QkIIIZ5ejCilXn/a17XWK4BGQBHiVy0s4WnLRNOJXwEZppSyV0rZEX/E+BTwo4Xub5W8sntR\nzaMat+/fZu7huabjPLc8efKwdetWlixZQq1atUzHsRhHR0fat2/PkSNH+PnnnylRogRnzpzhk08+\noWDBgowaNYqoqCjTMYUQwqo9a2Wk77MuoLWOAT4G7j/rvU/zYHVlP1D/wef7HnwkPIR/sDfFh/hH\nQ0eAvwAX4veL3H6Z+z/O39//X5t7nqTtorZUmlIpRY7ddizTEYAf9/5ok8dL3d3dqVevnukYycLe\n3p5WrVpx6NAhFi1axBtvvMGFCxfo3bs3Hh4eDB06lMjISNMxhRDiX4KDg/H39zcd4+kbWJVSG7TW\nNZN0oed4r7V7ng2sRcYVIfRqKMe6HaNojqLJmis6Jpr8o/PzT9Q/bPHdQtUCVZP1fuLFaa1ZtWoV\n33zzDXv2xPfXy5o1Kz169ODTTz8la9ashhMKIcT/s/YNrNWVUpuUUn5KqUpKqaedX0yeNqRWLrlb\nwj/KycGJzmU6AzB2z9hkv19K2bhxI3/++eez32hDlFLUq1ePXbt2sW7dOqpUqcK1a9fw9/fHw8MD\nPz8/Ll++bDqmEEJYhWcVI5HEt3z/BtgGXFNKrVZKfa6UKqv+3fDiP6dq0oKULEYAOpftjIOdA4uP\nLSbiesSzv8HK7d27lwYNGuDt7c2+fftMx7E4pRQ+Pj5s3bqVLVu2ULNmTSIjIwkKCsLDw4O+ffsm\n9C0RQoi06lnFyDqtdQkgF/ABMB/wIn7T6B7gqlJqiVKqJ4kfuU313DM+mE8TlTILQ7ky5aJFiRbE\n6Tgm/D4hRe6ZnEqUKEG1atW4cuUK1atXZ9eu5BgIbR2qVq3K+vXr2blzJ3Xr1iUqKoqRI0dSsGBB\nPv30U86ePWs6ohBCGPGsYmQYxA+j01rPfdBWvRDxhcfHwAqgDDAKeDs5g1qrlF4ZAehRvgcQPzzv\n9n2L7ttNcc7OzixevJhmzZpx48YNfHx82Lx5s+lYyapixYqsXLmSP/74gyZNmnD37l3GjRtHoUKF\n6Ny5M+Hh4aYjCiFEinrhDqz/uohSxYAlWutiLx/JvIft4JPSBv7ElROcu3kOr+xe5M6UO2UCAhUm\nV2DPuT38VP+nhFM2tiwmJoaPPvqIWbNm4erqSlhYGK6urqZjpYjDhw8TGBjI/Pnz0Vrj4ODABx98\nwIABAyhSpIjpeEKIVOxhW3ibbQf/nwspNU9r3coiFzPM2trBJ2bu4bm8v+h9Xn3lVQ53OYydeplm\nutYhLi6OTz/9lNq1a9OgQQPTcVJcSEgIQUFBzJkzh9jYWOzs7GjZsiUDBw6kZMmSpuMJIVIx06dp\nLFmMeGmtj1vkYobZQjFyL/YenmM8OXfzHMtaLaNB0bT3wzu1OnnyJEOHDmXGjBncvx/fvqdx48b4\n+flRtmxZw+mEEKmR6WLEYr9Op5ZCxFaks09H74q9ARiyfYhNNkETiStUqBCTJk3i5MmTdO/eHWdn\nZ5YsWUK5cuWoU6cO27ZtMx1RCCEsymIrI6mJLayMANy6d4sCowtw9c7VVN8E7dSpU3h6epqOYcTF\nixcZNWoUEyZM4NatW0D8yRw/Pz98fHz49wl7IYR4fqlmZUSkPJd0LnR/qzsQvzqSWq1du5bixYsT\nGBhoOooR7u7uDBs2jIiICAYNGoSrqytbt26ldu3alC9fnqVLlxIXF2c6phBCvDApRiygza9teO3H\n1/j75t8pfu/ub3Ung2MG1oSu4cAF25+Gm5jz589z//59vvjiCwYOHJhmH0lly5YNf39/IiIiGDp0\nKG5ubvz+++80btyY119/nXnz5hEbmyZ7DwohbJwUIxZw5NIRDl08xIVbKd9JM3uG7HR8M/5ob+C2\n1Lly4Ovry9y5c7G3t2fIkCH07NkzzRYkAJkzZ6Zfv36EhYUxZswY8ubNy+HDh2ndujXFixdn2rRp\nCRtfhRDCFkgxYgEmGp89qk+lPjjZO7HwyEL2n99vJENya9WqFb/++ivp0qVj7Nix9OjRw3Qk4zJk\nyMCnn35KaGgoEydOxNPTkxMnTvDRRx9RuHBhJkyYwN27d03HFEKIZ5Ji5An8/f0JDg5O0ntNFyN5\nMueha7muAHy5+UsjGVJCo0aNWL58OVmzZqVRo0am41gNJycnOnToQEhICLNnz6Z48eKcPn2abt26\nUbBgQUaOHJmw8VUIIR4VHByMv7+/6RhymiYxz3uapvfa3ozaPYrhNYfT9+2+yZjsyf6J+gfPMZ5E\n3Y9i50c7qZivopEcKSEyMpLMmTObjmG14uLiWLJkCYMHD2b//viVsmzZstGzZ0+6d++eZjrbCiGS\nTk7TpAKmV0YeZuhZoScAfpv8jOVICVKIPJ2dnR1NmzZl7969rFq1ikqVKnH16lW++uorChQowMCB\nA7l06ZLpmEIIkUBWRhLxvCsjZ26c4fyt83i4eiQUJiZcv3udgmMKcv3udda8v4bahWsby2JCXFwc\ndnZSXz9Oa82WLVsIDAxkw4YNAKRPn55OnTrRp08f8uTJYzihEMI0WRlJBfJlycdbed4yWogAuDq7\nMqDyAAA+W/cZ92PTzomKJUuW8M4773Dt2jXTUayOUgpvb2/Wr1/P7t27adCgAXfu3GH06NF4enrS\nuXNnwsLCTMcUQqRhsjKSCFvpwJqY6JhoSkwowclrJxlbZyzdy3c3HSnZ3b9/n1KlShESEsJrr73G\nunXrcHMzWxhau4MHDxIUFMSCBQvQWmNvb8/777/PgAEDKFYsVQzfFkI8B1kZERbl5ODEt7W+BWBQ\n8CCu3L5iOFHyc3R0ZP369Xh5eXHw4EHeeecdzp07ZzqWVXvttdf45ZdfOHLkCO3btwdg5syZvPrq\nq7Ro0YIDB1JnAz0hhHWSYiQVali0ITUK1uDa3Wt8tfkr03FSRL58+di6dSslS5bk2LFjVK1alfDw\ncNOxrF6xYsWYPn06J06coHPnzjg6OrJgwQLeeOMN6tevz65du0xHFEKkAVKMpEJKKUbXGY29sueH\nP35g99ndpiOlCHd3d4KDgylTpgwXLlzg779Tvj2/rSpYsCA//PADp06dolevXmTIkIGVK1dSqVIl\natSowebNm9N011shRPKSPSOJeJE9I60WtmLv+b2sarOKItmLJFOy59N/Q3+G7RhGiVdKsK/TPtLZ\npzMdKUXcuHGDo0ePUqFCBdNRbNalS5cYPXo033//PZGRkQBUrFgRPz8/6tatK5OChUhlZM9IKnEm\n8gyhV0ONzKd5kkHvDKJQ1kL8dekvhu8YbjpOismSJYsUIi/plVdeITAwkIiICL755huyZ8/Orl27\nqF+/PmXKlGHhwoUyKVgIYTFSjFiINTQ+e1x6x/RMajAJgG+2fsOxy8cMJxK2xtXVlS+++ILw8HBG\njhxJzpw52b9/P82bN6dkyZLMmjWLmJgY0zGFEDZOihELcctgfcUIQLWC1fjo9Y+4F3uPDxZ/kKZ6\njzzu559/ZuPGjaZj2CQXFxd69+5NWFgYEyZMIH/+/Bw9epR27dpRtGhRJk6cSHR0tOmYQggbJcWI\nhbi7uANwMeqi4ST/Nar2KApkKcAff/9BwJYA03GM+P333/nggw+oV68eK1euNB3HZjk7O9OlSxdC\nQ0OZNm0aRYoU4dSpU3Tq1IlChQoxZswYbt++bTqmEMLGSDFiIdb4mOahLM5ZmNVkFgrFkO1D2H56\nu+lIKa5MmTJ07NiR6OhoGjduzIIFC0xHsmmOjo74+vpy9OhR5s2bR6lSpTh37hw9e/bEw8ODoUOH\nJmx8FUKIZ5HTNIlQSulBgwbh7e2Nt7d3kr7nn6h/uBR1iTyZ8+DqbJ1TUQduHMiQ7UPwcPXgQKcD\nZHHOYjpSitJa8/nnnzNy5Ejs7OyYNm0a7dq1Mx0rVYiLi2PFihUEBgby22+/AfH7Tbp3706PHj3I\nnj274YRCiMQEBwcTHBxMQECA0dM0UowkwpbbwT/Nvdh7VJpSib3n99K0eFMWNl+Y5o5oaq35+uuv\n8ff3x9XVlZMnT5ItWzbTsVINrTUbNmwgMDCQLVu2AJAxY0a6dOlC7969yZkzp+GEQojEmD7aK8VI\nIlJrMQIQejWUMhPLEBkdyUifkfSu1Nt0JCO+++473nrrLd5++23TUVKt7du3ExgYyJo1awBwcnLi\n448/5vPPPyd//vyG0wkhHiXFiBVKzcUIwJJjS2jySxPslT2b22+mSoEqpiOJVGzv3r0EBgayePFi\nABwcHPjggw/o378/Xl5ehtMJIcB8MSIbWNOgxsUa83mlz4nVsbRc2NKqGrWJ1KdMmTIsWrSIw4cP\n06ZNG+Li4pg2bRrFixenVatWHDp0yHREIYRhaWJlRCnlBAQBCkgHlNZaV33K+1P1yghATFwMNWfW\nZEvEFt7O9zYb223EycHJdCzjNm3aRNWqVXFwcDAdJdUKDQ1l2LBhzJgxg/v34/veNGjQAD8/P8qX\nL284nRBpk6yMpIyhwCat9Wda60+Afslxk1YLW5H729wcuGD949cd7ByY99488mbOy44zO+i8snOa\nH4T2888/U7NmTVq0aCENvJJR4cKFmTRpEidPnuTTTz8lffr0LF++nAoVKlCzZk0ZyidEGmRVxYhS\nKpdSao1SymJDL5RSzsBHQGalVKBSagJw01LXf9TVO1c5f+u8zTz2yOmSk2WtlpHBMQPTD0xn1K5R\npiMZ5eHhQZYsWVi8eDENGjQgKirKdKRULV++fIwZM4bw8HAGDBhApkyZ2LhxI9WrV+ftt99m5cqV\nUpQIkUZYTTGilGoK7AQ8gSf+DaSUekUpNVspdUwpdVQptUAplecply4IZAIKaa39gHHAJqVUVkvm\nB+tufPYkb+R6gxmNZwDQd31fVp1YZTiRORUrViQ4OBg3NzfWr19PrVq1uH79uulYqZ6bmxtBQUGc\nPn36P0P53njjDebPn09sbKzpmEKIZGQ1xQjQF6gJ7HjSG5RSjsAGwBEoDrwKRAGblVIZHnlfd6XU\nCaXUcaA+8cXNYgCt9VHgHOBj6X8BWyxGAN579T0CvAPQaFotbMVf//xlOpIxr732Gtu2bSNfvnzs\n3LmTFi1amI6UZjw6lO/bb78lV65cHDx4kJYtW1KiRAmmT5+esMdECJG6WFMx8rbW+uQz3uMLlAQ+\n1w8Qv//DE+jy8E1a63Fa6yJaay9g9oOXH330cw+w+G5N94wP5tPcsr75NM/yZdUvaVmiJTfv3aTu\n3Lqcv3nedCRjvLy82L59O2XLlmXIkCGm46Q5Li4ufPbZZ5w6dYoffvgBDw8PQkJC+PDDDylSpAgT\nJkzg7t27pmMKISzIaooRrXVS9ok0BU5rrSMe+b6LwBGg2ROue5741ZYqAA8ez3gCW1428+MSVkZu\n29bKCMTvpJ7aaCoV8lbg9I3T1Jtbj1v3bpmOZUz+/Pn57bffKFOmjOkoaZazszOdO3fm+PHjzJgx\ng2LFihEREUG3bt0oWLAgI0eO5NattPtnVIjUxOqO9iqlpgHttNb2iXztHBCita7+2OtLgepa60xP\nuGZ+YAwQBuQFZmqtVzwlwwsd7b1+9zrX717HLaMbGRwzPPsbrNClqEtUmlqJ0KuhvFv4XZa1XoaD\nnRxzFebFxcWxaNEigoKC2L9/PwDZsmWjR48edO/enaxZLb4NTIg0w/TRXlsrRqKBNVrrRo+9Pgto\nA2TQWr/0mcy00GfkaUKvhlJxSkUu377Mx298zMQGE9PcDJun+eeff3BzczMdI83SWrNmzRoCAwPZ\nsSN+i5mLiwtdu3bls88+w93d3XBCIWyP6WLEah7TCOtROFthlrdejrODM5P3TyZwW6DpSFZj0qRJ\nFC5cmM2bN5uOkmYppXj33XfZtm0bwcHB+Pj4cOvWLYYPH46Hhwfdu3fn9OnTpmMKIZ6Dra2MvNBj\nmhfIoAcNGpTwube3N97e3pa4tE1ZfHQxzeY3Q6OZ0XgG7V5rZzqSUVprPvzwQ2bMmIGTkxMLFiyg\nQYMGpmMJ4PfffycwMJClS5cC8fNv2rVrR//+/SlSpIjhdEJYn+DgYIKDgxM+DwgIkMc0j3pGMbIa\nKKq19nzs9UPALa11JQtlSNOPaR41bs84Pl3zKQ52Dqxss5JahWqZjmRUXFwcn3zyCT/88AP29vbM\nnDmTNm3amI4lHjh8+DBDhgzhl19+IS4uDjs7O1q0aMHAgQMpVaqU6XhCWC15TPN8FgEFHmxIBUAp\n5U58z5GFxlKlYt3Ld6dPxT7ExMXQ9Jem/H7ud9ORjLKzs2P8+PH079+f2NhY2rZty4wZM0zHEg+U\nKlWKuXPnEhISwv/+9z/s7e2ZN28epUuXplGjRvz222+mIwohEmGNxcjTKrPpwGFgmFLKXillP+qg\ndwAAIABJREFUR/zcmVPAjymQ7ZlaLWyF61BXtkVsMx3FYob5DKNt6bZE3Y+i7ty6HL9y3HQko5RS\nDBkyhCFDhpA9e3bKlStnOpJ4TOHChZk8eTInT56ke/fuODs7s2zZMsqXL4+Pj4/MvxHCylhNMaKU\nGq6U2k98x1SUUvsefCScK9Va3ye+c2os8b1F/gJciN8vcttA7P+Ijo3mRvQNm+vC+jR2yo6pDadS\np3AdLt++TK1Ztfj75t+mYxnXv39/jh49yquvvmo6iniCfPnyMXbsWCIiIujfvz+ZMmViw4YNMv9G\nCCtjNcWI1vpzrfUbWuscWmt7rfWbDz5iHnvfJa11W611Ua11ca11c631OUvn8ff3/9fmnqRyy2Cb\nLeGfxdHekYXNF1I+T3kibkRQZ3Ydrt+VuS05cuQwHUEkgZubG0OGDCEiIoKvv/6abNmyJcy/efPN\nN1mwYIHMvxFpUnBwMP7+/qZjWN8GVmvwMhtYv9z0JYO3Dcb/HX8GeQ969jfYmCu3r1B5WmWOXT5G\nlfxVWNt2Lekd05uOZXViYmJwcJBmcdbq1q1b/PTTT4wcOZILF+KnbBctWpQBAwbQpk0bHB0dDScU\nImXJBtZUxt3lwXyaKNubT5MU2TNkZ23bteTJlIdtp7fR+tfWxMTFPPsb05BRo0ZRu3Ztbt68aTqK\neAIXFxd69+5NWFgYEyZMoECBAoSEhODr64uXlxc//PCDzL8RIgVJMWJhtjq593nkz5KftW3XktU5\nK0tDltJlRRd57v5AZGQkI0eOZNOmTVSrVo1//km9fw5SA2dnZ7p06cKJEycS5t+Eh4fTtWtXPD09\n+fbbb2X+jRApQB7TJOJlHtNE3Yvi9v3bZEufDXu7/7RKSVV2ntlJzZk1uRNzB78qfgyuPth0JKtw\n6tQpfHx8OHXqFF5eXqxbt44CBQqYjiWSIDY2lsWLFxMYGMiBAwcAmX8j0gbTj2mkGEmEND1LuhXH\nV9B4XmNidSxj64yle/nupiNZhQsXLlCnTh0OHjxI7ty5Wb9+vZy6sSFaa1avXk1gYCA7d+4EIFOm\nTHTt2pVevXrJ/BuR6pguRuQxzRO86GmatKa+V30mN5wMQI81PZj35zzDiaxDzpw52bJlC1WrViUu\nLg5nZ2fTkcRzUEpRt25dtm/fzubNm6lZsyY3b95k2LBheHh48Omnn3LmzBnTMYV4aXKaxorJysjz\nG75jOP029MPRzpGVbVbiU8jHdCSrcOfOHc6ePSvzUVKB3377jcDAQJYtWwaAo6NjwvybwoULG04n\nxMsxvTIixUgipBh5flpreq/rzXe7vyOjY0Y2td/EW3neMh1LCItLbP5Ny5YtGTBggMy/ETZLihEr\nZIliRGuNRmOn0s6TsDgdR/sl7Zl9aDbZ02dn24fbKP5KcdOxhEgWJ06cYNiwYcycOZP79+8D0LBh\nQ/z8/HjrLSnEhW0xXYyknZ+UKajVwlakG5yOdSfXmY6Soh62ja9XpB5X7lyh9uzanLkhz9UTExAQ\nwJdffilHom1YkSJFnjr/Jjg4WP7/FSKJpBhJBg52DsTExaTqXiNP4mjvyPzm83k739uciTxD7dm1\nuXL7iulYVuXYsWN88803DB48mA4dOhATI03jbNnD+Tfh4eH069cvYf5NtWrVqFy5MqtWrZKiRIhn\nkGIkGaSFxmdPk8ExA8tbL6eUWymOXj5K3bl1uXVPGkc9VKxYMRYtWkT69OmZMmUKjRo1ksZaqYC7\nuztDhw4lIiKCgIAAsmXLxs6dO6lXrx5lypRh4cKFxMXFmY4phFWSYiQZpPViBCBr+qysabsGD1cP\nfjv3G83mN+Ne7D3TsaxGw4YN2bRpE9mzZ2fVqlXSrTUVyZo1K1999RURERGMHDmSnDlzsn//fpo3\nb06JEiX+tcdECBFPipFk4J4xdc+nSarcmXKzru063DK6se7kOtotbkdsnExGfahChQrs3LmTggUL\ncunSJZkam8okNv/m2LFjtG/fXubfCPEYKUae4GWanj1cGbl+97oFE9mmItmLsOb9NWRKl4lf/vqF\nT1d/Ks/PH+Hl5cWuXbtYt24duXLlMh1HJINH599Mnz6dokWLyvwbYTWk6ZkVe9mjvdEx0cTExZAx\nXUYLprJtweHB1Jldh+jYaPzf8WeQ9yDTkYQwIjY2lkWLFhEYGMjBgwcByJ49e8L8G1dXV8MJRVpk\n+mivFCOJkKZnyWPJsSU0m9+MOB3H+Lrj6Vquq+lIVk1rjVLG/m4QyUxrzapVqwgMDGTXrl1A/Pyb\nbt260atXL9zc3AwnFGmJFCNWSIqR5DNl3xQ+Xv4xCsXcZnNpVbKV6UhWSWtNr169cHFx4ZtvvpGi\nJBXTWhMcHExQUBAbNmwAIH369HTo0IE+ffqQL18+wwlFWiDFiBWSYiR5Dds+jP4b++No58iKNiuo\nVaiW6UhW5/Dhw7zxxhvExsbywQcfMGnSJJycnEzHEslsz549BAUFyfwbkeKkGLFCUowkL601fdb1\nYdTuUWRwzMCmdpson7e86VhWZ+XKlbRs2ZKoqCiqVKnCokWLyJEjh+lYIgUcOnSIIUOGMH/+/H/N\nvxk4cCAlS5Y0HU+kQlKMWCFLFSMP+2qks0/30tdKbeJ0HB8u/ZCZB2eSLX02tn+4XebYJGL//v00\naNCAc+fOUahQIdauXUuhQoVMxxIp5MSJEwwdOpSZM2cmdOpt1KgRfn5+lCtXznA6kZqYLkbkaG8y\neX/R+zgNdmJ5yHLTUaySnbJjcoPJ1Peqz9U7V6k1uxanb5w2HcvqvPHGG+zZs4c33ngDBwcHsmXL\nZjqSSEFFihRhypQpnDx5kk8++QRnZ2eWLl3KW2+9Ra1atdiyZYsclRepghQjycTF0QVI211Yn8XR\n3pH5782ncv7KnI08S61Ztbh8+7LpWFYnT548bNu2jXXr1pE1a1bTcYQB+fPnZ9y4cQnzb1xcXFi/\nfj3e3t5UqVKF1atXS1EibJoUI8lEWsInTXrH9AlzbEKuhFB3Tl1uRt80HcvqZMyYkfz585uOIQx7\nOP/m9OnTCfNvduzYQd26dWX+jbBpUowkEylGks7V2ZW1bddS0LUgv//9O03nNyU6Jtp0LJtw7949\nbt++bTqGSGEP59+Eh4czYsQImX8jbJ4UI8nE3UXm0zyPXJlyse6DdbhndGfDqQ18sPgDmWPzDFpr\nOnfuTOXKlYmIiDAdRxiQKVMm+vTpQ1hYGOPHjyd//vz/mn/z448/yvwbYROkGEkmbhndUCiiY+U3\n/KQqnK0wa9quIbNTZhYcWcAnqz6R5+BPcfnyZbZt28b+/fspW7YsW7ZsMR1JGOLs7EzXrl0JDQ1l\n2rRpeHl5ER4eTpcuXfD09GTUqFFERUWZjinEE8nR3kRY4mjvw9/q7e3sLREpTdkSvoXas2sTHRvN\nV1W/IqBagOlIVuvatWu0atWKdevW4eDgwOjRo+natat0bE3jYmNj+fXXXwkKCvrX/JuePXvyySef\nyPwb8R+mj/ZKMZIIaXpm3tJjS2k6vylxOo6xdcbSvXx305GsVmxsLAMGDGDEiBEA/PTTT3Ts2NFw\nKmENZP6NSCopRqyQUkoPGjQIb29vvL29TcdJs6btn8ZHyz4CYG7TubQu1dpwIus2d+5cxo4dy4YN\nG3BxcTEdR1iRh/NvAgMD2bhxI/D/82/69u1L3rx5DScUpgQHBxMcHExAQIAUI9ZGVkasx/Adw+m3\noR8Odg4sb72cOoXrmI5k1R62DhfiSXbv3k1QUBDLl8c3ZHR0dKR9+/b0799fuvumYbIyYoWkGLEu\nfdf1ZeSukWRwzMDGdhupkLeC6UhC2LxDhw4RFBTE/Pnz0VpjZ2dHq1atGDhwICVKlDAdT6Qw08WI\n/AqVjLTWXLtzjdv3pQ/EyxjuMxzf1325ff829ebW48ilI6Yj2ZSoqCj69OlDZGSk6SjCipQuXZp5\n8+Zx7NgxPvzwQ+zs7Jg7dy4lS5akSZMm/PHHH6YjijREipFk1H5Je7INz8bCIwtNR7FpSikmNZhE\nw6IN4+fYzKpFxHXpq5FUn332Gd9++y1ly5bl0KFDpuMIK+Pl5cXUqVP/Nf9myZIllCtXjtq1a7N1\n61bTEUUakCaKEaXUbaXU3w8+ziulopRSzZL7vjkyxI97ly6sL8/BzoF5zeZRtUBVzt08R63ZtbgU\ndcl0LJvQp08fSpcuzYkTJyhfvjzTpk0zHUlYoYfzb8LCwvj8889xcXFh3bp1vPPOOzL/RiS7NFGM\nAJO01rkffOQC9gMrk/um0hLestI7pmdZq2W85v4ax68cp+5cmWOTFEWKFGH37t189NFH3L17l48+\n+ghfX9+EkfRCPCpnzpwMGzaMiIgI/P39yZo1K9u3b6du3bqULVuWX3/9VebfCIuzqmJEKZVLKbVG\nKWXRP+la6x6P3KM2sF1rnew9kqUYsbwszllY03YNnlk9+ePvP2j8S2OZY5ME6dOnZ8qUKUydOhVn\nZ2fs7OxwcHAwHUtYsWzZsjFo0CAiIiIYPnw47u7u7Nu3j/fee4+SJUsya9YsKWiFxVhNMaKUagrs\nBDyBJ64FKqVeUUrNVkodU0odVUotUErleY5bdQXGv2TcJHHPKPNpkkNOl5ysaxs/x2ZT2CbaLm4r\nc2yS6MMPP2Tfvn2MHTvWdBRhIzJlykTfvn0JCwvj+++/J3/+/Bw9epR27drh5eXFTz/9RHS0/EIg\nXo7VFCNAX6AmsONJb1BKOQIbAEegOPAqEAVsVkpleOR93ZVSJ5RSx5VS1R95vTBwV2t9Jpn+Hf7F\nLaMb6R3SY6+kJbylFcpWiLVt15LFKQsLjyyk26pu8jw7iYoXLy5N0cRzS58+Pd26dSM0NJSpU6fi\n5eVFWFgYnTt3xtPTk++++07m34gXZjV9RpRSdlrrOKXUNKCd1vo/P8GVUh2AHwFPrXXEg9fcgXNA\nP631t8+4xxhgodZ62zPeZ5E+Iw+vIXNCks+2iG3Uml2LuzF3+aLKF3xT/RvTkWzWwYMHCQ0NpVmz\nZN/bLVKBh/NvAgMDE05p5ciRg549e9KtWzeZf2NjpM/IA1rrpOwTaQqcfliIPPi+i8AR4Kl/gyql\nXIDXn1WIWJJSSgqRZFalQBXmvzcfe2XP4G2DGbtHHj+8iDt37tC6dWvee+892rRpw5UrV0xHElbO\n3t6eFi1acODAAZYvX06FChW4fPkyX3zxBQUKFGDgwIFcuiQn3kTSWE0xkkSlgbBEXg8DSj3je32B\n6RbOI6xAg6INmNxwMgA91vRgzqE5hhPZHicnJ7p160aGDBn4+eefKVGiBEuXLjUdS9gApRT169dn\n586dbNy4kerVqxMZGcmQIUMoUKAAPXv25OzZs6ZjCitna8VIDiCxs5yRQAallNOTvlFr/b3WWhos\npFK+r/sywid+aq3vUl9Wn1htOJFtsbOzo1u3bhw6dIiqVaty8eJFGjduTPfuMi1ZJI1SiurVq7Nx\n40Z27dpF/fr1uXPnDmPGjMHT05OOHTty8uRJ0zGFlbK1YkSIJ+pTqQ+fV/qcmLgYms1vxq4zu0xH\nsjmFChVi8+bNjBkzhvTp01O+fHnTkYQNqlChAsuXL+fAgQO0aNGCmJgYJk2ahJeXF23btuWvv/4y\nHVFYGavZwPrQMzawngNCtNbVH3t9KVBda53JQhn0oEGDEj739vbG29v7ha51L/Yel6Iu4ZLOhSzO\nWSwRTzyF1pqPl33M1ANTyeqclW0fbqOEmwz9ehHnzp0jd+7csu9JvLSQkBCGDRv2r94kTZo0wc/P\njzJlyhhOlzYFBwcTHByc8HlAQIBM7X3UM4qR1UBRrbXnY68fAm5prStZKIPFpvZ2Wt6JifsmMqHu\nBLqU62KRa4qni4mLofmC5iw5toTcmXKz46MdeLh6mI6Vajz8YSJN08TzioiIYMSIEUyePDmhN0nt\n2rXx8/OjSpUqhtOlbXKa5vksAgoopfI/fOHB0d7igFVOo5MurCnPwc6Bn5v9zDsF3uHvm39Ta1Yt\n+e9vQePGjaNcuXL8/vvvpqMIG1OgQAG+//57wsPD6du3Ly4uLqxdu5aqVatSpUoV1qxZI/2C0ihr\nLEaeVplNBw4Dw5RS9kopO2AocIr4/iNWR4oRM5wdnFnaaimv53ydE1dP8O6cd4mMjjQdy+bFxsYy\ndepUDhw4QPny5enQoQP//CN/tsXzyZkzJ8OHDyciIoJBgwYlzL959913KVeuHIsWLZL5N2mM1RQj\nSqnhSqn9QP0Hn+978JGwFqy1vg/4ALHE9xb5C3Ahfr/IbUvm8ff3/9fztBeVUIzclr+wU1oW5yys\neX8NhbIWYt/5fTSe15i7Mck+kihVs7e3Z/fu3Xz++efY29szefJkihQpwrfffktsrLTkF88nW7Zs\n+Pv7ExERwbBhw3B3d2fv3r00a9aMUqVKMXv2bJl/k8yCg4Px9/c3HcP69oxYA0vuGQkOD6bajGpU\nyV+FrR9utcg1xfM5de0Ub099mwu3LtC0eNP4Jml20qL/ZYWEhNC7d29WrlzJO++8w+bNm2Wzq3gp\nd+7cYcqUKQwfPpwzZ+Kndnh6etKvXz/at2+Pk9MTuzeIl2R6z4gUI4mwZDFy9NJRqk6vSvk85VnR\nZoVFrime36GLh6g6rSo3om/Q4c0O/FT/J/nBaSGrV68mb968lCr1rL6DQiTNvXv3mDNnDkOGDOHE\niRMA5MmThz59+tChQwcyZsxoOGHqI8WIFbJkMSKsx/bT2/GZ5cPdmLsMrDyQwBqBpiOlejdv3iRT\nJoucuBdpUGxsLAsXLiQoKOhf82969epFt27dyJJF2iVYiulixGr2jAiR3Crnr8yC5guwV/YEbQ9i\n9O7RpiOlaufOnSNfvnz06NFDZpSIF2Jvb0/Lli05cOAAy5Yto3z58ly+fBk/Pz/y58+Pn5+f/NlK\nJaQYEWlKfa/6TG00FYBea3sx6+Asw4lSr3Xr1hEZGcnYsWMpVKgQX3/9NTdu3DAdS9ggpRQNGjRg\n165dbNiwgWrVqhEZGUlQUBAeHh706tWLc+fOmY4pXoI8pkmEPKZJ/UbtGkXvdb2xV/YsbbWUel71\nTEdKlQ4dOsSAAQNYtWoVAK6urkyfPp1GjRoZTiZs3a5duwgMDGTlypUApEuXDl9fX/r164enp+cz\nvls8Th7TWClLHe0V1umzip/R/+3+xOpYmi9ozo7TO0xHSpVKly7NypUrCQ4O5p133iEyMpLixYub\njiVSgYoVK7JixQr2799P8+bNuX//PhMnTsTLy4sPPviAI0eOmI5oE+RorxWz9MrIzeibnL91Hldn\n14S+I8I8rTUdl3dk8v7JuDq7stV3K6Xc5URIcjp27BjFihUzHUOkQiEhIQwdOvRfvUmaNm3KwIED\nZf5NEsjKSBowKHgQRb8vysyDM01HEY9QSvFD/R9oUqwJ1+9ep/bs2oRdCzMdK1V7UiFy8OBB+vTp\nw9mzZ1M4kUgtihYtyrRp0zhx4gRdu3bFycmJRYsWUbZsWerUqcO2bdtMRxRPIcVICpCW8NbLwc6B\nuc3m4u3hzflb56k1uxYXb100HSvNGT58ON9++y2enp58+OGHssQuXpiHhwfjx48nLCyMPn36kDFj\nxoT5N1WrVmXt2rUy/8YKSTGSAqQYsW4P59i8kfMNQq+GyhwbAz777DNatGhBbGws06dPp0SJEjRs\n2JBTp06ZjiZsVK5cuRgxYgQRERF89dVXuLq6sm3bNurUqUO5cuVYvHixzL+xIlKMpAApRqxfZqfM\nrH5/NYWzFWb/hf00mtdI5tikoDJlyvDLL79w/PhxunTpgrOzMxs2bJCGaeKlZc+enYCAACIiIhg6\ndChubm7s3buXpk2bUrp0aebMmSPzb6yAFCMp4GExcjFKlv+tmbuLO+s/WE8ul1wEhwfT5tc2xMTJ\nX1IpqVChQkyYMIGIiAjmz5/PK6+8YjqSSCUyZ85Mv379CA8PZ+zYseTLl4+//vqLtm3bUrRoUSZO\nnEh0dLTpmGmWFCMpwD2jO/mz5Cd3ptymo4hn8HD1YG3btbg6u7L42GI6r+gsz5cNcHNzo379+ol+\nbe3atXTs2JE///wzhVOJ1CB9+vR0796d0NBQpkyZQpEiRTh16hSdOnWiUKFCjB49mqioKNMx0xw5\n2psIpZQeNGgQ3t7eeHt7m44jDNhxegc+s3y4E3OH/m/3Z0jNIaYjiQdq167NunXrAKhRowY9evSg\nbt262NvLJGbx/GJjY1mwYAFBQUEcPnwYSFvzb4KDgwkODiYgIEAG5Vkb6cAqAFYeX0mjeY2I1bF8\nW+tbPqv4melIgvheJePGjWPGjBkJv8F6eHiwYsUKSpQoYTidsFVxcXGsWLGCwMBAfvvtNwCyZMnC\nJ598Qs+ePcmRI4fhhMnLdJ8RKUYSIcWIeGj2odl8sPgDAGY0nkG719oZTiQeun79OlOmTGH8+PFE\nRkZy9uxZnJ2dTccSNk5rzcaNGwkMDEzowp0hQwY6depE7969yZMnj9mAyUSKESskxYh41Ojdo+m1\nthf2yp4lrZZQ3yvxvQzCjLi4OEJDQ/Hy8vrP1+7evUt0dHSqX2oXyWPnzp0EBQWlifk3posR2cAq\nxDP0rNCTgZUHJsyx2X56u+lI4hF2dnaJFiIAc+bMIXfu3HTo0IF9+/alcDJh6ypVqiTzb1KIFCMp\n5Oqdq+w/v58zN86YjiJewODqg+nwZgfuxtyl/tz6HLp4yHQkkQQHDhzg9u3bTJ48mTJlylCuXDkm\nTpxIZKQ0tRNJ9/rrrzN//nyOHDlC+/btAZg9ezYlSpSgWbNm7N2713BC2yfFSAoZtWsUb058k6n7\np5qOIl6AUoof6v1A0+JNuRF9g9qza3PqmnQHtXbjxo3j6NGj9OzZE1dXV/744w86derE7t27TUcT\nNqhYsWJMnz6d0NBQmX9jYVKMpBDpwmr77O3smdN0DtU8qnHh1gVqzZI5NragWLFifPfdd/z999/M\nmjWLli1bUrNmTdOxhA2T+TeWJ8VICkkoRm5LMWLLnB2cWdJqCW/mepOT105SZ04dbty9YTqWSIL0\n6dPTtm1b5s2bh53df//qu3DhAi1atGDt2rUys0Qkicy/sRwpRlKIrIykHg/n2Hhl9+LAhQMyxyaV\nmDlzJgsWLKBOnTp4enoSEBDA6dOnTccSNkDm37w8KUZSiHtGdwBZ1k8l3DK6sbbtWnJnys2WiC20\nWthK5tjYuPfff5/BgwdTsGBBIiIi8Pf3x8PDg3HjxpmOJmyEzL95cVKMpBB3F3defeVViuYoajqK\nsJCHc2yyOmdlachSOi3vJM+JbViePHnw8/MjNDSUDRs20KpVKxwdHalYsaLpaMLGyPyb5ydNzxIh\nTc/E89h5Zic1Z9bkTswdPq/0OcN8hpmOJCzk2rVruLq6otR/e0HNnTuXmjVr4ubmZiCZsCW2MP/G\ndNMzKUYSIcWIeF6rT6ym4byGxMTFMMJnBH0q9TEdSSSjEydO4OXlhYODA/Xq1cPX15d69erh6Oho\nOpqwYonNv8mcOXPC/JtXXnnFWDbTxYg8phHCAt4t8i7TG00HoO/6vsw4MMNsIJGsoqOjadCgAVpr\nli5dSpMmTciTJw9Dhw41HU1YMTs7Oxo2bMju3bvZsGED1apVIzIykqCgIDw8POjVqxfnzp0zHdMI\nKUaewN/fP2FIkhBJ8X7p9xlTZwwA/1v2P5aHLDecSCSXkiVLsmzZMs6ePcuIESN49dVXuXTpEjdv\n3jQdTdgApRQ1atRg06ZN7Ny5k3r16nH79m1Gjx6Np6cnnTp14tSplGmqGBwcjL+/f4rc62nkMU0i\n5DGNeBlfbvqSwdsG4+zgzNq2a6laoKrpSCKZaa35448/yJUrF3nz5v3P1yMiIsiTJw8ODg4G0glb\ncODAAYKCgli4cCFaa+zt7WndujX9+/enRIkSyX5/049ppBhJRHIVI/9E/cOJKydwy+hGkexFLH59\nYR201nRZ2YWf9v5EZqfMbPXdyms5XzMdSxiitaZ06dJcuXKFdu3a4evrS7FixUzHElYqJCSEoUOH\nMnv27ITeJE2aNMHPz48yZcok231NFyPymCYFTdk3hcrTKjNp3yTTUUQyUkoxvu543nv1PSKjI2WO\nTRp3+fJloqOjOX/+PMOGDaN48eJUqFCBH3/8URphif8oWrQo06ZN48SJEwnzbxYvXpzq599IMZKC\npAtr2mFvZ8/sJrOpUbAGF6Mu4jPLhwu3LpiOJQx45ZVXCAkJYfv27fzvf//DxcWFPXv2MGbMGOzt\n7U3HE1bqafNvqlSpwpo1a1JVXyN5TJOI5HpMszxkOQ3nNeTdwu+y6v1VFr++sD43o29SfWZ1/vj7\nD15zf41g32BcnV1NxxIGRUVFsXjxYtKlS0eLFi3+8/WYmBjZWyL+48qVK4wdO5axY8dy/fp1AMqU\nKcPAgQNp3LhxovOWnoc8pkkBSql3lVJrlFIjlFKLlFJtTeSQlZG0J5NTJla1WUXR7EU5ePEgDX9u\nyJ37d0zHEgZlzJiRtm3bJlqIAHz99de8/vrrjB49mn/+kb8rRLxH598MGzYsYf5Ns2bNKFWq1L/2\nmNiiNLEyopS6ALTVWm9QSrkBfwPuWusrT3h/sqyMhF8Pp+CYguTNnJczvc5Y/PrCekVcj+DtqW9z\n7uY5GhZtyK8tfsXBTn77Ff/15ptvsn//fgAcHBx49913ad++PfXr18fJyclwOmEt7ty5w5QpUxg+\nfDhnzsT/PPH09KRfv360b9/+uf+syMrII5RSuR6sYFh65vJZINeDf84FxAAp/h/dLaMb5XKXo1zu\ncil9a2FYAdcCCXNsloUso8PyDqnqea+wnF27drFgwQLq16+P1prly5fz3nvvERYWZjqasCLp06fn\nk08+SXT+jaenJ999951Nzb+xmpURpVRT4FvgPlBIa53ozi6l1CvAd0BZQAN/Aj211k+4QU3eAAAS\nCUlEQVRsW6eUKgT8AuwH3gL8tdaLn/J+6TMiksXus7upMbMGt+/fpm+lvgz3GW46krBiFy9eZO7c\nuRw6dIhp06aZjiOs2JPm3/Ts2ZNu3brh6vr0vWqmV0bQWlvFB7ALKARMA2Kf8B5H4CDxhYV68DEd\nOA5keOR93YETD16v9+Cfqz74Wn5gH5D1KVm0EMll9YnV2uFrB40/evj24abjCBu2b98+7ePjo2fP\nnq2joqJMxxFWIC4uTi9btkyXL19eE/8Lu86cObMeOHCg/ueff574fQ9+7pmrAUze/F9BwO7B/z6t\nGOkAxAIFHnnNnfjHLr2f8D1lgNuPvbYb8H1Klif+HyaEJcw9NFcrf6XxR0/dN9V0HGGjevTokfAD\nJ1OmTPqjjz7SW7Zs0XFxcaajCcPi4uL0hg0bdLVq1RL+jGTIkEH37NlTnzlz5j/vN12MWM1jmoeU\nUtOAdjqRxzRKqdVAMa11wcdePwTc0lpXSuR7cgCngVe11uFKqfRABNBKa73pCRm0tf13EanPuD3j\n+HTNp9gpOxa1WESjYo1MRxI25urVq/zyyy/MmDGDPXv2JLw+duxYunfvbjCZsCa7du0iMDCQlStX\nAuDo6Iivry/9+vWjUKFCgPnHNLZWjJwDQrTW1R97fSlQXWud6QnXbAJ8DBwDvICdWushT8kgxYhI\nEV9t/opvtn6Dk70Ta9uu5R2Pd0xHEjbq2LFjzJw5kzlz5rBjx45EZ+SItO3x+Td2dna0bt2aAQMG\nULJkSSlGHvWMYiQaWKO1bvTY67OANsTvG4m2QIZkK0b+vvk3f/7zJzldclLavXSy3EPYDq01XVd2\n5ce9P5LZKTNbfLfwes7XTccSNkxrjVL//ZkSFxfH/fv35XiwSHT+DWC0GJFGB0/w6Ehlb29vvL29\nLXLdG3dvEKfjyOKUxSLXE7ZNKcX3db/n6t2rzP9rPr/8+YsUI+KlJFaIANy9e5cMGTKkcBphjYoW\nLUr79u3JkiULO3bsYP/+/cTGxhrNZGsrIy/0mOYFMshjGpGiomOi+fnPn2n/Wvsn/jAR4mXExcW9\ndMtwkTqdP3+e3Llzy2OaRyVhA2tRrbXnY68/cQPrC2aQYkQIIUSaYXoDq62VyYuAAkqp/A9fUEq5\nA8WBhcZSCSGEEOKFWWMx8rTKbDpwGBimlLJXStkBQ4FTwI8pkE0IIYQQFmY1xYhSarhSaj9Q/8Hn\n+x58JGyy1VrfB3yIb3x2BPgLcCF+v8htS+bx9/cnODjYkpcUQgghrEpwcPC/DmyYYnV7RqyB7BkR\nQgiRlsieESGEEEKkaVKMCCGEEMIoKUaEEEIIYZQUI0IIIYQwSoqRJ5DTNEIIIVI7OU1jxeQ0jRBC\niLRETtMIIYQQIk2TYkQIIYQQRkkxIoQQQgijpBgRQgghhFFSjAghhBDCKClGhBBCCGGUFCNCCPF/\n7d19sF1Vecfx7y8vDQQbKiGGRGgMKJiUhBqpMLHFJFDUWl6MpchAsYipOENr1IDWFkgt9QUVhxLH\nlGGATqkWCbQdtUqnheiIjDSEEELeUAKJSYAkjhIEEkie/rHWjZuTfW7uvbn3rnPu+X1m9py711l7\nnXX3k33Pk/2ylpkV5WSkCQ96ZmZmQ50HPWthHvTMzMw6iQc9MzMzs47mZMTMzMyKcjJiZmZmRTkZ\nMTMzs6KcjJiZmVlRTkbMzMysKCcjZmZmVpSTETMzMyvKyYiZmZkV5WTEzMzMinIyYmZmZkU5GTEz\nM7OinIyYmZlZUU5Gmli4cCFLly4t3Q0zM7MBs3TpUhYuXFi6GygiSveh5UgK7xczM+sUkogIlfp8\nnxkxMzOzopyMmJmZWVFORszMzKyoEaU7MBgknQl8GHgCOBr4YkQsL9srMzMzgw64gVXSa4FNwJSI\n2CTpGOABYHJEvNxkG9/AamZmHcM3sFZImiDpe5L29mOzxwLDImITQH49DJjdj59hZmZmfdQyyYik\nucCPSMlD09MSksZJul3SWklrJN0p6fXdNL0GeE7SzLz9W4DDgd/ux+6bmZlZH7VMMgJcAZwB3N+s\ngqSRwP8AI4EpwFTgV8B9kkZX6v2lpMclrQdOBU4HLpZ0LfAHwArguYH6RczMzKznWuaeEUnDImKv\npFuBiyNieE2decBi4NiIeCqXjQc2A5+MiC/38LO2AGdGxKom7/ueETMz6xi+ZySLiJ7cJzIX2NiV\niOTtngFWA+9rtpGkGyQNyz+/E1jVLBGx9uYh/Nub49feHD/rq5ZJRnpoOrChpnwDMK2b7Q4H7pK0\nCLgAuGgA+mYtwH8M25vj194cP+urdhtn5EhgWU35c8BoSaMiYlfjmxHx5wPdMTMzM+ubdjszYmZm\nZkNMy9zA2uUAN7BuBtZFxJyG8v8E5kTEb/ZTH1prp5iZmQ2wkjewtttlmpXACTXlk4FH++tDSgbE\nzMys07TbZZq7gUmS9g1Ylh/tnQIsKdYrMzMz67NWTEa6OytxG+kMyBckDc+P636eNAHe4kHom5mZ\nmfWzlklGJF0n6WHgj/P68rzsu5SUJ7b7Q2APaWyRx4DXkO4XeaFAt8060gDNI2WDwLGzVtQyyUhE\nXBkRb4mIIyNieETMyMsrDfW2RcRFEXFCREyJiPMiYvPBfn4f5ryxfiRpkqSdlSR0uaSH8+uYSr3D\nJC3KcVqV/6hOrWlvhKS/z7FcKemHkt7e5LPnS3pM0gpJyySdM5C/a7vrxTxSxWIl6cJcZ0Vu80N9\n/HWHlF7EbneTY/HNDfUcu0Ei6SRJN0laLemRfEzdIOnIhnrtedxFRMcvpLluHgHuIF0mEumS0Hpg\ndOn+dcICTALu7UG97wI/AEbl9c8AzwITGuotBtYCR+T1S0nzGE1vqPepvP0b8voZwG7gnaX3Sasu\nwAPAccCtwJ5WixXwfuAl4K15fRrwPDCv9L4rvfQidk/0sD3HbvBitxa4Ezgkr08gTQS7tusYy+Vt\nedwV38GtsADzSJd+JlXKxgOvAJ8o3b9OWOhBMkK6RLcXeEelbCSwA7ixUnZ8jucHGrZfBXyrsn54\nPliuaaj3beDR0vukVRdgWH5t+oVWKlak/0hsBG5tqLcI2AaMLL3/Wj12+f0DJiOO3aDHbjUwuaHs\ngzkG783rbXvctcxlmsL6NOeNDbr3kTLyfTM7R7qP6H5eHae5+XVpw/b3Amfq1zM8vxs4tEm9qZKO\n75deDzHRs3mkSsXqbcDRTeodAczuQd+HrB7Grqccu8E1PSIap0PZQkoEXpvX2/a4czKS9HXOG+tf\nR0n6F0k/ztc7/1XSiZX3pwFbouE+IlKcxleunU4j/e9gY029EcDUSr2u8sZ6kP5dWN+UitV00r0Q\ndfWEY9pTh0n6qqQfSVon6T8k/X5DHcduENUcS5DG3doLfD+vt+1x52QkORLYWVO+b86bQe5PJ9oD\nvAx8OSJOAU7O6z+W9NZcp7s4AYyt1Hsh8nnCA9Sjps3nSAfQWKyvSsWqu3rV9qx7zwNLImIm6Ytk\nNbBU0lmVOo5dQUpDW3wQuDkifpqL2/a4czJiLSEifhYRJ0XEirz+PHAZ6YaqzxbtnFmHiYjjIuK+\n/POuiPg06WbJL5XtmVVcTbok87HSHekPTkaS7UDdvDZjSNnjfjMB28CLiJdIg9ydmou6ixOkm7S6\n6o2W1DiAXl09atpsrGe9VypWjunAeRB4o6Su+xMcu0IkXQL8CfCuiHix8lbbHndORpKVwBtqyvt1\nzhtrTtIYSSNr3toDdE2auBKYqMpAeNlk4JmI2F6pNww4pqbeK6RTzl31YP/YTyZd/1yJ9VWpWK0k\nnT4+UD1rIo9TcUjNW3vya/V4dOwGmaQ/I50NmR0RjV/ybXvcORlJPOdNeTfQ8ORSTk6mAQ/lortJ\nj6nNbKjzdl4dp3/Pr7MaPmM2cE/8erTe7wEv1tSbA6yOiPV9+D0sKRWrB4GfNan3c/a/29/2twCY\nX1N+MrC58oXm2A0ySRcBVwCnR8S2XPYeSfNylfY97ko/O90KSw7eCuAbpKx/GOk5/HV40LPBisGt\nwDLgqLw+nJSgvEwa7r+r3n+R7hw/NK//HfAM+w/o8zXSNe6xef0S0v0n0xrqfTJvPzmvnwHsAs4s\nvU9afSENDNjdWBVFYgWcD7wAzMjr00g31n2o9D5rlaW72AHXAE8Cx1XKFlA/LoVjN3gxuzDvm4/n\nn7uWxcDVlXptedwV38GtsgDjgNtJCcga0kh3ry/dr05ZgN8B/hF4mJQYbgL+Gzitod5o4MYcp8eA\ne4ApNe0NJ408uJZ0ivB+YGaTz/6r3NYK0lmYs0rvj1ZegOtynLbnL6jleRnRKrECLiCNqrwit3lp\n6f3WCktPYkcagPBzef8uJz3++UPgXMeuaOx25JjVLdVkpC2PO+UGzMzMzIrwPSNmZmZWlJMRMzMz\nK8rJiJmZmRXlZMTMzMyKcjJiZmZmRTkZMTMzs6KcjJiZmVlRTkbMzMysKCcjZh1A0nhJWyRdU7ov\nB0PSbZJWN5lU0czalJMRs84wijS99xFdBZImSdor6epy3dqfpCcl3dvk7bHAbwGNs5KaWRvzAW3W\nASJio6RxEfFS6b70QNM5KiLiLEm/ERG7B7NDZjawfGbErEPUJCIq0pGD5ETEbOhxMmI2xEk6W9JW\nSbskPZHLrgIeJJ2FWJDf3yLpryvbjZF0vaSNkrZLekrSIknVSz1X5W33SrpF0rmSlkn6ZS67WNIw\nSR+V9H1JmyT9XNIjki5r6OfpkrYCRwMzc3+2SvpWfn+DpF/kdk9r2HZk7su6vM0mSTdLmtDYftd+\nkHSipPskbZP0E0nza/bdJEm350tHWyStkbRY0ox+CY6ZJaWnRfbixcvgLMB9wBOV9UnAXuCqmrqH\nkKaPXw28KZcdn9dXAYc21N+by7+Wtx0NrAEuBg7L738M9s0Ufj7wMrCg5rM3APc2+R0+QJoy/bRK\nmYDvAE8DM3LZ60jT3m8EXlezH54F7gDG5rJP5T6eU6k3gjQN+y3AqFx2AvAT4JbS8fTiZSgtPjNi\nZnWXaxYAJwGXR8TjABGxHvgEMBX4cM0244H5EfFSRLwAzCedfdkDfDsivhIRkdu6A1iSP+dgXQi8\nG/hCRCzP7T8LfJR0luXzNduMBa6NiB15/XrgFeDcSp2pwJuAuyJiV253HXAtsLUf+m1mmZMRM6vz\np8Bu4AcN5cvy67tqtvm/ri9tgIi4JyLW5uTk7Jr6jwPjJI07yL6eT7rc9J1qYUQ8REoa5kpqTLhe\njIhHK3V3A9uAiZU6O0iJ1GcknVKpe1tE/M1B9tnMKvw0jZnVeSPp78PG/b/HeZ7KI8IVzzZrTNIs\n0mWaE0mXcIL0qDHAoQfZ1+Pya93Zii3ADNJZm6cr5dtr6u4G9o1fEhGbJV0OfAl4QNJPSWdzboqI\nDQfZZzOr8JkRM2tmZ0RMrFnGRMTbaurvrWtE0tnA/wK/AE6JiAkRMZH0JV9KbV8bRcQ/kS71XAZs\nBq4E1kg6bwD7ZtZxnIyYWZ31wOGS9jtrIekESdN70dYl+fXjEVF3RuJgPZ5fJ9a8NxHYCTzTl4Yl\nDYuIX0bETRExC/i93F7JRMpsyHEyYta5fpVfRwBIerOkri/ZO/Lr3OoG+d6LJcDpvficXU3KJ3XT\nr64+DZd0o6Sju2n/m6SbcN/T0NeTgQnAkq4bZ3tD0juAldWyfIPsUtIosGbWT5yMmHWOV938kc9S\nbCM9NQJwHjA7/3w98BDwD11jakh6DbCI9Hfj5l587p359YuSRue2/gi4oK5fpEeCj5U0CpgJzANe\nrNRtrP914LvAFTkBQdJRwFdIj/Z+uqF+bwZ7myLpLyQNy+3+LjAL+EYv2jCzA3AyYjbEdQ16BpwK\nHJMH77oovz0PmCbpadJZkI8A5KdiZgP/BizJ268gfZHPiYidue3L8nsBnJ/bftVjvxFxF3ApcDKw\nVdKjwPuBf85VHpS0sLLJ3wJPkhKJ24CPRMQOSRtICUYAd3cNhpbPepwDfBX4eu7PMlJSc2p+zBdJ\n02v2wxxJs2oGW5tCGmflStLYJk9J2gzcDlwHXN7bOJhZc+rD2UszMzOzfuMzI2ZmZlaUkxEzMzMr\nysmImZmZFeVkxMzMzIpyMmJmZmZFORkxMzOzopyMmJmZWVFORszMzKwoJyNmZmZWlJMRMzMzK+r/\nAZb8umDqgJ0UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cfc4cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.figure(figsize=(8,8))\n",
    "pyplot.xlabel(r'iterations', fontsize=18)\n",
    "pyplot.ylabel(r'$L_2$-norm', fontsize=18)\n",
    "pyplot.semilogy(numpy.trim_zeros(l2_diffJ,'b'),\n",
    "                'k-', lw=2, label='Jacobi')\n",
    "pyplot.semilogy(numpy.trim_zeros(l2_diffGS,'b'), \n",
    "                'k--', lw=2, label='Gauss-Seidel')\n",
    "pyplot.semilogy(numpy.trim_zeros(l2_diffSOR,'b'), \n",
    "                'g-', lw=2, label='SOR')\n",
    "pyplot.semilogy(numpy.trim_zeros(l2_diffSORopt,'b'), \n",
    "                'g--', lw=2, label='Optimized SOR')\n",
    "pyplot.legend(fontsize=16);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Jacobi method starts out with very fast convergence, but then it settles into a slower rate. Gauss-Seidel shows a faster rate in the first few thousand iterations, but it seems to be slowing down towards the end. SOR is a lot faster to converge, though, and optimized SOR just plunges down!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.  [Gonsalves, Richard J.  Computational Physics I. State University of New York, Buffalo: (2011): Section 3.1 ](http://www.physics.buffalo.edu/phy410-505/2011/index.html)\n",
    "\n",
    "2. Moin, Parviz, \"Fundamentals of Engineering Numerical Analysis,\" Cambridge University Press, 2nd edition (2010).\n",
    "\n",
    "3. Young, David M. \"A bound for the optimum relaxation factor for the successive overrelaxation method.\" Numerische Mathematik 16.5 (1971): 408-413."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Source+Code+Pro' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 750px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "    margin-top:0.8em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Alegreya Sans' sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 125%;\n",
       "    font-weight: 400;\n",
       "    width:600px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: 400;    \n",
       "    font-size: 45pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 30pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.3em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    margin-top:16px;\n",
       "    font-size: 22pt;\n",
       "    font-weight: 600;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: rgb(102,102,0);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-size: 14pt;\n",
       "    text-align: center;\n",
       "    margin-top: 0em;\n",
       "    margin-bottom: 2em;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Nixie One', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 16pt;\n",
       "    color: rgb(163,0,0);\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.8em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "    font-family: \"Source Code Pro\";\n",
       "    font-size: 90%;\n",
       "}\n",
       "\n",
       ".alert-box {\n",
       "    padding:10px 10px 10px 36px;\n",
       "    margin:5px;\n",
       "}\n",
       "\n",
       ".success {\n",
       "    color:#666600;\n",
       "    background:rgb(240,242,229);\n",
       "}\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"],\n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "css_file = '../../styles/numericalmoocstyle.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  }
 ],
 "metadata": {
  "hide_input": true,
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}