cds-numerical-methods/Week 5/9 Ordinary Differential Equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
     "grade": false,
     "grade_id": "cell-98f724ece1aacb67",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# CDS: Numerical Methods Assignments\n",
    "\n",
    "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
    "\n",
    "- Solutions must be submitted via the Jupyter Hub.\n",
    "\n",
    "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
    "\n",
    "## Submission\n",
    "\n",
    "1. Name all team members in the the cell below\n",
    "2. make sure everything runs as expected\n",
    "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
    "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
    "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
    "6. submit your solutions  **in time (before the deadline)**"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "team_members = \"Koen Vendrig, Kees van Kempen\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "7317367b26066eae942edeea4e7d4280",
     "grade": false,
     "grade_id": "cell-e9dc85dd9ad77a5e",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Ordinary Differential Equations - Initial Value Problems\n",
    "\n",
    "In the following you will implement your own easy-to-extend ordinary differential equation (ODE) solver, which can be used to solve first order ODE systems of the form\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = \\vec{f}(x, \\vec{y}),\n",
    "\\end{align*}\n",
    "\n",
    "for $x = x_0, x_1, \\dots, x_n$ with $x_i = i h$, step size $h$, and an initial condition of the form $\\vec{y}(x=0)$.\n",
    "The solver will be capable of using single-step as well as multi-step approaches.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "910692c2f4b93f251a3a128caf2b0c37",
     "grade": true,
     "grade_id": "cell-acb87410eaef9fe1",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import gridspec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "79e708efbac889db28e537bcf01e26c6",
     "grade": false,
     "grade_id": "cell-a27c54a734bf2232",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 1\n",
    "\n",
    "Implement the Euler method in a general Python function $\\text{integrator(x, y0, f, phi)}$, which takes as input a one-dimensional array $x$ of size $n+1$, the initial condition $y_0$, the *callable* function $f(x, y)$, and the integration scheme $\\Phi(x, y, f, i)$. It should return the approximated function $\\tilde{y}(x)$. \n",
    "\n",
    "The integration scheme $\\text{phi}$ is supposed to be a *callable* function as returned from another Python function $\\text{phi_euler(x, y, f, i)}$, where $i$ is the step number. In this way we will be able to easily extend the ODE solver to different methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "8396264309b023abf89a62200fa55379",
     "grade": true,
     "grade_id": "cell-1f9655333f2dc3cb",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def integrator(x, y0, f, phi):\n",
    "    \"\"\"\n",
    "    Numerically solves the initial value problem (IVP) given by ordinary differential\n",
    "    equation (ODE) f(x, y) = y' with initial value y0 using the integration scheme\n",
    "    provided by phi.\n",
    "\n",
    "    Args:\n",
    "        x:   size n + 1 numerical array of x values\n",
    "        y0:  an initial value to the function f\n",
    "        f:   a callable function with signature (x, y), with x and y the current state\n",
    "             of the system, and f such that f(x, y) = y' in the ODE\n",
    "        phi: a callable function with signature (x, y, f, i), with x and y the current state\n",
    "             of the system, i the step number, and f as above, representing the integration\n",
    "             scheme to use\n",
    "\n",
    "    Returns:\n",
    "        An n + 1 numerical array representing an approximate solution to y in y' = f(x, y)\n",
    "        given initial value y0, steps from x, and using integration scheme phi.\n",
    "    \"\"\"\n",
    "    \n",
    "    eta    = np.zeros(len(x))\n",
    "    eta[0] = y0\n",
    "    \n",
    "    for i in range(1, len(eta)):\n",
    "        h      = x[i] - x[i - 1]\n",
    "        eta[i] = eta[i - 1] + h*phi(x, eta, f, i)\n",
    "    \n",
    "    return eta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "35d6f653b74f2c8e2ff509890754df72",
     "grade": true,
     "grade_id": "cell-04ac3a0537092706",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def phi_euler(x, y, f, i):\n",
    "    \"\"\"\n",
    "    Returns the integrator phi = f(x, y) for the Euler method to solve the ODE\n",
    "    y' = f(x, y).\n",
    "\n",
    "    Args:\n",
    "        x: size n + 1 numerical array of x values\n",
    "        y: size n + 1 numerical array of y values with estimates up to index i - 1\n",
    "        f: a callable function with signature (x, y), with x and y the current state\n",
    "           variables of the system, and f such that f(x, y) = y' in the ODE\n",
    "        i: the current step number to calculate phi for\n",
    "\n",
    "    Returns:\n",
    "        The integrator phi(x, y, f, i) = f(x, y).\n",
    "    \"\"\"\n",
    "    \n",
    "    return f(x[i - 1], y[i - 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "a300f75f62faab4732ebcac7bde5e18e",
     "grade": false,
     "grade_id": "cell-047f2dfab5489a88",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 2 \n",
    "\n",
    "Debug your implementation by applying it to the following ODE\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = y - x^2 + 1.0,\n",
    "\\end{align*}\n",
    "        \n",
    "with initial condition $y(x=0) = 0.5$. To this end, define the right-hand side of the ODE as a Python function $\\text{ODEF(x, y)}$, which in this case returns $f(x,y) = y - x^2 + 1.0$. You can then hand over the function $\\text{ODEF}$ to the argument $\\text{f}$ of your $\\text{integrator}$ function.\n",
    "\n",
    "Plot the solution you found with the Euler method together with the exact solution $y(x) = (x+1)^2 - 0.5 e^x$.\n",
    "\n",
    "Then implement a unit test for your $\\text{integrator}$ function using this system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "692bfda705027a7fa4452a7e03c8fa7d",
     "grade": true,
     "grade_id": "cell-8b8cb282dfe95a03",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAHmCAYAAAAhuHt9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABu9ElEQVR4nO3deXxU5dn/8c81kw3CJmEnIFARRQiLYVFcQQWrgm2tu6K02s2ltbXFblrb/h5bfarVupTHDeuGu7jUXYpaF0BRdkFkCWsACQTIOtfvj5nEISQhQDJnknzfr9e8Zs59lvs7mRCu3LnPOebuiIiIiIhIcggFHUBERERERL6mAl1EREREJImoQBcRERERSSIq0EVEREREkogKdBERERGRJKICXUREREQkiahAFxERERFJIirQRUSkWTKz4Wb2vpnNNLPHzCw16EwiIqACXUREmq/VwGh3Pw5YAUwINo6ISJQKdBGRfWRmK8zspAPYf4GZnVB/iYJlZqebWWbQOfaVu69z912xxRIgEmQeEZEKKtBFpNkys2PM7L9mVmBmW8zsPTMbVs997FHMu/sR7j6jPvvZ1wz17DYgpQGP36DM7GDgFOCFejjWFWY228yKzezBOmzf3syeNbMdZrbSzM4/0Awi0vg12h+oIiIHwszaAC8CPwKeANKAY4HiIHM1NmZ2GPCOuxcEnWV/xL4P/gVc4u6le9n2BgB3v6GWzdYCfwLGAi3qEOFOoqP3nYHBwEtm9qm7L6jDviLSRGkEXUSaq0MB3P0xdy93913u/pq7fwZgZoeb2Qwz2xqbkjK+pgOZmZvZIXHLD5rZn8zsX0BP4AUzKzSzX8bWV45o762f2La/MLPPYiP908wsI7buV2a2xsy2m9kSMxtTTbaaMtTp/ZlZKzMrN7OucW0DzGydmbUGRgP/rMsXvJav3976OJBj/9XMnotbvtnM3jSzNDNLAR4H/uDuSw6knwru/oy7PwdsrkO2TOA7wO/cvdDd3wWmAxfVRxYRabxUoItIc/U5UG5mU83sVDM7qGJF7GoeLwCvAZ2AK4FHzKzfvnTg7hcBq4Az3L2Vu/81fv0+9HM2MA7oDeQAl8S2uQIY5u6tiY7YrqhLhn15f+5eCCwGhsY13wT8P3ffDrR29w+q7mdmL8aK/+oeL+5jHwfiL8CJZjbEzH5I9Ov4bXcvAc4DRgC/i/2ycs4B9rWvDgXK3P3zuLZPgSMSnENEkowKdBFpltx9G3AM4MD/AflmNt3MOgMjgVbATe5e4u5vEZ0Oc149x6hrP7e7+1p330K0sB4MlAPpQH8zS3X3Fe7+RT33W2EWseLZzI4D+hMbNXf3v1S3g7uf7u7tanicvi99HAh33wzcCkwFrgO+WTEdx93/5e5Z7n5C7DHtQPvbR62AbVXaCoAD+quBiDR+KtBFpNly90Xufom7ZwMDgG5ET3jsBqx29/ireqwEutdzhLr2sz7u9U6glbsvA34K3ABsNLPHzaxbPfdbobJ4Bv5KdEpGSR37qqs69REb6fYaHu/WcOxPgIHAde6+el9Cxf8lAJgMTK7pLwH7oRBoU6WtDXCgfzUQkUZOBbqICODui4EHiRbqa4EeZhb/M7InsKaG3XcCLeOWu8QfupZu97WfqpkfdfdjgINj/VQ7ml1Nhn3tdxYw1My+A2QAj+4tm5n9OzbnvbrHv/e3j9hIt9XwOKaaHAOBu4mOoE/aW+5q+qv8SwDRaTc37eUvAfvicyDFzPrGtQ0CdIKoSDOnAl1EmiUzO8zMfm5m2bHlHkSneHwAfEi06P6lmaVa9JrlZxA9obA6c4HzzSxsZuOA4+PWbQD61LDfvvYTn7+fmY02s3SgCNhFzdfxrpphX/v9lOgvHf9LdBS6tl86AHD3U2Nz3qt7nFoffeyNmXUnOiXoh8CPgYHWwNefN7OU2Em8YSBsZhmxk1H34O47gGeAG80s08xGEb1Z0r8aMqOIJD8V6CLSXG0neoLgh2a2g2hhPh/4eWxqxRnAqcAm4C7g4tgoe3Wujm2/FbgAeC5u3f8Av41NifhF/E770U+8dKIjupuIToHpRHSOdXV2y7Cv/bp7MTAPWOHu1Y1+H7D67sOil098Gfibu093953AzcCfD/TYe/Fbor8sTQYujL3+bVyuf5vZr+O2/zHRyzFuBB4DfqRLLIqI1cMghYiINGFmlgYsA86u7ootjaUPEZHGQiPoIiKyN9cD7zVw4ZyIPkREGgUV6CIiUi0zG2pmBcBxRK+V3ij7EBFpbDTFRUREREQkiWgEXUREREQkiahAFxERERFJIirQRURERESSSLU3T2jOOnTo4L169Qo6hoiIiIg0cXPmzNnk7h2rtqtAr6JXr17Mnj076BgiIiIi0sSZ2crq2jXFRUREREQkiahAFxERERFJIirQRURERESSSJMv0M1snJktMbNlZjY56DwiIiIiIrVp0gW6mYWBO4FTgf7AeWbWP9hUIiIiIiI1a+pXcRkOLHP35QBm9jgwAVgYaKoqIuXlRCLlQcdIKLMm/bvhbsws6AgJ05zeK4CFms/3sYiIJE5TL9C7A6vjlvOAEQFlqdGcF//JsE+uCzqGiOyHcjcihHDACRHB8NgjguFW8ToEFW27bRNrt/j2UOV+FctfbxP9pcAttMf2YJXtVPYdwi1MxMLRdZYSfQ6lVC4TW/76OYyHwmBhLJQSex2CUAoWCkMoBULRdfHPoXDFciqhcAhLSSeckoalpJKSkk4oNZ1QahopqemkpKSRkpYefZ2aTkpaGqlp6aSmphMKhwP7PEVEkkFTL9DrxMwuBy4H6NmzZ8L7z/rGkby/5YcJ71cSwD3oBAnUnN4r4BHMHcfBI7EHmJdHX8e2waMleeU2lftGAAf3yte7tVHDNsRv70A0h8Wev96mHNwJeST6a4CXRx9ECFFOyL9+DhNtr3hO8YrlCKmW+L/ulXqYMsKUkkKZpVAW91xuqZRbCmWWSnkojbJQOuWhdMrDGXg4jUhKBh5Oh5QWkJKBpWZAagtCqRmE0loQSmtBOK0F4bSWpKZnkJLektT0lqS1yKRFZltaZLYmLS0t4e9ZRCReUy/Q1wA94pazY227cfcpwBSA3NzchFcZfQaMoM+ApBvYFxEBwCPllJeVUVZWSiRSRllZOV5eSllZGV5eRll5GZGyMjxSRnl5tK28vIxIaQllpSVEyoopLyvGy0opL4sue1kJXlaKl5dEX5eXQlkJREqhvATKS7FICZSXEYpE20ORUkKRsthzCWEvoUVZAamRYtK8mFQvIY0S0ryUDEoI2f79OC/2VHZaBkVkUBRqQUnsURZuQVlKJpHUlkRSMiEtE0/LJJTeCkvPJCWjNWmZB5Heqh0t2mTRsk17Wrdtr4JfRPZZUy/QZwF9zaw30cL8XOD8YCOJiDQuFgqTkhYmJS096Ch15pEIpaXFFO3aSUnxDkqKdlFatJPS4h2UFRdRVryD8pIiykt2EikpIlKyk0jxDrxkB5TsIFQafYTLdpJSvovU8p1klhSQXrSLDC8iw4vItKI6ZdnhGRRaS3aGWlEUyqQ4pTWlqa0pT2tNJL0tpLfBWh5EaqsOpLftRMt2nWjdvjNt23ckvRF9zUWk/jTpAt3dy8zsCuBVIAzc7+4LAo4lIiINzEIhUtNbkJreAshqkD48Uk7xrh3sLCygaMd2incWUFz4FcWFX1G24yvKdxXguwqgeBuh4m2klGwjtayQVmVbyCheTabvoJXvqHUaUYFnsi3UhsJwW4pS21GSdhBlGe2hRXusVUfS23UhM6s7bTpmk9WpG+karRdpEpp0gQ7g7i8DLwedQ0REmhYLhcnIbENGZpv9P4g7RTu3U7h1E4VfbWTn1g2UbMuntHATkR2bsZ1bSCneQlrJVtqWbCSzaCnttm4j3Ur3OFS5G/nWlq2h9uxIzaIooyPlLTtC6y6ktu1CRlY2B3XpQ8euPclIVyEvksyafIEuIiKStMwqi/wO3fvUbR93indtZ9umdWzLX8OOr9ZQ8tU6ItvWE9qxkbSifFqVbKL7ti84aOtWUiyy2+6lHmatteerlE4UZnSltFVXrG02aVkH07pzL7K69SErqxOhUPO6bKpIMlGBLiIi0piYkd6yDR17tqFjz361buqRcrZ/tYGvNqymMH8VRZtWUb51NSnb19Bi1zp67viMDtvfJnX97tNsCjyT9eGuFLTIprh1L0JZvWnZpS9ZPQ6ja3YvUlN0KUyRhqQCXUREpImyUJjWWd1ondWNmm4D4uVlFGxay5a1y9m+cQUlm1diX60go3AV2bsW06lwJinrIxA7g2uXp7Eq1IWv0ruzs00frNPhtMk+gm6H5NAhK6vZ3bBMpCGoQBcREWnGLJxC2849adu5+vuAeFkJm9d8waa8xexYv4zIpi9I376KzrtW0WXjR6RuLIf50W3XkcX6tF4Utu4DHQ8jM7s/XQ4ZStfOnVW4i+wD82Z1I5W9y83N9dmzZwcdQ0REJOl5WQmb8z4nf/mn7FqzgNCWpbTd/gVdy1aTQUnldnl0Ym1GX3a2P5y07oPo1G84B/fqq6ky0uyZ2Rx3z63arhF0ERER2S+WkkaHXgPo0GvA7isiEQrWfcGG5Z+yY9WnhDfOp3vhYrqu+S+htQ6z4CtvxarUPmxteziWnUvn/sfQ5xv9VLSLoBH0PWgEXUREpGGU7Sxg/dKP2fzFHCLrPqPt1kV0L/2SdKKXjdzo7ViRcTiFHQfTsvcIDh54DF06dtD0GGmyahpBV4FehQp0ERGRxPGyYjYsncPGRe9B3mw6FHxGt/K1QPTa7sutJ2vaDoVeo8gefBJ9Du6lS0BKk6ECvY5UoIuIiASrZNsm8ua/w7ZlH5CxfjYH75xHC4oBWE53VrcZSqTn0XTNGUPfQw4lrIJdGikV6HWkAl1ERCS5eFkJG5d8yIZ5b5Ka9196Fn5GJrsA+ILurGx3FOG+Yzh0xFi6dsgKOK1I3alAryMV6CIiIkmuvIz8ZbNZ/9kbpK38DwcXziWDEoo9lfkpR7Cl67G0yzmVgYNHkJGm62FI8lKBXkcq0EVERBoXL9lJ3qdvsvWzf3PQunfJLlsJwBrvwOK2xxLufzoDjxpHVttWAScV2Z0K9DpSgS4iItK4FW1aycqPXsCX/JveBR+RTglbPZPPWoygpO+pHHr0mfTs2inomCIq0OtKBbqIiEjT4cWFrJr1EoWfPk+PTTNp49sp9lRmpQ2jsO8EBpzwXbI7ad66BEMFeh2pQBcREWmiysvYuGAG+R89Sbe1r3JQ5Cu2ewvmtDiK0sO+xeATvk3HdpoGI4mjAr2OVKCLiIg0A5FyNs57g03vP0KPDW/Q2newxVsxp9WJpORezFGjRusEU2lwKtDrSAW6iIhIM1NWzNrZL7Jt9uP03jSDdEpYwsEs7X4mfU6cxOHfOFh3M5UGoQK9jlSgi4iINF+RnV+xYsZDpH72CD2KllDsKXyQNpJdAy7gqJPOom1mWtARpQlRgV5HKtBFREQEYPvKT1jz1hS6rXqBNr6d5d6Ned3P4fBxP+DQnl2DjidNgAr0OlKBLiIiIrspLSLvvUeJfPhPeu5azHZvwbutTiHzmB9x9PARpIRDQSeURkoFeh2pQBcREZGabFv6Putf/zu9N75OKmX8N3QkW4b8hBNPHk9mRmrQ8aSRUYFeRyrQRUREZG/KCtbx5at30nnRQ7TxAubSj5WHX86xp11I+1YZQceTRkIFeh2pQBcREZE6K9nJqjen0HLOXXQo28BSz2Zer0sZMf5yume1CTqdJLmaCvSkmzRlZjeb2WIz+8zMnjWzdnHrrjOzZWa2xMzGxrWPi7UtM7PJce29zezDWPs0M9Op1yIiIlJ/0lrS89Sf0uG6Bawf83daZ6Ty7ZV/pOz2XJ66/39Z/9WOoBNKI5R0BTrwOjDA3XOAz4HrAMysP3AucAQwDrjLzMJmFgbuBE4F+gPnxbYF+Atwq7sfAnwFfC+h70RERESah3AqXY69hC6/+pjNZzxIWovWnLXqRgpvG8ZTD93Bpu27gk4ojUjSFeju/pq7l8UWPwCyY68nAI+7e7G7fwksA4bHHsvcfbm7lwCPAxMsekeB0cBTsf2nAmcm6G2IiIhIcxQKkXXkt+j6y1lsOnUKrTNSOGv5b8m/ZQRPPnovBTtKgk4ojUDSFehVTAL+HXvdHVgdty4v1lZTexawNa7Yr2gXERERaVihEB1GnEPnX33CxpNup0NaKd/9/OcsufkEnv/3vyktjwSdUJJYIAW6mb1hZvOreUyI2+Y3QBnwSALyXG5ms81sdn5+fkN3JyIiIs1FKEynYybScfJnrBv1Jw4PreaMD87jjf85i//M/hRdrEOqkxJEp+5+Um3rzewS4HRgjH/9nbsG6BG3WXasjRraNwPtzCwlNooev33VPFOAKRC9iss+vRkRERGRvQmn0vXkK/FjLmT1839izOIHKXvhZJ6ccTYDzv4d/Xt2DjqhJJGkm+JiZuOAXwLj3X1n3KrpwLlmlm5mvYG+wEfALKBv7IotaURPJJ0eK+zfBs6K7T8ReD5R70NERESkKmtxED3P/V/sills7HI8Zxf+i1b3jeKRf/0f24tKg44nSSLpCnTgH0Br4HUzm2tm9wC4+wLgCWAh8ArwE3cvj42OXwG8CiwCnohtC/Ar4BozW0Z0Tvp9iX0rIiIiIntK7dCHXj96ksJznyMjI5MLvvgFs/5yGq+/P0fTXkQ3KqpKNyoSERGRhCorYd0rN9N+9t8pc3iu3USOvuC39O7UNuhk0sAazY2KRERERJqVlDS6nv4bUq78kC2dRnBBwRR23HkCT738GuURDaQ2RyrQRURERJJAOKs3PX48nYLT76VneAvjPzyXx//3Kpat+yroaJJgKtBFREREkoUZbXO/S+trZpOffTIX7HiIontO5MmXXtFoejOiAl1EREQkyVirjnS/bBoFZ9xLz/BXTPjofB7+289Zvbkw6GiSACrQRURERJJU2yO/S+ufz2FztxOYWHgfq27/Jq9++GnQsaSBqUAXERERSWKW2YGulz/NlhP/Qq4t5siXT+fe++6hsLgs6GjSQFSgi4iIiCQ7M9of/0PCP5yBZ3bm+6t/xSs3T2T+yvygk0kDUIEuIiIi0kikdOlPx5+9y/rDL+Wsshcpve9Unv/PR0HHknqmAl1ERESkMUnNoMs5t7F9/H0cFs5j1Fvf4Z8PPkBRaXnQyaSeqEAXERERaYRaDz2LtB/9B2uZxfe//BmP/u0aVm3aEXQsqQcq0EVEREQaqXCnfmT99F029RzHpF0Psvgf3+a/i1YHHUsOkAp0ERERkcYsvRWdJz3GV6Ou5yQ+pPVjZ/D0DM1Lb8xUoIuIiIg0dmYcdPI1FJ31MIeE1zPq7bP557RnKCuPBJ1M9oMKdBEREZEmouWA00m7/HUy0tO5eOEPuPueW9lWVBp0LNlHKtBFREREmpBw14G0u+oddhx0GFfm38ijt01mfUFR0LFkH6hAFxEREWlqWnWiw09eZ1OPsfyw6F5ev/2HLNuwPehUUkcq0EVERESaotQMOlz6GJsPv5CLyp9l4d0X8MmKjUGnkjpQgS4iIiLSVIXCZJ39D7aO+AXj+Q8FD5zDfxasDDqV7IUKdBEREZGmzIx2p/6O7SfdzLE2l8xp3+Xfsz8POpXUQgW6iIiISDPQ+pjLKTnzXgaHltFl+nm8+NHCoCNJDVSgi4iIiDQTLQZ/h/KzHmRgaAUHv3gez703L+hIUg0V6CIiIiLNSPqA8UTOeYR+oTX0e/V8npz5SdCRpAoV6CIiIiLNTNrh4+C8x+kT3sCgNy5g2ttzgo4kcZK2QDezn5uZm1mH2LKZ2e1mtszMPjOzoXHbTjSzpbHHxLj2I81sXmyf283MgngvIiIiIskmrd9JhC54il7hTQx8+1Ke1nSXpJGUBbqZ9QBOAVbFNZ8K9I09Lgfujm3bHrgeGAEMB643s4Ni+9wNXBa337hE5BcRERFpDFIPOQ7OfYS+obX0efUSXpqlq7skg6Qs0IFbgV8CHtc2AXjIoz4A2plZV2As8Lq7b3H3r4DXgXGxdW3c/QN3d+Ah4MyEvgsRERGRJJfW72T8rAfICS0n64WLef3TL4OO1OwlXYFuZhOANe7+aZVV3YHVcct5sbba2vOqaRcRERGROGkDzqB0/D0MDy0m/emJzFyUt/edpMEEUqCb2RtmNr+axwTg18DvE5zncjObbWaz8/PzE9m1iIiISFLIGHoOReNu5bjQp2x/7Pt8vHJz0JGarUAKdHc/yd0HVH0Ay4HewKdmtgLIBj42sy7AGqBH3GGyY221tWdX015dninunuvuuR07dqyfNykiIiLSyLQceSk7jvsdp4XeZ/4DV7M8vzDoSM1SUk1xcfd57t7J3Xu5ey+i01KGuvt6YDpwcexqLiOBAndfB7wKnGJmB8VODj0FeDW2bpuZjYxdveVi4PlA3piIiIhII5F54s/ZljOJi3mBl6b8lo3bi4KO1OwkVYG+Fy8THWFfBvwf8GMAd98C/BGYFXvcGGsjts29sX2+AP6d4MwiIiIijYsZbc68ha29TuXK0gd44J//S2FxWdCpmhWLXuBEKuTm5vrs2bODjiEiIiISrNIitk45jRYb53JL55v45Q++R2q4MY3tJj8zm+PuuVXb9VUWERERkT2lZtBu0lMUte7JTzZcz9+ffBUN7CaGCnQRERERqV6Lg2g76RnSUlOYsPDnPP7O/KATNQsq0EVERESkZu17k3H+I/QOrafbGz/hnSXrg07U5KlAFxEREZFahfocS9mpt3B86FO+fOwaXX6xgalAFxEREZG9yhgxie2DvsfFvMTT9/4PBTtLg47UZKlAFxEREZE6aT3+rxR0O5ari+7h7w89RiSik0Ybggp0EREREambcAptL3yIkpad+N66G7jvNV2auiGoQBcRERGRumvZnswLH6VTaDv9/vszZi7WSaP1TQW6iIiIiOwT6z6EyKl/5bjQPBY//mvWbN0VdKQmRQW6iIiIiOyz9GGXsP3wc7icp7nv/nsoLisPOlKToQJdRERERPadGa2//Xe2tT2Mqwpu5h/Pzgg6UZNRa4FuZhlmdpaZ/d3MnjSzh8zsl2Z2RKICioiIiEiSSm1Bm4sfpUXYOWber3l9/tqgEzUJNRboZvYH4D3gKOBD4J/AE0AZcJOZvW5mOQlJKSIiIiLJKesbhE7/X0aEFrP06RtYX1AUdKJGL6WWdR+5+/U1rPubmXUCejZAJhERERFpRFKHns/2xa/zgyVPcONDI/j9TyYRDlnQsRqtGkfQ3f0liE5zqbrOzDq4+0Z318UvRURERITW3/47uzK7c9mm/8f9b84NOk6jVpeTRGeZ2ciKBTP7DvDfhoskIiIiIo1ORhsyz59KV/uKbjMn8/HKLUEnarTqUqCfD9xhZjeb2SPAZcDoho0lIiIiIo2NZedSevyvOS38AS89ejs7S8qCjtQo7bVAd/d5wJ+BHwInAle4e15DBxMRERGRxifj+J+xveORXFX0T+6a/k7QcRqlvRboZnYf8FMgB7gUeNHMftLAuURERESkMQqFaX3u/9EyVM7QT2/gvaX5QSdqdOoyxWUecKK7f+nurwIjgKENG0tEREREGq2sb8DJf2B0eC4zp93K9qLSoBM1KnWZ4nKbu3vccoG7f69hY4mIiIhIY5Y68gds73oUV5Tezz+eeTvoOI1KbTcqesHMzjCz1GrW9TGzG81sUsPGExEREZFGKRSi9dn/JC1sHLfoBt5etD7oRI1GbSPolwHHAovMbJaZvWxmb5nZcqJ3FZ3j7vcnJKWIiIiIND4HHUxo3J8ZFV7A+0/dpqkudVTbjYrWu/svgduBy4E/AtcAA9z9ZHd/vqFCmdmVZrbYzBaY2V/j2q8zs2VmtsTMxsa1j4u1LTOzyXHtvc3sw1j7NDNLa6jMIiIiIrKn1GGXsr3LSH5SNpV7XtStdOqiLieJdgKeBH4GdAF2NWQgMzsRmAAMcvcjgFti7f2Bc4EjgHHAXWYWNrMwcCdwKtAfOC+2LcBfgFvd/RDgK0Bz50VEREQSyYzWZ91JZqiMwz79f3y86qugEyW9upwk+lugL3AfcAmw1Mz+n5l9o4Ey/Qi4yd2LY/1vjLVPAB5392J3/xJYBgyPPZa5+3J3LwEeByaYmRG9odJTsf2nAmc2UGYRERERqUmHQyg/9uecEf6AZx+/n9LySNCJklpdRtCJXcVlfexRBhwEPBU//aQeHQocG5ua8h8zGxZr7w6sjtsuL9ZWU3sWsNXdy6q0i4iIiEiCpR93DYVtDuEHO+7igbfnBx0nqdXlRkVXm9kc4K/Ae8BAd/8RcCTwnf3p1MzeMLP51TwmAClAe2AkcC3wRGw0vMGY2eVmNtvMZufn62L6IiIiIvUuJY1WZ91Ftm0ideb/8OWmHUEnSlp1GUFvD3zb3ce6+5PuXgrg7hHg9P3p1N1PcvcB1TyeJzrS/YxHfQREgA7AGqBH3GGyY201tW8G2plZSpX26vJMcfdcd8/t2LHj/rwlEREREdmbniPYOegSLg79m3ufeIa4W+1InLrMQb/e3VfWsG5R/UfiOeBEADM7FEgDNgHTgXPNLN3MehOdF/8RMAvoG7tiSxrRE0mnx6blvA2cFTvuRKDBrjwjIiIiInvXctwfKEk7iG+v/zsvf7Y26DhJqU5z0BPsfqCPmc0nesLnxNho+gLgCWAh8ArwE3cvj80xvwJ4FVgEPBHbFuBXwDVmtozonPT7EvxeRERERCRei3akj/sjR4aWMueFu9lZUrb3fZoZ058Wdpebm+uzZ88OOoaIiIhI0xWJUHj3aHZtXM6jw5/h6tNyg04UCDOb4+57vPlkHEEXERERkaYsFKLVt24ly7bR5oP/1QmjVahAFxEREZHE6zaE4oEXclHoFf7v6Zd0wmgcFegiIiIiEogW4/5AWWorTs+7lTcXbgg6TtJQgS4iIiIiwcjMIvXk33F0eCEzXphKSZnuMAoq0EVEREQkQOHcSexo8w0u3fkAj76/LOg4SUEFuoiIiIgEJ5xCy9P+H98IrWP9W3dTsLM06ESBU4EuIiIiIoGyQ8eyo9soLo88wf+9/nHQcQKnAl1EREREgmVG5hk30c520Hb2HazavDPoRIFSgS4iIiIiweuaQ3H/s5kY+jf3vvh20GkCpQJdRERERJJCi3E3YOEww5bdzpyVXwUdJzAq0EVEREQkObTpho/8CWeEP2DaCy8225sXqUAXERERkaSRduzVFKe0YdyGe3ln6aag4wRCBbqIiIiIJI8W7Qgf9zNGh+fy0ovPNMtRdBXoIiIiIpJUUkb+gF3pHfj21vt5df66oOMknAp0EREREUkuaZmkj/4VI0KLmfHvaZRHmtcougp0EREREUk6oSMvYWfL7pxfOJXnPs4LOk5CqUAXERERkeSTkkbGSb8mJ/Qlc157iJKySNCJEkYFuoiIiIgkpdCgc9nRpg8Tix7liVkrg46TMCrQRURERCQ5hVNoedJ19AvlseitR5vNKLoKdBERERFJWjbgO+xs3ZsLiqfx7Merg46TECrQRURERCR5hcK0GH0t/UMr+eSNxygtb/qj6CrQRURERCSpWc7Z7MzswXlF03j+kzVBx2lwKtBFREREJLmFU2kx+loGhZYz640nmvx10ZOuQDezwWb2gZnNNbPZZjY81m5mdruZLTOzz8xsaNw+E81saewxMa79SDObF9vndjOzIN6TiIiIiBwYG3Qeu1p245ydj/Lip017FD3pCnTgr8Af3H0w8PvYMsCpQN/Y43LgbgAzaw9cD4wAhgPXm9lBsX3uBi6L229cYt6CiIiIiNSrlDTST/gFQ0PLeO/1p5v0KHoyFugOtIm9bgusjb2eADzkUR8A7cysKzAWeN3dt7j7V8DrwLjYujbu/oG7O/AQcGYi34iIiIiI1J/Q0AvZldGJCdun8dqC9UHHaTDJWKD/FLjZzFYDtwDXxdq7A/HX1smLtdXWnldNu4iIiIg0RinppB9zJaPCC3jjzVeIjsE2PYEU6Gb2hpnNr+YxAfgR8DN37wH8DLgvAXkuj813n52fn9/Q3YmIiIjIfgrlXkJJSitO3Pw4H325Jeg4DSKQAt3dT3L3AdU8ngcmAs/ENn2S6LxygDVAj7jDZMfaamvPrqa9ujxT3D3X3XM7dux4oG9PRERERBpKRhtCw77HqeGPeObNd4NO0yCScYrLWuD42OvRwNLY6+nAxbGruYwECtx9HfAqcIqZHRQ7OfQU4NXYum1mNjJ29ZaLgecT+k5EREREpN6lHPUjsDBHrPwXS9ZvDzpOvUvGAv0y4H/N7FPg/xG9YgvAy8ByYBnwf8CPAdx9C/BHYFbscWOsjdg298b2+QL4d4Leg4iIiIg0lDZdKRtwNt8N/4dH3vo46DT1LiXoAFW5+7vAkdW0O/CTGva5H7i/mvbZwID6zigiIiIiwUo/7qcw7xE6LJrKuoLhdG3bIuhI9SYZR9BFRERERGrX8VB29hnLhaFX+dd/Fgadpl6pQBcRERGRRqnlCT+nvRVSMvththWVBh2n3qhAFxEREZHGqecIdnYawnn8myc+Whl0mnqjAl1EREREGq2Wx/yEb4TWseS95yiPNI0bF6lAFxEREZHGq/8EijI6ctrO6by1eGPQaeqFCnQRERERabxS0kgdeTknhD/l1Rn/CTpNvVCBLiIiIiKNWjj3UsotlYFrn2Dx+m1BxzlgKtBFREREpHFr1ZGyI77DWeGZTHtnXtBpDpgKdBERERFp9NJH/ZhMKybts8f4akdJ0HEOiAp0EREREWn8ug5iZ9cRXGCv8PhHXwad5oCoQBcRERGRJqHlMT+mZyif5e89Q1l5JOg4+00FuoiIiIg0DYedRlFGJ04tfoUZS/KDTrPfUoIO0BiUlpaSl5dHUVFR0FGkkcrIyCA7O5vU1NSgo4iIiDRd4VRSh13CCe/czC/f+5CT+o8POtF+UYFeB3l5ebRu3ZpevXphZkHHkUbG3dm8eTN5eXn07t076DgiIiJNWjh3IpF3bqH3yqdYveUkerRvGXSkfaYpLnVQVFREVlaWinPZL2ZGVlaW/gIjIiKSCG2zKe5zMmeHZ/DEh18EnWa/qECvIxXnciD0/SMiIpI4LUZ+n45WwMZZz1JS1vhOFlWB3kiEw2EGDx5c+bjppptq3f7BBx/kiiuuSFC6Pc2YMYP//ve/lcuXXHIJTz311H4f70D3h7p9TVasWMGjjz5auTx79myuuuqqA+pXREREEuyQMezK7M740ld4feGGoNPsM81BbyRatGjB3LlzG+z4ZWVlpKTU37fDjBkzaNWqFUcffXS9HTMRKgr0888/H4Dc3Fxyc3MDTiUiIiL7JBQmffgkRr39R6559z1Oyzkr6ET7RCPojVyvXr3YtGkTEB3tPeGEE/bYJj8/n+985zsMGzaMYcOG8d577wFwww03cNFFFzFq1Cguuuii3faZMWMGxx9/PBMmTKBPnz5MnjyZRx55hOHDhzNw4EC++OKLGo+9YsUK7rnnHm699VYGDx7MO++8A8DMmTM5+uij6dOnT+VouLtz7bXXMmDAAAYOHMi0adMq26+44gr69evHSSedxMaNG6t9/7fffjv9+/cnJyeHc889F4AtW7Zw5plnkpOTw8iRI/nss8/22K/qiHyrVq0AmDx5Mu+88w6DBw/m1ltvZcaMGZx++um1HveGG25g0qRJnHDCCfTp04fbb799bx+biIiINLDQ0IsotzCHrX2GL/ILg46zTzSCvo/+8MICFq7dVq/H7N+tDdefcUSt2+zatYvBgwdXLl933XWcc845dTr+1Vdfzc9+9jOOOeYYVq1axdixY1m0aBEACxcu5N1336VFixZ77Pfpp5+yaNEi2rdvT58+ffj+97/PRx99xN///nfuuOMObrvtthqP/cMf/pBWrVrxi1/8AoD77ruPdevW8e6777J48WLGjx/PWWedxTPPPMPcuXP59NNP2bRpE8OGDeO4447j/fffZ8mSJSxcuJANGzbQv39/Jk2atEfGm266iS+//JL09HS2bt0KwPXXX8+QIUN47rnneOutt7j44ovr/NeHm266iVtuuYUXX3wRiP6iUqG24y5evJi3336b7du3069fP370ox/pkooiIiJBat2Z0r6n8d0lb3H3+0v59fghQSeqMxXojcSBTHF54403WLhwYeXytm3bKCyM/iY5fvz4aotzgGHDhtG1a1cAvvGNb3DKKacAMHDgQN5+++29HruqM888k1AoRP/+/dmwITof7N133+W8884jHA7TuXNnjj/+eGbNmsXMmTMr27t168bo0aOrPWZOTg4XXHABZ555JmeeeWblMZ9++mkARo8ezebNm9m27cB/qartuKeddhrp6emkp6fTqVMnNmzYQHZ29gH3KSIiIvsvY8QkMj6fzrZPnqPkm4NIS2kck0dUoO+jvY10J1pKSgqRSPTs5Jou4xeJRPjggw/IyMjYY11mZmaNx05PT698HQqFKpdDoRBlZWV7PXZtx3P3vW5fFy+99BIzZ87khRde4M9//jPz5s2r037xX7dIJEJJSckB5Yh/b+FwuPLrIyIiIgHqfTxFmd355rY3eXPRDzl1YNegE9VJ4/g1QmrUq1cv5syZA1A5ulvVKaecwh133FG5XJ8nm9Z07NatW7N9+/a97n/ssccybdo0ysvLyc/PZ+bMmQwfPpzjjjuusn3dunWVI/bxIpEIq1ev5sQTT+Qvf/kLBQUFFBYWcuyxx/LII48A0SkqHTp0oE2bNrvtG/91mz59OqWlpXvNXZfjioiISBIJhUg78kKOCc/njQ/mBJ2mzgIp0M3su2a2wMwiZpZbZd11ZrbMzJaY2di49nGxtmVmNjmuvbeZfRhrn2ZmabH29Njystj6Xgl7gw2gYg56xWPy5OiX4Prrr+fqq68mNzeXcDhc7b633347s2fPJicnh/79+3PPPffUW66ajn3GGWfw7LPP7naSaHW+9a1vkZOTw6BBgxg9ejR//etf6dKlC9/61rfo27cv/fv35+KLL+aoo47aY9/y8nIuvPBCBg4cyJAhQ7jqqqto164dN9xwA3PmzCEnJ4fJkyczderUPfa97LLL+M9//sOgQYN4//33K/+SkJOTQzgcZtCgQdx666277VOX44qIiEhyCQ05nxBOt5XPsWFb47hpoNXXVIN96tTscCAC/BP4hbvPjrX3Bx4DhgPdgDeAQ2O7fQ6cDOQBs4Dz3H2hmT0BPOPuj5vZPcCn7n63mf0YyHH3H5rZucC33H2vZ1Xm5ub67Nmzd2tbtGgRhx9++IG/cWnW9H0kIiISjF33nsbGVZ/z0okv8eMTD937DgliZnPcfY/rOQcygu7ui9x9STWrJgCPu3uxu38JLCNarA8Hlrn7cncvAR4HJlj09oyjgYrr5U0Fzow7VsUQ51PAGNPtHEVERESanRbDLubg0EY+/+i1ejsPriEl2xz07sDquOW8WFtN7VnAVncvq9K+27Fi6wti24uIiIhIc3L4GZSmtOLYwleYs/KroNPsVYMV6Gb2hpnNr+YxoaH63F9mdrmZzTaz2fn5+UHHEREREZH6lNYSBnyHU0Mf8fyHi4NOs1cNdplFdz9pP3ZbA/SIW86OtVFD+2agnZmlxEbJ47evOFaemaUAbWPbV5d1CjAFonPQ9yO3iIiIiCSx1CMvInXuVGzBs+woHk5mevJebTzZprhMB86NXYGlN9AX+IjoSaF9Y1dsSQPOBaZ7dBLR28BZsf0nAs/HHWti7PVZwFveGCYdiYiIiEj9y85lV9tDGM8MXp63Lug0tQrqMovfMrM84CjgJTN7FcDdFwBPAAuBV4CfuHt5bHT8CuBVYBHwRGxbgF8B15jZMqJzzO+Ltd8HZMXarwEqL80oIiIiIs2MGRnDLiY39DnvffRh0GlqFdRVXJ5192x3T3f3zu4+Nm7dn939G+7ez93/Hdf+srsfGlv357j25e4+3N0PcffvuntxrL0otnxIbP3yxL7LxuO5555j4cKFB3ycGTNmcPrpp9e6zdatW7nrrrsql9euXctZZ51Vyx4iIiIi9cMGnoVjHLz2ZdZu3RV0nBol2xQXCUB9Feh1UbVA79atG0899VQte4iIiIjUk7bdKc4+mgmh93j+kzV73z4gKtAbiYcffpjhw4czePBgfvCDH1BeXs6sWbPIycmhqKiIHTt2cMQRRzB//nwKCwsZM2YMQ4cOZeDAgTz//POVx3nooYcq79x50UUX8d///pfp06dz7bXXMnjwYL744ovd+n3yyScZMGAAgwYN4rjjjgOgqKiISy+9tPIOnm+//fYeeW+44QZuueWWyuUBAwawYsUKJk+ezBdffMHgwYO59tprWbFiBQMGDKj1uA8++CDf/va3GTduHH379uWXv/xlvX99RUREpHnIGHoefULrWTj77aS9Jnrynr6arP49GdbPq99jdhkIp95U4+pFixYxbdo03nvvPVJTU/nxj3/MI488wsUXX8z48eP57W9/y65du7jwwgsZMGAAZWVlPPvss7Rp04ZNmzYxcuRIxo8fz8KFC/nTn/7Ef//7Xzp06MCWLVto374948eP5/TTT692qsmNN97Iq6++Svfu3dm6dSsAd955J2bGvHnzWLx4Maeccgqff/55nd7qTTfdxPz585k7dy4AK1asqFxX23Hnzp3LJ598Qnp6Ov369ePKK6+kR48e1fQgIiIiUov+4yl/8RqGFrzOonXn0L9bm6AT7UEFeiPw5ptvMmfOHIYNGwbArl276NSpEwC///3vGTZsGBkZGdx+++0AuDu//vWvmTlzJqFQiDVr1rBhwwbeeustvvvd79KhQwcA2rdvv9e+R40axSWXXMLZZ5/Nt7/9bQDeffddrrzySgAOO+wwDj744DoX6LWp7bhjxoyhbdu2APTv35+VK1eqQBcREZF9l9GW8kPGcsaS//B/n6ykf7eBQSfagwr0fVXLSHdDcXcmTpzI//zP/+yxbvPmzRQWFlJaWkpRURGZmZk88sgj5OfnM2fOHFJTU+nVqxdFRUX71fc999zDhx9+yEsvvcSRRx7JnDlz6rRfSkoKkUikcnl/+6+Qnp5e+TocDlNWVlbL1iIiIiI1SxtyLh0+f4H1n7xC+akDCIcs6Ei70Rz0RmDMmDE89dRTbNy4EYAtW7awcuVKAH7wgx/wxz/+kQsuuIBf/epXABQUFNCpUydSU1N5++23K7cdPXo0Tz75JJs3b648DkDr1q3Zvn17tX1/8cUXjBgxghtvvJGOHTuyevVqjj32WB555BEAPv/8c1atWkW/fv12269Xr158/PHHAHz88cd8+eWXe+2rLscVEREROWB9T6YktQ3HF8/gg+XV3scyUCrQG4H+/fvzpz/9iVNOOYWcnBxOPvlk1q1bx0MPPURqairnn38+kydPZtasWbz11ltccMEFzJ49m4EDB/LQQw9x2GGHAXDEEUfwm9/8huOPP55BgwZxzTXXAHDuuedy8803M2TIkD1OEr322msZOHAgAwYM4Oijj2bQoEH8+Mc/JhKJMHDgQM455xwefPDB3Ua4Ab7zne+wZcsWjjjiCP7xj39w6KGHApCVlcWoUaMYMGAA11577W771OW4IiIiIgcsJZ3QgG8xLjSLl+YsCzrNHixZz14NSm5urs+ePXu3tkWLFnH44YcHlEiaCn0fiYiIJJEV78GD3+T6lJ/y+1/fEMg0FzOb4+65Vds1B11EREREmp+eR1HeOpvfd1qgOegiIiIiIoELhQgPOY9wOAUi5UGn2Y1G0EVERESkeTrxN2DJNXoOGkGvM83VlwOh7x8REZEklITFOahAr5OMjAw2b96sIkv2i7uzefNmMjIygo4iIiIijYCmuNRBdnY2eXl55OfnBx1FGqmMjAyys7ODjiEiIiKNgAr0OkhNTaV3795BxxARERGRZkBTXEREREREkogKdBERERGRJKICXUREREQkiZiuTLI7M8sHVgbQdQdgUwD9SmLpc24e9Dk3D/qcmz59xs1DkJ/zwe7esWqjCvQkYWaz3T036BzSsPQ5Nw/6nJsHfc5Nnz7j5iEZP2dNcRERERERSSIq0EVEREREkogK9OQxJegAkhD6nJsHfc7Ngz7npk+fcfOQdJ+z5qCLiIiIiCQRjaCLiIiIiCQRFehJwMzGmdkSM1tmZpODziP1z8x6mNnbZrbQzBaY2dVBZ5KGYWZhM/vEzF4MOos0DDNrZ2ZPmdliM1tkZkcFnUnqn5n9LPbzer6ZPWZmGUFnkgNnZveb2UYzmx/X1t7MXjezpbHng4LMCCrQA2dmYeBO4FSgP3CemfUPNpU0gDLg5+7eHxgJ/ESfc5N1NbAo6BDSoP4OvOLuhwGD0Ofd5JhZd+AqINfdBwBh4NxgU0k9eRAYV6VtMvCmu/cF3owtB0oFevCGA8vcfbm7lwCPAxMCziT1zN3XufvHsdfbif6H3j3YVFLfzCwbOA24N+gs0jDMrC1wHHAfgLuXuPvWQENJQ0kBWphZCtASWBtwHqkH7j4T2FKleQIwNfZ6KnBmIjNVRwV68LoDq+OW81Dh1qSZWS9gCPBhwFGk/t0G/BKIBJxDGk5vIB94IDaV6V4zyww6lNQvd18D3AKsAtYBBe7+WrCppAF1dvd1sdfrgc5BhgEV6CIJZWatgKeBn7r7tqDzSP0xs9OBje4+J+gs0qBSgKHA3e4+BNhBEvw5XOpXbA7yBKK/kHUDMs3swmBTSSJ49PKGgV/iUAV68NYAPeKWs2Nt0sSYWSrR4vwRd38m6DxS70YB481sBdGpaqPN7OFgI0kDyAPy3L3iL2BPES3YpWk5CfjS3fPdvRR4Bjg64EzScDaYWVeA2PPGgPOoQE8Cs4C+ZtbbzNKInoQyPeBMUs/MzIjOWV3k7n8LOo/UP3e/zt2z3b0X0X/Hb7m7RtyaGHdfD6w2s36xpjHAwgAjScNYBYw0s5axn99j0MnATdl0YGLs9UTg+QCzANE/1UmA3L3MzK4AXiV6lvj97r4g4FhS/0YBFwHzzGxurO3X7v5ycJFEZD9dCTwSG1RZDlwacB6pZ+7+oZk9BXxM9Cpcn5CEd5uUfWdmjwEnAB3MLA+4HrgJeMLMvgesBM4OLmGU7iQqIiIiIpJENMVFRERERCSJqEAXEREREUkiKtBFRERERJKICnQRERERkSSiAl1EREREJImoQBcRERERSSIq0EVEREREkogKdBER2S9mNszMPjOzDDPLNLMFZjYg6FwiIo2dblQkIiL7zcz+BGQALYA8d/+fgCOJiDR6KtCr6NChg/fq1SvoGCIiIiLSxM2ZM2eTu3es2p4SRJhk1qtXL2bPnh10DBERERFp4sxsZXXtmoMuIiIiIpJEVKCLiIiIiCQRFegiInXUa/JLXP/8/KBjiIhIE6c56CIi+2Dq+yv5wwRdSVCkMSotLSUvL4+ioqKgo0gzk5GRQXZ2NqmpqXXaXgW6iIiINAt5eXm0bt2aXr16YWZBx5Fmwt3ZvHkzeXl59O7du077aIqLiIiINAtFRUVkZWWpOJeEMjOysrL26S83KtBFRESk2VBxLkHY1+87FegiIiIiIklEBbqIiIhIkmnVqlWt67du3cpdd91Vubx27VrOOuusBsuzYsUKBgw48BPkb7jhBm655ZZat5k7dy4vv/zyAfe1L775zW+ydevWPdrrkrchqEAXERERSTB3JxKJ7Pf+VQv0bt268dRTT9VHtMAdaIFeVla2z/u8/PLLtGvXbr/7rG8q0EVEREQSYMWKFfTr14+LL76YAQMGsHr1am6++WaGDRtGTk4O119//R77FBYWMmbMGIYOHcrAgQN5/vnnAZg8eTJffPEFgwcP5tprr91thHvkyJEsWLCg8hgnnHACs2fPZseOHUyaNInhw4czZMiQymPVpT+IFr4XXHABhx9+OGeddRY7d+6szNK/f39ycnL4xS9+UfleR48eTU5ODmPGjGHVqlV79FWRC2DTpk306tWLkpISfv/73zNt2jQGDx7MtGnT6pR7xowZHHvssYwfP57+/ftTXl7OtddeW/m1/ec//wnAunXrOO644xg8eDADBgzgnXfeAaBXr15s2rQJgD//+c8ceuihHHPMMSxZsqTWvECNfR0IXWZRREREmp0/vLCAhWu31esx+3drw/VnHFHrNkuXLmXq1KmMHDmS1157jaVLl/LRRx/h7owfP56ZM2dy3HHHVW6fkZHBs88+S5s2bdi0aRMjR45k/Pjx3HTTTcyfP5+5c+cC0YK4wjnnnMMTTzzBH/7wB9atW8e6devIzc3l17/+NaNHj+b+++9n69atDB8+nJNOOonMzMy99gewZMkS7rvvPkaNGsWkSZO46667uPTSS3n22WdZvHgxZlY5TeTKK69k4sSJTJw4kfvvv5+rrrqK5557bq9fw7S0NG688UZmz57NP/7xD4A65Qb4+OOPmT9/Pr1792bKlCm0bduWWbNmUVxczKhRozjllFN45plnGDt2LL/5zW8oLy+v/CWjwpw5c3j88ceZO3cuZWVlDB06lCOPPLLWzPfdd1+1fdX1korV0Qi6iIiISIIcfPDBjBw5EoDXXnuN1157jSFDhjB06FAWL17M0qVLd9ve3fn1r39NTk4OJ510EmvWrGHDhg219nH22WdXTnd54oknKuemv/baa9x0000MHjyYE044gaKioj1Gtmvrr0ePHowaNQqACy+8kHfffZe2bduSkZHB9773PZ555hlatmwJwPvvv8/5558PwEUXXcS7776731+zuuQGGD58eGVR/Nprr/HQQw8xePBgRowYwebNm1m6dCnDhg3jgQce4IYbbmDevHm0bt16t2O88847fOtb36Jly5a0adOm8peTveWrrq8DkfARdDPrWcdNt7p7/f5qKyIiIgJ7HeluKPGjvu7Oddddxw9+8IMat3/kkUfIz89nzpw5pKam0qtXr71eT7t79+5kZWXx2WefMW3aNO65557K/p5++mn69eu3X/1VvVSgmZGSksJHH33Em2++yVNPPcU//vEP3nrrrb1+HQBSUlIq5+HX9p7qkhv2/NrecccdjB07do/tZs6cyUsvvcQll1zCNddcw8UXX3xAeWvra38FMYI+FXgw9lzT40HgzACyiYiIiCTE2LFjuf/++yksLARgzZo1bNy4cbdtCgoK6NSpE6mpqbz99tusXLkSgNatW7N9+/Yaj33OOefw17/+lYKCAnJycir7u+OOO3B3AD755JM99qupP4BVq1bx/vvvA/Doo49yzDHHUFhYSEFBAd/85je59dZb+fTTTwE4+uijefzxx4Fo0X/sscfu0VevXr2YM2cOwG4nuFZ9b3XJXdXYsWO5++67KS0tBeDzzz9nx44drFy5ks6dO3PZZZfx/e9/n48//ni3/Y477jiee+45du3axfbt23nhhRf2mremvg5EwkfQ3f3ERPcpIiIikmxOOeUUFi1axFFHHQVEL6348MMP06lTp8ptLrjgAs444wwGDhxIbm4uhx12GABZWVmMGjWKAQMGcOqpp/KTn/xkt2OfddZZXH311fzud7+rbPvd737HT3/6U3JycohEIvTu3ZsXX3xxt/1q6g+gX79+3HnnnUyaNIn+/fvzox/9iIKCAiZMmEBRURHuzt/+9jcA7rjjDi699FJuvvlmOnbsyAMPPLDH+//FL37B2WefzZQpUzjttNMq20888cTKKS3XXXddnXJX9f3vf58VK1YwdOhQ3J2OHTvy3HPPMWPGDG6++WZSU1Np1aoVDz300G77DR06lHPOOYdBgwbRqVMnhg0btte8NfV1IKzitxGJys3N9YozdEVE4vWa/BIAK246bS9bikgyWrRoEYcffnjQMaSZqu77z8zmuHtu1W2T6iRRM7s06AwiIiIiIkFKqgId+EPQAUREREREghTEVVw+q2kV0DmRWURERKR5cfc9rkYi0tD2dUp5EDcq6gyMBb6q0m7AfxMfR0RERJqDjIwMNm/eTFZWlop0SRh3Z/PmzWRkZNR5nyAK9BeBVu4+t+oKM5uR8DQiIiLSLGRnZ5OXl0d+fn7QUaSZycjIIDs7u87bB3GZxe/Vsu78RGYRERGR5iM1NfWAbr8ukijJdpKoiIiIiEiz1qgLdDMbZ2ZLzGyZmU2uZv01ZrbQzD4zszfN7OAgcoqIiIiI1FXgBbqZnbGf+4WBO4FTgf7AeWbWv8pmnwC57p4DPAX89UCyioiIiIg0tMALdODP+7nfcGCZuy939xLgcWBC/Abu/ra774wtfgDUfXa+iIiIiEgAkqFA39/rHHUHVsct58XaavI94N/72ZeIiIiISEIEcZnFqvbtyu37wcwuBHKB42tYfzlwOUDPnj0bOo6IiIiISI2SYQR9f60BesQtZ8fadmNmJwG/Aca7e3F1B3L3Ke6e6+65HTt2bJCwItK47etd4ERERPZXYy7QZwF9zay3maUB5wLT4zcwsyHAP4kW5xsDyCgiTYTqcxERSZRkKNA37M9O7l4GXAG8CiwCnnD3BWZ2o5mNj212M9AKeNLM5prZ9BoOJyIiIiKSFAKfg+7uJx/Avi8DL1dp+33c65MOIJqISCUNoIuISKIkwwi6iEjS0xx0ERFJFBXoIiIiIiJJJPAC3cwyY3cFFRFJWho/FxGRREl4gW5mITM738xeMrONwGJgnZktNLObzeyQRGcSEdkbzXAREZFECWIE/W3gG8B1QBd37+HunYBjgA+Av8RuLCQikjRcY+giIpIgQVzF5SR3L63a6O5bgKeBp80sNfGxRERERESCl/AR9Iri3Mz+bmZW2zYiIslCU1xERCRRgjxJdDsw3cwyAcxsrJm9F2AeEREREZHABXajInf/rZmdD8wwsxKgEJgcVB4RERERkWQQWIFuZmOAy4AdQFdgkrsvCSqPiEhtNMVFREQSJcgpLr8BfufuJwBnAdPMbHSAeUREaqSruIiISKIEOcVldNzreWZ2KtGruBwdVCYRERERkaAFcaOimq7csg4YU9s2IiJB0RQXERFJlEBuVGRmV5pZz/hGM0sDjjKzqcDEAHKJiNRI9bmIiCRKEFNcxgGTgMfMrDewFcgAwsBrwG3u/kkAuUREauQaQhcRkQRJeIHu7kXAXcBdsTuGdgB2ufvWRGcREREREUk2gV3FxcxeB/q7+zoV5yKS7DR+LiIiiRLkZRZ/BdxmZg+YWdcAc4iI7JVmuIiISKIEVqC7+8fufiLwIvCKmV1vZi2CyiMiIiIikgyCHEGvuJziEuBu4EpgqZldFGQmEZFqaQRdREQSJMg56O8Ba4Bbge7AJcAJwHAzmxJULhGR6uhOoiIikigJv4qLmR0FfABcDiz0Pa9ddqWZLUp0LhGR2mgOuoiIJEoQI+gXA3OA3wETzaxLNducVpcDmdk4M1tiZsvMbHI1648zs4/NrMzMzjqw2CIiIiIiDS+I66D/CMDMDgNOBR40s7bA28ArwHvuvnxvxzGzMHAncDKQB8wys+nuvjBus1VEp878ol7fhIg0OxpAFxGRRAnyKi6L3f1Wdx8HjAbeBb4LfFjHQwwHlrn7cncvAR4HJlTpY4W7fwZE6jG6iDRDupOoiIgkSsJH0Ksys0ygyN1fBl7eh127A6vjlvOAEfWZTUREREQk0RI+gm5mITM738xeMrONRC+zuN7MFprZzWZ2SACZLjez2WY2Oz8/P9Hdi0gjoPFzERFJlCCmuLwNfAO4Duji7tnu3hE4hujVXf5iZhfW4ThrgB5xy9mxtn3m7lPcPdfdczt27Lg/hxCRJk4zXEREJFGCmOJykruXVm109y3A08DTZpZah+PMAvqaWW+ihfm5wPn1mlREJEbXQRcRkURJ+Ah6RXFuZn+P3Um0xm32cpwy4ArgVWAR8IS7LzCzG81sfKyPYWaWR/Tk03+a2YL6eh8iIiIiIg0hyJNEtwPTzexcd99hZmOB37v7qLoeoLoTS93993GvZxGd+iIicmA0gC4iIgkSWIHu7r81s/OBGWZWAhQCe9xsSEQkGag+FxGRRAmsQDezMcBlwA6gKzDJ3ZcElUdEREREJBkEdqMi4DfA79z9BOAsYJqZjQ4wj4hIjXQVFxERSZQgp7iMjns9z8xOJXoVl6ODyiQiUhNdxUVERBIliBsV1XTllnXAmNq2EREJikbQRUQkUQK5UZGZXWlmPeMbzSwNOMrMpgITA8glIiIiIhK4IKa4jAMmAY+ZWR/gK6AF0V8WXgNuc/dPAsglIlIjDaCLiEiiJLxAd/ci4C7grtgdQzsAu9x9a6KziIjUlWuOi4iIJEiQl1m8BrgY2ALMM7NPgU+B+e5eHFQuEREREZEgBXkn0SuBk4FyIAcYDJwBDDCzYncfEGA2EZHdaABdREQSJcgCfQnwhUf/bvwl8HzFCjNrG1gqEREREZEABXmjoo3A/WbWu+oKdy8III+IiIiISOCCLNDnEx3Bf9bM8szsdTP73wDziIjUqLQ8EnQEERFpJoK8k+hfK16bWQrQDxgYVB4RkdoUlapAFxGRxEh4gV71BkVxtgP/jVu/1d23JSiWiEitisvKg44gIiLNRBAj6FOJ3vPDatnGgQeBhxIRSERkbzSCLiIiiRLEjYpOTHSfIiIHqkgj6CIikiBBniQqItJoFGsEXUREEkQFuohIHWwrKq18XaYruoiISANSgS4iUgebCosrX2/ZURJgEhERaepUoIuI1MHarbsqX+fHFesiIiL1TQW6iEgdzMsrIByKXnwqf7sKdBERaTiB3aioPpjZOODvQBi4191vqrI+neilGo8ENgPnuPuKROcUkcbtw+Wb+TSvgAtG9OSRD1exrqAo6Egi0oDcHXeIuBOJPX+9HG3zuHUV673K8u77x20fqfsx49fVeszK/as7/l62j8te2zaVx4zU7ZjxGao+V/+ea3sPNW9f+9e9lv0j0fXdD2rBKz89Luhvu9002gLdzMLAncDJQB4wy8ymu/vCuM2+B3zl7oeY2bnAX4BzEp9WRBqDSMTZUVLG9qIyNm4v5stNhXz05Vc883EevbJa8stxh/Hqgg28u3QT5w2v6Z5r0pCqFilONYUR4JHouvj/pD1un6r/sVc+E7ddNcVXdcd0vi5YaswTKzicfcyz2/qvC6KKPqPHq7JvNX3skTfuPdb8HuL23e097L7v131W/x4iDrBncQTVbFflPUTb9iy22O091L0oq/r1qKlwi3UhQMggZEbIDKt8zdfLIatss7h1Feut1v3jt69+/3DISA1ZjdsbseVQxf5V8u21P6N9ZmrQX+Y9NNoCHRgOLHP35QBm9jgwAYgv0CcAN8RePwX8w8zMvXH/04tEnHJ3yiN7/oCr+h9T/H8A8f9pVPwwjt+/6g/12M9Foq+o/IFV+cyeP8Ti1+2+XLHe415X7lXDttHtq7bFH5/q9qml7/hlan1v1e/jXnvemt571W+5Grf1r7ev6Krys4t7fxXLVPl67v65xeWv5ThfZ9p7P3u0x3/N47PX0E9N31cVbcR9n1aXN76Pr7+2X3/vRiJOWST6byP6HKF8t2WnrDz6XPFvqCwSobTM2VFcRmFJ2R7f0y3TwpyW05VfjTuMti1SmTC4Gw+89yVPzcnj20O6YwZf7Sxl684Stu4qZduuUorLIpTEHtHX5RSXRSiPZSyPQLlHC62KLBXFUXkkWjCUV4xQVayv8m8zvjBx9iycqhZE8Z9N1cKvalFU9XupusIo/mdJ9fvv+fOnukIt/vulLrnkwMUXL7sVOBXLFr2TYEXhZXxdOFUUPJVtITB23zcUd0yzavatXB+/zir7TNmtuIsrrojvI+74WCxr7UXZvheKsbbQ3rentv1366+aQrTWwvLr7Pt0zN2K571vb7VsI8FozAV6d2B13HIeMKKmbdy9zMwKgCxgU0IS1tEbCzdw86tLKItEKIsrICoKjLLy+KIjEvtPTiSxKv7TrviPtOI/xuhKdmurui3xy1XXxfZjt/32PE5lhrh1sUNXHjclZIRDodjz14/UcIiMVNutPSUUij1Hl1tlpNA6PYXWGam0ykghKzONPh0z6dk+k7SUr0/XufqkvsxdvZVfPPkpf3ppISVlEXaW7N9NjMwgbNHiIhxXDIRjy2ZGOETl6+qKofhCKFT5daq+mNp9ufaiqGKb6gqjaEEU+/x3K8jic8UVFrHtoeZcVmXd7setqZjc/X1VZPu6mKmmSGTPgif6Nfh6390Klt2+nl8fs2q+aJEW1xYrjCrec9XPqLIwqylPZe49C6j474HaC+O4Yi3uMxeR5NeYC/R6Y2aXA5cD9OyZ+D9bt8pIoXeHTMLhaLGQUlFghL8uHlJCRko4rvAw+3p0I+4/zd3/I4r95k/1P8Dj/zPZ7Yd/lTagSiH1dUtFgQBxRVSVbb8upL6uqPa6TzXHp5p1u+9bXaY981bdp7a+48Wvr+m97Xn8r/uvMe9ejl+1IK58f1WX2XNbKguk6o9TXb97bKv/0HfTJiOVaZeP5NUFG5j5eT6Z6SlkH9SC9plptG2ZSpuMVNJTQmSkhkgLh0lLCZGeEiI1JfrvN2TRf8ManRIRkZo05gJ9DdAjbjk71lbdNnlmlgK0JXqy6G7cfQowBSA3Nzfh49Mj+2Qxsk9WorsVkf2UEg5xWk5XTsvpGnQUERFpghrzZRZnAX3NrLeZpQHnAtOrbDMdmBh7fRbwVmOffy4iIiIiTVujHUGPzSm/AniV6GUW73f3BWZ2IzDb3acD9wH/MrNlwBaiRbyIiIiISNIyDSjvzszygZUBdN2BJDt5VRqEPufmQZ9z86DPuenTZ9w8BPk5H+zuHas2qkBPEmY2291zg84hDUufc/Ogz7l50Ofc9Okzbh6S8XNuzHPQRURERESaHBXoIiIiIiJJRAV68pgSdABJCH3OzYM+5+ZBn3PTp8+4eUi6z1lz0EVEREREkohG0EVEREREkogK9CRgZuPMbImZLTOzyUHnkfpnZj3M7G0zW2hmC8zs6qAzScMws7CZfWJmLwadRRqGmbUzs6fMbLGZLTKzo4LOJPXPzH4W+3k938weM7OMoDPJgTOz+81so5nNj2trb2avm9nS2PNBQWYEFeiBM7MwcCdwKtAfOM/M+gebShpAGfBzd+8PjAR+os+5yboaWBR0CGlQfwdecffDgEHo825yzKw7cBWQ6+4DiN4QUTc7bBoeBMZVaZsMvOnufYE3Y8uBUoEevOHAMndf7u4lwOPAhIAzST1z93Xu/nHs9Xai/6F3DzaV1DczywZOA+4NOos0DDNrCxxH9E7VuHuJu28NNJQ0lBSghZmlAC2BtQHnkXrg7jOJ3l0+3gRgauz1VODMRGaqjgr04HUHVsct56HCrUkzs17AEODDgKNI/bsN+CUQCTiHNJzeQD7wQGwq071mlhl0KKlf7r4GuAVYBawDCtz9tWBTSQPq7O7rYq/XA52DDAMq0EUSysxaAU8DP3X3bUHnkfpjZqcDG919TtBZpEGlAEOBu919CLCDJPhzuNSv2BzkCUR/IesGZJrZhcGmkkTw6OUNA7/EoQr04K0BesQtZ8fapIkxs1Sixfkj7v5M0Hmk3o0CxpvZCqJT1Uab2cPBRpIGkAfkuXvFX8CeIlqwS9NyEvClu+e7eynwDHB0wJmk4Wwws64AseeNAedRgZ4EZgF9zay3maURPQllesCZpJ6ZmRGds7rI3f8WdB6pf+5+nbtnu3svov+O33J3jbg1Me6+HlhtZv1iTWOAhQFGkoaxChhpZi1jP7/HoJOBm7LpwMTY64nA8wFmAaJ/qpMAuXuZmV0BvEr0LPH73X1BwLGk/o0CLgLmmdncWNuv3f3l4CKJyH66EngkNqiyHLg04DxSz9z9QzN7CviY6FW4PiEJ7zYp+87MHgNOADqYWR5wPXAT8ISZfQ9YCZwdXMIo3UlURERERCSJaIqLiIiIiEgSUYEuIiIiIpJEVKCLiIiIiCQRFegiIiIiIklEBbqIiIiISBJRgS4iIiIikkRUoIuIiIiIJBEV6CIiIiIiSUQFuoiIiIhIElGBLiIiIiKSRFKCDpBsOnTo4L169Qo6hoiIiIg0cXPmzNnk7h2rtqtAr6JXr17Mnj076BgiIiIi0sSZ2crq2jXFRUREREQkiahAFxEREZFm6c1FG/jX+ysoj3jQUXajAl1EREREmp3C4jJ++9x8Hv1oNe7JVaBrDnodlJaWkpeXR1FRUdBRpB5lZGSQnZ1Nampq0FFEREQkwW57/XPWbyvizguGkhJOrjFrFeh1kJeXR+vWrenVqxdmFnQcqQfuzubNm8nLy6N3795BxxEREZEEWrC2gAf+u4LzhvdkaM+Dgo6zh+T6dSFJFRUVkZWVpeK8CTEzsrKy9FcRERGRZiYScX7z7HzatUjlV2MPCzpOtRpVgW5m95vZRjObH9fW3sxeN7OlseeDYu1mZreb2TIz+8zMhh5g3wcaX5KMPlMREZHm57FZq5i7eiu/Pf1w2rZMzmmujapABx4ExlVpmwy86e59gTdjywCnAn1jj8uBuxOUMaFatWp1wMd48MEHueKKK2rdZsWKFTz66KMH3Ne++P73v8/ChQv3aK9LXhEREZGq1hcUcdO/F3NUnyzOHNw96Dg1alQFurvPBLZUaZ4ATI29ngqcGdf+kEd9ALQzs64JCdoEHWiB7u5EIpF92ufee++lf//++92niIiISAV35zfPzqO0PML/fHtgUv8lvVEV6DXo7O7rYq/XA51jr7sDq+O2y4u1NUpnnnkmRx55JEcccQRTpkzZbd3PfvYzjjjiCMaMGUN+fj4At99+O/379ycnJ4dzzz0XgC1btnDmmWeSk5PDyJEj+eyzz/bo55JLLuGpp56qXK4YoZ88eTLvvPMOgwcP5tZbb6W8vJxrr72WYcOGkZOTwz//+c89jrVixQr69evHxRdfzIABA1i9ejU333xz5T7XX389ADt27OC0005j0KBBDBgwgGnTpgFwwgknVN7V9YEHHuDQQw9l+PDhvPfee3vNC1Tbl4iIiDRPz89dy5uLN/KLU/rRq0Nm0HFq1aSu4uLubmb7fCFLM7uc6DQYevbsWeu2f3hhAQvXbtu/gDXo360N159xRK3b3H///bRv355du3YxbNgwvvOd75CVlcWOHTvIzc3l1ltv5cYbb+QPf/gD//jHP7jpppv48ssvSU9PZ+vWrQBcf/31DBkyhOeee4633nqLiy++mLlz59Yp40033cQtt9zCiy++CMCUKVNo27Yts2bNori4mFGjRnHKKafscUWUpUuXMnXqVEaOHMlrr73G0qVL+eijj3B3xo8fz8yZM8nPz6dbt2689NJLABQUFOx2jHXr1nH99dczZ84c2rZty4knnsiQIUNqzVtTX8cdd1yd3q+IiIg0HRu3F3HDCwsY2rMdl45K/qu3NYUR9A0VU1dizxtj7WuAHnHbZcfa9uDuU9w9191zO3bs2KBh99ftt9/OoEGDGDlyJKtXr2bp0qUAhEIhzjnnHAAuvPBC3n33XQBycnK44IILePjhh0lJif4e9u6773LRRRcBMHr0aDZv3sy2bfv3y8Zrr73GQw89xODBgxkxYgSbN2+uzBTv4IMPZuTIkZX7vPbaawwZMoShQ4eyePFili5dysCBA3n99df51a9+xTvvvEPbtm13O8aHH37ICSecQMeOHUlLS6t8v3vLV11fIiIi0vxc//wCdpaU89ezBhEOJe/UlgpNYQR9OjARuCn2/Hxc+xVm9jgwAiiImwqz3/Y20t0QZsyYwRtvvMH7779Py5YtOeGEE2q8PGDFfKqXXnqJmTNn8sILL/DnP/+ZefPm1amvlJSUyrnikUiEkpKSardzd+644w7Gjh1b6/EyM7/+E5K7c9111/GDH/xgj+0+/vhjXn75ZX77298yZswYfv/73x9Q3tr6EhERkebjxc/W8u/56/nVuMM4pNOBX1wjERrVCLqZPQa8D/Qzszwz+x7RwvxkM1sKnBRbBngZWA4sA/4P+HEAketFQUEBBx10EC1btmTx4sV88MEHlesikUjlHOxHH32UY445hkgkwurVqznxxBP5y1/+QkFBAYWFhRx77LE88sgjQLTo79ChA23atNmtr169ejFnzhwApk+fTmlpKQCtW7dm+/btlduNHTuWu+++u3L9559/zo4dO2p9H2PHjuX++++nsLAQgDVr1rBx40bWrl1Ly5YtufDCC7n22mv5+OOPd9tvxIgR/Oc//2Hz5s2Ulpby5JNP7jVvTX2JiIhI87GuYBe/eXY+g7LbctmxyT+1pUKjGkF39/NqWDWmmm0d+EnDJkqMcePGcc8993D44YfTr1+/yikjEB2h/uijj/jTn/5Ep06dmDZtGuXl5Vx44YUUFBTg7lx11VW0a9eOG264gUmTJpGTk0PLli2ZOnXqHn1ddtllTJgwgUGDBjFu3LjKEfCcnBzC4TCDBg3ikksu4eqrr2bFihUMHToUd6djx44899xztb6PU045hUWLFnHUUUcB0RM6H374YZYtW8a1115LKBQiNTWVu+/e/YqYXbt25YYbbuCoo46iXbt2DB48eK95a+qrU6dO+/z1FxERkcYnEnF+8eSnlJRFuPWcwaSEG8+4tEXrWKmQm5vrFVcOqbBo0SIOP/zwgBJJQ9JnKyIi0jT938zl/PnlRfzPtwdy3vDaLwISFDOb4+65Vdsbz68SIiIiIiJ1sHDtNm5+dQkn9+/MucN67H2HJKMCXURERESajKLScq5+/BPatkzlL9/JSeobEtUkIXPQzayuf1fY6u71e5FxEREREWk2/vjiQpZuLGTqpOG0z0wLOs5+SdRJolMBB2r7FcaBB4GHEhFoX7l7o/wNTGqm8y9ERESalufnruGRD1fxg+P7cPyhyXlvm7pISIHu7icmop+GkpGRwebNm8nKylKR3kS4O5s3byYjIyPoKCIiIlIPlm3cznXPzGNYr4P4xSn9go5zQAK5zKKZZQJF7l4eRP/7Kjs7m7y8PPLz84OOIvUoIyOD7OzsoGOIiIjIAdpZUsaPHv6YFqlh7jhvKKmN6JKK1UnUHPQQcC5wATAMKAHSzSwfeAn4p7svS0SW/ZGamkrv3o3n4vYiIiIizYW789tn57Msv5B/TRpBl7aN/6/jifr14m3gG8B1QBd3z3b3jsAxwAfAX8zswgRlEREREZEm4uEPVvLMJ2u4ekxfjunbIeg49SJRU1xOcvfSqo3uvgV4GnjazFITlEVEREREmoD3v9jMH15YyOjDOnHl6L5Bx6k3CRlBryjOzezvVsNZltUV8CIiIiIi1Vm9ZSc/fmQOB2e15LZzBxMONZ0LeSR6Bv12YHrsJFHMbKyZvZfgDCIiIiLSiO0oLuOyh2ZTHnHunTiMNhlNayJGQq/i4u6/NbPzgRlmVgIUApMTmUFEREREGq9IxPn5E5/y+YbtPHjpcHp3yAw6Ur1LaIFuZmOAy4AdQFdgkrsvSWQGEREREWm8bnplMa8sWM9vTzuc4xrxzYhqk+gpLr8BfufuJwBnAdPMbHSCM4iIiIhIIzT1vyuYMnM5Fx91MN87puleAjvRU1xGx72eZ2anEr2Ky9GJzCEiIiIijcurC9ZzwwsLOOnwzlx/xhFN+u7uCRlBr+XKLeuAMbVtIyIiIiLN28ervuKqxz4hJ7sdd5w3pEldsaU6CbtRkZldaWY94xvNLA04ysymAhMTlEVEREREGomlG7bzvQdn0blNBvdNzKVFWjjoSA0uUVNcxgGTgMfMrDewFcgAwsBrwG3u/kmCsoiIiIhII7Bq804uvO9DUsIhHpo0nA6t0oOOlBAJKdDdvQi4C7grdsfQDsAud9+aiP5FREREpHFZX1DE+fd+QHFZhGmXH0WvJng5xZok9CouZvY60N/d16k4FxEREZHqbCos5oJ7P2DrzlIemjScfl1aBx0poRJ9mcVfAbeZ2QNm1jXBfYuIiIhIkvtqRwkX3/cRa7bu4v5LhpGT3S7oSAmX0ALd3T929xOBF4FXzOx6M2uRyAwiIiIikpw2FxZz3v99wLL8Qv55US7De7cPOlIgEj2CXnE5xSXA3cCVwFIzuyjROUREREQkeeRvjxbnX27awX0Tczm+id4ltC4SPQf9PWANcCvQHbgEOAEYbmZTEplFRERERJLDhm1FnDvlfVZv2cUDlw7j2L7NtziHBN9JFLgcWOjuXqX9SjNblOAsIiIiIhKw1Vt2ctF9H5K/vZipk4Y322kt8RJaoLv7glpWn5awICIiIiISuEXrtjHx/o8oKi3noe+N4MiDDwo6UlJI9Ah6jdx9edAZRERERCQxPli+mcumzqZVRgpP/ehoDu3cvC6lWJuEnyQKYGZnBNGviIiIiATvlfnruPj+j+jcNoOnVZzvIZACHfhzQP2KiIiISEDcnfve/ZIfPfIxA7q14akfHkW3drridlVBTXGxgPoVERERkQCUlEW4fvp8HvtoNWOP6Mxt5wyhRVo46FhJKagCvepVXA6Yma0AtgPlQJm755pZe2Aa0AtYAZzt7l/Vd98iIiIiUrOvdpTww4fn8OGXW7jixEO45uRDCYU0XluToKa4NJQT3X2wu+fGlicDb7p7X+DN2LKIiIiIJMjSDduZcOd7fLJ6K7edM5hfjO2n4nwvmlqBXtUEYGrs9VTgzOCiiIiIiDQvz89dw4Q732NnSTmPXz6SM4d0DzpSoxDUFJcNDXBMB14zMwf+6e5TgM7uvi62fj3QuQH6FREREZE4xWXl/OnFRfzrg5XkHnwQ/zh/KF3aZgQdq9EIpEB395Mb4LDHuPsaM+sEvG5mi6v06bHifQ9mdjnRu5zSs2fPBogmIiIi0jzkfbWTnzzyMZ/mFXDZsb355bjDSA039Ukb9StpblR0oNx9Tex5o5k9CwwHNphZV3dfZ2ZdgY017DsFmAKQm5tb7yewioiIiDQH/563jsnPzCMSce65cCjjBnQNOlKj1CR+nTGzTDNrXfEaOAWYD0wHJsY2mwg8H0xCERERkaarsLiMXzz5KT965GMOzmrJ9CuPUXF+AAIZQY8V0UXuXl5Ph+wMPGtmEH1Pj7r7K2Y2C3jCzL4HrATOrqf+RERERASYs/IrfjZtLnlf7eSKEw/h6pP6akrLAUpIgW5mIeBc4AJgGFAMpJvZJuAloid1Ltvf47v7cmBQNe2bgTH7e1wRERERqV5RaTl3vLWUe/6znK5tM5j2g6MY1qt90LGahESNoL8NvAFcB8x39whA7EZCJwJ/MbNn3f3hBOURERERkf00a8UWfvX0ZyzP38F3j8zm92f0p3VGatCxmoxEFegnuXtp1UZ33wI8DTxtZvpURURERJLY9qJS/vrKEv71wUqyD2rBQ5OGc9yhHYOO1eQkpECvKM7N7O/AT919jyulVFfAi4iIiEjw3J1X5q/njy8uZN22IiaN6s3PTzmUzPQmc0HApJLor+p2YLqZnevuO8xsLPB7dx+V4BwiIiIiUgdLN2znhhcW8N6yzRzWpTX/uGAoQ3seFHSsJi2hBbq7/9bMzgdmmFkJUAhMTmQGEREREdm7bUWl3Pb6Uqa+v4JW6SncOOEIzh/ekxRdoaXBJbRAN7MxwGXADqArMMndlyQyg4iIiIjUrKQswuOzVnH7m0vZvKOEc4f15Nqx/WifmRZ0tGYj0VNcfgP8zt3fNbOBwDQzu8bd30pwDhERERGJE4k4L85bxy2vLmHVlp2M6N2eBy7pz8DstkFHa3YSPcVldNzreWZ2KtGruBydyBwiIiIiEuXuvLN0E399dTHz12zjsC6teeDSYZxwaEdiN4GUBEvUjYqshiu3rItNe6lxGxERERGpf+7OjCX53P7WUj5ZtZXu7Vpw6zmDmDCoO6GQCvMgJexGRWb2NPC8u6+qaDSzNOAoM5tI9GZGDyYoj4iIiEizFIk4ry/awD/eWsa8NQV0b9eCP505gO/mZpOeEg46npC4An0cMAl4zMz6AF8BLYAQ8Bpwm7t/kqAsIiIiIs1OcVk50+eu5b53v2Tx+u0cnNWSv34nh28N7U6qrsySVBJ1o6Ii4C7grtgdQzsAu9x9ayL6FxEREWmuNhcW8/AHq/jXByvZVFhMv86t+dvZgxg/qJsumZikEn2ZxWuAi4EtwDwz+xT4FJjv7sWJzCIiIiLSlM1fU8C/3l/Js3PXUFIW4cR+HfneMX0YdUiWTv5Mcom+zOKVwMlAOZADDAbOAAaYWbG7D0hwHhEREZEmo7C4jOlz1/LYR6uYt6aAjNQQZx2ZzaRRvTmkU6ug40kdJbpAXwJ8Ebtay5fA8xUrzEwX2RQRERHZR+7OZ3kFPD5rFc/PXcvOknIO69KaGyccwYTB3WnbIjXoiLKPEl2gbwTuN7Mb3f3L+BXuXpDgLCIiIiKN1srNO3juk7U8P3cNyzftoEVqmDMGdeW84T0Z3KOdprE0Yoku0OcDA4FnzawDsAj4zN1/nuAcIiIiIo3OpsJiXvpsHc9+soa5q7diBiN6t+fy4/rwzZyutMnQaHlTkOg7if614rWZpQD9iBbsIiIiIlKN1Vt28uqC9by2YAOzV24h4nBYl9Zcd+phnDGoG93atQg6otSzRN1JtGcNq7YD/41bv9XdtyUik4iIiEgycneWbNjOaws28Mr89SxcFy2NDuvSmitG9+W0gV3p16V1wCmlISVqBH0q4EBtk6Gc6J1EH0pEIBEREZFkUbCrlPeWbeI/S/KZuTSfdQVFmMHQngfx628extgjunBwVmbQMSVBEnWjohMT0Y+IiIhIY1BaHuGzvIJoUf55PnNXb6U84rTOSOGYQzpw1ZiOjDm8E51aZwQdVQKQ6JNERURERJqd4rJy5q7aykdfbuHDL7cwZ+VX7Cotxwxyurflxyd8g+MP7cjgHu10d09RgS4iIiJS3zZuL2Luqq18mreV2Su+4pPVWykpi2AG/Tq35pxhPRjRuz0j+mTRPjMt6LiSZFSgi4iIiByAnSVlzMsr4NO8rcxdvZVPVxewZusuAFJCxuFd23DxyIMZ0SeLYb0Ool1LFeRSOxXoIiIiInXg7qwrKGLRum2xx3YWrdvGl5t34B7dpkf7Fgzp2Y5LR/ViSM92HNGtLRmp4WCDS6OjAl1EREQkjruTv72YZfmFfJG/gy82FrJo3TYWr99Owa7Syu16tG/B4V3acMagbuRkt2VQj3Z0aJUeYHJpKlSgi4iISLO0s6SM1Vt28eWmrwvxL/ILWZ6/g+3FZZXbtUwLc2jn1nxzYFf6d23N4V3b0K9La1rrrp3SQFSgi4iISJNUWh5h7dZdrN6yi9Vf7WTVlp2s3rKT1V/tIm/LTjbvKNlt+y5tMvhGp0y+NbQ73+jYim90bEWfjpl0aZNBKFTbrVxE6pcKdBEREWl0dpaUsb6giPXbitiwrYh1BUVsiC2v31bM+oJd5G8vJuJf75MSMrq1a0GP9i04uX9nerRvSfZBLejToRW9O2bSKl1lkSQHfSeKiIhI4NydbbvK2LSjmM2FJWwuLGbTjhK2FJawOda2qbCYzTtK2LitiG1FZXsco3VGCl3aZNClbQZ9O3Wka9sMehzUkuz2LejZviVd2mToGuPSKKhAFxERkXpTXFZOwa5Stu0qpWBXKVt3Rp+rPrbFvd66s5QtO0ooix/ujtOuZSpZmWlkZabTt1Mrjv5GFl3aZkSL8TYZdI69ztQIuDQRTf472czGAX8HwsC97n5TwJFEREQCV1Yeoagswq6ScopKyykuK6eoNEJRaTk7S8rZUVxGYXEZO4rL2BFbjraVs7Mkbl1xOTtKvn5dUh6ptd9W6Sm0bZFKmxaptG2RQu8OmbRrkUZWqzSyWqXToVUa7WPFeIdWaRyUmUaqRr2lmWnSBbqZhYE7gZOBPGCWmU1394XBJhMRkabO3SmPOOUVz3GPsohTUhahtDxCablTWh6hpDxCaVmV5YpHme++XB6/f2y5PEJxaYSisnKKS78utneVRgvwotLIbkV4TaPVNQmHjMy0MK3SU8iMPVqlp9ChVfpuba0zUmLF956PNhkpmmIiUgdNukAHhgPL3H05gJk9DkwAkqpAX7t1F/PWFOx1O6/zz9K6bVjX49W127ofL6h8+/af0d6PV8ft6vv9Bvh1qedvwQC/F+p6vGDy1fWA9f5+G6DfSOyYEXciXqUtEmujYp1XrvfYPh5rj/iey1/vV9FWsX/sWLv1G9sXJxL5uo+KY1UtoMs9WkRHIvHPESJO9DkSfS6PQHkkstt+Fa/3sf7dL6lhIzUcqnykhY2MtDAZKWEyUkNkpIZp0yI1+jolvMe6yueUMOmpIVqkhslIDdMyLUzLtJRY4R0mMz2F9JQQZrqSiUgiNPUCvTuwOm45DxhRdSMzuxy4HKBnz56JSRbng+WbueaJTxPer4hIMjGDkBlG7Nm+bqtoN4NQqOo2RsjAiD3H2kO2+3LFPnvsZxAOhQgbpIRChEKQFgoTChnhinWhinVGSih6jJSQVS6Hqz5sz7bK/cLR5WhBXVFcG6kpIdLDIVJTvm6rXJ9SZTm2XgWzSNPU1Av0OnH3KcAUgNzc3ASMeexu9GGdeOmqY+q0bfS/pTpsV8ef2XXerr77rdtmdT5eXY9Y//nq2G+dj1fH7er5/e6Lun9vNY2vTT1/Cwb4bymgzyOuSA7FimbYvXiueBYRkaimXqCvAXrELWfH2pJKu5ZptGuZFnQMEREREUkCTf1MjVlAXzPrbWZpwLnA9IAziYiIiIjUqEmPoLt7mZldAbxK9DKL97v7goBjiYiIiIjUyOr7yhaNnZnlAysD6LoDsCmAfiWx9Dk3D/qcmwd9zk2fPuPmIcjP+WB371i1UQV6kjCz2e6eG3QOaVj6nJsHfc7Ngz7npk+fcfOQjJ9zU5+DLiIiIiLSqKhAFxERERFJIirQk8eUoANIQuhzbh70OTcP+pybPn3GzUPSfc6agy4iIiIikkQ0gi4iIiIikkRUoCcBMxtnZkvMbJmZTQ46j9Q/M+thZm+b2UIzW2BmVwedSRqGmYXN7BMzezHoLNIwzKydmT1lZovNbJGZHRV0Jql/Zvaz2M/r+Wb2mJllBJ1JDpyZ3W9mG81sflxbezN73cyWxp4PCjIjqEAPnJmFgTuBU4H+wHlm1j/YVNIAyoCfu3t/YCTwE33OTdbVwKKgQ0iD+jvwirsfBgxCn3eTY2bdgauAXHcfQPRmh+cGm0rqyYPAuCptk4E33b0v8GZsOVAq0IM3HFjm7svdvQR4HJgQcCapZ+6+zt0/jr3eTvQ/9O7BppL6ZmbZwGnAvUFnkYZhZm2B44D7ANy9xN23BhpKGkoK0MLMUoCWwNqA80g9cPeZwJYqzROAqbHXU4EzE5mpOirQg9cdWB23nIcKtybNzHoBQ4APA44i9e824JdAJOAc0nB6A/nAA7GpTPeaWWbQoaR+ufsa4BZgFbAOKHD314JNJQ2os7uvi71eD3QOMgyoQBdJKDNrBTwN/NTdtwWdR+qPmZ0ObHT3OUFnkQaVAgwF7nb3IcAOkuDP4VK/YnOQJxD9hawbkGlmFwabShLBo5c3DPwShyrQg7cG6BG3nB1rkybGzFKJFuePuPszQeeRejcKGG9mK4hOVRttZg8HG0kaQB6Q5+4VfwF7imjBLk3LScCX7p7v7qXAM8DRAWeShrPBzLoCxJ43BpxHBXoSmAX0NbPeZpZG9CSU6QFnknpmZkZ0zuoid/9b0Hmk/rn7de6e7e69iP47fsvdNeLWxLj7emC1mfWLNY0BFgYYSRrGKmCkmbWM/fweg04GbsqmAxNjrycCzweYBYj+qU4C5O5lZnYF8CrRs8Tvd/cFAceS+jcKuAiYZ2ZzY22/dveXg4skIvvpSuCR2KDKcuDSgPNIPXP3D83sKeBjolfh+oQkvNuk7Dszeww4AehgZnnA9cBNwBNm9j1gJXB2cAmjdCdREREREZEkoikuIiIiIiJJRAW6iIiIiEgSUYEuIiIiIpJEVKCLiIiIiCQRFegiIiIiIklEBbqIiIiISBJRgS4iIiIikkRUoIuIiIiIJJH/D62bgWYG68rYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    For the given initial value problem y' = (x + 1)^2 + 1.0, taking x from 0 to 10 in 3000 steps,\n",
      "    we find a large deviation between exact and estimate solutions at the point that y' = 0.\n",
      "    We could fix this by taking some other measure, like the absolute difference, but as\n",
      "    the solutions decrease so hard, that will suck later on.\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "# This here is for debugging. We thought the keep it here,\n",
    "# as debugging was part of the task.\n",
    "\n",
    "def ODEF(x, y):\n",
    "    return y - x**2 + 1.0\n",
    "\n",
    "x = np.linspace(0, 10, 3000)\n",
    "y0 = .5\n",
    "\n",
    "eta = integrator(x, y0, ODEF, phi_euler)\n",
    "y_exact = (x + 1)**2 - 0.5*np.exp(x)\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "gs = gridspec.GridSpec(3, 1, height_ratios=[2, 1, 1])\n",
    "\n",
    "# Above plots\n",
    "ax0 = plt.subplot(gs[0])\n",
    "ax0.plot(x, eta, label=\"Euler method solution\")\n",
    "ax0.plot(x, y_exact, label=\"exact solution\")\n",
    "ax0.legend()\n",
    "plt.title(\"Solutions to $y' = y - x^2 + 1.0$\")\n",
    "plt.xlabel(\"x\")\n",
    "ax0.set_ylabel(\"y(x)\")\n",
    "\n",
    "# Relative absolute residue plots\n",
    "ax1 = plt.subplot(gs[1])\n",
    "ax1.plot(x, np.abs(eta/y_exact - 1), label=\"relative absolute residue\")\n",
    "ax1.legend()\n",
    "ax1.set_ylabel(\"$|\\eta(x)/y(x) - 1|$\")\n",
    "\n",
    "# Absolute residue plots\n",
    "ax2 = plt.subplot(gs[2])\n",
    "ax2.plot(x, np.abs(eta - y_exact), label=\"absolute residue\")\n",
    "ax2.legend()\n",
    "ax2.set_ylabel(\"$|\\eta(x) - y(x)|$\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "print(\"\"\"\n",
    "    For the given initial value problem y' = (x + 1)^2 + 1.0, taking x from 0 to 10 in 3000 steps,\n",
    "    we find a large deviation between exact and estimate solutions at the point that y' = 0.\n",
    "    We could fix this by taking some other measure, like the absolute difference, but as\n",
    "    the solutions decrease so hard, that will suck later on.\n",
    "    \"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "aea0c6713e9bd3871c2e93e65e93d903",
     "grade": true,
     "grade_id": "cell-715c304f22415d45",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def test_integrator():\n",
    "    def ODEF(x, y):\n",
    "        return y - x**2 + 1.0\n",
    "    \n",
    "    x = np.linspace(0, 3, 100000)\n",
    "    y0 = .5\n",
    "\n",
    "    eta = integrator(x, y0, ODEF, phi_euler)\n",
    "    y_exact = (x + 1)**2 - 0.5*np.exp(x)\n",
    "    \n",
    "    # Assert the relative error of the approximation is lower than TOL.\n",
    "    TOL = 1e-4\n",
    "    assert np.all(np.abs(eta/y_exact - 1) < TOL)\n",
    "    \n",
    "test_integrator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "8048460a3a4d5e94b9e9f086d90ff0d0",
     "grade": false,
     "grade_id": "cell-99399435c164c96d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 3\n",
    "\n",
    "Extend the set of possible single-step integration schemes to the modified Euler (Collatz), Heun, and 4th-order Runge-Kutta approaches by implementing the functions $\\text{phi_euler_modified(x, y, f, i)}$, $\\text{phi_heun(x, y, f, i)}$, and $\\text{phi_rk4(x, y, f, i)}$. These can be used in your $\\text{integrator}$ function instead of $\\text{phi_euler}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "dc82ac823cbd262952a6990776e26d09",
     "grade": true,
     "grade_id": "cell-ea865d2efbc574d6",
     "locked": false,
     "points": 9,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def phi_euler_modified(x, y, f, i):\n",
    "    h = x[i]-x[i-1]\n",
    "    phi = f(x[i-1]+0.5*h, y[i-1]+0.5*h*f(x[i-1],y[i-1]))\n",
    "    return phi\n",
    "    \n",
    "def phi_heun(x, y, f, i):\n",
    "    h = x[i]-x[i-1]\n",
    "    phi = 0.5*(f(x[i-1],y[i-1])+f(x[i-1]+h,y[i-1]+h*f(x[i-1],y[i-1])))\n",
    "    return phi\n",
    "    \n",
    "def phi_rk4(x, y, f, i):\n",
    "    h = x[i]-x[i-1]\n",
    "    k1 = f(x[i-1],y[i-1])\n",
    "    k2 = f(x[i-1]+0.5*h,y[i-1]+0.5*h*k1)\n",
    "    k3 = f(x[i-1]+0.5*h,y[i-1]+0.5*h*k2)\n",
    "    k4 = f(x[i-1]+h,y[i-1]+h*k3)\n",
    "    phi = (1/6)*(k1+2*k2+2*k3+k4)\n",
    "    return phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "996fa4ae58be6c15aed465d8d7e5618e",
     "grade": false,
     "grade_id": "cell-a783908a0db32a24",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 4\n",
    "\n",
    "Add the possibility to also handle the following multi-step integration schemes:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\Phi_{AB3}(x, y, f, i) = \n",
    "    \\frac{1}{12} \\left[ \n",
    "        23 f(x_i, y_i) \n",
    "        - 16 f(x_{i-1}, y_{i-1})\n",
    "        + 5 f(x_{i-2}, y_{i-2})\n",
    "    \\right]\n",
    "\\end{align*}\n",
    "\n",
    "and\n",
    "\n",
    "\\begin{align*}\n",
    "\\Phi_{AB4}(x, y, f, i) = \n",
    "    \\frac{1}{24} \\left[ \n",
    "        55 f(x_i, y_i) \n",
    "        - 59 f(x_{i-1}, y_{i-1})\n",
    "        + 37 f(x_{i-2}, y_{i-2})\n",
    "        - 9 f(x_{i-3}, y_{i-3})\n",
    "    \\right]\n",
    "\\end{align*}\n",
    "        \n",
    "In these cases the initial condition $\\text{y0}$ must be an array consisting of several initial values corresponding to $y_0$, $y_1$, $y_2$ (and $y_3$). Use the Runga-Kutta method to calculate all of the initial values before you start the AB3 and AB4 integrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "d301f6514b2bc50a77d929518f6928df",
     "grade": true,
     "grade_id": "cell-0bd05d93759bd652",
     "locked": false,
     "points": 6,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def phi_ab3(x, y, f, i):\n",
    "    # The first three y values are to be calculated using the Runga-Kutta method.\n",
    "    if i < 3:\n",
    "        return phi_rk4(x, y, f, i)\n",
    "    \n",
    "    phi = (1/12)*( 23*f(x[i - 1], y[i - 1]) - 16*f(x[i - 2], y[i - 2]) + 5*f(x[i - 3], y[i - 3]) )\n",
    "    return phi\n",
    "\n",
    "# TODO: AB4 looks a bit off in the plots in task 5.\n",
    "def phi_ab4(x, y, f, i):\n",
    "    # The first three y values are to be calculated using the Runga-Kutta method.\n",
    "    if i < 4:\n",
    "        return phi_rk4(x, y, f, i)\n",
    "    \n",
    "    phi = (1/24)*( 55*f(x[i - 1], y[i - 1]) - 59*f(x[i - 2], y[i - 2]) + 37*f(x[i - 3], y[i - 3]) - 9*f(x[i - 4], x[i - 4]) )\n",
    "    return phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "c621385b0e71eb6fd2e4d6041268bd15",
     "grade": false,
     "grade_id": "cell-3fac52fc46f0a497",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 5\n",
    "\n",
    "Plot the absolute errors $\\delta(x) = |\\tilde{y}(x) - y(x)|$ with a $y$-log scale for all approaches with $0 \\leq x \\leq 2$ and a step size of $h=0.02$ for the ODE from task 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "22f375f78eeb46befb192c08db0f48e9",
     "grade": true,
     "grade_id": "cell-e61d969c23914a5f",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAH0CAYAAAD/taEiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACicklEQVR4nOzdd3yV9fn/8dcn55zkZE/2CpuwRxCU4aAyFHBvUWtdVVvqt1ptawu21trqz7rrqBYXuIuIGwUVUCDBoECQvVf2HifnfH5/nBDDCISQcDLez8cjD3Lu+z73fd137oQrn1z39THWWkRERERE5PgFBToAEREREZGmSsm0iIiIiEgdKZkWEREREakjJdMiIiIiInWkZFpEREREpI6UTIuIiIiI1JGSaRERERGROlIyLSIijZ4x5hRjzDfGmK+MMXOMMa5AxyQiAkqmRUSkadgBnGWtHQtsBc4LbDgiIn5KpkWkWTPGbDXG/OwE3r/GGHNG/UUUWMaYycaY8EDHcbystXustSWVL8sBXyDjERE5QMm0iDQJxpjRxpilxpg8Y0y2MWaJMWZ4PR/jsMTbWtvPWruoPo9zvDHUs0cBZwPuv0EZY7oA44H362FftxtjUowxZcaYWbXYPs4Y8z9jTJExZpsx5soTjUFEmr4m+wNVRFoOY0wUMB/4JfAmEAyMAcoCGVdTY4zpA3xtrc0LdCx1UXkfvAJcZ631HGPbmQDW2plH2Ww3cD8wAQitRQhP4R8VbwMMBj4wxqyy1q6pxXtFpJnSyLSINAW9AKy1c6y1XmttibX2U2vt9wDGmCRjzCJjTG5lWcbUmnZkjLHGmB7VXs8yxtxvjHkF6Ay8b4wpNMb8rnJ91UjxsY5Tue2dxpjvK0fQ3zDGuCvX3W2M2WWMKTDG/GiMGXeE2GqKoVbnZ4yJMMZ4jTHtqi3rb4zZY4yJBM4Cnq3NBT/K9TvWMU5k3/80xsyt9vohY8znxphgY4wTeB24z1r744kc5wBr7bvW2rlAVi1iCwcuAv5krS201i4G5gHT6iMWEWm6lEyLSFOwHvAaY14yxkwyxsQeWFHZ1eF94FOgNfAr4DVjTO/jOYC1dhqwHZhirY2w1v6z+vrjOM6lwESgKzAQuK5ym9uB4dbaSPwjoVtrE8PxnJ+1thBYBwyttvhB4AFrbQEQaa399tD3GWPmVybqR/qYf5zHOBH/AM40xgwxxtyC/zpeaK0tB64ARgB/qvzF4rITPNbx6gVUWGvXV1u2Cuh3kuMQkUZGybSINHrW2nxgNGCB54EMY8w8Y0wbYCQQATxorS231n6BvyTkinoOo7bHedxau9tam40/CR4MeIEQoK8xxmWt3Wqt3VTPxz1gBZWJrjFmLNCXytFoa+0/jvQGa+1ka21MDR+Tj+cYJ8JamwX8C3gJ+D1wzoGSFGvtK9baeGvtGZUfb5zo8Y5TBJB/yLI84IRG40Wk6VMyLSJNgrU23Vp7nbW2I9AfaI//Ybr2wA5rbfXuDtuADvUcQm2Ps7fa58VAhLV2I/AbYCaw3xjzujGmfT0f94CqRBf4J/6yhPJaHqu2anWMyhFkW8PH4hr2/R0wAPi9tXbH8QRVfYQduAe4p6YR9jooBKIOWRYFnOhovIg0cUqmRaTJsdauA2bhT6p3A52MMdV/nnUGdtXw9mIgrNrrttV3fZTDHu9xDo15trV2NNCl8jhHHCU+QgzHe9wVwFBjzEWAG5h9rNiMMR9V1mgf6eOjuh6jcgTZ1PAx+ghxDAD+jX9k+vpjxX2E41WNsOMvPXnwGCPsx2M94DTG9Ky2bBCghw9FWjgl0yLS6Blj+hhjfmuM6Vj5uhP+ModvgWX4E+TfGWNcxt8Tegr+h9WOJA240hjjMMZMBE6vtm4f0K2G9x3vcarH39sYc5YxJgQoBUqouU/yoTEc73FX4f8F4f/hH9092i8IAFhrJ1XWaB/pY1J9HONYjDEd8JfF3ALcCgwwDdzf2xjjrHxA1AE4jDHuygcdD2OtLQLeBf5ijAk3xozCP3HMKw0Zo4g0fkqmRaQpKMD/8NkyY0wR/iR6NfDbyvKCKcAkIBN4GrimcvT6SKZXbp8LXAXMrbbu78C9lWUBd1Z/Ux2OU10I/pHSTPxlIK3x1wQfyUExHO9xrbVlwA/AVmvtkUaVT1h9H8P4W959CDxirZ1nrS0GHgL+dqL7PoZ78f9icw9wdeXn91aL6yNjzB+qbX8r/hZ6+4E5wC/VFk9ETD0MKIiISCNhjAkGNgKXHqlzR1M5hohIU6GRaRGR5mUGsKSBk9yTcQwRkSZBybSISDNgjBlqjMkDxuLvRd0kjyEi0tSozENEREREpI40Mi0iIiIiUkdKpkVERERE6qjRJNPGmG7GmBeMMW8HOhYRERERkdpo0JppY8yLwGRgv7W2f7XlE4HH8DfK/4+19sFq69621l5cm/0nJCTYxMTE+g1aREREROQQqampmdbaVocuP+JMT/VoFvAk8PKBBcYYB/AUcDawE1hhjJlnrV17vDtPTEwkJSWlnkIVERERETkyY8y2Iy1v0DIPa+1XQPYhi08BNlprN1fO7PU6/ilZRURERESalEDUTHcAdlR7vRPoYIyJN8Y8AwwxxtQ0zS7GmJuMMSnGmJSMjIyGjlVEREREpEYNXeZRa9baLOCWWmz3HPAcQHJysppki4iIiEjABCKZ3gV0qva6Y+WyeuHxeNi5cyelpaX1tUuRKm63m44dO+JyuQIdioiIiDQCgUimVwA9jTFd8SfRlwNX1tfOd+7cSWRkJImJiRhj6mu3IlhrycrKYufOnXTt2jXQ4YiIiEgj0KA108aYOcA3QG9jzE5jzC+stRXA7cAnQDrwprV2TX0ds7S0lPj4eCXSUu+MMcTHx+uvHiIiIlKlQUemrbVX1LD8Q+DDhjquEmlpKLq3REREpLpGMwOi1M7cuXNZu/a4W3IfZtGiRUyePPmo2+Tm5vL0009Xvd69ezcXX1yr+XREREREWgQl001MfSXTtXFoMt2+fXvefluzvYuIiEhgNOTM3XWlZLoBvPrqq5xyyikMHjyYm2++Ga/Xy4oVKxg4cCClpaUUFRXRr18/Vq9eTWFhIePGjWPo0KEMGDCA9957r2o/L7/8MgMHDmTQoEFMmzaNpUuXMm/ePO666y4GDx7Mpk2bDjruW2+9Rf/+/Rk0aBBjx44F/DXkP//5zxkwYABDhgxh4cKFh8U7c+ZMHn744arX/fv3Z+vWrdxzzz1s2rSJwYMHc9ddd7F161b69+9/1P3OmjWLCy+8kIkTJ9KzZ09+97vf1fv1FRERkZbB+izluwsp+HonmbPWsP+J7wId0mEaTZ/p42GMmQJM6dGjx1G3u+/9NazdnV+vx+7bPooZU/rVuD49PZ033niDJUuW4HK5uPXWW3nttde45pprmDp1Kvfeey8lJSVcffXV9O/fn4qKCv73v/8RFRVFZmYmI0eOZOrUqaxdu5b777+fpUuXkpCQQHZ2NnFxcUydOpXJkycfsdziL3/5C5988gkdOnQgNzcXgKeeegpjDD/88APr1q1j/PjxrF+/vlbn+uCDD7J69WrS0tIA2Lp1a9W6o+03LS2N7777jpCQEHr37s2vfvUrOnXqdIQjiIiIiPzEWkvF/mLKNuVRtimXsi15+IorAHAmhBLSPRpb4cM4G894cJNMpq217wPvJycn3xjoWA71+eefk5qayvDhwwEoKSmhdevWAPz5z39m+PDhuN1uHn/8ccB/0/zhD3/gq6++IigoiF27drFv3z6++OILLrnkEhISEgCIi4s75rFHjRrFddddx6WXXsqFF14IwOLFi/nVr34FQJ8+fejSpUutk+mjOdp+x40bR3R0NAB9+/Zl27ZtSqZFRETkMNZavFmllG7KpWyzP4H2FXoAcMSE4E6KJ6R7NCHdY3BGhwQ42iNrksl0bR1tBLmhWGu59tpr+fvf/37YuqysLAoLC/F4PJSWlhIeHs5rr71GRkYGqampuFwuEhMT69x67ZlnnmHZsmV88MEHDBs2jNTU1Fq9z+l04vP5ql6faOu3kJCfbnaHw0FFRcUJ7U9ERESaj4rc0p9Gnjfl4c0rAyAoMhh3jxhCuscQ0i0aZ3xogCOtncYzRt5MjBs3jrfffpv9+/cDkJ2dzbZt2wC4+eab+etf/8pVV13F3XffDUBeXh6tW7fG5XKxcOHCqm3POuss3nrrLbKysqr2AxAZGUlBQcERj71p0yZGjBjBX/7yF1q1asWOHTsYM2YMr732GgDr169n+/bt9O7d+6D3JSYmsnLlSgBWrlzJli1bjnms2uxXRERExFtQTnHafnLe3cCeh1aw98EV5Ly1ntIfswnuHEnM+d1p89thtPvDKcRd3ofw4W2bTCINzXxkOhD69u3L/fffz/jx4/H5fLhcLp566im+/PJLXC4XV155JV6vl9NOO40vvviCq666iilTpjBgwACSk5Pp06cPAP369eOPf/wjp59+Og6HgyFDhjBr1iwuv/xybrzxRh5//HHefvttunfvXnXsu+66iw0bNmCtZdy4cQwaNIg+ffrwy1/+kgEDBuB0Opk1a9ZBI8cAF110ES+//DL9+vVjxIgR9OrVC4D4+HhGjRpF//79mTRpErfddlvVe2699dZj7ldERERaHl+xh7LNef7SjU15VOwvBsC4HYR0jSbi1PaEdI/B1SYME9T0528wjbHFSG0lJyfblJSUg5alp6eTlJQUoIikJdA9JiIi8hNfWQVlW/OryjY8uwvBgnEFEdw1mpBu0bh7xOBqH9Gkk2djTKq1NvnQ5RqZFhEREZFasx4fZdt/Sp7LdxSAz4LDENIliqifdSGkezTBHSMbVdeNhqJkWkRERERqZL0+yncWVibPuZRty4cKC0EQ3DGSyNM7+jtudInCuByBDvekUzItIiIiIlWsz+LZW+RPnDfmUrYlH1vuBcDVLpyIke0J6RFDSGIUQW6lkroCIiIiIi2YtZaKzJKqso2yTbkHTZQSNrS1f+S5WwyOcFeAo218mmQyXdsZEEVERETkcBW5ZT+NPG/KxZtfDoAjOhh3nzj/yHMjniilMWmSyXRjngFRREREpLHxFnl+qnnelEdFZgkAQeFOQrrFENIjBnf3GBzxboxpuh03AqFJJtONncPhYMCAAVWvL7/8cu65554at581axYpKSk8+eSTJyM8ERERaeZ8ZV7KtuRVjTx79hQBYEL8vZ7DR7QjpEfz6fUcSEqmG0BoaChpaWkNtv+KigqcTn3pRERExM9W+Cjfnk/pxkPa1TkNIZ2jiBrfhZDuMf52dQ4lz/VJGdlJlJiYSEpKCgkJCaSkpHDnnXeyaNGig7bJyMjglltuYfv27QA8+uijjBo1ipkzZ7Jp0yY2b95M586dmTNnTgDOQERERBoD67N4dvvb1ZVuzKV8az7W4wMDro6RRI6tbFeX2DLb1Z1MzTuZ/uge2PtD/e6z7QCY9OBRNykpKWHw4MFVr3//+99z2WWX1Wr306dP54477mD06NFs376dCRMmkJ6eDsDatWtZvHgxoaFNZ756EREROXEHddzYkEvp5jxsSWXHjdZhhA9vS0j3GEK6RRMU2rzTu8ZGV7sBnEiZx4IFC1i7dm3V6/z8fAoLCwGYOnWqEmkREZEWwptfTumBjhsbc/DmHei4EUJo33jclR03HFHBAY60ZWveyfQxRpBPNqfTic/nA6C0tPSI2/h8Pr799lvcbvdh68LDwxs0PhEREQkcX2lFVZ/n0o25VOwvBiAozOkfde7u77rhVMeNRqV5J9ONTGJiIqmpqUyaNIl33nnniNuMHz+eJ554grvuuguAtLS0g0pGREREpHmwFT7KtuVXjjznUr6zACwYVxDBXaMJH9bG33GjXbg6bjRiSqYbwKE10xMnTuTBBx9kxowZ/OIXv+BPf/oTZ5xxxhHf+/jjj3PbbbcxcOBAKioqGDt2LM8888zJCVxEREQajPVZPHuKKNuYc/BDg0EQ3DGSyDM74e4RQ3DnKIwzKNDhSi0Za22gY6iz5ORkm5KSctCy9PR0kpKSAhSRtAS6x0REpLYqskr87eoq+z1XTdPdOsxf89yj8qFBt8Y3GztjTKq1NvnQ5U3yK6fpxEVERKQxqpppcKO/7tmb7X9GyhH10zTd7h4xOKI0TXdz0SSTaU0nLiIiIo2B9Xgp21I5WcrGHDy7K2cadDsI6RZD5OgO/ocGW4XqocFmqkkm0yIiIiKBYH0Wz67CquS5bFs+VFhwGII7RxF1dhdCesYQ3EEzDbYUSqZFREREjqKq7nlDDqWbfposxdU2nIiR7XH3jCG4azRBwZppsCVSMi0iIiJSja/YU/XQ4EF1z9HB/slSelZOlhKpyVJEybSIiIi0cFX9njfkUroxB8+uQn+/5xAHId2i/XXPPWNwJqjuWQ6nZFpERERaFGstFfuKKd2QQ+mGXMq35P3U77lTFFHjOhPSM5bgjqp7lmNTMt3IJSYmkpKSQkJCAqeddhpLly4F4K677uLDDz/knHPOoXv37oSFhXHNNdfUer8REREUFhYettzhcDBgwICq15dffjn33HNPjfuZNWsWKSkpPPnkk8dxVgfbs2cPN954I/Pnzwdg+fLl3Hnnnezbt4+wsDCGDRvG448/TlhY2DFjmDlzJhEREdx55501Hm/u3Ln06tWLvn371rjNk08+SVhYGNdff32dz0tERBoPb345pRtyKks3cvAVeABwtgolLLkN7p6x6vcsdaI7pgk5kEgDPPfcc2RnZ+Nw1O/DDqGhoaSlpdXrPqurqKjA6Tz4tnvkkUe48UZ/l8N9+/ZxySWX8Prrr3PqqacC8Pbbb1NQUFBjMn285s6dy+TJk4+aTF9//fWMGjVKybSISBPlK/dStiXPX7qxIYeKfcUABIU7CekR658wpWcszhj1e5YT06yT6X8s/wfrstfV6z77xPXh7lPurnH91q1bmThxIiNHjmTp0qUMHz6cn//858yYMYP9+/fz2muvccopp5Cdnc3111/P5s2bCQsL47nnnmPgwIFkZWVxxRVXsGvXLk499VSqz1B5YDR56tSpFBYWMmzYMH7/+9+Tnp5eNRq7adMmbrvtNjIyMggLC+P555+nT58+bNmyhSuvvJLCwkLOO++84z7v6iPkKSkp3HnnnSxatOigbTIyMrjlllvYvn07AI8++iijRo1i5syZbNq0ic2bN9O5c2fmzJlz0Pveeecd7r//fgCeeuoprr322qpEGuDiiy8GqPGa1eT555/nueeeo7y8nB49evDKK6+QlpbGvHnz+PLLL7n//vt56qmnuO2226re88MPP7B582a6dOlCYmIiy5cv55RTTjnu6yUiIifXgam6SzfkULYhh7Kt+eC14DSEJEYTPrQ1IT1icbULxwSpdEPqjyZ+bwAbN27kt7/9LevWrWPdunXMnj2bxYsX8/DDD/PAAw8AMGPGDIYMGcL333/PAw88UFWicd999zF69GjWrFnDBRdcUJWYVjdv3ryqEeTLLrvsoHU33XQTTzzxBKmpqTz88MPceuutAEyfPp1f/vKX/PDDD7Rr167G2EtKShg8eHDVxxtvvFHr854+fTp33HEHK1as4J133uGGG26oWrd27VoWLFhwWCK9ZcsWYmNjCQnxjwysXr2aYcOGHXH/NV2zmlx44YWsWLGCVatWkZSUxAsvvMBpp53G1KlTeeihh0hLS2PUqFGkpaWRlpbGjTfeyEUXXUSXLl0ASE5O5uuvv671+YuIyMnlzSujKGUfWXPWsedv37L/ie/I/3grvqIKIk5rT8L1/ekw41Ra3TCAyNM7EdwhQom01LtmPTJ9tBHkhtS1a9equuN+/foxbtw4jDEMGDCArVu3ArB48WLeeecdAM466yyysrLIz8/nq6++4t133wXg3HPPJTY2ttbHLSwsZOnSpVxyySVVy8rKygBYsmRJ1fGmTZvG3Xcf+dqcSJnHggULWLt2bdXr/Pz8qrrsqVOnEhoaeth79uzZQ6tWrWq1/5quWU1Wr17NvffeS25uLoWFhUyYMKHGbZcsWcLzzz/P4sWLq5a1bt2adevq9y8bIiJSd1WlG+v9Dw5W7K8s3Yhw4e4VR0gvf/mGWtbJydQkk2ljzBRgSo8ePQIdyhEdGGUFCAoKqnodFBRERUVFgx3X5/MRExNTYzJ8Iu18nE4nPp8PgNLS0hqP/+233+J2uw9bFx4efsT3hIaGHrS/fv36kZqaWqdSlENdd911zJ07l0GDBjFr1qzDylIO2LNnD7/4xS+YN28eERERVctLS0uP+AuAiIicHNZn8ewt8k+Wsr566UYQIV2jCE9ug7tXLM42YWpZJwHTJMs8rLXvW2tvio6ODnQodTZmzBhee+01ABYtWkRCQgJRUVGMHTuW2bNnA/DRRx+Rk5NT631GRUXRtWtX3nrrLcDf+mfVqlUAjBo1itdffx2g6rjHIzExkdTUVICq0eFDjR8/nieeeKLqdW1GuHv16lU1Wg9w++2389JLL7Fs2bKqZe+++y779u2r8ZrVpKCggHbt2uHxeA4658jISAoKCgDweDxccskl/OMf/6BXr14HvX/9+vX079//mOcgIiL1x1tQTtHKfWS/vo49Dyxj/+PfkfdRtdKNX/Snw4yRtPrFACLHdsTVNlyJtARUk0ymm4OZM2eSmprKwIEDueeee3jppZcAf13wV199Rb9+/Xj33Xfp3Lnzce33tdde44UXXmDQoEH069eP9957D4DHHnuMp556igEDBrBr164a339ozfSBtngzZsxg+vTpJCcn19hB5PHHHyclJYWBAwfSt29fnnnmmWPGGx4eTvfu3dm4cSMAbdq04fXXX+fOO++kd+/eJCUl8cknnxAZGVnjNavJX//6V0aMGMGoUaPo06dP1fLLL7+chx56iCFDhrB06VJSUlKYMWNG1Tnv3r0b8Jd+nH322cc8BxERqTvr8VG6MYfcD7ew79GV7PnbMnLeXE/phlzcPWKIvaQX7f4wgja/GUrMud1w94zFuDRttzQepnq3iKYmOTnZpqSkHLQsPT2dpKSkAEUkdfG///2P1NTUqo4ejcF3333HI488wiuvvHLYOt1jIiJ1Z62lIqPEX7axIYeyzZUTpjgMIV2i/HXPPdV1QxofY0yqtTb50OVNsmZampcLLriArKysQIdxkMzMTP76178GOgwRkWbBV+yhdFMuZetzKV2fgzfP/3C8MyGU8OFtCekVS0jXaIJCNOIsTY+SaWkUqrfRawxU3iEiUnfWZynfWeDvurE+h/IdBWDBuB24u8cQclYn3D1jccYd/sC6SFOjZFpEREROmDevjNLKrhulG3KxJRVgwNUxksgzO+HuFUtwpyiMQ6Ub0rwomRYREZHjZit8lG3N89c+r8/Bs7ey53NkMKF943H3iiGkRyyOcFeAIxVpWEqmRUREpFYqMv0PDpauz6FsU+5PDw4mRhE9KZGQXnG42qrns7QsSqYbQERERNXMfwCzZs0iJSWFJ598MoBRiYiIHB9fuZeyzXmU/pjtf3Awyz/JliPeTVjlhCkh3WL04KC0aEqmRUREBKhsW7e/2D/6/GMOZVvywGsxriBCuscQOaqDf8bBBM0OKydHRUUF6at+YO2yFLK3bsdbXMrtTz4U6LAOomT6JMvIyOCWW25h+/btADz66KOMGjWKmTNnEhERwZ133glA//79mT9/PgCTJk1i9OjRLF26lA4dOvDee+9pmmsREakXvtIKyjblUvqjv3zDm1vZtq51GBGntfePPidGY1ya503qn8fjYeW3y9mUmkbu9l1U5BZgy8qwPg/WluGzxUBFtXc48Xg8uFyNpxa/WSfTex94gLL0dfW6z5CkPrT9wx+Ous2BWQQPyM7OZurUqQBMnz6dO+64g9GjR7N9+3YmTJhAenr6Ufe3YcMG5syZw/PPP8+ll17KO++8w9VXX33C5yIiIi2PtRbP3mJKf8ymbH0OZVvzwWcxIQ5CesQQeVYn3L3icMaEBDpUaQY8Hg+rv1vFumUp5GzZgSc3Hw5KlosA7yHvCiHIhOIwwYQEBRPsKCfCVUhCaBaJCR5oZBMONutkOlBCQ0NJS0uren2gZhpgwYIFrF27tmpdfn7+QfXVR9K1a9eq5HzYsGFs3bq1vkMWEZFmzFdaQemG3KoE2ptfDoCrXTiRYzrg7h1LcJcojEOjz3L8dm7bzoovvmT/j5soy8rBlpRivdWT5YpD3hGCozJZdjtcBAeVExVcQOuwLHq28dC+Y1uIbg0xnSC600//RncEZ+P7Ja9ZJ9PHGkEOBJ/Px7fffovbfXCjeqfTic/nq3pdWlpa9XlIyE83jsPhoKSkpOEDFRGRJqv66HPpj9mUbyvwjz67Hbh7xuLuFYu7dyyOqMaXmEjjU1ZaRuripWxIWUnBzr14C4qwnnLwleOjBGsPzUtcBJkwHMZFSFAMwQ5/stwqNJNebSro0LktxLSFmC4Q0/mnj8h2ENT0HmZt1sl0YzR+/HieeOIJ7rrrLgDS0tIYPHgwiYmJVTXSK1euZMuWLYEMU0REmhhfWQVlGyprn3/MPnj0eWxH/+hzZ02aIkeWmZHBtwsWsmf1Okr2Z1WOLpfXUIphMCYcByEEm3CCnSGEBxcRH5JD94RiEjvH4IpvDbFd/AnzgX8j24Gj+aWeTfKMjDFTgCk9evQIdCjH7fHHH+e2225j4MCBVFRUMHbsWJ555hkuuugiXn75Zfr168eIESPo1atXoEMVEZFGrKrzRmXyXLY13995I8SBu2cM7t5xGn2Wg+zctp3lCxay78eNeDJzKx/0K8dnS7C26JCtnQSZcH8pRpALt7OU6JB8OkRm06+zi4h20RDbCWITISax8t9OjbIMo6EZ28iKuI9HcnKyPVCLfEB6ejpJSUkBikhaAt1jIhIovnLvT5031mVXdd5wtQ0jpHdcZecN1T63ZDu3bWfZZ1+wf91GPNmVCbO3HB/FRyjHcFfWLjtxBVUQ5iomwZ1Nt/hiuneNwZXQA+ISIbar/yOuK4TGQgudlMcYk2qtTT50eZMcmRYREWkpKrJKKF2XTcmPOZRtzoUKiwkOIqRHLJFndsLdW503WprMjAyWfPwZe1evozQzB0pLKxPmEqwtPmhbY0IJIpRgE0Gw0014cAFtQrPp2bacLl07QlzbykS5mz9ZjukCLncNR5YjUTItIiLSiNgKH2Vb8ihdl03pjzlUZPpHE52tQokY2R53n8q+z06NPjdnxSUlfLNgIVtWfEfx3gxscYk/Ya4qyaheWeDGYcIINuG4nCFEVCbMfTpW0LFLJ4jr6E+U47r5P6I7Ncva5UDRlRQREQkwb34ZpetyKPkxm7INudhyLzgNId1iiDi1He4+cTjjNVlXc5T+w1pWLviC3E3bqSgowJaXY20pPlvIwQ/9uXCYcFwmFJfDRbirmDZhmfRuV06Xbp0gvoM/UY7vDnHd/W3kmmBnjKZIybSIiMhJZn2W8p0F/tHnddl4dvsf/nJEBxM2pBXu3nGE9IghKFjJUHNQVFzMko8/Y1vqKkr2ZWJLKkeZKcLa0mpbBlU+9BdCiCOGUEcJ8aE59GhdSM+uCbjatIH4Hv6PuG7+kgyNMAecvgIiIiInga+0gtL1OVXlG74iDxgI7hxF1IRE3H3icLUNw7TQh7uag13btrP4g0/J/HEjFbn52PKyGkaZQ3CYcIJNBCGuEKJD8ukYlc2A7qFEtu/gH11O6PlT0uzSXyUaMyXTIiIiDcSTUexPntOzq6btDgpzEtIrltA+cYT0jMUR7gp0mHKcVq9MY+Xni8jbvANfQWFlLXPxIe3lDEEmAocJwe2IIcxZTEJoFn3al9Gte0dMqy6VCXNP/79h8S22S0ZTp2S6gXi9XpKTk+nQoUPVZCyPPvooN910E2FhYQBEREQccyrxukpMTCQlJYWEhIQ6vf+6667jyy+/JDo6GmstjzzyCOPGjavnKGvnwHTsTz75JD6fj5///Oc4HA5eeOGFI47g5ObmMnv2bG699VYAtm7dytKlS7nyyitPdugi0sLYCh9lW/MoTfeXb1Rk+f+E72wTRuTYDrj7xBHcSROnNAUVFRWkfrOMtYuWULBjN7aoGOstO0JphpMgE4HLhBHsDCYyOJ8OkVkM6h5MTKeO/kQ5oZf/IzYRnMGBOiVpIEqmG8hjjz1GUlIS+fn5VcseffRRrr766qpkur5UVFTgdJ7Yl9Lr9eJwHFyb99BDD3HxxRezcOFCbrrpJjZs2HBCxzhR1lpuueUWPB4P//3vf2v8U2hubi5PP/30Qcn07NmzlUyLSIPwFpZX9X0uXZ+DLfM/POjuHkPEaH8C7YxVq7HGyuPxkPL1UtZ+tZSiXXuxxf6k2WuLgLJqWwbjMBEEm0hCXCHEhuTQJT6PQb1iCW7X7qeEuVVviGijUeYWRMl0A9i5cycffPABf/zjH3nkkUcA/8yHu3fv5swzzyQhIYGFCxcC8Mc//pH58+cTGhrKe++9R5s2bQ7aV3Z2Ntdffz2bN28mLCyM5557joEDBzJz5kw2bdrE5s2b6dy5M08++SRXXHEFu3bt4tRTT6X6ZDyvvvoqjz/+OOXl5YwYMYKnn34ah8NBREQEN998MwsWLOCpp55i9OjRRzyfU089lV27dgEHjxIDTJ48mTvvvJMzzjiDiIgIpk+fftj5bNq0iauuuoqioiLOO+88Hn300aoR+Yceeog333yTsrIyLrjgAu67774ar+uvf/1rsrKyeOONNwgKCmLmzJlERERw5513AtC/f3/mz5/PPffcw6ZNmxg8eDBnn302X3/9Nenp6QwePJhrr72WCy64gGnTplFU5P9z3JNPPslpp5123F9nEWmZrLVU7CumJD2b0vQsyncUgIWgyGDCBrXC3UcPDzZG/pHmb1mzcAmFO3ZXS5oLgfJqW4bgNOG4gyJxO4OJc2fRo00RfXq1xdUuERJ6+xPmhJ7gjg7Q2Uhj0qyT6a/fXE/mjvoto0joFMGYS48+1fdvfvMb/vnPf1JQUFC17Ne//jWPPPIICxcurCq9KCoqYuTIkfztb3/jd7/7Hc8//zz33nvvQfuaMWMGQ4YMYe7cuXzxxRdcc801pKWlAbB27VoWL15MaGgov/71rxk9ejR//vOf+eCDD3jhhRcA/2x9b7zxBkuWLMHlcnHrrbfy2muvcc0111BUVMSIESP4f//v/x31fD7++GPOP//8Y16bms5n+vTpTJ8+nSuuuIJnnnmmavtPP/2UDRs2sHz5cqy1TJ06la+++oqxY8cetu/Zs2eTlJTEokWLjjkK/+CDD7J69eqq67Ro0SIefvjhqnKb4uJiPvvsM9xuNxs2bOCKK67g0Jk0RUSqsxU+yjbnUZKe5Z95MKdy5sEOEUSN6+x/eLB9BCZIo5GNwdq0Vaz4+HPytuzAFhZVS5qrjzQfSJqjcDvLiA/LolfbUnr2bI+rXU9o1eenpFkPAMpRNOtkOhDmz59P69atGTZsGIsWLTrqtsHBwUyePBmAYcOG8dlnnx22zeLFi3nnnXcAOOuss8jKyqoqHZk6dSqhof5v8K+++op3330XgHPPPZfY2FgAPv/8c1JTUxk+fDgAJSUltG7dGgCHw8FFF11UY3x33XUXf/jDH9i5cyfffPPNMc+9pvP55ptvmDt3LgBXXnll1Ujyp59+yqeffsqQIUMAKCwsZMOGDUdMpocOHcq6detYvnw5o0aNOmYsR+PxeLj99ttJS0vD4XCwfv36E9qfiDRP3iJP5cODWZSu9/d+Nq4gQnrEEHlmJ0L7xOGI0syDgbR7xy6+fm8+Ges24s0rwFaU4bNFh8wC6C/PCAmKItRZSnxoFr3aldCzZwdc7XtVJs19/EmzU19POX7NOpk+1ghyQ1iyZAnz5s3jww8/pLS0lPz8fK6++mpeffXVw7Z1uVxVdb8Oh4OKiorjOlZ4ePgxt7HWcu211/L3v//9sHVut/uwOunqDtRMP/HEE1x//fWkpqbidDrx+XxV25SW/vQQxvGej7WW3//+99x8880HLX/qqad4/vnnAfjwww8B6NOnD3/5y1+49NJL+eSTT+jXr99RYzmaf/3rX7Rp04ZVq1bh8/lwu1XLKCKV5RsZJZSmZ1GyNpvy7fn+8o2oyt7PfeJw94jBuFS+cbIVFRfz5fsfsW1FGuUZWdjy0sqZAAuqbeXAYSIJNuGEOF3EunPo3rqA/n3a4GrfDVr3rUyae2m6bKlXzTqZDoS///3vVYnrgfKCA4l0ZGQkBQUFx9VhY8yYMbz22mv86U9/YtGiRSQkJBAVFXXYdmPHjmX27Nnce++9fPTRR+Tk5AAwbtw4zjvvPO644w5at25NdnY2BQUFdOnSpdYx3H777bz44ot88sknJCYm8vTTT+Pz+di1axfLly8/5vtHjhzJO++8w2WXXcbrr79etXzChAn86U9/4qqrriIiIoJdu3bhcrm47bbbuO222w7bz2mnnca///1vJk+ezJdffkliYmJV6cbKlSvZsmUL8NN1PuDQ13l5eXTs2JGgoCBeeuklvF4vItIyWa/9qftGelZV9w1Xhwgiz+pMaFIcrg4R6v18Eq1OTWPFR59RsG0ntrgYn7cUny0AfhqgCTKROI2bYIeT6JB8OsXkMrRvNGEdO0ObvtA6yV/bHBIRuBORFkPJ9El00003MXHiRNq3b1/1AOKxzJw5k+uvv56BAwcSFhbGSy+9dMTtZsyYwRVXXEG/fv047bTT6Ny5MwB9+/bl/vvvZ/z48fh8PlwuF0899dRxJdPGGO69917++c9/smDBArp27Urfvn1JSkpi6NChx3z/gS4mf/vb35g4cSLR0f4HNsaPH096ejqnnnoq4G8V+Oqrr1aVoRzJlClTyMzMZOLEiXz22We8/PLL9OvXjxEjRtCrl/8vEfHx8YwaNYr+/fszadIkHnjgARwOB4MGDeK6667j1ltv5aKLLuLll19m4sSJtRrhF5Hmo2rylLVZlPyYgy2pAIchpHsMEWM64k6KwxmtP/c3tIL8Ar743zx2pa3Gk5WL9ZQeVqJhjBsH4YQ4ogl3FdM2IoshPR207t4dWverTJz7QlhcAM9EWjpTvetDU5OcnGwPfXAsPT2dpKSkAEUkR1JcXExoaCjGGF5//XXmzJnDe++9F+iw6kz3mEjTU5FbSml6NiVrsyjbnAdeS1C4E3efeEKT/JOnBIWofKOhpH+/hm/nf0TBlh3YomJ8vpLK0eYDfxkMqhxtDsbtLCfOnU1Sh2L6JCUS1K4/tOnvT5pjOqvlnASMMSbVWpt86HKNTEuDS01N5fbbb8daS0xMDC+++GKgQxKRZs5ai2d3ESVrsyhNz8Kz298K09kqlIhRHQjtG0dw5yh136hnHo+HLz/8hA1LllG+LwvKS/HaooNmBjww2ux2xhDhKqRjdDbDkiKISuwCbfpB2/7+Eg3VNUsT0SSTaWPMFGBKjx49Ah2K1MKYMWNYtWpVoMMQkWbOeivb163NonRtNt68MjAQ3DmK6EldcfeNw9WqfifNasmyM7P5/J257PthLd7cAmxFKV5bwE89mw1BJhKXCSPE6SIuNJuk9iX0SkrA1XEQtBngT5w1wYk0cU0ymbbWvg+8n5ycfGOgYxERkcDxlVZQ+mOOP4H+MRtbWtm+rmcsUWf7+z87IjR984navHETi9+dR+6GbdiiosoyjXzgQEclJw4ThdsRRZirmHYRWQztHUzrHh39CXPbAdAqSaPN0iw1yWRaRERaLm9+OSXpWZSsyaJsU25V/XNovwRC+8YT0lOzD56IH1K/Y/n8TyjatgtbWozXV3xQC7qfyjRiiQrJp0tcDsMGtCa0cx9oO9CfPMd00WiztBhKpkVEpNHz7C/2jz6vqZy+G3DEu4k4rT2hfeMJ7qL657pYsfgb0j5eQMmuvdjSkiPUN4fjNKGEOJ3EunP8DwX2bY2r05DKxHmAOmlIi6dkWkREGh3rs5TvLPC3r1uTRUVGCQCujhFEje9CaL94nK3D1P/5OCz78mtWfbKQ0t378JWV4LOFB7Wh89c3hxLidJEQnsXARA89+nWDdoP8H637QrBqzkUOpWS6gcydO5cLLriA9PR0+vTpw9atW0lKSqJ3795YawkPD+e///0vvXv3Zvny5dx0002A/wn0mTNncsEFFwT4DERETq6qBwjXZFGyNgtffjkEQUi3GCJOa487KR5njPo/18aKxd+Q9tFnlOyqKXGOwmXCCXU5aRWexdCehk5JnX5KnFv1BocrgGcg0nQomW4gc+bMYfTo0cyZM4f77rsPgO7du5OWlgbAs88+ywMPPMBLL71E//79SUlJwel0smfPHgYNGsSUKVNwOvXlEZHmzVfupWx9jj+BTs/GllZgXEG4e8Xi7p9AaO9YgsKU1B3N6tQ0vp33AcXb9/hrnI+SOLcJz2RYbyftk7pAu8HQfjDE94Ag1ZiL1JWytQZQWFjI4sWLWbhwIVOmTKlKpqvLz88nNjYWgLCwn/5sVlpaqj9bikiz5iv2UJKe7X+AcEMO1uMjKMxJaN84Qvsn4O4Zg3EpuTuSzRs38eWb/6Ng0zZsSTFeXxHWFlat95dq+BPn1hFZJPdy0D4p0Z80txtcmTgHBSp8kWapWSfTC2c9x/5tm+t1n627dOPM62466jbvvfceEydOpFevXsTHx5Oamkp8fDybNm1i8ODBFBQUUFxczLJly6res2zZMq6//nq2bdvGK6+8olFpEWlWvPnllKzNomRNJmWb8sBncUQHE5bchtB+CYR0jcY4NJBQXWZGBp/NeYuMNRuwBYX4fEWV7ej8jInAady4nU5ah2cxpBd06dsF2g/xf8R314izyEmgjK0BzJkzh+nTpwNw+eWXM2fOHG6//faDyjzeeOMNbrrpJj7++GMARowYwZo1a0hPT+faa69l0qRJuN3qxykiTVdFdiklazIpWZ1F+fZ8sOBMCCVybAdC+yXg6hihv8RVKistY8G777H1m1QqcnLxVRTjs3kc6ONsjBuniSDMGUecO5tBiWX0GtwLOgz1J84JvcGh/9JFAqFZf+cdawS5IWRnZ/PFF1/www8/YIzB6/VijOG22247aLupU6fy85///LD3JyUlERERwerVq0lOPmz6dxGRRs2zv5iSHzIpWZNZNYW3q304UT/rQmh/deA4IHXpMlLe/5jS3Xux5SV4ffn8NHOgC6eJwu2MJTokj97tCxk4uAOuLsOh/VD/lNtOPYgp0lg062Q6EN5++22mTZvGs88+W7Xs9NNPZ8eOHQdtt3jxYrp37w7Ali1b6NSpE06nk23btrFu3ToSExNPZtgiInVircWzu4iS1ZmUrM6samEX3DmS6HO6+lvYxYcGOMrA2rd3L5+99ibZ6ZuwxYWH1DkbHCYatyOKCFcRXRKyOWVILGFdh0GHYf7OGiERAY1fRI5OyXQ9mzNnDnffffdByy666CL+/ve/V9VMW2sJDg7mP//5D+BPrB988EFcLhdBQUE8/fTTJCQkBCJ8EZFjstZSvqOgMoHOwptdCgZCukb7J1HpF48jqmWOnHo8Hha+/zEbv1qCJzPnCOUa/klQQl0O2kZkMWJAMK379PEnzh2GQUTrwJ6AiBw3Y60NdAx1lpycbFNSUg5alp6eTlJSUoAikpZA95i0RNZnKd+e7y/hWJ2JN68cHIaQ7jGE9U/A3TcOR0RwoMM86bZt3sIXs9+mYPM2fCXF+Gw+1pZWrnXiNFEEOwyxoTn071xK/2EDoEMydEyGuO7qrCHShBhjUq21h9XgamRaRESOyPosZVvyqkagfQXl4DS4e8YSNT6B0KS4FtUD2uPx8MV7H7Dxq2+oyMrB5z0w6uwflAoyUQSbSMJDHHSKzeHUoTGE9zjFnzi3GwSull3uItJcKZkWEZEq1msp25zrT6DXZOEr9PgnUekdS+iABNx94ggKaRn/dezZtZvPXn2d3B+3YEuK8B406ux/SDDUGUtcWDbDevnoPqg7dBzu/4hqF9DYReTkaRk/EUVEpEbW66NsU+UI9JpMfEUVmOAg3H3i/Al07ziCgpt/v+JvF37Jd/M/o3xfxmG1zkEmsmrUuUt8NiOHtyKs20h/4tx2gKbeFmnBlEyLiLRABxLo4u8zKF2bha+4AhPswJ0UR9iABNy9Y5v1LITFJSV89Mob7EpdhS0owOsrwNqiyrUOnCYatzOGuNBchvSooNfQHtDpFOh4CkS2CWjsItK4KJkWEWkhqkagK/tA+4orMCEOQpMqR6B7Nd8Eevf2HXzyyhvkb9yCLS3G68vjQF9nY0JxmnDCXA46xGQyalgronqPgE4joO1AcLa8BytFpPaUTIuINGNVNdCVXTgOjECH9o0jdECrygS6+XWU+G7ZCpa9+z5lu/djPcV4bS4/lWz4+zpHBhfQu0MBQ4d3x9VtlD95jukMmlRGRI6DkukGMnfuXC644ALS09Pp06cPW7duJSkpid69e2OtJTw8nP/+97/07t276j3bt2+nb9++zJw5kzvvvDOA0YtIU1bVheP7DH8CXVRZwtG3soSjmY1AV1RU8OUHH7NuwVdUZOfg8xbhs/mVa/0lG2HOWOLDshmeZOk6uC90HuGvd3ZHBzR2EWn6lEw3kDlz5jB69GjmzJnDfffdB0D37t1JS0sD4Nlnn+WBBx7gpZdeqnrP//3f/zFp0qRAhCsiTZz1Wcq35VP8fQYlP2T+1IUjKY6wga2aVQ10eWk5H73xJtu+ScXm51fWOxdXrg3GFRRFuDOG9lGZjEqOJDZpBHQe6S/Z0IOCIlLPlEw3gMLCQhYvXszChQuZMmVKVTJdXX5+PrGxsVWv586dS9euXQkPDz+ZoYpIE1Y1E+EqfwLtzS/3J9B94ggd2Hy6cBTkF/DBS7PZv2oNtqiQCl8eUAaAMWG4TBjhriASW+UwemR7gnuMgS6nQVw3lWyISINr1sl07vubKN9ddOwNj0Nw+3BipnQ/6jbvvfceEydOpFevXsTHx5Oamkp8fHzVdOIFBQUUFxezbNkywJ98/+Mf/+Czzz7j4Ycfrtd4RaR5sdbi2V3kH4H+PgNvThk4DO7ecUQPTMCdFE9QSNNOoLOzsvlw1qtkr1mPr6QIry8X8ABgTCQhQVFEhhTSu0Mew0b0wtV9DHQ+FSLbBjRuEWmZmmQybYyZAkzp0aNHoEM5ojlz5jB9+nQALr/8cubMmcPtt99+UJnHG2+8wU033cTHH3/MzJkzueOOO4iIiAhg1CLSmHn2F1O8KoOSVRlUZJZAkMHdM4aon3UhtF88Qe4m+eMcgMyMDD588RVy1m3Clh5InisA/6yCbkc00e58BnYtZeDIXtBllL9NXWhMIMMWEQHAWGsDHUOdJScn25SUlIOWpaenk5SUFKCIIDs7m44dO9KqVSuMMXi9XowxfPnll0yZMoXVq1cDUFJSQnx8PMXFxYwZM4YdO3YAkJubS1BQEH/5y1+4/fbbA3YeUrNA32PSclRkl1Yl0J69RWAgpFs0oQNbEdo/AUd406z/zczI4IP/vEzej5vxlRXh9eUAXsDfaSMkyEVMaC7Dennpfcpp/pKNjsMhOCywgYtIi2aMSbXWJh+6vOkOZTRSb7/9NtOmTePZZ5+tWnb66adXJcsHLF68mO7d/eUiX3/9ddXymTNnEhERoURapIXyFpT7SzjSMijfUQBAcOdIoqd0I2xAKxxRTa/ncXZmNvNffInc9I34Sg9NnmMIdcQSG5bD8CRDj+GDocto6DAUnCEBjVtEpDaUTNezOXPmcPfddx+07KKLLuLvf/97Vc20tZbg4GD+85//BChKEWlMfMUeStZkUbwqg7JNuWDB1S6cqImJhA1shTPOHegQj0tBXj7zXniJrB9+xFdaeFDZhqMyeY4Ly+GUfkF0Sx4CiWOg3WBNjiIiTZLKPESOk+4xqQ/W46UkPZvitAxKf8wGr8UR7yZsUCvCBrfG1brplDQUFRczf9ar7Fu5Gl9RwUEPDB4o24gNy2V4EvQ4ZQwkjoX2g9WmTkSaFJV5iIgEmPVayjblUpy2n5I1WdgyL0GRwUSMbOdPoDtGYJpAKzePx8PHb7zDlq+XYfPzDmpVd+CBwRh3PsP6VNBn5ODK5HmIRp5FpFlSMi0i0oCqekGnZVD8fYZ/MhW3g9ABCYQNbk1It2hMUONPoL/6+DO+n/8p3uwcvL78qklSjAknJCiKKHc+g7uXMXDUAOg6xv/AoGqeRaQFUDItItIAPBnFFKdlUJK2n4qsUnAaQvvEETa4Ne4+cRhnUKBDPKq1363iq9feomzPfnzeAnzW/zCkMW5cJpKIEENSpwKGjeqJq+eZ/j7P6rYhIi1Qs0ymrbVN4k+l0vQ05WcMpOEd6MRR/N1+PDsL/a3suscQeWYnQvsnNOpe0Pv37uOD518if8MWfOWF+GxO5RoXrqBoIl2Grq1zGDO2M8G9x/l7PavPs4hI80um3W43WVlZxMfHK6GWemWtJSsrC7e7aXVWkIblK/dSujaL4u/2U7ohB3zgah9O9DldCRvUCkd04yx1KCst4/2XZ7N72Up8xQXVOm4E4TQxhLniaB+TxRmnJRDZ/2zoOlYzDIqIHEGzS6Y7duzIzp07ycjICHQo0gy53W46duwY6DAkwKzPUrYxl+Lv9lOyJhNb7sMRE0Lk2E6EDWmFq014oEM8om+++JKV735ARXY2Xl8u1pYC/ocGQx0xxIfnMGpoKB2HjYFuZ0J8d9CghIjIUTW7ZNrlctG1a9dAhyEizVD5niKKv9tH8XcZ+ArKMW4HYYNbEza4NcGJUY3uQcI9u3bzwbP/pWjLdnyeQnw2D/ip7jnKDUN6ljFwzBDodha0GwSOZvffgohIg9JPTRGRo/Dml1Octp/ilfv9U3oHGdy9Ywkb2prQPvEYV+N5kNDj8fD+y7PZ+U1KZb/nAzMNOnCaGCKCo0lsncPpp3cluM/ZkDgKghvnKLqISFOhZFpE5BAH6qCLVu6nbEOOf0bCTpHEnNed0IGtcIQ3nslG1nyXxpevvIln7/6DWtYdKN1oHZHDmJExtBk6AbqdAVHtAxuwiEgzo2RaRAR/HXT51jyKVu6n5IdMbJnXXwd9RifChrbG1apxtH0rKi7mvWdfJHNVOr7SArw2B7BACMFBUUS6YUD3EoadMQy6j4O2AyGo8Yyei4g0N0qmRaRFq8gqoWjlfopX7sObU4YJrpxQZWhrQro2jglVli9eyoo35lKRlUWFN5cDsw06TCzhzhg6xudw5hkdCe83ERJHQ0hkQOMVEWlJlEyLSIvjK62g5IdMilL3Ub41398PukcM0eMTcfeLJyjYEdD4ioqK+N8zL5D1/Tp8ZQVVPZ+NCSUkKIro0HxGDAii16ix0GMcxHULaLwiIi2ZkmkRaRGsz1K2OZfi1P2UrM7Eenw4W4USNTGRsCGtcQa4H/SqFSksfvUtKjKyqPDm4B99NjhMLBGuWBJbZXPmuB4E953on23QGRzQeEVExE/JtIg0axVZJRSl7qN45X68uWX+dnZDWxM2rA3BnSIDNrmTx+Ph3ednsW95WmXtczbgb1sXEhRFTGg+pw1y0m30WdDjZxCt/uYiIo2RkmkRaXZ8Zd7KMo69lG+pLOPoGUv0pERC+yYErJ3dlg2b+OT5lyjbtQevN6+q84Z/9DmOTgnZ/OxnPQnuNwk6jdTos4hIE6BkWkSaBWst5dvyKUrZR8n3mdhyL854N1ETuhA2tE3Ayjg+f28+a+Z/hq8wD68vG/ABwQQHRRPlhlP6+Ug6fSz0HA+xXQISo4iI1J2SaRFp0rz55RSt3Edxyj4qMkswwUGEDmxFeHIbgrtEnfQyjqKiIt596nmyV/+IrzwPn80HKvs+O2NpF53NuDM7EjVoMnQdA67QkxqfiIjUr6Mm08YYNzAZGAO0B0qA1cAH1to1DR+eiMjhrNdH6bpsilbso3R9NvggODGK2DM6EjqgFUEhJ7cbx6Yf1/Ppcy9TvmdvtYcHg3CaWCKDo0hKLGTUhOHQayK06gMBqtMWEZH6V2MybYy5D38ivQhYBuwH3EAv4MHKRPu31trvT0KcIiJ4MoopWrGP4pX78BV6CIoMJnJsR8KGtTnpk6osnP8xP8z7GF9BbrXyDf/EKbGh+YweHkriqHP9Dw+GxZ3U2ERE5OQ52sj0cmvtjBrWPWKMaQ10boCYRESq+Mq9lHyfSdGKvZRvy4cgcPeJJ3x4G9y94jCOkzPKe6D7xt5l3+Ery8dncwF/+UaYM5b2sdn8bFw3wgdP8beuczSeKcdFRKTh1JhMW2s/AH+ph7W2tPo6Y0yCtXY//tFqEZF6V76rkKLleyhOy8CWeXEmhBI9KZGwoW1wRJ6cLhfZWdn874lnKdy4BW9FHtYWAQaniSUiOJo+XQoZPWE4ps+5kNBT5RsiIi1QbR5AXGGMudFa+y2AMeYi4O/4yz1EROqNr7SC4rT9FC3fi2d3ETiDCBuQQPgpbQlOPDkPE27esJFPnvkv5Xv2UeHNBsoBJ66gGKLclpEDDX3OmujvvhGe0ODxiIhI41abZPpK4EVjzCL8DyHGA2c1ZFAi0nJYaynfXkDR8r2UfJ+B9fhwtQsn5rzuhA1uTVBowzcdSl36LUtffQtfTjYVvmzAC7gJCYomPiKPs8a0ps3IC6Dr6eByN3g8IiLSdBzzfylr7Q/GmL8BrwAFwFhr7c4Gj0xEmjVfsYfi7/ZTuHwvFfuKMcEOwoa0JvyUtrg6RDT4KPTn8+azZt6n/v7PNgsAYyIIdcTQPi6Hsyf0InzwBdBhGAQFZpIXERFp/I6ZTBtjXgC6AwPxl3bMN8Y8Ya19qqGDE5HmpWpileV7Kf4+Eyp8uDpGEHNhD8IGtSIopOFGoSsqKvjwtTfYunAp3tK8ag8QRhPujKV7+zzOOGcYrv7nQStVsYmISO3U5n+uH4AbrLUW2GKMGQE80rBhiUhz4iup8I9CL9vjH4UOcRA+rDXhp7QjuENEgx3X4/Hw7nP/9XfgqDaBisPEERUcTd9uxYw69wzocy5EtW+wOEREpPmqTZnHo4e8zgN+0VABiUjzYK3Fs7OQwmV7KFlVWQvdMYLYC3sSOqjhJlYpKynl7aeeIzNtTWUHjkL8HTjiiHZHkNzPx8Dxk6DXBPV/FhGRE3a0SVveB54DPrbWeg5Z1w24DthqrX2xQSMUkSbFV+71d+T4dg+e3UWY4KCqWujgjpENcszikhLeeuzf5K5eh7ciF2uLgSBcQXFEuS2nDgmm98/Ohx5nQ0jDjYSLiEjLc7SR6RuB/wP+ZYzJATLwz4CYCGwCnrTWvtfgEYpIk+DZW0Thsj0Ur9yPLfPiahtOzPmVHTnc9V8LXVRUxJv/eor8dRvxVuRgbQngwBUUS0yoZczISLqecQl0O1MdOEREpMEcbdKWvcDvjDE7ga/xJ9IlwHrrH/YJGGPMFGBKjx49AhmGSItnK3yUrM6k8Ns9lG/NB6chbEArwke2I7hzZL135KhKoNM34vUeSKD9PaDjwixnjImn45hLIHEsOE/OxC4iItKy1Wa4qDXwFrASeBFY1aAR1YK19n3g/eTk5BsDHYtIS1SRW0rRsr0UrdiLr9CDI95N9DldCRvWBkd4/U6jXVRczNuP/ZvcNT9WG4H2J9Dx4ZazTm9Fu9GXQ5dR4Gj4ntQiIiLV1eYBxHuNMX8CxgM/B540xrwJvGCt3dTQAYpI42B9lrJNuRR+s4fSdH9fZnefOCJObU9IjxhMUP2NQpeXlvPmk/8mK21NtRron0agx53Z+qcEOqhhHmQUERGpjVoN41hrrTFmL7AXqABigbeNMZ9Za3/XkAGKSGD5SisoStlH0bd7qMgsISjcReTpnQg/pS3OuPqrRfZ4PLz97xfYvyINrycHa4s4UAMdG2Y56/QEOoy9Qgm0iIg0KrWZtGU6cA2QCfwHuMta6zHGBAEbACXTIs2QZ28Rhd/spvi7/dhyH8GdI4m9rDdhAxIwzvqbEXDurFfZtnAJ3rJcrC3gQBeOmFDLGaNj6HzmldBltEo4RESkUarN/05xwIXW2m3VF1prfcaYyQ0TlogEgvVaStZmUfTNbso25/kfKBzUmohT29VrW7tP353Huvc/wVuSi8/mUdUHOjScsadG0m3c5dD1dHDUb/21iIhIfatNzfSMo6xLr99wRCQQvEUeilbspeibPXjzynDEhBA9KZGw5Lb19kDh8q+X8O3Lb+IryMZrcwBwmjii3JGcMsTFgEmXQ/ezwBlSL8cTERE5GfR3U5EWrHx3IYVLd1OclgEVPkK6RxMztRvupPh6eaBw07p1fPzkC3iyMvD6sgCLw8QQ6YphSF8Pw6eeB70mQnD4iZ+MiIhIACiZFmlhDpRyFC7dRfmWfIwriPBhrYk4tT2utiee1GZnZvPWQ49Run0nFb5MwIsxEYQ7Y+ndpZAzL5kIfSZDaMwJH0tERCTQlEyLtBC+kgqKVuylcOluvLllOGJDiD6nK+HJbQgKO7FSjrLSMt547Glyvk+noiILKMMYN25HLF1a5zHholG4Bl8CkW3r52REREQaCSXTIs2cJ6OYwiW7KV65z9+Vo2s0MZO74e574qUcH8x+k40fL8RbnoO1hYCT4KBYWkfnMWnqQKJGXAXx3evnRERERBohJdMizZC1lrKNuRQu3kXpjzngMIQNbk3EqPYEt484oX2vWr6CL//zGt78LHw2hwOdOGLCffzszHZ0OOsaaD8U6nkqcRERkcZIybRIM2I9PorT9lO4ZBeevcUERbiI+llnwke2wxERXOf9Zu7fz1v/fJzyXbsr66B9BJloooKjGTHEwcAp06DbGeoFLSIiLY7+5xNpBryF5RR9u4fCb/fgK/TgahtO7MW9CBvcqs4TrHg8Ht564lkyUr+vVgcdRpgzlh6dijj78smQNBlC6q//tIiISFOjZFqkCfPsL6bw610UfbcPKizuPnFEjO5ASPdoTB3LLBZ/uoCVr8/FW5yNz+YDDoKD4mgbk8+5Fw4n7JQrIbpD/Z6IiIhIE6VkWqSJsdZStimPwq93+uuhnUGED2tDxOgOuFqF1Wmfe3fv4d1/Pkb5vr14fZkAOEwccaERnHV6K7pMuB7aDVYdtIiIyCGUTIs0Edbro+T7TAq+3olndxFB4SdWD+3xeHjrqWfJWPE9FRWZQDnGhBPujGNA73JGXXoZ9DgbnHWvtRYREWnulEyLNHK+sgqKlu+jcMkuvLllOFuFEnthT8KGtMa4jr8eOnXptyx5YTbeoszKMg5/O7t2cflMuWw0IcOuhPD4+j8RERGRZkjJtEgj5c0vo3DJbgqX7cGWegnuGkXM1O64+8Qdd3/ogvwCXn/wEYq3bKvqxuEwccSGRnLmGQl0nXgTtB3QMCciIiLSjCmZFmlkPPuLKfhqJ8Xf7QefJbR/AhFjOhDSOeq49/XJO3NZ994neMuzsLYYY9yEOmPp3bWEcVddAj0nqIxDRETkBCiZFmkkyrbnU7BoJ6XpWeAIInx4WyLHdMAZH3pc+9m3dy/vPPhotYcJDa6geOIjLedeOISYUT+HqHYNcxIiIiItjJJpkQCy1lL6Yw4FX+6gfEs+JtRJ5JmdiDit/XE/VDh31qtsW/A1FZ5M/D2hI4hwxXLKEAdDLr4ZOo9UNw4REZF6pmRaJACs11LyQwYFi3bi2VuEIzqY6HO7EX5KW4JCHLXez+7tO/jfPx/Hk7kPr80GgggOiqdtXAFTrxpHyNBLwR3dcCciIiLSwimZFjmJrMdH0cp9FHy5E292Kc7WocRe0ouwQcc3U+H/XnyZ7V8srhyFLifIRBEVEs2Y0yLpc94vod2ghjsJERERqaJkWuQk8JVVUPTtXgoW78RX4MHVMYKYc5NwJ8XXujPHnl27effBR/Fk7D1oFLpDq0KmTJuIa9DFEFy3SVtERESkbpRMizQgX7GHgiW7KVy6G1tSQUiPGCIv60hI95haT/f9wZy32PjBgqpaaP8odAxjRkfT57zboE3fhj0JERERqZGSaZEG4C0op2DxLoq+2YMt9+LuG0/UmZ0I7hRZq/fn5WQz+/7/R9nuXZUdOfyj0O1bFTD1mkn+UWjX8XX5EBERkfqnZFqkHlXkllH41U4Kl+8Fr4/Qga2IOrMTrrbhtXr/0s++YMWrb+Ety8TaEowJJ9IVy4gRYQy65FfQtn8Dn4GIiIgcDyXTIvWgIruUgkU7KErdBxbChrYm8oxOuBKOPXrs8Xh49cFHyF/7IxW+DMDiDEqgTQxMnfYzwpKvgJCIBj8HEREROX5KpkVOQEVmCfkLd1D83T4wxj/Ryhkdcca4j/neTevW8eEj/6YiPwOfzQdCCHPGMag/nDbtdugwTH2hRUREGjkl0yJ14MkopuCLHRSn7QdHEBGntidybEcc0SHHfO8Hs99g04ef4/FkAB4cJpb48EjOuaA/rc+6BcLjG/4EREREpF4omRY5Dp79xRR8sZ3iVRkYZxARozv4k+jIo89WWFRUxGv3P0Tx1q0HPVCY2LGYKTdOgx4/g6Da95kWERGRxkHJtEgtHJZEj+1I5JgOx5zy+8fVa/j0X8/gKcrA2kKMCSPCFcvoUdH0u/QOiO9+ks5AREREGoKSaZGjqGsS/cHsN9n04QI8nv1ABQ4TR3xUOFOvPp3oU6+F4Np19xAREZHGTcm0yBF4Mkso+Hw7xWn7a51EezweXvrLPyncuBGvLwMIIiQonu6dS5h0y42QOFoPFIqIiDQzSqZFqqnILiX/8+0Uf7cP4wgiYkxlTfRRkug9O3fxzt/+H56cvfhsPsaEEuGK47RR0Qy4/LcQm3jyTkBEREROKiXTIvgnWylYuJ2iFfsgCH93jjM6HfXBwmWLvmbZi7PxlO3HP813DPHhkZx/1anEjL5evaFFRERaACXT0qJ5C8opWLiDwmV7AAg/pS1RZ3Y6aou7t5/7L7sWLaHCux/w4QpKoGNrF1NuvgpXnwnqyiEiItKCKJmWFslX7KHgq10ULtmF9foIG9qGqHGdccYeebIVj8fDS3/9J4UbNlS2tnPidsQxsH8QY35xJ7Tpe3JPQERERBoFJdPSovjKvBQu2UXBVzuxpV5CB7Ui6uwuNU77nZWRyZz7HsSTuQefzfPXQwfHctb4DvS84LcQ0eokn4GIiIg0JkqmpUWwFT4Kv91DwcId+Io8uJPiiBqfSHC7I7eo27B2LR89/DQVxfuxtpggE01ceCTnXz2G2DHXg+vY04WLiIhI86dkWpo167MUf7ef/M+24c0tI6RbNFETEwnpHHXE7ZcuWEjKS6/jKd8PeHCaeNq2Mlz4yytxJU1SPbSIiIgcRMm0NEvWWkrTs8n7ZCsV+4pxdYgg9qKehPSIwRyh1/O8l15jyycLKx8qtAQHJdA90cM5v7oD2g8+6fGLiIhI06BkWpqdsq155H20lfJt+TgTQom7sg+h/RMwQYcn0a898gSZK76jwrefAw8VDhkawmm/+L36Q4uIiMgxKZmWZsOzv5i8j7ZQmp5NUFQwMRf0IDy5DcZxcGmGf6bCf1C4cQNeXxYQQrgrnrHj2tH30rsgPD4wJyAiIiJNjpJpafK8+WXkL9hO0Yq9mGAHURMSiRjVnqBgx0HbFRUV8fKf/kbp7u34bC7GhBEZEsOkiwfSaeLtEBwWoDMQERGRpkrJtDRZvtIKCr7cSeHiXVifJeK09kSe1RlHuOug7fJysnnl3r/hydqNzxYQZCKJDYvm/OvOJG70teBw1XAEERERkaNTMi1NjvVailbsIf+z7fiKPIQOakX0+C444w/uFb1vzx7emvkPyvP2YG0RQSaa1lFRXHjrpYQPmqrOHCIiInLClExLk2GtpXRdNnkfbqEio4TgrlHEnNOP4E6RB223c+s2/vfXh/EU7cXaEhwmljbxDi6+4xe4epwJR+jmISIiIlIXSqalSSjfVUjeB5sp25yHMyGU+Gl9cfeNO6jN3eYNG5n/wL/wFO8DSnGaeDq0dXHeb27F1fXUwAUvIiIizZaSaWnUvHll5H2yleLv9hMU6iRmanfCR7Q9qEPHpnXr+PDBxykv2QeU4TQJdO7k5oLf3gVt+wcueBEREWn2lExLo+Qr91L49S4KFu3wP1w4piNRZ3YiKPSnW3bd6jV8+vCTeEr2AeW4ghJI7BLO1N/+AVr1ClzwIiIi0mIomZZGxVpLyaoM8j7aijevjND+8URP6nrQw4Ub1qbz0T8fPyiJ7t4Nzv2/P0F898AFLyIiIi1Oo0mmjTHhwNNAObDIWvtagEOSk6x8RwG572+ifHsBrvbhxF3Wi5BuMVXrN6/fwAcP/KuqnMMVlED37oZzfztDsxWKiIhIQDRoMm2MeRGYDOy31vavtnwi8BjgAP5jrX0QuBB421r7vjHmDUDJdAvhzS8n7+MtFK/cT1Cki9iLexI2tE3V9N/bNm7mvb89XPlgYWUS3S2Sc+/6C8R0DmzwIiIi0qI19Mj0LOBJ4OUDC4wxDuAp4GxgJ7DCGDMP6Aj8ULmZt4HjkkbAVvgoXLKL/M93YL0+Is/oSOSZnQgK8d+We3bu4u0Zf6e8cC9QijMogW6JkUz53UyNRIuIiEij0KDJtLX2K2NM4iGLTwE2Wms3AxhjXgfOw59YdwTSAM2m0cyVrMsmb/5mKjJLcCfFEXNuN5wJ/rro7MxsZv/hPsrzd2NtCU4TT2KXMM67ewbEdQ1w5CIiIiI/CUTNdAdgR7XXO4ERwOPAk8aYc4H3a3qzMeYm4CaAzp31J/6mxpNZQt77myj9MQdnq1ASft4Pd+84AAry8nn5D/dRlrULawtxmDg6dnBz8R//rAcLRUREpFFqNA8gWmuLgJ/XYrvngOcAkpOTbUPHJfXDV+6lYOEOCr7aiXEGEX1OVyJOa49xBlFWWsaLf5hJ6e5t+Gw+DhNDm9YxXHzPnbg6DAx06CIiIiI1CkQyvQvoVO11x8pl0gxZayldk0Xu/M14c8sIG9Ka6EldcUQFU1FRwUt//jsF69fhtTkEmShax8Zw+R9+g6tzcqBDFxERETmmQCTTK4Cexpiu+JPoy4ErAxCHNDBPZgm58zZRtj4HV9sw4m4aSEi3aABmP/Ik+1eswOvLwphwYiNiuOx3NxLe+/QARy0iIiJSew3dGm8OcAaQYIzZCcyw1r5gjLkd+AR/a7wXrbVrGjIOObmsx0v+wh0UfFlZ0jG5GxGntsc4DO+/MofNH35GhW8/xriJcsdy4e0XEJ98ARgT6NBFREREjktDd/O4ooblHwIfNuSxJTBKf8wm571NeLNL/SUd53TFERnMskVfs+z5l/BU7AOchLvimXT1aXQZfyMEqXmLiIiINE2N5gFEadq8eWXkzt9MyQ+Z/i4dNwzA3SOGzes3MP+2/4enbC/gw+1I4PTJPel/2V3gcAU6bBEREZETomRaToj1Wgq/2U3+p9uwPkvU+C5Eju1IfkEu/7np15W9oksJDmrFwFNiOP22+yE4PNBhi4iIiNQLJdNSZ+U7C8h5dwOe3UW4e8cSM7U7NsrJ83/4M0XbN+GzBf4JV7pHc94f/wHhCYEOWURERKReNclk2hgzBZjSo0ePQIfSIvnKvOR/to3CJbsIiggm7qo+hPZP4M2nnmXP4m/w2iyCTBRt4mO5bOa9uFr3DHTIIiIiIg3CWNt05z1JTk62KSkpgQ6jRSlZl03u3I14c8sIH9GW6IldWfzlF6S9/AYe7z6McRMZGsplv7ueqKQzAx2uiIiISL0wxqRaaw+bCKNJjkzLyectKCf3/U2UfJ+Js3UYrW4ZSAa5vHLzbZUPF0KoM4EJV4+m+8RfqM2diIiItAhKpuWorLUUp+4nd/5mbIWXqLO7EDwigRfu+ROlmTuwtpjgoFYMGd2O0TfPBGdwoEMWEREROWmUTEuNKrJLyfnfBso25BKcGEXsRT158+XnyHg+Fa/NwWFi6dw1igtnPAShsYEOV0REROSkUzIth7G+ynZ3n2wFDDHndSctbw0rpz9eWRcdSmxkHFfN/B0hHfsHOlwRERGRgFEyLQfx7C8m550NlG/LJ6RXLIyN5b9/+yvlxXsAS5gzgXNumkCX0484uaWIiIhIi6JkWgD/5CsFX+8kf8E2goIdxF7aizlvPknBpxvx2QJcQQkMGdWRMb+cCQ7dNiIiIiKgZFqoHI1+az3lOwpw94vnO5PO2oeeo8KXQZCJon27VlzxwP+DsLhAhyoiIiLSqCiZbsGsz1L49S7yPttKULADxiUw5+VHKS/fAziICInn0t/fRGzSqECHKiIiItIoNclkWjMgnjhPRuVo9PYCQpJi+WD122S/kI61RQQHteK0yf0YduVv1S9aRERE5Cg0A2ILc6BTR95HWzGuIHa0y2HZwrcqSzqi6dg5jEvufxSCwwMdqoiIiEijoRkQhYrcMnLeXk/ZxlxM13A++PZlCtZuAYKIcCdw+Z9vJ7r7YfeIiIiIiNRAyXQLYK2l+Lv95M7bBD7LWuePrFn4OT5bQHBQK06Z1JcR0+5USYeIiIjIcVIy3cx5izzk/m8DJauz8MQZPv/hLfLKNxFkImnbrg1X/eNRCIkMdJgiIiIiTZKS6Was5Mdsct5aj6+kgrUlP7A69Qss5YS6Erjodz+nzcDTAx2iiIiISJOmZLoZsh4feR9toXDpbkpCylm8fT7Z5RtwmDiShvVgwm8fgKCgQIcpIiIi0uQpmW5mPHuLyJqzjop9xWwsWsd3Wz/EZ4OIjkhg2oN/IaRV50CHKCIiItJsKJluJqy1FC7ZTd7HWyj3lrNs36fsLl5NcFArxlx6KoMvvCnQIYqIiIg0O0qmmwFvQTnZb62nbH0Oe0p2sGz/XMqtoXWbtkx7+HEIDgt0iCIiIiLNkpLpJq50Qw7Zb/yIp7CUtKyv2Zi/jBBHK87/1SV0O+2cQIcnIiIi0qw1yWRa04mD9frI/2wb+Yt2UlCew9L971JQUU7n7l245K+PgaNJfmlFREREmhRNJ94EVeSUkj1nHeXbC9hUsJrvsj4l2BnDxff+ktZ9Tgl0eCIiIiLNjqYTbyZKVmeS+UY6FWXlrMj8iF3Fe+jatxvn/+lhzWAoIiIicpIpmW4ibIWPvA/9vaNzyvaxdP97+IyTafffTUKPgYEOT0RERKRFUjLdBFTklLL7PysJyvLyY94Kfsj5jsR+bTn/3n9qNFpEREQkgJRMN3IlP2az979pWJ+PFRkfkeXJ5pq/301c1/6BDk1ERESkxVMy3UhZn2X3u6uxKbkUleXwTcZnxHYP49aZL2o0WkRERKSRUDLdCHkLykl/6DNiyqPYXLCKNbmrOO+eG2g/YESgQxMRERGRapRMNzI5a/eQ+d/viTChLM9eQElMJr984lkICgp0aCIiIiJyCCXTjciKp+fTZlsEXl85SzM+ZcQN4+gzVrMYioiIiDRWSqYbgfLiMlLufZ3Owd3YU7KV9PJvuPa/j4HDFejQREREROQolEwH2I9LvqPs7S10DulGeu53uEc7uebqpwMdloiIiIjUQpNMpo0xU4ApPXr0CHQoJ2T+zMfpU9SbSFcUy7M/Z8KDvyQ8rk2gwxIRERGRWmqST7VZa9+31t4UHR0d6FDqpDC/gM9ufYKBJYMo95axNvwbLnxuphJpERERkSamSY5MN2XL536Ia1EhSVGD2V64kXbX9ubc4X8KdFgiIiIiUgdKpk+i9373MP3sAJzBCXxfuJhJj92FcehLICIiItJUKZM7CbL2ZvDDX99hWOQIcsszyei6h3N+9ftAhyUiIiIiJ0jJdAP76tW3iUk19Ijsx+aCdHpNH0v/XhcEOiwRERERqQdKphvQvLseoh+DMU5IK/2KyU/+AYwJdFgiIiIiUk+UTDeA/OwcUv40h6GRI8kpyyCnbw6Tb/xjoMMSERERkXqmZLqerZj3Ec4vCugVOYAtBT/S8zdjGdCzV6DDEhEREZEGoGS6Hr3/x/9Hn/J+uIITSCtawuQn71ZZh4iIiEgz1iQnbWlsysvK+OTWRxlUkYzHV87GdulMfuIeJdIiIiIizZxGpk/QzvR0dj69gn5Rw9hZtIXWP+/Hz4ZODXRYIiIiInISKJk+AYv+8zKt1kbRPqwLawtSOOORXxIc4g50WCIiIiJykiiZrqP3f/cP+tmhWIdllfmWKU/dHeiQREREROQkUzJ9nEqLivj6rhcYHHkquZ4scgcVM+U6JdIiIiIiLZGS6eO04A9PMjjqNLYXbqLLr05jQG+1vRMRERFpqZpkNw9jzBRjzHN5eXkn/dg/e+B2vitfSvL/u4IOSqRFREREWjRjrQ10DHWWnJxsU1JSAh2GiIiIiDRzxphUa23yocub5Mi0iIiIiEhjoGRaRERERKSOlEyLiIiIiNSRkmkRERERkTpSMi0iIiIiUkdKpkVERERE6kjJtIiIiIhIHSmZFhERERGpIyXTIiIiIiJ11KRnQDTGZADbAnDoBCAzAMdtynTNjp+u2fHTNTt+umbHT9fs+OmaHR9dr+N3Mq5ZF2ttq0MXNulkOlCMMSlHmk5SaqZrdvx0zY6frtnx0zU7frpmx0/X7Pjoeh2/QF4zlXmIiIiIiNSRkmkRERERkTpSMl03zwU6gCZI1+z46ZodP12z46drdvx0zY6frtnx0fU6fgG7ZqqZFhERERGpI41Mi4iIiIjUkZJpEREREZE6UjJ9CGPMRGPMj8aYjcaYe46wPsQY80bl+mXGmMRq635fufxHY8yEkxp4ANXimv2fMWatMeZ7Y8znxpgu1dZ5jTFplR/zTm7kgVOLa3adMSaj2rW5odq6a40xGyo/rj25kQdOLa7Zv6pdr/XGmNxq61rcfWaMedEYs98Ys7qG9cYY83jl9fzeGDO02rqWeo8d65pdVXmtfjDGLDXGDKq2bmvl8jRjTMrJizqwanHNzjDG5FX7/vtztXVH/Z5ujmpxve6qdq1WV/7siqtc11LvsU7GmIWVecQaY8z0I2wT2J9n1lp9VH4ADmAT0A0IBlYBfQ/Z5lbgmcrPLwfeqPy8b+X2IUDXyv04An1OjeSanQmEVX7+ywPXrPJ1YaDPoZFes+uAJ4/w3jhgc+W/sZWfxwb6nBrDNTtk+18BL1Z73RLvs7HAUGB1DevPAT4CDDASWFa5vEXeY7W8ZqcduBbApAPXrPL1ViAh0OfQCK/ZGcD8Iyw/ru/p5vJxrOt1yLZTgC+qvW6p91g7YGjl55HA+iP8nxnQn2camT7YKcBGa+1ma2058Dpw3iHbnAe8VPn528A4Y4ypXP66tbbMWrsF2Fi5v+bumNfMWrvQWltc+fJboONJjrGxqc19VpMJwGfW2mxrbQ7wGTCxgeJsTI73ml0BzDkpkTVS1tqvgOyjbHIe8LL1+xaIMca0o+XeY8e8ZtbapZXXBPSzDKjVfVaTE/k52GQd5/Vq8T/HAKy1e6y1Kys/LwDSgQ6HbBbQn2dKpg/WAdhR7fVODv+CVW1jra0A8oD4Wr63OTre8/4F/t8eD3AbY1KMMd8aY85vgPgao9pes4sq/1z1tjGm03G+t7mp9XlXlhF1Bb6otrgl3mfHUtM1ban32PE69GeZBT41xqQaY24KUEyN1anGmFXGmI+MMf0ql+k+OwpjTBj+pO+daotb/D1m/KW1Q4Blh6wK6M8zZ33vUKQmxpirgWTg9GqLu1hrdxljugFfGGN+sNZuCkyEjcr7wBxrbZkx5mb8fw05K8AxNRWXA29ba73Vluk+k3pjjDkTfzI9utri0ZX3WGvgM2PMuspRyJZuJf7vv0JjzDnAXKBnYENqEqYAS6y11UexW/Q9ZoyJwP/LxW+stfmBjqc6jUwfbBfQqdrrjpXLjriNMcYJRANZtXxvc1Sr8zbG/Az4IzDVWlt2YLm1dlflv5uBRfh/42zujnnNrLVZ1a7Tf4BhtX1vM3U85305h/xptIXeZ8dS0zVtqfdYrRhjBuL/njzPWpt1YHm1e2w/8D9aRpnfMVlr8621hZWffwi4jDEJ6D47lqP9HGtx95gxxoU/kX7NWvvuETYJ6M8zJdMHWwH0NMZ0NcYE47+ZD33yfx5w4GnQi/E/HGArl19u/N0+uuL/zXv5SYo7kI55zYwxQ4Bn8SfS+6stjzXGhFR+ngCMAtaetMgDpzbXrF21l1Px14gBfAKMr7x2scD4ymXNXW2+NzHG9MH/kMk31Za11PvsWOYB11Q+BT8SyLPW7qHl3mPHZIzpDLwLTLPWrq+2PNwYE3ngc/zX7IjdGloaY0zbyueKMMacgj/vyKKW39MtkTEmGv9fcN+rtqzF3mOV988LQLq19pEaNgvozzOVeVRjra0wxtyO/0I78HcDWGOM+QuQYq2dh/8L+ooxZiP+hwgur3zvGmPMm/j/k64Abjvkz8zNUi2v2UNABPBW5c/U7dbaqUAS8Kwxxof/B+yD1tpmn+TU8pr92hgzFf+9lI2/uwfW2mxjzF/x/0cE8JdD/gzYLNXymoH/+/H1yl9wD2iR95kxZg7+TgoJxpidwAzABWCtfQb4EP8T8BuBYuDnleta5D0Gtbpmf8b/jMzTlT/LKqy1yUAb4H+Vy5zAbGvtxyf9BAKgFtfsYuCXxpgKoAS4vPL784jf0wE4hZOqFtcL4ALgU2ttUbW3tth7DP8AyDTgB2NMWuWyPwCdoXH8PNN04iIiIiIidaQyDxERERGROlIyLSIiIiJSR0qmRURERETqSMm0iIiIiEgdKZkWEREREakjJdMiIiIiInWkZFpEREREpI6UTIuIiIiI1JGSaRERERGROlIyLSIiIiJSR0qmRURERETqSMm0iIiIiEgdKZkWEREREakjJdMiIiIiInWkZFpEREREpI6UTIuIiIiI1JGSaRERERGROlIyLSIiIiJSR0qmRURERETqSMm0iIiIiEgdKZkWEREREakjJdMiIiIiInWkZFpEREREpI6cgQ7gRCQkJNjExMRAhyEiIiIizVxqamqmtbbVocubdDKdmJhISkpKoMMQERERkWbOGLPtSMtV5iEiIiIiUkdKpkVERERE6kjJtIiIiIg0CdbaQIdwmCZdM30kHo+HnTt3UlpaGuhQpBlyu9107NgRl8sV6FBERERaDJ/Xy8qP5rH9hzQuuHsGJqjxjAc3u2R6586dREZGkpiYiDEm0OFIM2KtJSsri507d9K1a9dAhyMiItIi7Nuyic+ee4J9mzfSbehwyktLCQkLC3RYVZpdMl1aWqpEWhqEMYb4+HgyMjICHYqIiEiz5ykrZelbs0n9YC6hkVFM/s3d9Bo5utHleM0umQYa3UWW5kP3loiISMPbumolC/7zFHn79zHgrPGMvep63BERgQ7riJplMh1oDoeDAQMGVL2+/PLLueeee2rcftasWaSkpPDkk0+ejPBEREREGqWi3By+fPVF0r9eSGy7Dlw64+906jvg2G8MICXTDSA0NJS0tLQG239FRQVOp750IiIi0jxYn4/vP/+Er+fMwlNaxsgLL2PEBZfhDA4OdGjH1HgehWwBEhMTyczMBCAlJYUzzjjjsG0yMjK46KKLGD58OMOHD2fJkiUAzJw5k2nTpjFq1CimTZt2MsMWERERaTD7t25mzp/vYsF/nqJ1l25c89ATjLpsWpNIpEEj0w2ipKSEwYMHV73+/e9/z2WXXVar906fPp077riD0aNHs337diZMmEB6ejoAa9euZfHixYSGhjZE2CIiIiInTXlpCUvfms3KD9/DHRHJpNv+j6QxZza555OadTJ93/trWLs7v1732bd9FDOm9DvqNidS5rFgwQLWrl1b9To/P5/CwkIApk6dqkRaREREmjRrLeu/XcKiV/5DYVYmA8dNZPSV1xIaERno0OqkWSfTjY3T6cTn8wHUOKmMz+fj22+/xe12H7YuPDy8QeMTERERaUhZu3bwxX+fZfsPabRK7Mbk6XfToXdSoMM6Ic06mT7WCPLJlpiYSGpqKpMmTeKdd9454jbjx4/niSee4K677gIgLS3toJIRERERkaamvLSEb999g9T5c3GFhHDWz29m0NnnEORwBDq0E6YHEBvAgZrpAx8H2uLNmDGD6dOnk5ycjKOGm+fxxx8nJSWFgQMH0rdvX5555pmTGbqIiIhIvbHW8uM3i/nv//2SFe+9TdLoM7j+0WcZMnFKs0ikAYy1NtAx1FlycrJNSUk5aFl6ejpJSU37zwXSuOkeExERObaMbVtYOOs5dqz9gVaJ3Rh3/S+bdEmHMSbVWpt86PJmXeYhIiIiIidXSUE+S958je8/+4iQ8HDG/eJWBo6b0GxGog+lZFpERERETpjP62XVgo9Y+uZrlBUXMWj8OZx26VVNtktHbSmZFhEREZETsu2HNBa9/B8yt2+lc/+BnHntTSR0Tgx0WCeFkmkRERERqZPs3Tv58pUX2LxyBVGt2jDl/35Pz1NOa3ITr5yIE06mjTGda7lprrW2fmdQEREREZGTrqSwgG/ens2qTz/EGRzMmCuvY+ikqU1mCvD6VB8j0y8BFjjaryAWmAW8XA/HExEREZEA8FZ4WPXph3zz9hzKiosZcNZ4Trv0KsJjYgMdWsCccDJtrT2zPgIRERERkcbJWsuGZUv4es5L5O7dQ+cBgznjmhto1ULqoo+m3idtMcaEG2OaZ++TAEhMTCQzMxOA0047rWr5XXfdRb9+/bjrrrt45plnePnl4xv0j4iIOOJyh8Nx0IQzDz744FH3M2vWLG6//fbjOvah9uzZw+TJk6teL1++nLFjx9K7d2+GDBnCDTfcQHFxca1imDlzJg8//PBRjzd37lzWrl171G2efPJJXnzxxeM4CxERkeZp14/pzPnzXbz/rwdxOF1ccPcMLv7jX5VIV6qPmukg4HLgKmA4UA6EGGMygA+AZ621G2u5LzfwFRBSGdvb1toZJxpjc7F06dKqz5977jmys7NrnEmxrkJDQ0lLS6vXfVZXUVGB03nwbffII49w4403ArBv3z4uueQSXn/9dU499VQA3n77bQoKCggLC6uXGObOncvkyZPp27dvjdtcf/31jBo1iuuvv75ejikiItLU5OzZxdezX2LD8qWEx8Zx9k2/ov8ZP2u2/aLrqj5GphcC3YHfA22ttR2tta2A0cC3wD+MMVfXcl9lwFnW2kHAYGCiMWZkPcR40mzdupU+ffpw3XXX0atXL6666ioWLFjAqFGj6NmzJ8uXLwcgOzub888/n4EDBzJy5Ei+//57ALKyshg/fjz9+vXjhhtuoPoMlQdGk6dOnUphYSHDhg3jjTfeOGg0dtOmTUycOJFhw4YxZswY1q1bB8CWLVs49dRTGTBgAPfee+9xn1f1EfKUlBTOOOOMw7bJyMjgoosuYvjw4QwfPpwlS5YA/tHiadOmMWrUKKZNm3bY+9555x0mTpwIwFNPPcW1115blUgDXHzxxbRp06bGa1aT559/nuHDhzNo0CAuuugiiouLWbp0KfPmzeOuu+5i8ODBLFmy5KCReIfDwbZt2wgLCyMxMbHq6yUiItJSFOXmsOCFfzPrt7ey9fvvOO3Sq/jFo88164lXTkR9PID4M2ut59CF1tps4B3gHWOMqzY7sv7MsbDypavyo+7znX90D+z9oc5vP6K2A2DS0UsfNm7cyFtvvcWLL77I8OHDmT17NosXL2bevHk88MADzJ07lxkzZjBkyBDmzp3LF198wTXXXENaWhr33Xcfo0eP5s9//jMffPABL7zwwmH7nzdvHhEREVUjyDNnzqxad9NNN/HMM8/Qs2dPli1bxq233soXX3zB9OnT+eUvf8k111zDU089VWPsJSUlDB48uOr173//ey677LJaXZrp06dzxx13MHr0aLZv386ECRNIT08HYO3atSxevJjQ0NCD3rNlyxZiY2MJCQkBYPXq1Vx77bVH3H9N16wmF154YdWI97333ssLL7zAr371K6ZOncrkyZO5+OKLAar28dRTT/Hll1/SpUsXAJKTk/n666855ZRTanX+IiIiTVlZcREp779Lygdz8VVUMOCsCZx68RUt+uHC2qiPBxA9AMaYx4Df2OpDqYdsUxuV9dapQA/gKWvtskPW3wTcBNC5c2278p1cXbt2ZcCAAQD069ePcePGYYxhwIABbN26FYDFixfzzjvvAHDWWWeRlZVFfn4+X331Fe+++y4A5557LrGxtb+BCwsLWbp0KZdccknVsrKyMgCWLFlSdbxp06Zx9913H3EfJ1LmsWDBgoNqkfPz8yks9P9uNHXq1MMSafDXS7dq1apW+6/pmtVk9erV3HvvveTm5lJYWMiECRNq3HbJkiU8//zzLF68uGpZ69atq0b2RUREmquK8nLSPv2AZXPforQgn96njmHU5dOIbds+0KE1CfU5aUsBMM8Yc7m1tsgYMwH4s7V21PHsxFrrBQYbY2KA/xlj+ltrV1db/xzwHEBycvLRR62PMYLcUA6MsgIEBQVVvQ4KCqKioqLBjuvz+YiJiakxGT6RBupOpxOfzwdAaWlpjcf/9ttvcbvdh60LDw8/4ntCQ0MP2l+/fv1ITU3lvPPOq3OsB1x33XXMnTuXQYMGMWvWLBYtWnTE7fbs2cMvfvGLqhH/A0pLS4/4C4CIiEhz4PN6WfvVFyx9ezYFmRl0GTiEMVdcS5tuPQIdWpNSb908rLX3AnOARcaYJcD/AfecwP5y8ddjT6yXABuZMWPG8NprrwGwaNEiEhISiIqKYuzYscyePRuAjz76iJycnFrvMyoqiq5du/LWW28B/jY2q1atAmDUqFG8/vrrAFXHPR6JiYmkpqYCVI0OH2r8+PE88cQTVa9rM8Ldq1evqtF6gNtvv52XXnqJZct++oPEu+++y759+2q8ZjUpKCigXbt2eDyeg845MjKSgoICADweD5dccgn/+Mc/6NWr10HvX79+Pf379z/mOYiIiDQl1udj3ZIvmfXbW/nkmccIj47hkj/9jYv/+Fcl0nVQb8m0MWYccCNQBCQAv7bWfn2c+2hVOSKNMSYUOBtoln9nnzlzJqmpqQwcOJB77rmHl156CfDXBX/11Vf069ePd99997hLWV577TVeeOEFBg0aRL9+/XjvvfcAeOyxx3jqqacYMGAAu3btqvH9B2qmD3zcc889VXFNnz6d5OTkGjuIPP7446SkpDBw4ED69u3LM888c8x4w8PD6d69Oxs3+hu+tGnThtdff50777yT3r17k5SUxCeffEJkZGSN16wmf/3rXxkxYgSjRo2iT58+Vcsvv/xyHnroIYYMGcLSpUtJSUlhxowZVee8e/duwF/6cfbZZx/zHERERJoCay0bU5bxyt2/5oPHH8LhdHLenfdy5d8eoXP/QYEOr8kyRyhxrtuOjPkCf1nHYmPMAOAV4P+stV8cxz4G4p9R0YE/0X/TWvuXmrZPTk62KSkpBy1LT08nKSmpLqcgAfK///2P1NRU7r///kCHUuW7777jkUce4ZVXXjlsne4xERFpSqy1bPv+O5a8+Sp7N64ntl17Tr3kKvqcOgYTVO9TjjRbxphUa23yocvrrWbaWntWtc9/MMZMwt/N47Sa33XYPr4HhtRXTNI0XHDBBWRlZQU6jINkZmby17/+NdBhiIiI1Jm1lu2rV7H0rdns/nEtkfGtGH/zr+l3+ji1uKtH9TFpi6mhg8eeytKPGrcROeCGG24IdAgHUXmHiIg0ZdtXf8/St15j17o1RMTFM+76X9L/rPE4XbXqVizHoT5GphcaY94B3rPWbj+w0BgTDJxqjLkW/4OEs+rhWCIiIiJSgx1rf2DpW6+xc+1qImLjOOvnNzPgrAk4g4MDHVqzVR/J9ETgemCOMaYrkAu48dc9fwo8aq39rh6OIyIiIiKHsNay/YdVfPvu6+xMX014TCxnXncTA8dNVBJ9EtTHpC2lwNPA05UzHSYAJZWt7URERESkAVhr2ZKWwrfvvM6eDT8SERfPmdfdxIBxE3AFhxx7B1Iv6u0BRGPMZ8Cd1tpV9bXPpioiIqJq5j+AWbNmkZKSwpNPPhnAqERERKQ5sD4fm1KX8+27r7Nv80YiE1rxsxtupd8ZZ6smOgDqcwbEu4FHjTFbgT9Ya/fU475FREREWjSf18uPS79i2dy3yNq5neg2bRl/86/pO/ZMHE4l0YFSn63xVgJnGmMuAj42xrwL/NNaW1Jfx2gOMjIyuOWWW9i+3f+s5qOPPsqoUaOYOXMmERER3HnnnQD079+f+fPnAzBp0iRGjx7N0qVL6dChA++9956muRYREWkhKsrLWfPlAlbMe4e8/fuI79iZc27/Lb1PG6sWd41AfY5MY4wxwI/Av4H7gRuNMb+31h4+80UzdmAWwQOys7OZOnUqANOnT+eOO+5g9OjRbN++nQkTJpCenn7U/W3YsIE5c+bw/PPPc+mll/LOO+9w9dVXN+QpiIiISICVFRfz/YKPSP1gLkW5ObTt0YszrrmR7sNO0WQrjUh91kwvAboCa4BvgevwTwU+3Rgzxlp7U30dq7b+sfwfrMuu39nI+8T14e5T7j7qNqGhoaSlpVW9PlAzDbBgwQLWrl1btS4/P/+g+uoj6dq1a1VyPmzYMLZu3Vqn2EVERKTxK8zJZuVH81j16YeUlxTTecBgzvnVnXTqNxD/uKU0JvU5Mn0TsPYIk7P8yhhz9KHXFsTn8/Htt9/idrsPWu50OvH5fFWvS0tLqz4PCfnpiVyHw0FJiSpnREREmpvs3TtZMe9d0r/+Ap/XR88RpzF8yoW07dEr0KHJUdRnzfSao6w+t76OczyONYIcCOPHj+eJJ57grrvuAiAtLY3BgweTmJhYVSO9cuVKtmzZEsgwRURE5CSw1rL7x3RS5r/LxpRlOJ0u+p85nmGTzye2bftAhye1UK810zWx1m4+GcdpCh5//HFuu+02Bg4cSEVFBWPHjuWZZ57hoosu4uWXX6Zfv36MGDGCXr30W6iIiEhz5fN62bB8KSnz/8fejetxh0cw8sLLGDJhMmHRMYEOT46DObwq4wR3aMwUa+379brTGiQnJ9sDtcgHpKenk5SUdDIOLy2U7jEREamr8pJifvjiM1Z+9B75GfuJaduOYeecT7/Tx+E6pARUGhdjTKq1NvnQ5Q0xMv034KQk0yIiIiJNQd7+fXz3yXxWf/EpZcVFdOjTjzOu9XfmCApSe7umrCGSaT1mKiIiIi2etZad6atZ+eE8NqUsAwO9Roxi2OTzadejd6DDk3rSEMl0/daNiIiIiDQhFeXlrFv6FSs/mkfG1s24IyIZft5FDB5/LpHxCYEOT+rZSXkAUURERKS5y8/cz6rPPuKHzz+hpCCf+I6dOfum20kafQauENVDN1dKpkVERETqyFrLjjXf893H8/2lHED35FMYPGEynfsP0iQrLUBDJNP7GmCfIiIiIo1GWXExa7/+grRPPiB71w7ckVEMn3ohg84+h6hWrQMdnpxE9T6xu7X27PreZ1Pk9XoZMmQIkydPrlr26KOPUlxcXPU6IiKiwY6fmJhIZmZmnd9/3XXXVU1jPmjQID7//PN6jO74zJo1i9tvvx3wzyB57bXXcv3111NTW8fc3Fyefvrpqtdbt25l9uzZJyVWERFp3vZv3cxnzz3Js7dcwxcvPoMrJISJt97BzU/PYsyV1ymRboFU5tFAHnvsMZKSksjPz69a9uijj3L11VcTFhZWr8eqqKjA6TyxL6XX68XhOLg1z0MPPcTFF1/MwoULuemmm9iwYcMJHeNEWWu55ZZb8Hg8/Pe//63xT2cHkulbb70V+CmZvvLKK09muCIi0kx4ystY/81iVn36IXs2/ojTFUzvUWMZfPY5mupb6n9kWmDnzp188MEH3HDDDVXLHn/8cXbv3s2ZZ57JmWeeWbX8j3/8I4MGDWLkyJHs23d4hUx2djbnn38+AwcOZOTIkXz//fcAzJw5k2nTpjFq1CimTZtGVlYW48ePp1+/ftxwww0Hjdq++uqrnHLKKQwePJibb74Zr9cL+EfGf/vb3zJo0CC++eabGs/n1FNPZdeuXcDBo8QAkydPZtGiRVX7O9L5bNq0iZEjRzJgwADuvffeg0bkH3roIYYPH87AgQOZMWPGUa/rr3/9a7Kysnj55ZcJCgpi5syZPPzww1Xr+/fvz9atW7nnnnvYtGkTgwcP5q677uKee+7h66+/ZvDgwfzrX/9i69atjBkzhqFDhzJ06FCWLl161OOKiEjLlLljGwtnPcdzt1zLx0//i9LiIs689kZufuZlJv7yN0qkBWiAZNoYE26MadHdx3/zm9/wz3/+k6Cgny7vr3/9a9q3b8/ChQtZuHAhAEVFRYwcOZJVq1YxduxYnn/++cP2NWPGDIYMGcL333/PAw88wDXXXFO1bu3atSxYsIA5c+Zw3333MXr0aNasWcMFF1zA9u3bAf9sfW+88QZLliwhLS0Nh8PBa6+9VnX8ESNGsGrVKkaPHl3j+Xz88cecf/75xzzvms5n+vTpTJ8+nR9++IGOHTtWbf/pp5+yYcMGli9fTlpaGqmpqXz11VdH3Pfs2bNZuXIlr7/++jFH4R988EG6d+9OWloaDz30EA8++CBjxowhLS2NO+64g9atW/PZZ5+xcuVK3njjDX79618f89xERKRl8JSVsubLz5nz59/x0p23kfbph3QZOIRL/vQAP3/k3ww95zzcDVimKU3PCZd5GGOCgMuBq4DhQBkQYozJBD4AnrXWbqzFfjoBLwNt8Peqfs5a+9iJxLb3gQcoS193Irs4TEhSH9r+4Q81rp8/fz6tW7dm2LBhVSO2NQkODq6qqR42bBifffbZYdssXryYd955B4CzzjqLrKysqtKRqVOn8v/bu/PwOKoz3+PfV2trl6zNtmxZssG2bINjYzZjMNsFO4lhiOcmzsIAvjPOM8FzuZObfWDi5CaQZPJMZu5MyJiHgUAWwVxIgMkQlsRAYgI2Fhi8yLstb/KmfVer+9w/qtW0ZNla3HZr+X309FNVp05VvXW61P2qdKoqJSUFgD/84Q/86le/AuBjH/sYOTk5APz+97+noqKCyy+/HIC2tjYKCrz+XPHx8SxfvvyM8X35y1/mG9/4BocPHz7rmev+9uett97iueeeA+Azn/kMX/rSlwAvmX7llVeYN28eAM3NzezevZvrrrvutHXPnz+fHTt2sHHjRq655pp+Yzkbv9/P6tWrw39c7Nq165zWJyIiI9+JA/vYsu4VKv/4Gh2tLeRMKGLx51Yya/FNpGZmxTo8Gcai0Wf6NeB3wNeBrc65IICZjQNuAL5vZr92zv28n/V0Af/bOfeumWUAFWb2qnNuexRivGDefPNNXnjhBV588UXa29tpbGzkc5/7HD//+em7n5iYGO73Gx8fT1dX16C2lZaW1m8d5xx33XUXDz300GnzfD7faf2kI3X3mf6Xf/kXVq5cSUVFBQkJCQSDwXCd9vb28Phg98c5x9e//nU+//nP9yj/8Y9/HD6r/eKLLwIwc+ZMvv3tb/PJT36Sl19+mdmzZ581lrP50Y9+RGFhIe+//z7BYBCfT/f+FBEZi9qbm6l883W2rnuVEwf2Ep+QwMVXXsOlNy9hUtkc3dZOBiQayfTNzjl/70LnXC3wLPCsmSX2txLnXDVQHRpvMrNKoAgYcjJ9tjPI58tDDz0UTlxff/11fvjDH4YT6YyMDJqamsjLG/jTj6699lp+8Ytf8MADD/D666+Tl5dHZmbmafWuu+46fvnLX3L//ffz29/+lrq6OgBuuukmbr/99nD3htraWpqampgyZcqAY1i9ejWPPfYYL7/8MiUlJTz88MMEg0GOHDnCxo0b+13+qquu4tlnn+VTn/oUTz31VLj81ltv5YEHHuCzn/0s6enpHDlyhMTERO69917uvffe09azcOFCfvKTn/Dxj3+cN954g5KSEn7zm98A8O6777J//37gw3bu1nu6oaGBSZMmERcXxxNPPBHuQy4iIqOfCwY5uO0Dtr72Krs3/omA309+yVRuvOfzzFx0PSnpGbEOUUaYc06muxNpM/tn4H+5Pu5X1leyfTZmVgLMAzaca3zDyapVq1iyZEm47/RArFmzhpUrV3LppZeSmprKE0880We9b37zm3z6059m9uzZLFy4kOLiYgBmzZrFd77zHW655RaCwSCJiYn8+Mc/HlQybWbcf//9/OAHP+B3v/sdpaWlzJo1i7KyMubPn9/v8t13Mfnud7/LkiVLyMry/l12yy23UFlZydVXXw14FzD+/Oc/D3dD6cuyZcs4deoUS5Ys4dVXX+XJJ59k9uzZXHnllUyf7l0IkpubyzXXXMOcOXNYunQpDz74IPHx8cydO5e7776bL3zhCyxfvpwnn3ySJUuWDOgMv4iIjGx1x46y/Q/r2PbG72k6dZLktDQuufEW5txwC4Wl02IdnoxgdqZ79Q56RWbfAeYCK5xzLWZ2K/D3zrlBdXA1s3TgDeC7zrlf9TF/FbAKoLi4+LKqqqoe8ysrKykrKxviXsj50NraSkpKCmbGU089RXl5Oc8//3yswxoyHWMiIiNDR2sru95ez7Y3fseRHdvBjJJL5zFr8U1cdPlVJCYlxzpEGUHMrMI5t6B3edTuM+2cu9/MPgO8bmadQDPwtcGsI9Qd5FngF30l0qHtPAI8ArBgwYLo/CUg51VFRQWrV6/GOUd2djaPPfZYrEMSEZFRKhgIULVlM9v/sI4977xNV2cHORMnsejTdzHr2hvIyB14V0uRgYhaMm1mNwF/BbQAE4CVzrmdg1jegH8HKp1z/xituCT2rr32Wt5///1YhyEiIqOUc44T+/ey/Y+vsePNN2htqMeXls7sxTcye/HNjL9oui4mlPMmmk9A/DvgAefcejO7BHjazL7onFs3wOWvAe4EtpjZ5lDZN5xzL0YxRhERERklGk4cZ8ebb1C5/nVqDh8kLj6BqfMvZ9Z1N1A673ISEvu9/4HIOYtmN48bI8a3mNlSvC4bCwe4/HpAfzaKiIjIGbU21LPz7fXsWP8GR3dVAjBxxixu/st7mX71It2NQy64aDy0xc5wB4/qUNePM9YRERER6U9Hawt73nmbHX/6A1UfvIcLBskrLmHRp+9i5sLryCoojHWIMoZF5aEtZvYs8Lxz7mB3oZklAVeb2V14D3b5aRS2JSIiImOAv72dvRUb2PnWH9m/uYKA309GXj6XL/sEMxddT35xSaxDFAGik0wvAVYC5WZWCtQDPiAeeAX4J+fce1HYzojy3HPPcccdd1BZWcnMmTM5cOAAZWVlzJgxA+ccaWlpPP7448yYMYONGzeyatUqwLuIYs2aNdxxxx0x3gMREZELy9/Rzv7NFex8az37KjbS1dlBWs445t68lBkLr2XCxTN1IaEMO+eUTJtZvHOuHXgYeDh0a7s8oM05Vx+F+Eas8vJyFi1aRHl5Od/61rcAmDZtGps3bwZg7dq1PPjggzzxxBPMmTOHTZs2kZCQQHV1NXPnzmXZsmUkJETz+lAREZHhx9/ezv7Nm9j59pvsf/cd/B3tpGRmMXvxTcxYeC1FM2cRFxcf6zBFzuhcs7UfmFmyc261mb0KfMk5N+bvgdbc3Mz69et57bXXWLZsWTiZjtTY2EhOTg4Aqamp4fL29nb91S0iIqNaZ3sb+9/bxK631rNv8ya6OjpIycyi7NrrmX7VIibPuoS4eCXQMjKcazJdB3SExr8K/JOZHcC7pV31Oa57xHr++edZsmQJ06dPJzc3l4qKCnJzc9m7dy8f+chHaGpqorW1lQ0bPnxa+oYNG1i5ciVVVVX87Gc/01lpEREZVdqam9hXsZHdG//EgfffJeD3k5qVzezrbmL6VYuYNGu2zkDLiHSuGdsVwKsAzrl3gRvMbDnwkpn9CviBc67tHLcxZH/8j12cOtQc1XXmTU7n2k9OP2ud8vJy7rvvPgBWrFhBeXk5q1ev7tHN4+mnn2bVqlW89NJLAFx55ZVs27aNyspK7rrrLpYuXYrP54tq7CIiIhdSU+0p9m7yEuhD2z7ABYNk5OYz9+alXHzFQibOLFMCLSPeOSXTzrnbzCynezr0FMOdwE+A7wB/ZWZfd8797NzCHDlqa2tZt24dW7ZswcwIBAKYGffee2+Perfddhv33HPPacuXlZWRnp7O1q1bWbDgtMe/i4iIDFvOOWoOH2Tvpg3seectju3dDUDOhCIuv205F1+xkMKpF6k7o4wq59yXwDlXB2BmbwKlwDbgbeBuYAdwn5ld65xbda7bGqz+ziCfD8888wx33nkna9euDZctXryYQ4cO9ai3fv16pk2bBsD+/fuZPHkyCQkJVFVVsWPHDkpKSi5k2CIiIkMSDAQ4urOSve9uZM87b1F/zOvlOf6i6Sxa8RdcdPlVjCuarARaRq1odsxdBWzv4+Esf2NmlVHczrBWXl7OV7/61R5ly5cv56GHHgr3mXbOkZSUxKOPPgp4ifX3vvc9EhMTiYuL4+GHHyYvLy8W4YuIiPSro7WFA++/y96Kjex/bxPtzU3ExSdQPOdSFnz8DqZddiXp43JjHabIBWHn+mBCMyseQLV4oMY513hOG+tlwYIFbtOmTT3KKisrKSsri+ZmRHrQMSYiY1Fd9RH2vbuJfe9u5HDlVoKBAL6MTKbOW8C0y65gyqXzSY64O5XIaGNmFc650/rgRuPM9BOAA872/xuH9wTEJ6OwPRERETnPuvx+DlduZf+777DvvXfC3TfGFU3mso/9GVMvu4KJ02fqAkIZ86LRZ/qGaAQiIiIisdVw4jgH3q9g/+YKDm55H39HOwmJSUyefQnzP3o7U+ctIKtgfKzDFBlWdDNjERGRMaqrs9M7+7y5ggObK6g9ehiAzPxCZl13I1PnX87k2ZeQmKxbtYqciZJpERGRMcI5R82hKg588B5VH7zH4e1b6fJ3Ep+YyORZl3DpzUspnXcZOROKdPcNkQFSMi0iIjKKtdTXcXDr+1R98B4HPniPlrpawOv7fOnNS5gydx6TZ+nss8hQRTWZNrMbnXPruofRXLeIiIj0r7O9jcOVWzm4ZTNVW97n1MEDAPjSM5hyyUeYMnceUy6ZR2ZefmwDFRklon1m+ofA/IjhmPXcc89xxx13UFlZycyZMzlw4ABlZWXMmDED5xxpaWk8/vjjzJgxI7zMwYMHmTVrFmvWrOFLX/pSDKMXEZGRosvvp3r3Dg5t+4CDWz+gevcOgoEA8YmJFM2czbWfuZviOXMpKJ2qO2+InAfnq5vHmO9oVV5ezqJFiygvL+db3/oWANOmTWPz5s0ArF27lgcffJAnnngivMwXv/hFli5dGotwRURkhAh0dXF83x4ved72AUd3VtLV2YFZHAWl01jw8TsovuQjTJxRRmJScqzDFRn11Gf6PGhubmb9+vW89tprLFu2LJxMR2psbCQnJyc8/dxzz1FaWkpaWtqFDFVERIa5cPK8fQuHt2/hyM5K/O1tAOQXl3DpTbcyec5cJpXNxpeWHuNoRcYeJdPnwfPPP8+SJUuYPn06ubm5VFRUkJubG36ceFNTE62trWzYsAHwku/vf//7vPrqq/zwhz+McfQiIhJLXZ2dHNu7iyM7tnNo+xaO7qzE39EOQO6kYmYvvpFJZZcwefYlpGZmxThaERnVyfRrP32EE1X7orrOgilTueHuVWetU15ezn333QfAihUrKC8vZ/Xq1T26eTz99NOsWrWKl156iTVr1vC3f/u3pKfrjIKIyFjT2dbK0V07OLJjG4crt1G9ZycBvx8IJc/X38TkWZcwqWwOqVnZsQ1WRE4T7WS6OTRsivJ6R4za2lrWrVvHli1bMDMCgQBmxr333tuj3m233cY999wDwIYNG3jmmWf4yle+Qn19PXFxcfh8PlavXh2LXRARkfOoqfYUR3dWcmTHdo7s2M7Jqv04F8Ti4igsncZHbvkYk8rmUDRzFikZmbEOV0T6EdVk2jl3XeRwsMzsMeDjwAnn3Jxzjae/M8jnwzPPPMOdd97J2rVrw2WLFy/m0KFDPeqtX7+eadOmAfDHP/4xXL5mzRrS09OVSIuIjALBYIBTB6s4umsHR3dVcnTndhpOHAcgITmZiRfP4MpPfJKiGbOYOH0mSSmpMY5YRAZruHXz+Cnwr8CTMY5jyMrLy/nqV7/ao2z58uU89NBD4T7TzjmSkpJ49NFHYxSliIicD23NTRzbvZOju3dwdGcl1Xt2hS8WTMvOYeKMMuYtuY2iGWXkl0wlPmG4fQ2LyGCZcy7WMfRgZiXAbwZyZnrBggVu06ZNPcoqKyspKys7T9GJ6BgTEU8wEODkwQNU795J9e4dVO/eSV31EQAsLo78KaVMnF7GxOkzmTi9jMz8Aj2iW2QEM7MK59yC3uVR/5PYzNKAdudcINrrFhERiQXnHI0nj1O9ZxfH9uykes9uTuzfS1dnBwCpWdlMuHgms6+/mYkXz6Bw2sUk+VJiHLWIXAjnnEybWRywAvgscDnQASSb2Sngv4C1zrk957qdiO2tAlYBFBcXR2u1IiIiYS31dRzbu5vj+3ZzbO9uju3ZRVtTIwAJiUkUlE7j0ptuZcLFM5hw8UyddRYZw6JxZvo14HfA14BtzrkggJmNA24Avm9mv3bO/TwK28I59wjwCHjdPKKxThERGbtaGxs4vm9P6LWbY/v20FxzCgCzOMYVTWLqZVcw4aLpjJ82nbziEvV1FpGwaHwa3Oyc85vZPzvn7usudM7VAs8Cz5pZYhS2M2DOOZ0hkPNiuF1jICKD01xXy4n9ezm+fw/H93nD7sQZIGdCEZNmzmb8tOkUTruIgpKp6q4hImd1zsm0c84fGm0ys/8EVjjnWszsVuDvnXPXRNQ5KzMrB64H8szsMPBN59y/DyYen89HTU0Nubm5Sqglqpxz1NTU4PP5Yh2KiPTDBYPUnzjGif37OHFgLycO7OPE/r20NtR7FczCiXNh6TQKp15EQek0klPTYhq3iIw8Ufs/lXPufjP7DPC6mXXiPcDla4Ncx6fPNY5JkyZx+PBhTp48ea6rEjmNz+dj0qRJsQ5DRCL4O9o5dbCKk1X7OXFgHyer9nPy4IHwLeni4uPJLZpM6Ucuo6BkKgUl0ygonap7OotIVEQtmTazm4C/AlqACcBK59zOaK1/oBITEyktLb3QmxURkfPMBYM0njrByaoDnDy4n1NVBzh5qIr66qOELtchKSWV/CklzF58I/lTplJYOo3cScUkJCXFOHoRGa2ieQXF3wEPOOfWm9klwNNm9kXn3LoobkNERMaA1sYGTh2s4tShKmoOVXHy0AFqDlXR2dYWrpNdOIG84hJmLryW/CmlFJRMJTO/UF38ROSCimY3jxsjxreY2VK8CxAXRmsbIiIyurQ1N1FzqIqaw4eoOXyQmsNVnDp08MO+zYAvPYO84inMuu5G8qeUkl9cSu7kYl0YKCLDQjTuM22uj1scOOeqQ10/zlhHRERGP+ccrQ311Bw+RO2RQ9QcOUTtkYPUHD5ES31duF6iL4XcSZOZOv9y8iZPIXfyFPKLS0jNytbZZhEZtqJyn2kzexZ43jl3sLvQzJKAq83sLrx7Uf80CtsSEZFhKhgI0HDiGLVHD1N75DC1R49Qe8RLoNtbmsP1klJSGFc0mZK5l5E7uZi8ScXkTi4mIzdfSbOIjDjRSKaXACuBcjMrBeqBFCAOeAX4J+fce1HYjoiIxJhzjramRuqOHqGu+gi11UfC43XVRwkGusJ1U7OyyZlQxPSrF5FbNJlxk4rJLZpM+jjdulRERo9o3Ge6HXgYeDj0cJY8oM05V3+u6xYRkdjoaG2l/thRL0k+dpS66qPUVx+l7thR2pubwvXi4hPILhxPzsQipl52BeMmFDGuaBI5EyeRkp4Rwz0QEbkwotFneiXwC+dcB/BRYDywBfjTua5bRETOD+cc7S3N1B87Sv3xY97wWLX3Ol7d4wJAgIzcfHImTGD6VdcwbuIkciYWkTOhiKz8QuLi42OzEyIiw0A0unnc55x7zMzWADcAbwGfMLMM4BPOuWNR2IaIiAxSMBCgqeYk9ceP0XDimDc8foz649U0HD9GR2tLj/rpuXnkFE5g2mVXkD1+IjnjJ5I9YSLZheNJTNaTP0VE+hKNZLozNPwocLVzLgBgZh/D6/7xiShsQ0REenHBIC31dTScPEHjyeM0nPBejSePUX/8OE01J3HBYLh+XHwCWQUFZBWMZ8LFM8kuHE9W4Xhyxk8kq3A8iUnJMdwbEZGRKRrJ9CEz+ylQgHfhYTOAc+6/zOw7UVi/iMiYFAwEaK6rofHkCRpPnQwNT3jD0HjA7++xTGpWNlkFhUycPpOsguvJKigkq6CQ7MIJpOfmEhenLhkiItEUjWT6bmA58CPgWTN7CdgGzOPDs9YiIhLBOUdHawvNNadoqjnlJcunTtB06iSNp07SVHOSpppTPc4sA6RkZJKZX0B+cQnTFlxJVn4hmQUF3jC/QN0xREQusGjczaMReBzAzP478Hm8BLsO+NS5rl9EZKTpkSjX1tBUc4rm2lM01dSEk+SmmlP429t6LBcXH0/6uDwy8/IpmjmbjNw8MvMKyMwv8IZ5+ST6lCyLiAwnUXucOIQT63+I5jpFRIaTQJeflvo6mmtraa6robm216uulqbaU3R1dPRc0Iy0rGwycvPILZpMyaXzyMjNIz03j4zcfDLz80nLzlE3DBGRESaqybSIyEjl7+ygpa6OlrpaWhpCw/o6mutqvfG6Wprramlrajxt2fiEBNJycknPGUf+lFKmzl9A+rg8L1nOySUjN4+0nBziExJjsGciInI+KZkWkVGry++ntaE+/GppqKO13hu21NfTWl9HS4M37H2bOACLiyMtK5u0nFwyCwqZOKOMtJxxpGWPIz1nHOnjckkfl0tKRqae6CciMkYpmRaREcMFg7Q1N9HW2EBrY4M3bGigtbG+17CB1oY6OlpOT5ABklJSScvOIS07h/wppaTNnUdaljftJcs5pOeM85LkuLgLvJciIjKSKJkWkZhwztHV0UFbcyNtjY20NTWGE+W2pkZv2NhIa1No2NhAe1MTzgX7XJ8vI5PUzCxSs7LILy4hNWsuqVnZpGXnkJqVQ1pWNqlZ2aRmZemOFyIiEjVKpkXknAW6/LQ3N4deTbS3NNHW1OSNNzfR1tRIe1MTbc1NtIeS5vamJrr8Z7h7phkp6RmkZGSSkplFzoQiimbMIjUri5RM75WakfXhdEYm8Qn6OBMRkQtP3z4iAnj9iztamulobaG9uZmOlmbaW1voaA6VtTR/WN7cFJ5ub2k+7RZvkSwuDl96BinpGfgyMsksKKSg9CJSMjO9ZDkjE19GKHFOzyQlMxNferruaiEiIiOCkmmRUSDQ1UVnWyudba10tLbS0drSczw07Ghp8Ybdr5YPh2c8SxySkJRMcloavrR0fOnpZOYXUFAyFV96Oslp6aSke0mwLz0j4pVOcmqaLs4TEZFRS8m0SIwEAwE629vobGsLJcJtdLa34Q8Nw2VtrXS0teEPDbuT5u5kubOtja7Ojn63F5+YSHJqmvdK84YZufnhcV9aemh+Ksnp6aHp9HBCnJCUdAFaRUREZGRRMi3Sj0BXF/6Odu/V3kFXZwf+9nb87W34Ozrwd7TT2d49v3vY5pV1DzvaQ0lyd1kbAb9/QNuPi08gKSWFpJRUb+hLITUzi+zCCSSnppGU6pUnp6SSlJoWGqaGE+OklFSS09JJSNQ9jkVERKJtWCXTZrYE+GcgHnjUOfe9GIckw5RzjmCgi67OztCrg67OTvwdHeHxrs5O/J2h6Y6Iso72cL3Iob+jPVQvlCx3duBv7yAY6BpUbAlJyST6fCT5fCT6Ukj0+UhM9pGWlU2SL4XElFSSfD5v3JcSSpRTQvO8YWTiHJ+YqG4SIiIiw9SwSabNLB74MfDfgMPAO2b2gnNue2wjk764YJCuLj8Bf+jV5aer0xt2l3X5O73pTr9Xt7OzR72uzk4C/k66/H5v2On3Et/wdKe3jtDQH5E4Bzr9Z7xFWn8SEpNISE4mISmJxORkEpKSSUhOJjHZR2pmlpcMJ/u8ecnJJCZ3T/tCiXGv6VCynOTzkZCcrAvnRGREcs7hHLjucQhNe+Uf1utZ1rs+EfMj50UuS3he9zp61Xcf1u9zfniei1hvz9hOi7OfuOgx7/Q26N7Z07YT2X7hOpExnB5z5Hb62r/T4u0d82n7cno8RMbfvU44+3b6en/dWbbTx3Fy5ve6n22csT0il3X4EuP59u1zGE6GTTINXAHscc7tAzCzp4DbgTGRTHsHSpBgV4BgoItAIECwq4tAVxfBQIBAl59gj7LQsMurG/D7vTO1/u4yb36gyx9epvsV7PKfYbq7bmja313mD493bycYCERlv+MSEohPTCI+MZG4hEQSQuPxSUnEJyQSn5xKQnoWqYlJXr2k0Pzu6cQk4pKSek6HhglJScQnJXuvBG/ZuMRELC6ux5dAd/v3nA6/M6fVDTpox9HeXc+BawPaunD4ey3fc33dH8g9vph6bT9y2z2X7bWuM8Tc/eHVe329t33aF0CosM8Ye6yv71jdILd5+rp7f+n1H19f+957W/TRTn3tQ+S+RS56pu313U6nt8NpMfW57bO0X19t38++R7bb2fYh8hjr6z3rsXw/x0D4yzIcQMQXYh8xnPnY6bnu3u9L7wSjr/ajV/mZ1tv7vTh9vT3b/Gzr7JEwnCXW3u/N6evs9b70tc7u6T7avs/3u8/2751IigyNGRhgZqEheGNAeJ5X1rsukdMR46FZWKhSnH24jrTk4XfCajgl00XAoYjpw8CVMYrljJ686yucbN8VxTU6IEjEV/95FId3KMaFxuMwixjHvHGL61mHZLAUIC78YwmR83uvK94rC68nvsf8nmXm7Xpn6AUEQq+zi1hA5DyyXkMRGUn0mzvadCQcAa6PdRg9DKdkekDMbBWwCqC4uPiCb98X5yM5fnLU1hdOYMNJrjdtWEQyGlHWXbdHYmofjltEsmy9l9GHykijd0xERORDqTb87iw1nJLpI0BkljopVNaDc+4R4BGABQsWXPB/Tn3y8W9f6E2KiIiIyDAVF+sAIrwDXGxmpWaWBKwAXohxTCIiIiIiZzRszkw757rMbDXwMl6H2secc9tiHJaIiIiIyBkNm2QawDn3IvBirOMQERERERmI4dTNQ0RERERkRLHe91UdSczsJFAVg03nAadisN2RTG02eGqzwVObDZ7abPDUZoOnNhsctdfgXYg2m+Kcy+9dOKKT6Vgxs03OuQWxjmMkUZsNntps8NRmg6c2Gzy12eCpzQZH7TV4sWwzdfMQERERERkiJdMiIiIiIkOkZHpoHol1ACOQ2mzw1GaDpzYbPLXZ4KnNBk9tNjhqr8GLWZupz7SIiIiIyBDpzLSIiIiIyBApme7FzJaY2U4z22NmX+tjfrKZPR2av8HMSiLmfT1UvtPMbr2ggcfQANrsi2a23cw+MLPfm9mUiHkBM9sceo2Zx8cPoM3uNrOTEW3zlxHz7jKz3aHXXRc28tgZQJv9KKK9dplZfcS8MXecmdljZnbCzLaeYb6Z2f8NtecHZjY/Yt5YPcb6a7PPhtpqi5n9yczmRsw7ECrfbGabLlzUsTWANrvezBoifv/+PmLeWX+nR6MBtNeXI9pqa+iza1xo3lg9xiab2WuhPGKbmd3XR53Yfp455/QKvfAeY74XmAokAe8Ds3rV+QLwb6HxFcDTofFZofrJQGloPfGx3qdh0mY3AKmh8b/ubrPQdHOs92GYttndwL/2sew4YF9omBMaz4n1Pg2HNutV/2+AxyKmx+Jxdh0wH9h6hvkfBX4LGHAVsCFUPiaPsQG22cLutgCWdrdZaPoAkBfrfRiGbXY98Js+ygf1Oz1aXv21V6+6y4B1EdNj9RibAMwPjWcAu/r4zozp55nOTPd0BbDHObfPOdcJPAXc3qvO7cATofFngJvMzELlTznnOpxz+4E9ofWNdv22mXPuNedca2jybWDSBY5xuBnIcXYmtwKvOudqnXN1wKvAkvMU53Ay2Db7NFB+QSIbppxzfwBqz1LlduBJ53kbyDazCYzdY6zfNnPO/SnUJqDPMmBAx9mZnMvn4Ig1yPYa859jAM65aufcu6HxJqASKOpVLaafZ0qmeyoCDkVMH+b0NyxcxznXBTQAuQNcdjQa7H7/D7y/Hrv5zGyTmb1tZn92HuIbjgbaZstD/656xswmD3LZ0WbA+x3qRlQKrIsoHovHWX/O1KZj9RgbrN6fZQ54xcwqzGxVjGIarq42s/fN7LdmNjtUpuPsLMwsFS/pezaieMwfY+Z1rZ0HbOg1K6afZwnRXqHImZjZ54AFwOKI4inOuSNmNhVYZ2ZbnHN7YxPhsPKfQLlzrsPMPo/335AbYxzTSLECeMY5F4go03EmUWNmN+Al04siiheFjrEC4FUz2xE6CznWvYv3+9dsZh8FngMujm1II8Iy4E3nXORZ7DF9jJlZOt4fF//LOdcY63gi6cx0T0eAyRHTk0JlfdYxswQgC6gZ4LKj0YD228xuBv4OuM0519Fd7pw7EhruA17H+4tztOu3zZxzNRHt9Chw2UCXHaUGs98r6PWv0TF6nPXnTG06Vo+xATGzS/F+J293ztV0l0ccYyeAXzM2uvn1yznX6JxrDo2/CCSaWR46zvpzts+xMXeMmVkiXiL9C+fcr/qoEtPPMyXTPb0DXGxmpWaWhHcw977y/wWg+2rQP8e7OMCFyleYd7ePUry/vDdeoLhjqd82M7N5wFq8RPpERHmOmSWHxvOAa4DtFyzy2BlIm02ImLwNr48YwMvALaG2ywFuCZWNdgP53cTMZuJdZPJWRNlYPc768wLwF6Gr4K8CGpxz1YzdY6xfZlYM/Aq40zm3K6I8zcwyusfx2qzPuzWMNWY2PnRdEWZ2BV7eUcMAf6fHIjPLwvsP7vMRZWP2GAsdP/8OVDrn/vEM1WL6eaZuHhGcc11mthqvoePx7gawzcy+DWxyzr2A94b+zMz24F1EsCK07DYz+w+8L+ku4N5e/2YelQbYZv8ApAP/L/SZetA5dxtQBqw1syDeB+z3nHOjPskZYJv9TzO7De9YqsW7uwfOuVoz+z94X0QA3+71b8BRaYBtBt7v41OhP3C7jcnjzMzK8e6kkGdmh4FvAokAzrl/A17EuwJ+D9AK3BOaNyaPMRhQm/093jUyD4c+y7qccwuAQuDXobIE4JfOuZcu+A7EwADa7M+BvzazLqANWBH6/ezzdzoGu3BBDaC9AO4AXnHOtUQsOmaPMbwTIHcCW8xsc6jsG0AxDI/PMz0BUURERERkiNTNQ0RERERkiJRMi4iIiIgMkZJpEREREZEhUjItIiIiIjJESqZFRERERIZIybSIiIiIyBApmRYRERERGSIl0yIiY4CZXW5mH5iZL/Q0tW1mNifWcYmIjHR6aIuIyBhhZt8BfEAKcNg591CMQxIRGfGUTIuIjBFmloT3WN12YKFzLhDjkERERjx18xARGTtygXQgA+8MtYiInCOdmRYRGSPM7AXgKaAUmOCcWx3jkERERryEWAcgIiLnn5n9BeB3zv3SzOKBP5nZjc65dbGOTURkJNOZaRERERGRIVKfaRERERGRIVIyLSIiIiIyREqmRURERESGSMm0iIiIiMgQKZkWERERERkiJdMiIiIiIkOkZFpEREREZIiUTIuIiIiIDNH/B4Cb54k2bhN4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "algos = [\n",
    "    (phi_euler, \"Euler\"),\n",
    "    (phi_euler_modified, \"modified Euler (Collatz)\"),\n",
    "    (phi_heun, \"Heun\"),\n",
    "    (phi_rk4, \"4th order Runge-Kutta\"),\n",
    "    (phi_ab3, \"AB3\"),\n",
    "    (phi_ab4, \"AB4\"),\n",
    "]\n",
    "\n",
    "def ODEF(x, y):\n",
    "    return y - x**2 + 1.0\n",
    "\n",
    "h   = .02\n",
    "x   = np.arange(0, 2, h)\n",
    "y   = (x + 1)**2 - 0.5*np.exp(x)\n",
    "eta = np.zeros(len(x))\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])\n",
    "\n",
    "ax0 = plt.subplot(gs[0])\n",
    "ax0.plot(x, y, label=\"exact solution\")\n",
    "ax0.set_yscale(\"log\")\n",
    "ax0.set_title(\"Solutions to $y' = y - x^2 + 1.0$\")\n",
    "\n",
    "\n",
    "ax1 = plt.subplot(gs[1])\n",
    "ax1.set_ylabel(\"$\\delta(x) = |\\~y(x) - y(x)|$\")\n",
    "\n",
    "\n",
    "for alg, alg_name in algos:\n",
    "    eta = integrator(x, y0, ODEF, alg)\n",
    "    ax0.plot(x, eta, label=alg_name)\n",
    "    ax1.plot(x, np.abs(eta - y), label=alg_name)\n",
    "\n",
    "ax0.legend()\n",
    "plt.xlabel(\"x\")\n",
    "ax0.set_ylabel(\"y(x)\")\n",
    "ax1.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "9fd7c159424b1d536c11d53764cce54c",
     "grade": false,
     "grade_id": "cell-3c99ac2c49457d60",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 6\n",
    "\n",
    "Study the accuracies of all approaches as a function of the step size $h$. To this end use your implementations to solve\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = y\n",
    "\\end{align*}\n",
    "        \n",
    "with $y(x=0) = 1.0$ for $0 \\leq x \\leq 2$. \n",
    "\n",
    "Plot $\\delta(2) = |\\tilde{y}(2) - y(2)|$ as a function of $h$ for each integration scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "0060c86da7630a17a6da4769bb0c728f",
     "grade": true,
     "grade_id": "cell-552ef9b066bc9b1c",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# We assume the definitions of the previous block as not to copy everything.\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])\n",
    "\n",
    "ax0 = plt.subplot(gs[0])\n",
    "ax0.plot(x, y, label=\"exact solution\")\n",
    "ax0.set_yscale(\"log\")\n",
    "ax0.set_title(\"Solutions to $y' = y - x^2 + 1.0$\")\n",
    "\n",
    "\n",
    "ax1 = plt.subplot(gs[1])\n",
    "ax1.set_ylabel(\"$\\delta(x) = |\\~y(x) - y(x)|$\")\n",
    "\n",
    "\n",
    "for alg, alg_name in algos:\n",
    "    eta = integrator(x, y0, ODEF, alg)\n",
    "    ax0.plot(x, eta, label=alg_name)\n",
    "    ax1.plot(x, np.abs(eta - y), label=alg_name)\n",
    "\n",
    "ax0.legend()\n",
    "plt.xlabel(\"x\")\n",
    "ax0.set_ylabel(\"y(x)\")\n",
    "ax1.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "0e3fcb92129c0b3ef4e2b8f99a38ddef",
     "grade": false,
     "grade_id": "cell-3070782f40a2599f",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 7\n",
    "\n",
    "Apply the Euler and the Runga-Kutta methods to solve the pendulum problem given by\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = \\vec{f}(x, \\vec{y})\n",
    "    \\quad \\leftrightarrow \\quad \n",
    "    \\begin{pmatrix}\n",
    "    y_0'(x) \\\\\n",
    "    y_1'(x)\n",
    "    \\end{pmatrix}\n",
    "    =\n",
    "    \\begin{pmatrix}\n",
    "        y_1(x) \\\\\n",
    "        -\\operatorname{sin}\\left[y_0(x)  \\right]\n",
    "    \\end{pmatrix}\n",
    "\\end{align*}\n",
    "        \n",
    "To this end the quantities $\\text{x}$ and $\\text{y0}$ in your $\\text{integrator}$ function must become vectors. Depending on your implementation in task 1 you might have to define a new $\\text{integrator}$ function that can handle this type of input.\n",
    "\n",
    "Plot $y_0(x)$ for several oscillation periods. Does the Euler method behave physically?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "9cb9ca1fb8e9c0722a7e6b5ec0ff0c2a",
     "grade": true,
     "grade_id": "cell-3a89944b55976666",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}