{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Methods\n",
    "\n",
    "\n",
    "# Lecture 4: Roots of Equations\n",
    "\n",
    "## Exercise solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some imports we will make at the start of every notebook\n",
    "# later notebooks may add to this with specific SciPy modules\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# scipy's optimization\n",
    "import scipy.optimize as sop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 4.1: </span>\n",
    "\n",
    "Complete the implementation of the method of successive approximations below to solve $x=\\mathrm{e}^{-x}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current iteration solution:  0.4065696597405991\n",
      "Current iteration solution:  0.6659307054401221\n",
      "Current iteration solution:  0.5137951132027094\n",
      "Current iteration solution:  0.5982209490817094\n",
      "Current iteration solution:  0.5497888689549504\n",
      "Current iteration solution:  0.5770716352569477\n",
      "Current iteration solution:  0.5615403562006408\n",
      "Current iteration solution:  0.570329875730725\n",
      "Current iteration solution:  0.5653389163484923\n",
      "Current iteration solution:  0.5681675528504282\n",
      "Current iteration solution:  0.566562684236262\n",
      "Current iteration solution:  0.5674726729169729\n",
      "Current iteration solution:  0.5669565140929677\n",
      "Current iteration solution:  0.5672492292377977\n",
      "Current iteration solution:  0.5670832110967041\n",
      "Current iteration solution:  0.5671773650126686\n",
      "Current iteration solution:  0.5671239655566297\n",
      "Current iteration solution:  0.5671542504764889\n",
      "Current iteration solution:  0.5671370745155531\n",
      "Current iteration solution:  0.5671468157234473\n",
      "Current iteration solution:  0.5671412910553173\n",
      "Current iteration solution:  0.5671444243313883\n",
      "Current iteration solution:  0.5671426473141187\n",
      "Current iteration solution:  0.5671436551372928\n",
      "Current iteration solution:  0.5671430835570621\n",
      "Current iteration solution:  0.5671434077249292\n",
      "Current iteration solution:  0.5671432238752903\n",
      "Current iteration solution:  0.5671433281443767\n",
      "Current iteration solution:  0.5671432690088631\n",
      "\n",
      "Picard used 29 function evaluations\n",
      "\n",
      "Solution from Picard:  0.5671432690088631\n",
      "\n",
      "Check this against the solution from SciPy: sop.newton(f, 0.9)= 0.5671432904097843\n"
     ]
    }
   ],
   "source": [
    "def picard(f, x, atol=1.0e-6):\n",
    "    \"\"\" Function implementing Picard's method.\n",
    "    \n",
    "    f here is the function g(.) described in the lecture and we are solving x = g(x).\n",
    "    \n",
    "    x is an initial guess.\n",
    "    \n",
    "    atol is a user-defined (absolute) error tolerance.\n",
    "    \"\"\"\n",
    "    # record number of function evaluations so we can later compare methods\n",
    "    fevals = 0\n",
    "    # initialise the previous x simply so that while loop argument is initially true\n",
    "    x_prev = x + 2*atol\n",
    "    while abs(x - x_prev) > atol:\n",
    "        x_prev = x\n",
    "        x = f(x_prev) # one function evaluation\n",
    "        fevals += 1\n",
    "        print('Current iteration solution: ',x)\n",
    "    print('\\nPicard used', fevals, 'function evaluations')\n",
    "    return x\n",
    "\n",
    "\n",
    "\n",
    "def g(x):\n",
    "    return np.exp(-x)\n",
    "\n",
    "\n",
    "print('\\nSolution from Picard: ', picard(g, 0.9, atol=1.0e-7))  # 0.9 is our initial guess\n",
    "\n",
    "\n",
    "# let's check our answer against a SciPy function: sop.newton.\n",
    "def f(x):\n",
    "    return x - np.exp(-x)\n",
    "\n",
    "print('\\nCheck this against the solution from SciPy: sop.newton(f, 0.9)=', sop.newton(f, 0.9))\n",
    "\n",
    "# NB. if we tighten up the atol toelrance in our call to picard the answer gets closer to SciPy \n",
    "# but of course takes more iterations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "## <span style=\"color:blue\">Exercise 4.2:</span>\n",
    "\n",
    "\n",
    "Let's consider an example:\n",
    "\n",
    "$$2\\,x + x \\, \\sin(x-3) = 5\\;\\;\\;\\; \\text{for}\\;\\;\\;\\; x \\in (-10,10).$$\n",
    "\n",
    "\n",
    "By means of visual inspection (i.e. do some plotting) find a subinterval $(a,b)$ such that \n",
    "\n",
    "\n",
    "1. there exists an $x^* \\in (a,b)$ such that $f(x^*) = 0$, and \n",
    "\n",
    "\n",
    "2. $f(x)$ is [monotonic](https://en.wikipedia.org/wiki/Monotonic_function).\n",
    "\n",
    "Where we define $f$ such that the solution to the above equation is a root - i.e. move all the terms on to one side and set equal to zero. Below we make the choice to define $f$ as $f(x) = 2\\,x + x \\, \\sin(x-3) - 5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAHJCAYAAADwyRpfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde5xN9f7H8dcX4zKuySVUxEGcVBrJ5SippHSqU9Kv26GIpCJ0oXLLJRE6pIOKikqdroSSogsN4xbR5VTklktEmYlhvr8/vpsz99kzs2evtfd+Px+P/dhmrbXX+uy1194fn/X9ru8y1lpERERERETEe8W8DkBEREREREQcFWgiIiIiIiI+oQJNRERERETEJ1SgiYiIiIiI+IQKNBEREREREZ9QgSYiIiIiIuITKtBERERERER8QgWaiIiIiIiIT6hA8wljzExjzLwi3sY8Y8zMotxGYDtF/l4kb8aYYsaYqcaYX40x1hjTNpdl7zTG/GiMORp4zUnGmF3GmHr52N5/jDH9QhJ8PgRzvIXr2BeJFcpZEmqxkrNEgmGstV7HEBMCSaZLNrOaWmvXGmMq4j6P34owhnnAXmtt16LaRmA7+X4vxpglwAZr7T1FFlgECOV+MMZcBbwFtAV+BPZZa49ks9yZwAbgBuBL4HdgCFDFWnt7PrbXBFgKnGGtPVDY+POx3TyPt3Ad+yLRQjkrz9csQTlLOUukiJTwOoAY8xFwW6ZpewGi6cchmt5LMIwxJbNLIj7wF2CntXZZHstdjUuwbwMYY+KB7sDf87Mxa+16Y8yPwK3AMwWIt0DCdbwV5efs42NIYptyVhTy8e9NTOSsYOX0Ofnp81NeLELWWj3C8ABmAvPymg9UBXYCg9PNOxv4E+iUbpoBHgR+AFKA9cCt6ebHB9b5B7ALGBRY/8w84lwC/Bt4GtgfeIwFiqVbphQwMbDeP3FnsP6W03sNrHMKMAqX3HcD446vM7C8zfSok0N8HYDPAnHtAz4AGqWbf2Egnj+AA0AicFYu7zfX9eWyj54NvIc9wMpg9ktey+RzPwTzGaRfz+Yc1vNdpuXeBjoBvxJoYQ8sdwNwGKidbtrTuOOverppg4HPC/lduQJ3RrRE4O/6gdieTbfMSGBRDsdbnsc+eXx/cvucs4k3qH0T5DGU1/G9hFy+S4FlygIvpXv/Awvy/vWI7Ufm71VO81HOUs6K8ZwVzDET5LGQ0+dUoM+P/OdS3+XFvPZZPvZ9xOVFTzYaiw+CTHaBf18OHAFaAmWAr4EZmZYfCXwbOHjPAG4GDgEdA/OnANsD6zoLeAM4SHDJ7ndgEnAm0BmXNPqlW+ZpXELuCDQCpgcO+hrZvdfAOg8Aw4EGgXUeBW4KzK8ILANeAE4JPIrnEN/1gUd93H8CXgf+C5TEtQjvD3wx6wXivznzlznY9QWxj54KbKNRMPslr2XyuR/y+gwqAsOArYH1VM1hPdUCx9EjgeUqBNb9YablDJAETA/8PQD3I1g/03IdcMdumUJ8V8oBqUCLwN934n6sv0m3zBfAIzkcb3ke++Tx/cntc84m3qD2TZDHUK7HI3l8lwLL/BvYAlwG/BV4LfCafL1/PWL7kfl7ldt8lLOUs2I4ZwVzzATz2eXyORXo8yP/udR3eTGvfZaPfR9xeTHsG4zVBy4BHA18eY4/FmSanz5BTMT1wZ4ROBjLpZtXFlfdt8m0jYnA/MCX8jBwS7p55YDfCC7ZfUfGM1GPAtvSbfsI8M9084vjzjiMyOG9LAGWZ9rOIuC5TMtMLsB+LQscA/4GVMadHbqoEJ/TifXlsY++yuZ1ee2XYJbJcz8Es57AtAHkcBYy3TKlcT/gbdNNewd4MZtl2weWfRj3I3p+NsucHfgM6hXy+5IIDAz8ezbu+oIU3H8K4gPvv3Xm4y2YY588vj+5fc65xJvnvsnrGArmeMzruxR4r0eA/8u0jv35ff96xPYD5SzlLOWs/HwOeR4zeX122X1Ohfn8AtOCyqVESF7M7njPa98ToXlRoziG16fAueke3XNZ9iECXz5c0voj3bzGuB+phcaYP44/gF64s3D1cGfnlh9/QeD164OM80sbODIDlgO1jDEVAuuOw511Ob7uY4FlGueyzq8y/b0DdyYsX4wx9YwxrxhjfjDGHMQ1VRcDTrfW7sMl2g+MMe8bY/oZY04r6PryCGVVpr+D2S8F3XeZhWo94M5UlwDWpptWBtddIgNr7Ye4LgcjgM7W2pXZrC8l3TqyMMbckv6YNca0ySGuJbgLxQEuAhYAKwLTWuN+9Fdk87pgjv28vj/pZf6csxXkvsksy7qDPB5z+y4dPzZO7Btr7SHcBfXH5ef9S2xTzvof5azslwlGLOQsyOOYCfKzyynnFOTzg+BzqS/zYj6O96jLixokJLySrbX/DXLZOsBpuDM7dXFnQY47Xlj/Hfg50+tSgZMKEWNeTODZZjMvu2nHpWazbEFOEMzFdYPpGXg+CmzEJXestbcbYybimqivBkYaY6611n5QkPXl4lCmv4PZLwXdd5mFaj3g/tO1xWYcvWwv2RxDxph2wDmB7e/KYX2VA897cpj/HhmP5e05LLcE6G2MaQyUx/1oLwEuDqx7mbU28zEF/9s3ucnr+5Ne5s85W0Hum8yyW3cwx2Nu36Xcjo3j8vP+JbYpZ2VcVjkr6zLBiIWcBXkfM8F8djnlnIJ8fhBkLjXG+DUvBnu8R11eVAuaDxlj4nBN0e/hmvyfNcakP1uwkf9dfPnfTI8tuO4lqUCLdOssizvzFIwLjDHp/6PbAthhrT0YWPcRXPeM4+sujrv2YGM+32p6R3DN8zkyxpyM65M8ylr7kbV2E+4HJ8OJBmvtOmvtGGttW9wPUZfCrC9IweyXYJbJcz8EuZ5gnUvGM5EAa8h0VtMYcw5u+ON7cd1JRuewvrNwx0q2P8bW2t8zHa8p2S2Huyi4FO6C3c8DZwaX4JJK28C/sxPMsZ/X9ydf8rFv8lpPKI7H4++/ebr1xlOE719EOSt7yln52law/JqzchXizw6C36fB5lLf5cUQ7rOIzItqQfOnx3FNs5fiLmLsALxsjLnYWptmrf3dGDMOGBdISp/i+ti2ANKstdOMMc8DY4wxe3BNvYPJ+0f0uJrARGPMFKAJ8ACuiRpr7SFjzLPAE8aYvcBPwP1AddxF3gW1GWhujKmDu9Zhn7U2LdMy+3Fnyu40xmwFauFG6zoKYIw5A3eW5T3cmZa6uP7lz+awzVzXlx/B7Jcg912e+yHEn8G5uKG00/sAd+ycbK391RhTG3edyHhr7QvGmBXAV8aYttbaJZle2wZYmM8YsrDW/mGMWY0b/vjhwOTluDP0Z+CSTU6vy/XYD+b7E2yc+dw3eSn08Rh4/y/g3v9e3AXkj+JOxtnAMiF7/yIBylnKWTnuh1jIWUEI2WcHwe/TYHOpT/NiSPZZxOZF68GFb7H4IPghiy8i6wWwp+BGvxmYbprBnZk4XvXvwV0UeVlgfvohRXcDj5G/IYsn4y7Q3o8bUad4umXSD+16mOCGLJ6c2/7AjbyzHEgm96F62+H6Df8ZeL488B674n6Y3sIlusO4Zuongbhc3m+O68tjH2W5KDqv/RLkvgt2PwSzrVwvuA4cQweBa7OZtxzojev+sQmYmmn+HLJelFsa95+zFiH6zjwR2AfNMu37Q+k/02yOpTyPffL4/uT2OaebH/S+yccxlOvxmN3rsnn/5YCXA/tpFy4pLybj0Mp5vn89YvuR+bjKaT7KWcpZylnBHjP5/n0v7OcXWC7YXOq7vBjM8R7kvo+4vGgCQYkAYIxZgrsB5D1exyLeMcZ0wA3h29i6LhHBvKY3cI21tn2RBif5YowphRteeKy19imv4xEJJeUsAeUsyZ9IyIvq4igiWVhrFxpjngFOxf2IBSMVd/ZJPGSMaYrrt78C11//ocDzHC/jEhEpKspZkptIzIsq0EQkW9baf+VzeV275B/9gIa4vvprgQuttdu8DUlEpOgoZ0keIiovqoujiIiIiIiIT2iYfREREREREZ+I2i6OVapUsXXq1PE6DBERKYRVq1bttdZW9ToOv1GOExGJfDnluKgt0OrUqUNSUpLXYYiISCEYY3Tz7Gwox4mIRL6ccpy6OIqIiIiIiPiECjQRERERERGfUIEmIiISJsaYF4wxu40xG9JNq2yMWWSM+T7wfJKXMYqIiLdUoImIiITPTKBDpmkPA4uttfWBxYG/RUQkRqlAExERCRNr7afAvkyTrwFeDPz7ReDasAYlIiK+ErWjOObl4MGD7N69m9TUVK9DEY/FxcVRrVo1KlSo4HUoIhKbqltrdwJYa3caY6plt5AxpgfQA+D000/PdYXKcXKccpxI5InJAu3gwYPs2rWLWrVqUaZMGYwxXockHrHWkpKSwvbt2wGUwETEt6y104BpAM2aNbM5LaccJ8cpx4lEppjs4rh7925q1apFfHy8EleMM8YQHx9PrVq12L17t9fhiEhs2mWMqQEQeC7Uj5FynBynHCcSmWKyQEtNTaVMmTJehyE+UqZMGXUFEhGvvAd0Cfy7C/BuYVamHCeZKceJRJaYLNAAnVWUDHQ8iEg4GGNeBZYDDY0x24wx3YAngMuMMd8DlwX+Lux2CrsKiSI6HkQiS0xegyYiIuIFa+1NOcy6JKyBiIiIb8VsC5qIiIiIiIjfqEATERERERHxCRVoIiIiIiIiPqECLQrt37+f6tWr88MPPwT9mk6dOjF+/PgijAqmT59O3bp1KVGiBD179vRljCIi4m/KcSIS7VSgRaFRo0Zx5ZVXUq9evaBfM2TIEEaMGMGBAweKJKZvvvmGXr168dRTT7F161bi4+N9F6OIhM7NN8Obb3odhUQj5TgR8Yq1sHo1PPAAjBtXdNtRgRZlkpOTee655+jWrVu+XtekSRPq1q3LrFmziiSu9957j7POOot//OMfVKxYkZkzZ/ouRhEJja1b4dVXYccOryORaKMcJyJe2LkThg+HRo0gIQEmToTvvy+67alAizCnnnpqli4Q69evp3Tp0mzcuJH58+dTrFgxWrdunWGZN954g1KlSrFly5YT0/r06UO9evXYtWsXAFdffTWvvvpqyGNu0KABDz30EOvWrcMYQ9myZX0Xo4iEzrJl7rlVK2/jkMijHKccJ+In69ZB165QuzYMHQqnnAJTp8Ivv7jnoqL7oAX07Qtr14Z3m+ee6yrw/GjZsiUrV67MMK1v3750796dxo0bM3XqVBISErLclLJTp06MGTOGESNGMH36dMaNG8err77KF198QfXq1QFo3rw5I0aMICUlhTJlymR4/ahRoxg1alSusS1YsIA2bdpkmf7555/Tpk0b/vnPf9KtWzcGDhzI9u3bQx6jiPjDF19AfDycc47XkchxynHKcSISvA0bYOBAmDfP5bOePaFPH/jLX8KzfRVoEaZly5ZMmTLlxN/vvPMOa9as4fXXXwdgy5Yt1KhRI8vrjDGMGjWKjh07Uq9ePUaOHMnHH39M/fr1TyxTs2ZNUlNT2bFjR5Z+83fddRedO3fONbZatWplO71ChQr8+OOPtG7dmlNOOYX9+/cXSYwi4g/LlsEFF0AJZRjJJ+U45TgRL23dCoMHw4svQoUKMHIk9OoFJ50U3jiUPgPye5bPKy1atKB///7s27ePsmXLMmDAAAYPHszJJ58MQEpKyokzcZm1b9+e888/n0cffZS5c+dy/vnnZ5h//GxdSkpKltdWrlyZypUrFyjmDRs2cPToUc4999wijVFEvPfHH66l5uGHvY5E0lOOU44TkZwdPep+JwcPhmPHoF8/GDQICvizUGi6Bi3CJCQkULJkSZKSkpg4cSIlSpSgd+/eJ+ZXqVKF/fv3Z/vajz/+mHXr1mGtzTZ57Nu3D4CqVatmmTdq1CjKlSuX6+Ozzz7Ldrtr166ldu3aVKpUqUhjFBHvrVzpklumy29EgqIcpxwnEm6rV7teHw88AJddBt9+60Zo9Ko4A7WgRZxSpUrRtGlT5s6dy4svvsgrr7xCXFzciflNmzZl5syZWV63bt06rrvuOiZNmsT777/PwIED+eCDDzIss2HDBmrWrJlt0ihM94+1a9eeOLNYlDGKiPeODxDSooW3cUhkUo5TjhMJl6NH4fHHXTfGqlXhjTfg+ush0+Wj3rDWRuUjISHB5mTjxo05zosEffv2tcYY2759+yzzvvrqK1usWDG7d+/eE9M2b95sa9asaYcNG2attXb9+vXWGGM/+eSTDK/t0qWLveOOO0Ieb+vWre2QIUN8HaO1kX9ciPjBFVdY+9e/hm59QJL1QU7x20M5zj/5QzlOJPL8/LO1bdpYC9bedpu1+/Z5E0dOOc7zJJMhGDgN+ATYBHwN9AlMrwwsAr4PPJ+U17qiOXnNnDnTFi9e3G7YsCHb+S1atLCTJ0+21lr766+/2jPPPNP26NEjwzKdO3e2LVq0OPF3SkqKrVChgl2+fHlIY01LS7Ply5e3b7/9tm9jPC7SjwsRrx07Zm2lStbeeWfo1qkCTTkuMz/lD+U4kcjz7rvWnnSSteXKWfvyy97GEikFWg3gvMC/ywPfAY2BJ4GHA9MfBsbkta5oTl6XXXaZvfvuu3Ocv2DBAtugQQN79OjRoNc5efJke9lll4UivKD4McZIPy5EvLZhg8sqM2eGbp0q0JTjMvNj/sjMjzFG+nEhUlhpadYOH+7y1HnnWfvdd15HlHOO89UgIdbandba1YF//45rSasFXAO8GFjsReBabyL0TlpaGrt27WLMmDGsX7+eESNG5Lhshw4d6N27N9u2bQt6/XFxcUyaNCkUoQYlEmIUkfzRDaqloJTjlONEilJyMtx4oxul8bbb3P06093hwneMK978xxhTB/gUOAv42VpbKd28/dbaLHckMMb0AHoAnH766QlbtmzJdt2bNm2iUaNGRRB10VmyZAnt2rWjYcOGPP/887TS/4BCLhKPCxE/6doV5s+HXbtCd5G1MWaVtbZZaNYWPZo1a2aTkpKynReJv2XKcUUvEo8LkVDYsQOuusrdAmbMGBgwwCcDgZBzjvPlKI7GmHLAm0Bfa+1BE+RetNZOA6aBS15FF2H4tW3blrS0NK/DEBHJ0bJlrvXML4lPIodynIgUhe++g8svh717Ye5c6NjR64iC46sujgDGmDhccTbbWvtWYPIuY0yNwPwawG6v4hMRkax274bvv1f3RhER8YdVq+Bvf4M//oBPPomc4gx8VqAZ11T2PLDJWjs+3az3gC6Bf3cB3g13bCIikrPly92zblAtIiJe++QTaNsW4uPd9WbNIqyjvN+6OLYGbgPWG2PWBqYNAp4AXjfGdAN+Bm7wKD4REcnG8uUQFwcJCV5HIiIisWzxYvj736FuXfjwQ6hZ0+uI8s9XBZq19nMgp6sXLglnLCIiErzERDjnHChd2utIREQkVh0vzurVg48/hqpVvY6oYHzVxVFERCLPsWOQlAQtWngdiYiIxKpoKc5ABZqIiBTSxo3uIuwLLvA6EhERiUVffBE9xRmoQBMRkUL68kv3rAJNRETC7auv3H3OTjvNtaJFenEGKtBERKSQEhOhcmX4y1+8jkRERGLJDz+4+5yVK+cGBKlWzeuIQsNXg4SIiEjkSUyE5s11g2oREQmfX36B9u0hNdV1a6xd2+uIQkctaCIiUmC//w5ff60BQkREJHwOHXLdGnftgvnzoVEjryMKLRVoUWj//v1Ur16dH374IejXdOrUifHjx+e9oIhIOitXgrW6/kzCRzlOJLYdOwa33AJr1sBrr7keHNFGBVoUGjVqFFdeeSX16tUL+jVDhgxhxIgRHDhwoEhimjJlCmeccQalS5cmISGBzz77rEi2IyLhlZjonqMxQYo/KceJxLYHH4R334WJE10rWjRSgRZlkpOTee655+jWrVu+XtekSRPq1q3LrFmzQh7TnDlz6NOnD4MGDWLNmjW0atWKK664gp9//jnk2xKR8EpMhPr13SAhIkVNOU4ktj37LIwfD/fe6x7RSgVahDn11FOzdNNYv349pUuXZuPGjcyfP59ixYrRunXrDMu88cYblCpVii1btpyY1qdPH+rVq8euXbsAuPrqq3n11VdDHvP48ePp2rUrd955J40aNWLSpEnUqFGDZ599NuTbEpHwsdYVaOreKKGiHCciOVm61BVlHTvChAleR1O0NIpjOm3bts0yrXPnztx9990kJydz5ZVXZpnftWtXunbtyt69e+nUqVOW+b169eLGG29k69at3HbbbRnmLVmyJN8xtmzZkpUrV2aY1rdvX7p3707jxo2ZOnUqCQkJmEzDqXXq1IkxY8YwYsQIpk+fzrhx43j11Vf54osvqF69OgDNmzdnxIgRpKSkUKZMmQyvHzVqFKNGjco1tgULFtCmTZsM044cOcKqVasYMGBAhunt27dn2bJl+XrvIuIvP//sRtHSACGRQTlOOU4kUv38M9xwg+ux8corULy41xEVLRVoEaZly5ZMmTLlxN/vvPMOa9as4fXXXwdgy5Yt1KhRI8vrjDGMGjWKjh07Uq9ePUaOHMnHH39M/fr1TyxTs2ZNUlNT2bFjR5a+/XfddRedO3fONbZatWplmbZ3716OHTt2IkEeV716dT766KO837CI+Nbx68/UgiahohwnIpmlpMA//gGHD8M770CFCl5HVPRUoKWT29m++Pj4XOdXqVIl1/mnnXZagc4mZtaiRQv69+/Pvn37KFu2LAMGDGDw4MGcfPLJAKSkpGRJFMe1b9+e888/n0cffZS5c+dy/vnnZ5h//IxiSkpKltdWrlyZyoW4yCTz2U5rbZZpIhJZEhOhVCk4+2yvI5FgKMcpx4lEGmuhRw83YuN770HDhl5HFB66Bi3CJCQkULJkSZKSkpg4cSIlSpSgd+/eJ+ZXqVKF/fv3Z/vajz/+mHXr1mGtzTbB7du3D4CqVatmmTdq1CjKlSuX6yO7UauqVKlC8eLF+eWXXzJM3717d45JVkQiQ2IinHcelCzpdSQSLZTjRCS9Z5+FWbNg2LDoHbExO2pBizClSpWiadOmzJ07lxdffJFXXnmFuLi4E/ObNm3KzJkzs7xu3bp1XHfddUyaNIn333+fgQMH8sEHH2RYZsOGDdSsWTPbpFLQ7h8lS5YkISGBRYsWccMNN5yYvmjRIq6//vq83q6I+FRqKqxeDT17eh2JRBPlOBE5LikJ7r/fDQryyCNeRxNm1tqofCQkJNicbNy4Mcd5kaBv377WGGPbt2+fZd5XX31lixUrZvfu3Xti2ubNm23NmjXtsGHDrLXWrl+/3hpj7CeffJLhtV26dLF33HFHyON97bXXbFxcnJ0+fbrduHGjve+++2zZsmXt5s2bQ76twoj040IknNassRasfeWVot0OkGR9kFP89lCOU47Lr0g/LiS27NtnbZ061p52mrXpvu5RJ6ccpxa0CHTuuedSrFixLEMRg7vXS/PmzXnttdfo3bs3+/bto0OHDlx11VUMHjwYgLPOOosbbriBgQMHsnz5cgD+/PNP3n777SxnHEPhxhtv5Ndff2XEiBHs3LmTs846i/nz51O7du2Qb0tEwkMDhEhRUY4TiW3Wwu23w7Zt8NlnELgENaYYV7xFn2bNmtmkpKRs523atIlGjRqFOaLQad++PfXr1+eZZ57Jdv7ChQvp06cPGzdupHiQ45A+88wzvPvuu3z44YehDDWiRPpxIRJO3brBu+/Cnj1QlGMhGGNWWWubFd0WIpNynHJcfkX6cSGxY8IE6NfPPfft63U0RSunHKcWtAiRlpbGnj17mDlzJuvXr2fOnDk5LtuhQwd69+7Ntm3bgj6DFxcXx6RJk0IVrohEuRUroHnzoi3OYo0x5n6gO2CB9cDt1to/vY0qPJTjRATcaI0PPQTXXAN9+ngdjXdUoEWITz/9lHbt2tGwYUPefPNNTjrppFyXv++++/K1/h49ehQmPBGJIb//Dl9/DRoDIXSMMbWA+4DG1toUY8zrwP8BMz0NLEyU40Tk0CG4+WaoWhWefz62TwCqQIsQbdu2JS0tzeswRERYvdpdI6Drz0KuBFDGGJMKxAM7PI4nbJTjRKRfP/j2W1i0KDavO0tP90ETEZF8WbHCPWe6D7AUgrV2OzAO+BnYCRyw1ma4YMoY08MYk2SMSdqzZ48XYYqIFIm334Zp0+DBB+GSS7yOxnsq0EREJF8SE6FuXahSxetIoocx5iTgGuAMoCZQ1hhza/plrLXTrLXNrLXNsrvZsohIJNqxA7p3h4QEGD7c62j8QQWaiIjky/EBQiSkLgV+stbusdamAm8BrTyOSUSkSFnrirOUFJg9G0qW9Doif9A1aCIiErSdO2HrVhVoReBnoIUxJh5IAS4Bsh9HX0QkSjz3HCxYAP/6FzRs6HU0/qEWNBERCdrKle5ZA4SElrU2EfgPsBo3xH4xYJqnQYmIFKGffnIDg7RrB717ex2Nv6gFTUREgrZiBRQvDk2beh1J9LHWDgGGeB2HiEhRO3YMunSBYsVgxgz3LP+jAk1ERIK2YgWcfTaUKeN1JCIiEqkmTYLPPnPF2emnex2N/6heFRGRoKSlaYAQEREpnB9+gEGDoGNH14omWalAiyBdu3blqquu8jqMsGrbti333HOP12GICPD993DggAo0KTrKcyLRLS0N7rwT4uJg6lQwxuuI/EkFWgR5+umnmTVrVtDL16lTh3HjxhVhRKEzc+ZMypUrl2X6W2+9xejRoz2ISEQyS0x0zxogRIqK8pxIdJs+HT75BJ56CmrV8joa/1KBVhizZ0OdOu7Kxjp13N9FqGLFilSqVKlIt5GdI0eOhH2bx1WuXJny5ct7tn0R+Z/ERChfHs480+tIJCzCnONAeU4kmm3dCg88AJdcAt26eR2Nv6lAK6jZs6FHD9iyxd1lb8sW93cRJrD0XT/atm3L3XffzaBBg6hSpQrVqlVjwIABpKWlnZi/ZcsWHnjgAYwxmHRtyMuWLeOiiy4iPj6eWrVq0atXLw4ePHhiftu2benVqxcDBgygatWqtG7dGoCpU6fSoEEDSpcuTdWqVbn88ss5evToidfNmDGDxo0bU7p0aRo0aMCECXrRKLAAACAASURBVBNOxANw8OBBevXqRY0aNShdujSNGjVizpw5LFmyhNtvv51Dhw6diHXo0KEnYknf9WP//v106dKFk046iTJlynDppZfy9ddfn5h//Azl4sWLOeussyhbtiwXX3wxP/30Uwg/CZHYlJgI55/vRnGUKOdBjgPlOVCek+hkLdx1lxu9cfp0dW3Miwq0gnrkEUhOzjgtOdlND5PZs2dTokQJli1bxuTJk5k4cSJz5swBXJeJU089lcGDB7Nz50527twJwPr162nfvj1XX30169at46233mLt2rXccccdGdY9a9YsrLV89tlnvPTSSyQlJdG7d2+GDBnCt99+y0cffUSHDh1OLD99+nQGDRrE8OHD2bRpE0899RRjxoxhypQpAFhrueKKK1i6dCkzZsxg48aNjB8/npIlS9KqVSsmTpxIfHz8iVgHDBiQ7Xvu2rUriYmJvPvuu6xYsYL4+Hg6dOhASkrKiWUOHz7M6NGjeeGFF1i+fDm//fYbd911V0j3vUisSUmBdevUvTFm+CDHgfKc8pxEi9dfh/nzYeRIOOMMr6OJANbaqHwkJCTYnGzcuDHHeUEzxlp3QiDjw5jCrzsHXbp0sR07drTWWnvRRRfZFi1aZJh/6aWX2m7dup34u3bt2nbs2LEZlrntttvsHXfckWHamjVrLGB37dp1Yt1NmjTJsMybb75pK1SoYA8ePJhtbKeddpp96aWXMkybMGGCbdSokbXW2g8//NAaY3Lc9zNmzLBly5bNMv2iiy6yvXv3ttZa+91331nALl269MT83377zVaoUMFOnz79xHoA+80335xYZtasWTYuLs4eO3Ys220fF5LjQiRKff65+4l7553wbhdIsj7IKX57RGOOs1Z5rijznHKceGXfPmurV7c2IcHao0e9jsZfcspxug9aQZ1+uuvykd30MDn77LMz/F2zZk12796d62tWrVrFf//73xNnIMEV6QA//PAD1apVAyAhISHD6y677DJq167NGWecweWXX0779u257rrrKF++PHv27GHr1q307NmTXr16nXjN0aNHT6x7zZo11KhRg0aNGhX4/W7atIlixYrRsmXLE9MqVqxIkyZN2Lhx44lppUqVomHDhif+rlmzJqmpqfz2229Urly5wNsXiWUaICTG+CDHgfIcKM9J5Hv4Ydi7FxYsUBf5YKlAK6iRI11//PRdQOLj3fQwiYuLy/C3MSZDX/jspKWl0b17d+6///4s82qlG06nbNmyGeaVL1+e1atX8+mnn7Jo0SJGjx7NoEGDWLlyJcUD37Z///vftGrVKtvtHk9ghZHbOtJfe1CiRIls5+W1b0QkZ4mJULs2nHKK15FIWPggx4HyXHrKcxKJPv8cpk2D/v2haVOvo4kcugatoG65xR1xtWu7Kx1r13Z/33KL15GdULJkSY4dO5Zh2nnnncfXX3/NX/7ylyyPMmXK5Lq+EiVK0K5dO0aPHs1XX33FoUOHmDdvHtWrV6dWrVr88MMP2a73+HZ37tzJpk2bgo41s8aNG5OWlsby5ctPTDt48CDr16+ncePGwewSESmgxES1nsWUCMhxoDwn4mdHjrjzPLVrw7BhXkcTWdSCVhi33OK7ZJVenTp1+Oyzz7j11lspVaoUVapU4aGHHqJFixbcdddd9OzZk/Lly/PNN98wd+5cpk6dmuO65s2bxw8//MCFF15I5cqV+eSTT/j9999PdOUYOnQo9957L5UqVeLKK68kNTWV1atXs337dgYOHMgll1zCBRdcwPXXX8+ECRNo0KAB//3vfzl06BDXXnstderU4c8//2TRokU0bdqU+Ph44uPjM8RQv359rrnmGnr27Mm0adOoVKkSjzzyCBUqVODmm28u0n0pEst++cX1drvvPq8jkbDyeY4D5TkRP3vqKdi0CebNg0wN1pIHtaBFseHDh7N161bq1atH1apVAdef/9NPP2Xz5s1cdNFFnHPOOQwcOJDq1avnuq5KlSrxzjvvcOmll3LmmWcybtw4nnvuOdq0aQNA9+7deeGFF3j55Zc555xzaNOmDdOmTeOMwFA9xYoVY8GCBbRu3Zpbb72VRo0a0adPnxP3nmnVqhV33XUXN910E1WrVuXJJ5/MNo4ZM2bQvHlzrr76apo3b05ycjILFy7M86yoiBScrj8Tv1KeE/GnzZvh8cfhH/+Ajh29jibymFD0mfajZs2a2aSkpGznbdq0qVAX8Up00nEhkr1Bg2DsWDh4EML9f0RjzCprbbPwbtX/lOMkv3RcSDhdey0sWuRa0MI8tlBEySnHqYujiIjkKjERzj47/MWZiIhEnnnz4N134YknVJwVlLo4iohIjo4dg5Ur1b1RRETylpLirldu1AiyGUhVgqQWNBERydE338Dvv0OLFl5HIiIifvfEE/DTT/Dxx1CypNfRRC61oImISI40QIiIiATjxx9hzBi46Sa4+GKvo4lsMVugRevgKFIwOh5Esvfll1CpEtSv73Ukkh/6TZP0dDxIOPTvDyVKuEGlpHBiskCLi4sjJSXF6zDER1JSUoiLi/M6DBHfWb7cdW8sFpPZIjIpx0lmynFS1D78EN55Bx55BGrV8jqayBeTKbdatWps376d5ORknVWKcdZakpOT2b59O9WqVfM6HBFf+e03+PpraNXK60gkP5Tj5DjlOAmHI0fcwCB/+Qv06+d1NNEhJgcJqVChAgA7duwgNTXV42jEa3FxcVSvXv3EcSEiTmIiWKsCLdIox0l6ynFS1CZNgm+/dcPrlyrldTTRISYLNHAJTD9WIiI5W7bMdW1s3tzrSCS/lONEJBx27YJhw6BjR/eQ0IjJLo4iIpK3ZcvgnHOgfHmvIxERET969FH4808YP97rSKKLCjQREcni2DE3gqO6N4qISHbWroXnn4d774UGDbyOJrqoQBMRkSw2bIA//lCBJiIiWVkLfftC5crw2GNeRxN9YvYaNBERydmyZe5ZBZqIiGT29tuwdClMmeLulSmhpRY0ERHJ4osvoEYNqF3b60hERMRPDh+GAQPgr3+FO+/0OpropBY0ERHJYtky13pmjNeRiIiInzz9NPz0k7s5dQlVEkVCLWgiIpLBzp0u+ap7o4iIpLd3L4wcCVdeCZdd5nU00UsFmoiIZLB8uXtWgSYiIukNGwaHDsHYsV5HEt1UoImISAbLlkGpUtC0qdeRiIiIX3z7Lfz73+66s8aNvY4muqlAExGRDJYtg2bNXJEmIiIC8NBDUKaMa0WTouW7As0Y84IxZrcxZkO6aZWNMYuMMd8Hnk/yMkYRkWj155+wahW0bOl1JCIi4hdLl8K778LDD0O1al5HE/18V6ABM4EOmaY9DCy21tYHFgf+FhGREPvySzhyBC680OtIRETED9LSoH9/OO00uP9+r6OJDb4r0Ky1nwL7Mk2+Bngx8O8XgWvDGpSISIxYutQNrd+mjdeRiIiIH8yZ43pWjBjhujhK0fNdgZaD6tbanQCB52wbV40xPYwxScaYpD179oQ1QBGRaLB0KZx7LlSq5HUkIiLitcOH4ZFH4Jxz4JZbvI4mdkRKgRYUa+00a20za22zqlWreh2OiEhEOXzYDbHftq3XkYiIiB88+6y7L+aYMVC8uNfRxI5Iuf/3LmNMDWvtTmNMDWC31wGJSPQ4cgTWrIEvvnAjGP7+OzRvDi1awAUXQJUqXkcYHitWuEFCLrrI60hERMRrBw64bo2XXgrt23sdTWyJlBa094AugX93Ad71MBYRiRLWwiuvQM2arhjr3x9Wr4Zdu2D0aLjqKqheHfr2hT/+8Draordkia4/85IxppIx5j/GmG+MMZuMMRpLU0Q8M2YM/PqrezbG62hii+8KNGPMq8ByoKExZpsxphvwBHCZMeZ74LLA3yIiBbZrF1x/vetTX78+vPEG7NgBP/4Ia9e6M4dLlrgbcj79NDRpAosWeR110Vq6FM4+GypX9jqSmPU0sNBaeyZwDrDJ43hEJEZt3w4TJrgced55XkcTe3zXxdFae1MOsy4JayAiErWWLIFOnVyr2NixbtjgzH3ry5Z1Xf0uusglqG7dXBePAQPgySej72zikSOue+edd3odSWwyxlQALgS6AlhrjwBHvIxJRGLX0KFw7Bg8/rjXkcQm3xVoIiJFKSkJ/v53dz+Xzz6DRo3yfk2bNrBuHfTrB+PGQbFi8MQT0VWkrVwJKSkaIMRDdYE9wAxjzDnAKqCPtfbQ8QWMMT2AHgCnn366J0GKSPT75ht44QW45x444wyvo4lNvuviKCJSVL79Fq64wg36sWhRcMXZcWXKwJQpcNddrgUt2s4qLl3qnnX9mWdKAOcBz1prmwKHgIfTL6CRikUkHB59FOLj3fD64g21oIlITNi6FS67zLV+ffgh1KqV/3UYA88841qahgxxCWzAgNDH6oWlS+Gss2JnxEof2gZss9YmBv7+D5kKNBGRorZiBbz5puviWC3buw5LOKgFTUSi3uHDrlvjgQOwcKEbFKSgihWD556Dzp3hgQdg3rzQxemV1FR3iwF1b/SOtfYXYKsxpmFg0iXARg9DEpEYYy08/DBUreq69It31IImIlFvyBB3DdncudC0aeHXV6IEvPii6zJ5++1u3TVrFn69Xlm1Cg4d0v3PfOBeYLYxpiTwI3C7x/GISAxZtAg++cSNXFy+vNfRxDa1oIlIVPviC3fN2J13uvuahUrp0vDaa5CcDLfd5ka7ilSLF7vnCy/0No5YZ61dG7jG7Gxr7bXW2v1exyQisSEtDQYOhDp1oGdPr6MRFWgiErX++AO6dHEJ56mnQr/+M8+Ef/0LPv7YFYGRauFCSEjQ9QYiIrHqrbdg9WoYNgxKlfI6GlGBJiJR68EH3Y2nZ84suu4ad9wBN94Ijz0GiYl5L+83v/0Gy5dDhw5eRyIiIl44etSN3Ni4sbvvp3hPBZqIRKWlS+HZZ91NqIuy654xMHWquwate3c34EYkWbzYdc9UgSYiEpteftldUz1iBBQv7nU0AirQRCQKpaW5wuy001zCKWoVK8KkSbBhg7u4OpIsXOjib9HC60hERCTcDh92Q+qffz5ce63X0chxKtBEJOq89BKsWQNPPOFuMB0O11zjhvIfMgR+/jk82ywsa2HBAnd/uBIa01dEJOZMnepy1qhRrkeI+IMKNBGJKocOwaBBcMEFcNNN4d32pEnuuU+f8G63oL7+GrZvV/dGEZFY9McfMHIkXHwxXHKJ19FIeirQRCSqPPkk7NwJ48eH/2xg7doweDC8846755rfLVzoni+/3Ns4REQk/CZNgt27XZGm1jN/UYEmIlFj2zYYO9aNqtiqlTcx9OsHf/0r3Hsv/PmnNzEEa+FCOOssOPVUryMREZFw+u03d0KzY0do2dLraCQzFWgiEjUGD3YDhDzxhHcxxMW5e6Nt2eLvAUP++AM++0zdG0VEYtH48a5Ie/xxryOR7KhAE5Go8MMPbnCQXr3cjam91K6dGzBk5EjXfcSPPvkEjhyBK67wOhIREQmnPXtgwgTo1AmaNvU6GsmOCjQRiQqjRrnWqwcf9DoSZ+xYSElxwxf70cKFULYstG7tdSQiIhJOTz4JyckwfLjXkUhOVKCJSMT76SfXetajB9So4XU0TsOGcNddbgjjr7/2OpqMrIV581xLX6lSXkcjIiLhsmMHTJ4Mt9wCjRp5HY3kRAWaiES8UaOgWDH/tJ4dN2QIlC8PDzzgdSQZJSa6+97ccIPXkYiISDiNGgVHj7r8JP6lAk1EItrmzTBzJtx5J9Sq5XU0GVWpAo895m4G/eGHXkfzP3PmQMmScPXVXkciIiLh8vPPMH063H471KvndTSSGxVoIhLRRo92rWcPP+x1JNm75x43aMmDD7oRJr2WlgZvvOFGb6xY0etoREQkXEaMcM+PPuptHJI3FWgiErG2bYMZM6BbN//ey6tUKdelZN06mD3b62hg2TLYvt3dK05ERGLDjz+6fHnnnXD66V5HI3lRgSYiEevpp12LkN+uPcvsxhuhWTN45BE3sqOXXn8dSpd2twEQEZHYMHw4lCgBgwZ5HYkEQwWaiESkAwfcCIk33OD9fc/yUqyYG9Z461aYNMm7OI4dc90br7zSDV4iIiLR77vv4OWX3X1Ca9b0OhoJhgo0EYlI06bB77/7b4TEnFx8MXTs6Lo7/vqrNzF8/jn88gt07uzN9kVEJPyGDXM9J/x6rbZkpQJNRCLOkSMwcSJccgmcd57X0QRvzBhXVD7+uDfbf/11KFMGrrrKm+2LiEh4bdwIr77qBqyqVs3raCRYKtBEJOK88oq72WaktJ4d99e/ugFNnnnGdTkJp6NH4T//ccVZ2bLh3baIiHhj2DD3mx9p+TLWqUATkYhiLYwbB2efDe3bex1N/j3+uGvF6t8/vNv98EPYvVvdG0VEYsX69e664/vuc/fllMihAk1EIsqCBfD11+5soDFeR5N/1au7e9DMmxfem1f/619Qo4ZuTi0iEiuGDYNy5cJ/QlAKTwWaiESU8eOhVq3Ivo9Xnz5Qty706+e6Hha1b76BDz5wI3iVLFn02xMREW+tXQtvvgn33w+VK3sdjeSXCjQRiRhffQWLF8O990JcnNfRFFypUjB2rGsJnD696Lc3aZIrzHr2LPptiYiI94YNg4oVXYEmkUcFmohEjIkTIT4eevTwOpLC+8c/oG1beOyxoh12/7ff4MUX4eabNYKXiEgsWL0a3nnH9dKoVMnraKQgVKCJSETYtQtmz4auXeGkk7yOpvCMgaefdjfcLsrrA154AQ4dcheJi4hI9Bs61BVmffp4HYkUlAo0EYkIzz7r7n8WTYXG2WfDgw+6Fq5Fi0K//mPHYPJkaNMGmjYN/fpFRMRfVq2CuXPdib+KFb2ORgpKBZqI+N6ff8KUKdCxIzRs6HU0ofXYY9Cggeu2eehQaNc9bx789FN0FbUiIpKzoUNdLxP97kc2FWgi4nuvvAJ79kTnxc6lS7uBQjZvhsGDQ7feo0dhyBCoXRuuvTZ06xUREX9audKdmOvfHypU8DoaKQwVaCLia9a6wUGaNIF27byOpmhceKEbYXHiREhMDM06J0+GdevcbQlKlAjNOkVExL+GDnVD6t97r9eRSGGpQBMRX1u8GNavd61nkXhj6mCNGePu79a5s2stLIzt213XySuucKNFiohIdFuxAubPhwED1HoWDVSgiYivTZzohoe/6SavIylaFSvC22/D7t2uSEtNLfi6+vd3XRwnT47uolZERJyhQ+Hkk+Gee7yOREJBBZqI+NZ338H770OvXu5arWiXkADTpsGSJe4saEEsWgRz5sCgQVC3bkjDExERH0pMhAULXN4oX97raCQUdGWCiPjW009DyZKuQIsVt90Ga9bAhAlw7rlw++3Bv3b3brev6td3w/eLiEj0GzbMtZ717u11JBIqKtBExJf274eZM+Hmm6F6da+jCa8nn4SvvoLu3WHfPujXL++uitu2waWXws6d8OGHUKpUeGIVERHvHG89Gz1arWfRRF0cRcSXpk+H5GTo29frSMKvRAl49103wMeAAdCli7sXXE5+/NHdjHrHDvjgA2jdOnyxioiId9R6Fp1UoImI76SmwqRJcPHFcM45XkfjjbJl4fXXXfJ9+WW46CJXtKW/mfWvv8KsWW6Y/oMH4eOP4W9/8y5mEREJH117Fr3UxVFEfOfNN12XvWee8ToSbxUr5m5e3aQJdOvmbjhdqpQrXJOT4fPPIS3N3Yx6yRK3nIiIxAa1nkUvtaCJiK9YC2PHQoMGcNVVXkfjD//4B+za5e4Jd/fd8NNPcOCAG6kxMdF1cVRxJiISO9R6Ft3UgiYivrJ0KaxeDVOnuhYkceLioF079xg/3utopKgYY4oDScB2a61OUYhIttR6Ft303x8R8ZWnnoKqVd1w8yIxqA+wyesgRMS/VqxQ61m0U4EmIr6xaRPMm+fOCJYp43U0IuFljDkV6Ag853UsIuJfaj2LfirQRMQ3xo+H0qXddVYiMWgi8CCQlt1MY0wPY0ySMSZpz5494Y1MRHxhxQqYP1+tZ9FOBZqI+MIvv8BLL0HXrq6Lo0gsMcZcBey21q7KaRlr7TRrbTNrbbOq+pKIxKThw6FyZbWeRTsVaCLiC8884+5/dv/9Xkci4onWwNXGmM3Aa0A7Y8wsb0MSET9ZuRLef1+tZ7FABZqIeO7gQZg82d3nq0EDr6MRCT9r7UBr7anW2jrA/wEfW2tv9TgsEfGR461n99zjdSRS1FSgiYjnpkyB336DRx7xOhIRERH/SUpyg2j176/Ws1ig+6CJiKeSk93gIJdfDgkJXkcj4j1r7RJgicdhiIiPDBum1rNYohY0EfHU9OmwZw88+qjXkYiIiPhP+tazChW8jkbCQQWaiHjm8GEYOxYuvBD+9jevoxEREfGf4cPhpJPUehZL1MVRRDzz0kuwfTu88ILXkYiIiPjPqlUwdy6MGKHWs1iiFjQR8URqKjzxBDRrBpdd5nU0IiIi/jNsmGs9u/deryORcFILmoh4YuZM+PFHmDABjPE6GhEREX9Zvdq1nj3+uFrPYk1EtaAZYzoYY741xvzXGPOw1/GISMGkpLizgi1bwt//7nU0IiIi/qPWs9gVMQWaMaY48AxwBdAYuMkY09jbqEQkX2bPhjp1KB1fjM+31+G5drPVeiYiIpLJ6tXw3nvQrx9UrOh1NBJukdTFsTnwX2vtjwDGmNeAa4CNnkYlIsGZPRt69IDkZAxQhy0woQc0Am65xevoREREfGPYMKhUSa1nsSqSCrRawNZ0f28DLshp4W+//Za2bdtmmNa5c2fuvvtukpOTufLKK7O8pmvXrnTt2pW9e/fSqVOnLPN79erFjTfeyNatW7ntttuyzO/fvz9///vf+fbbb+nZs2eW+Y8++iiXXnopa9eupW/fvlnmjxo1ilatWrFs2TIGDRqUZf7EiRM599xz+eijjxgxYkSW+VOnTqVhw4bMnTuXp556Ksv8l19+mdNOO405c+bw7LPPZpn/n//8hypVqjBz5kxmzpyZZf78+fOJj49nypQpvP7661nmL1myBIBx48Yxb968DPPKlCnDggULAHj88cdZvHhxhvknn3wyb775JgADBw5k+fLlGeafeuqpzJo1C4C+ffuydu3aDPMbNGjAtGnTAOjRowffffddhvnnnnsuEydOBODWW29l27ZtGea3bNmS0aNHA3D99dfz66+/Zph/ySWX8NhjjwFwxRVXkJKSkmH+VVddxYABAwCyHHegY+/ll1/mtEcecXelTi85mV+6daPE5Zfr2EPHXnbHnohIrDneejZ8uFrPYlXEdHEEsusIZTMsYEwPY0ySMSYpNTU1TGGJSFB+/jnbydUOHw5zICIiIv51vPXsvvu8jkS8Yqy1eS/lA8aYlsBQa+3lgb8HAlhrR2e3fLNmzWxSUlIYIxSRXNWpA1u2ZJ1euzZs3hzuaCRCGGNWWWubeR2H3yjHiUSn1ashIcG1ngU6T0gUyynHRVIL2kqgvjHmDGNMSeD/gPc8jklEgrT5zpEcIj7jxPh4GDnSm4BERER8Rq1nAhFUoFlrjwL3AB8Am4DXrbVfexuViAQjLQ1umX8L/ctP49iptd2Nz2rXhmnTNECIiIgIGrlR/ieSBgnBWjsfmO91HCKSPy+/DMuWQfcXbqH47SrIREREMlPrmRwXMS1oIhKZDhyABx+ECy6ALl28jkZERMR/1Hom6UVUC5qIRJ7Bg2HPHpg/H4rplJCIiEgWaj2T9PTfJREpMp9+CpMmwV13uVGpREREJKPjrWf9+6v1TBwVaCJSJA4cgH/+E+rWhSef9DoaERERfxo61LWe3Xuv15GIX6iLo4gUiT59YOtW+PxzKFfO62hERET8JykJ5s6Fxx9X65n8j1rQRCTk3nwTXnwRHnkEWrb0OhoRERF/GjoUKlfWtWeSkQo0EQmpbdugZ09o1gwee8zraERERPxpxQp4/3137VmFCl5HI36iAk1EQubQIbj6ajh82N37LC7O64hERET8aehQOPlkXXsmWekaNBEJibQ0uPVWWLfOjUZ15pleRyQiIuJPX34JCxbA6NFQvrzX0YjfqEATkZAYNAjeeQcmToSOHb2ORkRExL+GDoUqVeCee7yORPxIXRxFpNCmT4cxY9z9znShs4iISM6WLYMPPoAHHtAox5I9FWgiUiiTJ0OPHtChA/zrX2CM1xGJiIj415AhULUq9O7tdSTiVyrQRKTARo92Fzdfcw28/bYGBREREcnNZ5/BRx/Bww9D2bJeRyN+pQJNRPItLQ0eeshdd3bLLfDGG1C6tNdRiYiI+NuQIXDKKe6SAJGcaJAQEcmXnTvhtttg8WKXYJ55BorpVI+IiEiuPvnEPZ5+GuLjvY5G/EwFmogEbf586NLF3e/suefgjjt0zZmIiEherHWtZzVruuu2RXKj894ikqctW+Cf/3TD59esCatWQbduKs5ERESCsXixu/5s0CBdEiB5U4EmIjnav98NA9ywIbz+uruoOTERGjXyOjIREZHIYC089hicdhp07+51NBIJ1MVRRLJISoKpU+HVVyE52XVrHD7cJRcREREJ3oIF8OWXLq+WKuV1NBIJ8izQjDFnWmu/CUcwIuKd7793Q+XPmQOrV7sLmG+6Cfr0gSZNvI5OJPyU/0SksKyFwYPhjDPg9tu9jkYiRTAtaGuMMdOAodba/UUdkIiER2oqLF8OCxfCe+/B11+76eed524+feutULGitzGKeEz5T0QK5d133XXbM2boXqESvGAKtObABOB7Y8wwYIq19ljRhiUiRWHnTjcS4/vvuxtl/v47lCgBrVvDxIlw7bVQu7bXUYr4hvKfiBRYWpq79qxBA3fSUyRYeRZo1tr1wKXGmGuBsUAvY0x/a+2CIo9ORApt2zZ45RV3M+mkJDft1FNd98UOHaBdbGgy1AAAIABJREFUO7WUiWRH+U9ECuONN2DDBpg9250MFQlW0IeLtfYdY8x84H7gNWPMF0A/9c8X8Z/UVJcYnn/e3RTTWmjeHEaOhKuucteUaYh8keAo/4lIfh07BkOHQuPGcOONXkcjkSa/9Xw8sAp4EegNfGWM+TfwmLX2QKiDE5H8SU52Rdm4cfDzz1C3ruteceutUL++19GJRDTlPxEJ2qxZ8M038J//QPHiXkcjkSaYURz7AucHHvWAI8Ba4OnA8y3ARmPMddbaxCKMVURykJYGzz0HjzwCe/dCq1bwzDNw5ZVQTHc7FCkQ5T8RKYgjR2DYMGjaFK67zutoJBIF04LWH1gOPAt8Cayy1h5JN/8lY8xDwAvAX0MfoojkZtUquPtuWLECLrzQdWP829+8jkokKoQt/xljTgNeAk4B0oBp1tqnC7NOEfHGCy/ATz+5Abl0OYEURDCDhARza9oZwKjChyMiwUpNhUGDYPx4qFrVdae4+WYlA5FQCXP+Owr0t9auNsaUB1YZYxZZazeGYN0iEiYpKfD4464nyxVXeB2NRKpQjSmzB2gXonWJSB62b4fOnWHZMujRA8aMgUqVvI5KJCaFJP9Za3cCOwP//t0YswmoBahAE4kg//437NjhTprqhKkUVEgKNGutBZaGYl0ikruPPnItZcnJ8Oqr8H//53VEIrGrKPKfMaYO0BRIzDS9B9AD4PTTTw/lJkUkBP74A0aPhksugYsv9joaiWQaPkAkgsyYAZdf7ro0rlyp4kwk2hhjygFvAn2ttQfTz7PWTrPWNrPWNqtatao3AYpIjp5+GvbsgREjvI5EIp0KNJEIMWEC3HGHOzOXmAiNGnkdkYiEkjEmDleczbbWvuV1PCISvH37YOxYuPpqaNHC62gk0qlAE/E5a929zPr1g06dYO5cKFfO66hEJJSMMQZ4HthkrR3vdTwikj9jx8LBg2o9k9BQgSbicwMHuh/87t3htdegVCmvIxKRItAauA1oZ4xZG3hc6XVQIpK3nTtd98abb4YmTbyORqJBqEZxFJEiMG6cG6GxVy9342mNCCUSnay1nwP6hotEoJEj3a1vhg3zOhKJFmpBE/GpmTPhgQfgxhth8mQVZyIiIn7z008wbZrr5VKvntfRSLRQgSbiQ++9537s27eHl16CYvqmioiI+M7QoVC8ODz6qNeRSDTRf/tEfGbtWjd8fkICvPkmlCzpdUQiIiKS2YYN8PLLcO+9UKuW19FINFGBJuIje/bANdfAySe7VjSN1igiIuJPgwZBhQrw8MNeRyLRRoOEiPjEkSNuGP3du+Hzz/+/vfsOk7I6+zj+PaISWxTsCioWLEmUV1HBXlBEEUQRNTZs2BMLMXaNJTGWiCb22FBQETSiUoQoFsRCFNGoKIoFjBUNCooI5/3jLGGBBRaYmTO78/1c11y7O8/szI9nh5m5n3Oe+8Dqq+dOJEmSajJ8eFr25o9/hMaNc6dRfWOBJpWJ006DZ56BXr3S9EZJklR+YkyjZmusAb/5Te40qo8s0KQycNddcNNNcNZZaR0VSZJUngYMSDNdbroJllsudxrVR56DJmX273/DSSfBrrumqRKSJKk8TZ+ezj3bYAM45pjcaVRfOYImZTR5MnTpAiusAL17p1a9kiSpPPXuDaNHw333wVJL5U6j+soCTcrolFPgrbfgiSfSXHZJklSefvghrXe21Vbp4KpULBZoUiY9e6Zzzy64ANq0yZ1GkiTNzw03wEcfwZ13whKeJKQi8uklZTB2LJx8Muy8M1x0Ue40kiRpfr7+Gi6/HNq2hd12y51G9Z0FmlRi06bBoYfCkkvCPfd43pkkSeXuiivgm2/gz3/OnUSVwCmOUoldcgm89BL06QNNm+ZOI0mS5ufjj+G66+Cww2CLLXKnUSVwBE0qoWefTa30jzoKDjwwdxpJkrQgF16YFqe+9NLcSVQpLNCkEvnmm3T0rVmzdCROkiSVt9deg7vvht/+FtZdN3caVQqnOEolcsopMGECDB+e1j2TJEnlK0bo3h0aNUqLU0ulYoEmlcB990GvXvCHP8C22+ZOI0mSFmTwYBg6FK69FlZaKXcaVRKnOEpF9tFHcOKJ0Lq1R+AkSaoLpk+H3/0ONtgATjopdxpVGkfQpCKaPh2OPDJ9veee1FpfkiSVt7vvhjfeSB2Xl146dxpVGj8uSkV0zTUwbBjccUc6CidJksrbd9/B+edDq1bQuXPuNKpEFmhSkbz8Mpx3HhxwAHTtmjuNJEmqjauugv/8B/r1gxByp1El8hw0qQi+/RYOOQTWXBNuu80XeEmS6oLx41OBdtBB6dxxKQdH0KQiOPVUGDcuTW9s1Ch3GkmSVBvnngszZsAVV+ROokrmCJpUYPfdl04uPv982HHH3GkkSVJtjByZGnqdfjqst17uNKpkFmhSAY0ZA8cfD9tvDxdckDuNJEmqjRhTYbbaanDOObnTqNKVzRTHEMKBwMXApsA2McaR1badAxwDTAd+E2McnCWkCuaNN+Af/0jnak2dmi4bbAB77QW/+EXdPGdr8uTUEKRhwzSKZkt9SZLqhn794Lnn4Oab4ec/z51Gla6cPkK+AewP3FL9yhDCZsDBwC+AtYChIYTmMcbppY+oxfHDD2k9kVtugeefT9c1bJguSy0FX32VFoVce23Yd1846yxo1ixv5tqKEbp1gzffhCeegKZNcyeSJEm18f330L07bL45HHts7jRSGU1xjDG+FWMcU8OmjsD9McapMcZxwFhgm9Km0+IaOxa23DIt2vzll2l9sC+/TEXbf/+bvv/oI/j731PXpLvugubN03TBjz/OnX7BbroJeveGSy+FNm1yp5EkSbV1zTXw4YfQowc0aJA7jVRGBdp8rA1U/4g+vuq6uYQQuoUQRoYQRn7xxRclCacFe+op2HZb+OwzePRRePttOOMMWHnl2W/XtCkccww8+CC8914qzu68EzbcEC66CH78MU/+BRk+HE47DfbZx3nrkiTVJRMmwJ/+lE5R2HXX3GmkpKQFWghhaAjhjRouHef3azVcF2u6YYzx1hhjyxhjy1VXXbUwobVYbrkF9twT1lgDXnoJ2rev3flla60Ff/sbvPsudO4Ml1wCW28Nr75a/MwL4913oUOH1O2pZ09Yoi4c8pAkSQCcfTZMn57WPpPKRUk/TsYY28QYf1nD5ZH5/Np4oPoZPU2AT4qbVIVw//1wwgnQti2MGJGagCysddeFXr3gkUfg889TkXbxxfDTTwWPu9C++ALatUtF2YAB0Lhx7kSSJKm2RoyAe++FM8+sO+e8qzLUheP9/YGDQwgNQwjNgI2AlzJn0gL8+99puuIOO8DDDy9+R6QOHdJ9HnII/OEPaRrCRx8VJuui+P576NgxTY3o3z9Nw5QkSXXDjBnwm9/Ammt6eoLKT9kUaCGETiGE8UBr4PEQwmCAGOO/gT7Am8Ag4GQ7OJa3SZNg//1TUdanT+rQWAiNG6cFJO+9F0aNghYtUvFXalOnwsEHwwsvpCytW5c+gyRJWnR33JEWpr7qKlh++dxppNmVTYEWY3w4xtgkxtgwxrh6jLFttW2Xxxg3iDFuHGMcmDOn5i9GOOqo1OSjT590ZKrQDj00nYu2/vqpEDzpJJgypfCPU5MpU9LIWf/+6Ry5Aw4ozeNKkqTC+PrrNGq2ww7w61/nTiPNrWwKNNUPN94IDz2UjkjtuGPxHmfDDdNaameemVrcb701jB5dvMcD+O671KnxiSfScgAnnVTcx5MkSYV34YUwcSL89a+1a1wmlZoFmgrmq6/gggtgjz1S2/liW3ppuPpqGDw4PfY226S1TIrRQOSTT1I3ymefTdMajzmm8I8hSZKKa/TodDD5xBPTqRJSObJAU8FcfHFadPraa0t7RGrPPdML7p57Qvfuac21V14p3P0PHAhbbAGvvZbWaHM6hCRJdU+McOqp0KhRWr5HKlcWaCqIN99MUw1POAF+8YvSP/5qq6VW/A88kDorbrNNWgx7cdYrnzo1FXx7753WZRs5Ejp1KlxmSZJUOr17wzPPwB//6NI4Km8WaCqIM8+EFVZILfBzCQG6dIG33kpTEHv0SOuonX56Ktpq6/vv07z0DTdMUyZPOil1bNx00+Jll6QQwl4hhDEhhLEhhLNz55HqhV69YL31iEsswU5HrMf56/fi2GNzh5LmzwJNi23AABg0CC66CFZZJXeaNHXhllvSqF6XLqnYWn/91ODjr3+Fd99N0xyq+/rrdC7b+efDeuultVGaNYMhQ+CGG2CZZbL8UyRViBBCA+AGoB2wGXBICGGzvKmkOq5XL+jWDT78kBAjTWd8yMWfdGOJ+3rlTibNV4hzflKtJ1q2bBlHjhyZO0a9N2MG/OpXqTHH66+nxh3lZtw4uP56ePzxVJxBGu1bfnlYdtn0bxg3Ll0fQmpyct55sNNO+TJLSkII/4oxtsydo9hCCK2Bi2cuMRNCOAcgxvinmm6/wgorxK222mq267p06cJJJ53ElClT2Hvvvef6na5du9K1a1e+/PJLOnfuPNf2E088kYMOOoiPP/6Yww8/fK7tZ555Jvvuuy9jxozh+OOPn2v7+eefT5s2bRg1ahSn1dAp6o9//CPbbbcdzz//POeee+5c23v06EGLFi0YOnQol1122Vzbb7nlFjbeeGMeffRRrrnmmrm233PPPTRt2pQHHniAm266aa7tffv2ZZVVVuGuu+7irrvummv7gAEDWHbZZbnxxhvp06fPXNuHDRsGwNVXX81jjz0227ZlllmGgQPTKkCXXnop//znP2fbvvLKK9OvXz8AzjnnHEaMGDHb9iZNmnDvvfcCcNpppzFq1KjZtjdv3pxbb70VgG7duvHOO+/Mtr1Fixb06NEDgMMOO4zx48fPtr1169b86U/pqXTAAQfw1VdfzbZ9991354ILLgCgXbt2fP/997Ntb9++Pd27dwdgl112YU7l+ty7/4UXWGPq1Llu/2nDhhzcqtX/fva553Mv1+ve008/XeN73JJz3VJaCI8/nkaqevcuz+IM0kjYtdemy3vvpZGyMWPSmmZTpqTi8thjoVWr1K5/hRVyJ5ZUgdYGPq7283hg2+o3CCF0A7oBNGzYsHTJpDpqtRqKs/ldL5ULR9C0WHbeGT74IBU+S1ruSyqwChpBOxBoG2M8turnw4FtYoyn1nR73+OkWlhvPfjww7mvX3fd9OFFymxe73Geg6ZF9tJLqRvSaadZnEnSYhoPNK32cxPgk0xZpHrhm7MuZwrLzn7lssvC5ZfnCSTVkgWaFtk118CKK2I3JElafC8DG4UQmoUQlgYOBvpnziTVaSc+eygnNLiVaWutm04yX3dduPVWOPTQ3NGk+bJA0yIZNw769oXjj/ecLUlaXDHGn4BTgMHAW0CfGOO/86aS6q5Bg+D++2GDCw5lqQkfpI5gH3xgcaY6wYlpWiQ9esASS6R29JKkxRdjHAAMyJ1DquumTElrmG68MZztioKqgyzQtNC+/hpuvx1+/WtYe+3caSRJkma59NI002fYMLDhqeoipzhqofXsCZMnw+mn504iSZI0y6hRcNVVcNRRqdO0VBdZoGmhxAh33AEtW0KLFrnTSJIkJT/9BMccA6usAldfnTuNtOic4qiF8uqrMHo03HBD7iSSJEmzXHstvPIK9OkDjRvnTiMtOkfQtFDuvDPN5z7kkNxJJEmSkrFj4cILoWNH6Nw5dxpp8VigqdamToXevWG//aBRo9xpJEmS0ukX3brB0kunGT4h5E4kLR6nOKrW+veHiRPTibeSJEnl4NZb4amn4JZb7C6t+sERNNXanXdCkybQpk3uJJIkSfDhh9C9O+y2Gxx3XO40UmFYoKlWJkyAwYPhyCOhQYPcaSRJUqWLEY49Nn1/++1ObVT94RRH1UrPnjBjBnTtmjuJJElSmto4dCjcdBOst17uNFLhOIKmWundG7bfHjbcMHcSSZJU6apPbTz++NxppMKyQNMCvfUWvPEGHHRQ7iSSJKnSzZgBRx+dvndqo+ojpzhqgR58ML34HXBA7iSSJKnS/fWv8OSTaYqjUxtVHzmCpgV68EHYYQdYa63cSSRJUiV76y04+2xo335WgxCpvrFA03zNnN544IG5k0iSpEo2bRocfjgstxzcdptTG1V/OcVR8+X0RkmSVA4uuwz+9S/o2xfWWCN3Gql4HEHTfD34YOre6PRGSZKUy4gRcPnlaQTNg8aq7yzQNE9vv52mN3bpkjuJJEmqVJMmwaGHQtOmqUGIVN85xVHz5PRGSZKU28knw0cfwTPPwIor5k4jFZ8jaJqnPn2c3ihJkvLp1QvuvRcuvBC22y53Gqk0LNBUo3fftXujJEnKZ9w4OPHEdLD43HNzp5FKxwJNNerfP33t2DFvDkmSVHl+/BEOOgiWWCKNoC3pSTmqID7dVaP+/WHzzWHddXMnkSRJleacc+Dll1NL/fXWy51GKi1H0DSXr76C4cOhQ4fcSSRJUqV59FH4y19ScxAblakSWaBpLgMHwvTpsO++uZNIkqRK8tFHcOSR8H//B1dfnTuNlIcFmubSvz+ssQa0bJk7iSRJqhQ//giHHAI//ZQ6Sf/sZ7kTSXl4Dppm8+OPMGgQHHxwOjFXkiSpFH73O3j+eXjgAdhww9xppHz8CK7ZPP00fPut0xslSVLp3H8/XH89nHYadOmSO42UlwWaZvPoo7DMMrD77rmTSJKkSvDmm3DssWm9syuvzJ1Gys8CTf8TYzr/bI89YNllc6eRJEn13aRJsP/+sPzy6byzpZbKnUjKzwJN//P66/Dhh05vlCRJxTdjBhx+OIwdm847W2ut3Imk8mCTEP3PY4+lr/vskzeHJEmq/y6+OM3cuf562Hnn3Gmk8uEImv5nwIDUWn/NNXMnkSRJ9Vm/fnDppXD00XDKKbnTSOXFAk0ATJwII0ZAu3a5k0iSpPps9Oi0GHWrVnDjjRBC7kRSebFAEwBDhqS54BZokiSpWD77DDp0gJ//PI2iNWyYO5FUfjwHTQAMHAiNG8M22+ROIkmS6qPvv4f99oPPP4dnn7UpiDQvFmhixoxUoLVtCw0a5E4jSZLqmxjT+WYvvJBGzrbaKnciqXw5xVG8+mo6muX0RkmSVAwXXwz33w9XXJHWPZM0bxZoYuDA9LVt27w5JElS/XPnnXDJJWkE7ayzcqeRyp8FmhgwALbeGlZbLXcSSZJUnwwcCMcdB3vuCTffbMdGqTYs0CrcxInw4otOb5SkXEIIV4UQ3g4hjA4hPBxCWCl3JqkQXn4ZOneGzTeHvn1hqaVyJ5LqBgu0CvfEE7bXl6TMhgC/jDFuDrwDnJM5j7TY3nsP9tknzc4ZMABWWCF3IqnusECrcAMGwMorpymOkqTSizE+EWP8qerHF4AmOfNIi2vCBNhjj3QAeNAgWGON3ImkusUCrYLNmAGDB6d54bbXl6SycDQwsKYNIYRuIYSRIYSRX3zxRYljSbXz5Zfpc8WXX6bibOONcyeS6h7XQatgo0bZXl+SSiGEMBSoaRzhvBjjI1W3OQ/4CehV033EGG8FbgVo2bJlLFJUaZFNmpQ+U7z/firOWrbMnUiqmyzQKtigQenrnnvmzSFJ9V2Msc38tocQjgTaA7vHGC2+VOdMmQIdOqSDvw8/DDvvnDuRVHdZoFWwwYPh//4PVl89dxJJqlwhhL2A3wM7xxin5M4jLawpU2DffeHZZ+Hee6F9+9yJpLqtbM5Bm1+b4RDCOSGEsSGEMSEEl1MugEmT4PnnXZxaksrA34AVgCEhhFEhhJtzB5Jq6/vvYb/94Kmn4O674ZBDcieS6r5yGkEbApwTY/wphPBnUpvh34cQNgMOBn4BrAUMDSE0jzFOz5i1znvySfjpJ9hrr9xJJKmyxRg3zJ1BWhQ//AD77w9Dh8Kdd8Jhh+VOJNUPZTOCNp82wx2B+2OMU2OM44CxwDY5MtYngwfD8stD69a5k0iSpLpm8uQ0rXHQILjtNjjyyNyJpPqjbAq0OVRvM7w28HG1beOrrpuLLYhrJ8b0grr77rD00rnTSJKkumRmt8Ynn4S77oJjjsmdSKpfSlqghRCGhhDeqOHSsdpt5mwzHGq4qxo7XMUYb40xtowxtlx11VUL/w+oJ959Fz74wOmNkiRp4UycmBahHjEC7rvPkTOpGEp6DtoithkeDzStdrMmwCfFSVgZZrbXt0GIJEmqrU8+SQd3x4yBfv1SW31JhVc2UxyrtRnuMEeb4f7AwSGEhiGEZsBGwEs5MtYXgwdD8+bQrFnuJJIkqS4YMwa22w7GjYPHH7c4k4qpbAo05tFmOMb4b6AP8CYwCDi5FB0cp0yBsWOL/Sil98MPqRWuo2eSJKk2XnoJtt8+fTYaNgzazHc+lKTFVTZt9ufXZjjGeDlweQnj0LYtTJsGL7xQykctvueeS2uWWKBJkqQFefTRtLbZaqulGTgbbZQ7kVT/ldMIWllp1w5efBEmTMidpLAGDYKGDWGXXXInkSRJ5SpGuO466NgRNtkEhg+3OJNKxQJtHjp1Sl//8Y+8OQpt4EDYaSdYbrncSSRJUjn66Sc49VQ47TTYbz94+mlYc83cqaTKYYE2D5tuChtvDA8/nDtJ4Xz0Ebz5pu31JUlSzSZOhL33hhtugO7doW9fD+pKpWaBNh+dOqWTYSdOzJ2kMAYPTl8t0CRJ0pxefx223jqNmP3973DVVbCEnxSlkvO/3Xzsvz9Mnw6PPZY7SWEMGgRNm6bRQUmSpJkefBBat06NxJ5+Go45JnciqXJZoM1Hy5bQpEn9mOY4bRoMHZqan4SQO40kSSoHU6fCb38LXbrA5pvDv/4FrVrlTiVVNgu0+QghnRw7eDBMnpw7zeIZMQImTXJ6oyRJSsaNgx13hOuvTw1Bhg2zGYhUDizQFqBTpzTcP/P8rbpq0CBYcknYbbfcSSRJUm59+8KWW8I778BDD8G118LSS+dOJQks0BZop52gceO6P81x0CDYbjtYccXcSSRJUi6TJkHXrnDggdC8ObzyyqylhSSVBwu0BVhySejQITUKmTYtd5pF8+mn8OqrTm+UJKmSDR8OLVrAPffABRfAc8/B+uvnTiVpThZotdCpE3zzTZqbXRfZXl+SpMo1ZQqcfno63yxGeOYZuOQSWGqp3Mkk1cQCrRb22AOWXRb+8Y/cSRbNoEGw+uqwxRa5k0iSpFJ6+unUnbFHDzjxRBg9GrbfPncqSfNjgVYLyyyTRp/+8Q+YMSN3moXz00+pQGvXzsUmJUmqFF9+mdYy22WXNGr21FNwww2wwgq5k0laED+y19J++8Enn8DIkbmTLJznn0/TM9u3z51EkiQV24wZcOedsMkm0LMnnHVWGjXbZZfcySTVlgVaLbVvDw0a1L1ujo8/nuaY77FH7iSSJKmYXnwxTV88+mjYdNPUIOzPf4bllsudTNLCsECrpUaN0tGnunYe2mOPpaUCfv7z3EkkSVIxTJgAhx8OrVrBBx+kEbSnn4Zf/jJ3MkmLwgJtIXTqBG+/nS51wbhx8OabTm+UJKk++vprOPts2GgjePBBOPfctPB0166edy7VZf73XQgdO6avdWUU7fHH09d99smbQ5IkFc7kyWnq4vrrw5VXwgEHwFtvweWX2wREqg8s0BZCkyaw9dZ15zy0xx6D5s3TkTVJklS3ffddKsiaNUsjZ9tvD6NGpYWnmzXLnU5SoVigLaT99oOXXkrzvcvZd9+llrqOnkmSVLdNnAiXXQbrrQe//z1suSUMH54OxG6+ee50kgrNAm0hdeqUvj7ySN4cC/LPf8KPP3r+mSRJddUHH8BvfwvrrAMXXADbbgsjRqT1TbfbLnc6ScVigbaQNtkENt4YHnood5L5e/zxNA99hx1yJ5EkSbUVIwwbls4r22ADuPHG9P3o0em9vVWr3AklFZsF2kIKATp3Ti+eX3yRO03NYkwv4m3bwtJL504jSZIW5L//TcXYFlvArrumzxndu8P778Pdd8OvfpU7oaRSsUBbBJ07w/Tp5dvN8V//gk8+cXqjJEnlLEZ47rnUFn/NNeHkk2HJJeH222H8+NSpsWnT3CklldqSuQPURVtsARtuCH37wnHH5U4zt4ceggYNYN99cyeRJElzev/91HmxZ8/0/fLLp4Wmu3WDrbbKnU5SbhZoi2DmNMerroKvvoKVV86daJYYoV+/ND2icePcaSRJEqTuz336wAMPwIsvps8Su+0GF10E+++fijRJAqc4LrKZ0xzLrZvjm2/CO++kF3tJkpTPuHHwl7+khl1Nm8IZZ6QOy1dcAR9+CEOHwhFHWJxJmp0jaItoyy3TeiR9+8LRR+dOM0u/fumo3H775U4iSVJlmTEDRo5M65P17w+vvZaub9ECLrkEunSB5s3zZpRU/izQFlEIcOCB0KMHfP01NGqUO1Hy0ENpbZQ118ydRJKk+u/rr2HIEBg4MK1P9umnsMQS6b34mmvS+qnNmuVOKakusUBbDDPPQ3v00TRFIbf33ktH6665JncSSZLqp2nT0jlkQ4aky4svppGzlVaCPfdMDbratSuv89Ml1S0WaIth663TnPIHHyyPAm3m4tmefyZJUmHMmJEWiX7yyXR55hn49ts0StayJZx7birIttkmtciXpMXlS8limNnN8YYbymOa40MPzTo3TpIkLbwZM+CNN9JC0cOGwdNPw8SJaVvz5nDoobDHHqlbcu73fUn1kwXaYjr0ULj22tQ69/jj8+WYMAFeeAEuuyxfBkmS6prp09PpAU8/nS7PPjurIGvWDDp2TMXYrrtCkyZ5s0qqDBZoi2nLLWGzzdKCkzkLtIcfTl8POCBfBknSogshdAeuAlaNMX6ZO099NW1a6rT4zDOpIBs+HCZNSts22CB1Qd5553RZd928WSVVJgu0xRQCHH44nHNOatKxwQZ5ctxzD/zqV7DJJnl3ocmiAAAR6UlEQVQeX5K06EIITYE9gI9yZ6lvpk6Fl16aNV1xxAiYMiVt23RT+PWvYaedYMcdHSGTVB4s0Arg0EPTScL33gsXXVT6x3/77fTmc/XVpX9sSVJBXAucBTySO0hd9+OP8PLLqaHHU0+lguyHH9IB1c03h2OPTaNjO+wAq62WO60kzc0CrQCaNk1z0++5By68ML0JlFLPntCgQSoUJUl1SwihAzAhxvhamM8bSAihG9ANYJ111ilRuvI3Ywa8/jr8858wdGiaujh5cnov3mILOOGE9B6944429ZBUN1igFcjhh8NRR6VGHa1bl+5xp09PhWHbtrDGGqV7XElS7YUQhgI1vUqfB5wL7Lmg+4gx3grcCtCyZctY0IB1zGefwRNPpMuQIelnSNP8u3aF3XZLo2SuRSapLrJAK5ADDoCTTkqjWaUs0IYNg/Hjnd4oSeUsxtimputDCL8CmgEzR8+aAK+EELaJMX5awohlbcaM1Njj8cdhwID0PcCqq6aW93vuCW3awNpr580pSYVggVYgK6wAnTrBAw9Ajx7QsGFpHvfuu2HFFVMbYElS3RJjfB3435lQIYQPgJZ2cYTvv0/TFvv3h0cfhU8/TYtDt2qVlpRp1w5atEjXSVJ9YoFWQEccAb17w2OPlabd/bffQr9+cNhh8LOfFf/xJEkqpsmTYeBA6Ns3vZdOnpwOgLZrBx06wF57OW1RUv1ngVZAbdqkNVP+9rfSFGgPPZRaBR9xRPEfS5JUfDHG9XJnKLWpU2HQILjvvjRSNmVKmrp46KGw//6wyy6lm5UiSeXAAq2AGjSAU06B3/0ORo9O7XyL6e67YcMNYbvtivs4kiQVUozw/PPpfezBB+Gbb2CVVdIBx4MOSh0XGzTInVKS8nDmdoEdcwwsuyxcd11xH2f06LS+y9FHl76tvyRJi2LCBLj8cmjePK1D1rs37Ltvmtb4ySdw001pxMziTFIls0ArsEaN0hHAXr3giy+K9zh/+Qsstxwcf3zxHkOSpMU1fXrqvtixI6yzDpx/PjRpAnfdlRp/9OyZzi1baqncSSWpPFigFcGpp6Y59bfeWpz7nzAhHXU85hho3Lg4jyFJ0uKYOBGuuipNxW/fHl58Ec46C8aOTTNAjjwSll8+d0pJKj8WaEWw2WZpXZYbb4Rp0wp//9dfn45InnZa4e9bkqTFMWZMmt2x9tqpIFt3XejTBz7+GP70J9hgg9wJJam8WaAVyW9/m+bT9+1b2PudNAluvhk6d4ZmzQp735IkLapnn03TGDfZJDX/OOwweO01GDYMDjzQKYySVFsWaEXSrh1stFGa3jFjRuHu9+9/T0Xa735XuPuUJGlRxAhPPAE77ZQuw4fDhRfCRx/BbbcVv5uxJNVHFmhFssQS6U3q1Vfh3nsLc5/TpkGPHrDzztCyZWHuU5KkhRVjavyx7bbQti2MG5em33/0EfzhD7DaarkTSlLdZYFWRL/+dSqkzjkHJk9e/Pu7/fY0h79798W/L0mSFlaMMHRoWn+zfXv48ss0Uvbee6lB1rLL5k4oSXWfBVoRLbEEXHttOhftqqsW774+/RTOPht22w322acw+SRJqq2RI2H33VMTrAkTUqfiMWPg2GNh6aVzp5Ok+sMCrch22CGdHH3llTB+/KLfzxlnwPffp86QLkwtSSqVcePgkENg663hjTfSVMZ334XjjrPxhyQVgwVaCfz5z6kt/rnnLtrvP/EE3Hdf+v2NNy5sNkmSavLtt/D736eujI88Auedl9YwO/VUaNgwdzpJqr8s0EqgWTM4/XS45550UvXC+P57OPFEaN48TXGUJKmYZsyAnj3T+86VV6bRs3ffhcsug5//PHc6Sar/LNBK5PzzYcstoUuXNI+/NmJMDUbefz+tfeYRS0lSMY0eDTvuCEceCeusAy+8AHfdlRadliSVhgVaiSy/fBo9W3XV1ORj3Lj53z5GOPNMuO46OOUU2HXX0uSUJFWe775L62tuuSW88w7ccQeMGJHa6EuSSssCrYTWWAMGDkzrmbVrB199VfPtpk+H449PHSBPPTUVaZIkFcOgQfCLX8DVV8PRR8Pbb8NRR6VOxJKk0vPlt8Q23TSdbP3BB7DRRmnq46efpm1Tp8Izz8DBB6d1Zc49NxVnvklKkgpt4kTo2jUdMFxuOXjuudQ6f+WVcyeTpMrmR/8MdtwRhg9P0xb/+EdYd11o3RpWWgl23hkeegiuuAIuv9yW+pKkwnv0UdhsM+jVKx0ofPVV2H773KkkSQBL5g4wUwjhUqAjMAP4HOgaY/wkhBCA64C9gSlV17+SL2lhbLUV9OuXWhb/5S/wyitwwgmpaNtxR2jUKHdCSVJ99Z//wFprpemNLVrkTiNJqq5sCjTgqhjjBQAhhN8AFwInAO2Ajaou2wI3VX2tFzbcMC0+LUlSqRx3XDrPzIWmJan8lM0UxxjjpGo/LgfEqu87Aj1j8gKwUghhzZIHlCSpngjB4kySylU5jaARQrgcOAL4LzCzsfzawMfVbja+6rr/1PD73YBuAOuss05Rs0qSJElSoZV0BC2EMDSE8EYNl44AMcbzYoxNgV7AKTN/rYa7ijVcR4zx1hhjyxhjy1VXXbU4/whJkiRJKpKSjqDFGNvU8qa9gceBi0gjZk2rbWsCfFLgaJIkSZKUXdmcgxZC2Kjajx2At6u+7w8cEZJWwH9jjHNNb5QkSZKkuq6czkG7IoSwManN/oekDo4AA0gt9seS2uwflSeeJEmSJBVX2RRoMcYD5nF9BE4ucRxJkiRJKrmymeIoSZIkSZXOAk2SJEmSyoQFmiRJkiSVCQs0SZIkSSoTFmiSJEmSVCYs0CRJkiSpTFigSZIkSVKZsECTJEmSpDJhgSZJkiRJZSLEGHNnKIoQwhfAh4t5N6sAXxYgTrGZs7DMWVh1JSfUnayVlHPdGOOqhQhTn1TYe1ypuD9mcV/M4r6YnftjlqK9x9XbAq0QQggjY4wtc+dYEHMWljkLq67khLqT1ZwqBP8+s3N/zOK+mMV9MTv3xyzF3BdOcZQkSZKkMmGBJkmSJEllwgJt/m7NHaCWzFlY5iysupIT6k5Wc6oQ/PvMzv0xi/tiFvfF7NwfsxRtX3gOmiRJkiSVCUfQJEmSJKlMWKBJkiRJUpmo+AIthHBgCOHfIYQZIYSWc2w7J4QwNoQwJoTQdh6/3yyE8GII4d0QwgMhhKVLkPmBEMKoqssHIYRR87jdByGE16tuN7LYuWp4/ItDCBOqZd17Hrfbq2ofjw0hnJ0h51UhhLdDCKNDCA+HEFaax+2y7M8F7Z8QQsOq58TYqufieqXKVi1D0xDCUyGEt6r+P/22htvsEkL4b7Xnw4WlzlmVY75/x5BcX7U/R4cQtsyQceNq+2lUCGFSCOG0OW6TbX+GEO4IIXweQnij2nWNQwhDql4Lh4QQGs3jd4+sus27IYQjS5VZs8v9ultOano+V6LavI5XkhDCz0IIL4UQXqvaH3/InSm3EEKDEMKrIYTHcmfJreifCWOMFX0BNgU2BoYBLatdvxnwGtAQaAa8BzSo4ff7AAdXfX8zcGKJ818DXDiPbR8Aq2TctxcD3RdwmwZV+3Z9YOmqfb5ZiXPuCSxZ9f2fgT+Xy/6szf4BTgJurvr+YOCBDH/rNYEtq75fAXinhpy7AI+VOtvC/h2BvYGBQABaAS9mztsA+JS0mGVZ7E9gJ2BL4I1q110JnF31/dk1/T8CGgPvV31tVPV9o9zPiUq7lMPrbjldano+V+KlNq/jlXSpeg9Yvur7pYAXgVa5c2XeJ2cAvcvhvTz3pdifCSt+BC3G+FaMcUwNmzoC98cYp8YYxwFjgW2q3yCEEIDdgL5VV90N7FfMvDU8fhfgvlI9ZhFsA4yNMb4fY/wRuJ+070smxvhEjPGnqh9fAJqU8vEXoDb7pyPpuQfpubh71XOjZGKM/4kxvlL1/bfAW8DapcxQQB2BnjF5AVgphLBmxjy7A+/FGD/MmGE2McZngIlzXF39eTiv18K2wJAY48QY49fAEGCvogXVvGR/3S0n83g+V5x69jq+2KreA76r+nGpqkvFdtYLITQB9gH+njtLJaj4Am0+1gY+rvbzeOZ+oVoZ+Kbah/uablNMOwKfxRjfncf2CDwRQvhXCKFbCXNVd0rVNLE75jHlqTb7uZSOJo2e1CTH/qzN/vnfbaqei/8lPTezqJpi+X+ko41zal01XWRgCOEXJQ02y4L+juX2nDyYeR+EKYf9OdPqMcb/QPqgB6xWw23Kbd9WKv8Omq8FvI5XjKopfaOAz0kHlyp5f/QAzgJm5A5SJor6mXDJQt9hOQohDAXWqGHTeTHGR+b1azVcN+eRk9rcZpHUMvMhzH/0bPsY4ychhNWAISGEt6uOFBbM/HICNwGXkvbJpaTpmEfPeRc1/G7Bj1DVZn+GEM4DfgJ6zeNuir4/a5D1ebiwQgjLA/2A02KMk+bY/Appmt53IZ2P+A9go1JnZMF/x3Lan0sDHYBzathcLvtzYZTNvq1w/h00Twt4Ha8oMcbpQIuqc9MfDiH8MsZYcecqhhDaA5/HGP8VQtgld54yUdTPhBVRoMUY2yzCr40Hmlb7uQnwyRy3+ZI0/WnJqpGLmm6zSBaUOYSwJLA/sNV87uOTqq+fhxAeJk1rKWhBUdt9G0K4DajppNLa7OfFVov9eSTQHtg9Vk0uruE+ir4/a1Cb/TPzNuOrnhcrkmG6TghhKdKbeq8Y40Nzbq/+Rh9jHBBCuDGEsEqM8ctS5qzF37Ekz8laage8EmP8bM4N5bI/q/kshLBmjPE/VVNCP6/hNuNJ587N1IR0/q9Kq5ye4yojC3odr1Qxxm9CCMNIU7IrrkADtgc6VB0M/Bnw8xDCvTHGwzLnyqbYnwmd4jhv/YGDQ+qQ14x0ZPql6jeo+iD/FNC56qojgXmNyBVaG+DtGOP4mjaGEJYLIaww83tSI4ySvqjMcd5Op3k8/svARiF1w1yaNJ2rfynyzRRC2Av4PdAhxjhlHrfJtT9rs3/6k557kJ6LT86ryCyWqnPebgfeijH+ZR63WWPmuXEhhG1Irz9flS5lrf+O/YEjQtIK+O/MqXsZzHOUvBz25xyqPw/n9Vo4GNgzhNCoasrznlXXqbSyv+6q/NTmdbyShBBWrRo5I4SwDFWfu/KmyiPGeE6MsUmMcT3S68WTlVycleIzYcUXaCGETiGE8UBr4PEQwmCAGOO/SR0a3wQGASdXDXUTQhgQQlir6i5+D5wRQhhLOu/n9hJFn+u8lBDCWiGEAVU/rg48F0J4jVRYPh5jHFSibDNdGVIL0tHArsDpc+asGnk8hfQh7S2gT9W+L6W/kTpWDQmpXerNc+Yk0/6c1/4JIVwSQuhQdbPbgZWrnoNnkDroldr2wOHAbqHasgohhBNCCCdU3aYz8EbVPrye1P201NOqavw7zpFzAKm74FjgNlKXzJILISwL7AE8VO26stifIYT7gBHAxiGE8SGEY4ArgD1CCO9W5b6i6rYtQwh/B4gxTiRNd3656nJJ1XUqoTJ53S0b83g+V6IaX8dzh8poTeCpqs8wL5POQav49vICSvCZMJT+85EkSZIkqSYVP4ImSZIkSeXCAk2SJEmSyoQFmiRJkiSVCQs0SZIkSSoTFmiSJEmSVCYs0CRJkiSpTFigSZIkSVKZsECTJEmSpDJhgSbVISGEA0MIU0MI61a77roQwnshhNVzZpMkaXH4HiclIcaYO4OkWgohBOBl4NUY43EhhO7AWcD2McZ386aTJGnR+R4nJUvmDiCp9mKMMYRwLvB4COE94DxgN9+4JEl1ne9xUuIImlQHhRCeB7YB9o0xDsydR5KkQvE9TpXOc9CkOiaEsBuwBRCAzzLHkSSpYHyPkxxBk+qUEMIWwNPAGcA+wPIxxrZ5U0mStPh8j5MSCzSpjqjqavU8cEuM8ZIQwi+B0aT5+cOyhpMkaTH4HifNYoEm1QEhhMbAcOCZGOPx1a5/AFgnxtg6WzhJkhaD73HS7CzQJEmSJKlM2CREkiRJksqEBZokSZIklQkLNEmSJEkqExZokiRJklQmLNAkSZIkqUxYoEmSJElSmbBAkyRJkqQyYYEmSZIkSWXi/wEIwrOxcC9FRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# for this problem we can define f as 5 - (2*x + x*np.sin(x-3)), or (2*x + x*np.sin(x-3)) - 5.\n",
    "# it makes no difference to the root-finding, just the plots\n",
    "def f(x):\n",
    "    return 2*x + x*np.sin(x-3) - 5\n",
    "\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))\n",
    "fig.tight_layout(w_pad=4, h_pad=4)\n",
    "\n",
    "x = np.linspace(-10. ,10. , 100)\n",
    "ax1.plot(x, f(x),  'b', label='$y(x) = f(x)$')\n",
    "xlim = ax1.get_xlim()\n",
    "ax1.plot([xlim[0], xlim[1]], [0., 0.], 'k--', label='$y(x)=0$')\n",
    "ax1.plot(2.7903546180675676, 0., 'ro', label='intersection')\n",
    "ax1.set_xlim(xlim)\n",
    "ax1.legend(loc='best', fontsize=14)\n",
    "ax1.set_xlabel('$x$', fontsize=14)\n",
    "ax1.set_ylabel('$y$', fontsize=14)\n",
    "ax1.set_title('Fixed point as a root of $f(x)$ - wider $x$ range', fontsize=14)\n",
    "\n",
    "x = np.linspace(0., 5., 100)\n",
    "ax2.plot(x, f(x),  'b', label='$y(x) = f(x)$')\n",
    "xlim = ax2.get_xlim()\n",
    "ax2.plot([xlim[0], xlim[1]], [0., 0.], 'k--', label='$y(x)=0$')\n",
    "ax2.plot(2.7903546180675676, 0., 'ro', label='intersection')\n",
    "ax2.set_xlim(xlim)\n",
    "ax2.legend(loc='best', fontsize=14)\n",
    "ax2.set_xlabel('$x$', fontsize=14)\n",
    "ax2.set_ylabel('$y$', fontsize=14)\n",
    "ax2.set_title('Fixed point as a root of $f(x)$ - narrower $x$ range', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 4.3: </span>\n",
    "    \n",
    "For $2x + x \\mathrm{sin}(x-3) = 5$, use the subinterval $x \\in(a,b)$ you found in Exercise 4.2 and complete the code below to implement a root bracketing algorithm. Derive the concept from the Figure above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bracket =  (2.700000000000001, 2.800000000000001)\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    return 2*x + x*np.sin(x-3) - 5\n",
    "\n",
    "\n",
    "def root_bracketing(f, a, b, dx):\n",
    "    \"\"\" Function to perform root bracketing on the function f(.)\n",
    "    between a and b, with fixed interval size dx.\n",
    "    Returns the bracket of size dx that contains the root.\n",
    "    \"\"\" \n",
    "    # The sign function returns:  -1 if x < 0;  0 if x==0;  1 if x > 0.\n",
    "    sign = np.sign(f(a))\n",
    "    while sign == np.sign(f(a)):\n",
    "        a += dx\n",
    "        if a >= b:\n",
    "            raise RuntimeError('no root within [a,b]')\n",
    "    return (a-dx, a)\n",
    "\n",
    "\n",
    "a = 0.\n",
    "b = 5.\n",
    "dx = 0.1\n",
    "# print out the output from our root_bracketing function\n",
    "print('Bracket = ', root_bracketing(f, a, b, dx))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAFtCAYAAAB1DwLeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5xU1f3/8deHooCoqBQrrL1HJWtv+FWjxvpTE6OIYqLYS6KxYcFuookmMYooigpGjSX2EhVE7ItgrFiXoiigKEqR9vn9ce64s3fv7Mzszu7M7Lyfj8c8ZvfeO/eeuXNn7ueee87nmLsjIiIilaldsQsgIiIixaNAQEREpIIpEBAREalgCgREREQqmAIBERGRCqZAQEREpIIpEBBpw8xsiJm5mfUrdllyFZV3TJG2XRVtf0Qxti9SDAoERNJEJ4H0xxIz+8bMxpjZQDOzIpVrjJkp6UcBFDPQEClFHYpdAJESdUn03BFYD/h/wK5ANXBKsQolLe5zYGPgu2IXRKS1KBAQSeDuQ9L/N7MdgbHASWb2F3f/rCgFkxbl7ouAD4pdDpHWpFsDIjlw95cIJwgDfh6fb2btzOwEM3vDzH4ws7nR3yeaWeL3zMx2N7OnolsPC8zsQzO72sxWTFumKrolsGv0f/ptizH5vAczO9rMJpjZfDObYWa3mdmqCcuNida/jJldZGaTzOzH1H1zM1vRzP5oZs+b2TQzW2hmM83sETPbrpHtbxRtszZa3wwze9HMTsyx/H80s6Vm9pKZrRybd7iZjTaz2dG+fN/MLjCzZdOWGZh2e2XX2L4cEi2T2EbAzEZE06vM7HgzezvazldmNiz9M4u9bq+ovHOjz/k/0X74aX25vHeRlqQaAZHcpdoHLEqYdxdwBDAVuBVwwu2EG4GdgP71VmR2PHATMBf4NzAD6AecA+xvZju6+7fAt4TbFAOBPtTdsgCozaPsvwd+AdwLPBWV6Rign5lt6+4zE17zALA18CTwn6iMEKrOryDUkDwOzAZ6AwcA+5jZ/u7+VOz97hu9z2Wj7f8L6AZsAZwd7YtEUSB1PXAq8BBwhLsvSJs/HPgtMA14kLDPtgMuA3Y3sz3dfTEwkbD/LgYmAyPSNjMm0/Zj/gzsBTwKPAPsBhxHuH30f7FyHwbcDfwI3AdMB3YAXgHeynF7Ii3P3fXQQ4/oQTiBe8L0XYAlhB/11WLzDo9e9ybQNW36ckBNNO+ItOl9ovXMATaKrevGaPlhseljksqVw/sZEq1vIbBVbN510bzhSdsC/gd0T1jnihmmrwl8Abwfm96dcM99IbBr0usSPoMx0d+dCAGJA/8A2sWWHRjNexDonOG9n55p/QllqYrmj4hNHxFNnwL0TpvegRAQObBN2vTlCQHSj8AWsXVdnTrOgKpiH/N66KFbAyIJom53Q8zsCjO7F3iWUCNwlrtPjy3+2+j5XHf/ITXR3ecSrvABjk1b/khgGeAGd4/fjx4MfA8MSK/WLoC73H1CbNoQwgn6iAzbutDdZ8Unuvt3GaZPA+4HNjKz3mmzjgZWAG5y9xcyvK6BqPr/WULNyrnufqq7L40tdjqwGPitu8+PzbsM+JpYbUwzXeruU1L/eKhpuD36d5u05Q4k1HiMcvf41f/lhFoLkZKgWwMiyS6O/e/A79z99oRl+wJLSa5efoFQk7BVbHmA5+MLu/tsM5tAqIHYiMJVISedgL8zs4mE9gcbE6rO072eaWVR48nTge2BnoTAJt0ahKtnCNX0EG4x5KoX8BKwDjDA3UcllKEL4dbCLOCMDD07fyS8t0KpSZg2NXpeKW1a6vMeF1/Y3X+I9nu/ApZLpMkUCIgkcHcDMLPlCCe74cBQM5vs7vET+IrAN+6+MGE9i81sFuFkmb48hHvGSVLTuzW1/Am+yjD9y+g5qbHblwnTMLP/R7jyXwD8F/iE0NZhKeHktiuhLUBK6n18nkd5VyXUIkwDXsywzEqEWpoeNAzcWkrSlfzi6Ll92rTU/sy03zNNF2l1ujUg0gh3n+vuzwL7E37o74iuRNN9B6xsZh3jrzezDoR75HNiy0M42SVZLbZcIfTKMD1VhgbbcvdMCYwuI9zvr3b3g9z9THe/yEOXy0kJy6dOnmvkUd63CLcU1gDGmtk6CcukyjzB3a2xRx7bLZTU551pv2eaLtLqFAiI5MDd/wfcQmgQ9/vY7AmE79IuCS/dhRBAvBlbHhKqhs2sG7Al4Wr7/bRZS6L57eOvydGuCdtaMcO2slkPeM/d670mat2/U8Lyr0bP++SxDdx9JPAbYHVCMLBBbP4PwLvApvHuhFkspf7Ve0tIfcYN9oeZdSXsd5GSoEBAJHeXE06aZ5lZ+v3g26Lnq9JrC6K/r47+HZ62/EhCF8RTzWy92DYuI1SJj3T3H9Omfx0996ZpBpjZVrFpQwhV2P+KbSubWmB9M1s9NcHCDfqLgU0Slr+DcIV8opk1CJbMbM1MG3L3+4FDCLUqL5jZprFF/kpon3BbFETF172SmfWNTf4aWCvTNgvkYUKNRX8z2yI27wIKe9tHpFnURkAkR+7+uZndTGgkdzZwXjT9bjM7EPg18K6Z/YfQuPAgYG3gvvTGbu5ea2ZnAP8E3jSz+4CZhKv27QmJi86hvueAXwEPmtkTwHxgsrvflWPxnwReirY1nXCluhPhpH5uXjsidDscCkwwswcIQc2OhCDgUcJtlJ+4+ywzO4LQrmC0mT1J6Jq4AvAzwkl57Uwbc/dHo/37EDDGzPZItcR399vM7OfAScAnZvY0oZHiytE6dyG06j8hbZXPAb8xs0eB8YR7/GPdfWye+yEjd59jZicRgr6X0/b7DoQGji8QPu94LwiR1lfs/ot66FFKDzLkEUib34vQMG4u0CttejvCyagGmBc9xgMnE+v7nvaaXxCS0qT6m39MSFjTLWHZ9sCVwKeEE2/GvvCx1w2Jlu1H6HM/kRBEzCScIFdLeM2YxvZBtExqXXMJrfYfAjZP317CazYF7iQ0GlxIaDD3AjAo4TNo8N6i9/A98A2wdWzefsBjhKRHCwkNHV8n1OLEczX0JCT6+Ypwy8WBIdG8KhrPI1CVoVw/rSM2bx/g5eh4mE2oKdgoKqsnfdZ66NHaD3PXgGYiIq0laufxKbCsu2dqMCrSatRGQESkBZhZt3gPk6gtxQWEth4PFqVgIjGqERARaQFmtjdhbIdnCG0xuhKSK21JSEJU7e4zMq5ApJUoEBARaQFmtjahjcKOhKRHHQgJkh4DrnR3JRWSkqBAQEREpIKpjYCIiEgFUyAgIiJSwcoyEDCzO81sRjQgTLHK0M/M3MyG5Lj8wGj5gS1crlozqy3GtsudmY2I9lNVEbZdFW17RGtvu6mKub8qWSGPlaTfi1Kn465xZnammS0ys41yfU3ZBQJmVk0Yz/1qD+O9SwkysyHRl7VfscuSUoplkszayg++AvH86HvauBwuQm8kJNa6Ntd1lmOK4SsJectvKnZB8vQQYfCVTEPPttVti0hhfQ5sTGFGp9y9AOuQEuLu883sb8CfzGwHd38522vKqkYgGn1sD0Lu9vnFLk8+3P07d//A3Qs5tGzJb1tECsvdF0Xf52YH9u7+ibt/UohySUkZSUiffVJOSxc7x3E+D8JIbg7snjCvXzRvCLAN8DghJ3m9/OCEYWRvIKT4/JEwEtkjxHKXpy3fizBy3FeEHO0TCeOk/7S9HMs+MFp+YGx6bfToAlxDGDAllXf+HKIunrHXGHAKYQjWBYQrhBsII8nVArW5bDttf/wd+Cha1zeEHO0XxpbbDRgGvEeokZkPvEMYca5TwnvypEdsuS6EgXtSOet/AF4BDs/wno8m5G2fGZV1KvA0cFgO+z9rmUjLJw8cD7wdbeer6L2vmGHdeR1TGdZRFW17BCEX/X+iz2IuMA74RWPHFLA3YYyA72Lv6SDCj8KHaft4PHAamcdA6BIdezWE3P4/EIYp/jv1x1f4aX/FXr8F4ZicA+wZm7dR9Lqp0b76ipD3f8PYcomfFbFjO8s+PYwwwNA30edYC/yLkMgnfbllCQMv/Y8wJsAc4EXg11k+pyrgHsJYCwui/bVfbPkxjbyXqmiZ1YGLgJcIYyQsBL6I9svGjZUhNv2nz4Mcj1+y/F4QvvdjouNgDuF3tUGZotdtADxAGFNhLuG7ui+N/P6U+fd0A0LCqBmEwaP6pS23PvXH1fgi+n/9DOtdEbgKmBS9l9mE37Y9MnzGSY9+sWVHR+taIdv7KrdbA3sQopxXG1lme8LJZRxheNjuhA+CaDjSZwgjkz1NSPHZnfBjOc7M/p+7P5FakZmtQjiY14nWNw5YjTDy2jMFfF8do/WtThglbnFUpquBTsAlseWvJ/yQTycc+IuAA4FtCUOyLsxlo1F7i6cJ+2MsYX90IYwiN4QwJG7KOYQf8ZcJPwadCIlShgD9ohHhlqSV7yDC6Gp3EL7c8W13A54HtgLeJHxW7YC9gLvNbFN3vyDtJVcQPtfPgPsIJ7zVgK0Jo/Ldm+XtZi1Tmj9H5XiU8LnsBhwHrAf8X+x95HVM5WBtQjD0DnAz4T0eBjxpZke4e9L7PJQQCDxJODar0uZdTfiReo3wo7Ri9B7+Rth3A2LvZyXCD8gWhB+l2wjH07rAb6P3lzERjpntHi0zF9jF3Semzds7mteRsG8/Jvw4Hwzsa2a7ufub0eKXEPbhFlFZv42mp54zitL43k4IHGdF25wZbWu36H3VRMsuQ/jcdiWM+vhPwnfgUOBeM9vS3c9P2EwfQsD8KXAX4fM/DHg4+i6MjpYbEZX5QMKAQxPT1pF6L7sQApHRhBPpD4QTyaHAAWa2o0ejLeYo5+M3i/2icqeOq02AXwJbm9km7j4rtWDUMO0lwn54nBBUrUO4LZnP8V8u39N1Cd+pD4FRQGdCoISZbQ08CyxPCDTeI/x29gcONLPd3b0mrWzdCPtuE+CNaB90J4xm+oyZnejuN0eL/yd6PpowYNeYtDLVxsr4EuGCdRdCEqvMco2Eiv0AliOcIN/OML8fdZHR8QnzOxB+eBYAu8bmrU74kZxOGAgkNX1YtL7rYstXUzcC3JAcyz+QzDUCTviydE6b3pPwQ/Et0DFt+g7R8h8DK6dN70Q4gTS4akraNiFg+CyafkRCedeK/b8OybUTl0XrOCw2fQgZRqGLRbZnx6Z3Ap4inLy2TJv+NSErW5eEdXXP8TPItUxTgN6xY2dsNG+b5hxTjZStKu34vSbD8TabtOg+7XNdCuydYb3rJkxrR/iRdWDb2Ly7o+k3EasxIPywrZj2f2p/VUX/H0kIGt4D+sReu1JU/lnAJrF5mxJOfm9m+Dyqcvl80143KHrd68SuDgmjOK6W9v951H3/OsS+f7XRvB0yfE4Xx9a9V2pd2b5/sfk9geUTpm8R7ZcnMxwrI5pz/Ebzasn8e7GYWO0r4ao16Xv7XDT9xNj0fdL2V+L7L+Pv6ZUJ841Qe+ZA/9i8w6LpH5D23SIE/B49W9r09QkXPD9Sv1a7HzmcewhBnAN/zvq+8vmCFfNBqIZx4JkM81M7Z0KWnXJNhvmnR/N/Gf3fkXBVM4fkKrURuXwYacunvlwDY9Nro+nrJbwm9WO9Wdq0W6JpxzSyD2qzbRs4JJr2cDM/l1Wi9dwWm57xyxy9ZjHwRoZ1bhE/gAmBwGe5fFkbKWvGMsU+02MT5h0TzTulqcdUlrJVRct+S/JJIVW2oxM+14easC/6Rq+9KG1aT0KN2xfAcjmsI1WmKkKN0VJClfpKjeyLkzOs67po/iZJ68/zvb0dvW6rHJb9KCr3Rgnzfhc/ttM+p1qgfcJrJgOzYtMafP/yeC+PEE5gHRPKMKI5x280vZbMvxcjE9azdjTv/rRpa0XTPiLhdhPw33zeP+XxPf2ShN8iQi2pAy9neP2L0fxdov9T55nvSbuwS1s+daGV/j3tR26BwLbRcvdke1/ldGtgleh5dpblXs8wffvouU+GbhfrR88bE64ONiJUEb7oyY3sxhCqZwrhO3f/OGH61Oh5pbRpfaPnFxKWf5Fwgs3FdtHzk7ksHOVsOB34f4SgbHlC9JuyRo7bhVAl3R7I1AWmY/S8cdq0UcCpwLtm9m/C+38lw2fTXDUJ05I+i3yPqVy86e7fJ0wfQzjetiIEiOkyHfOp21t/JFTprkOoWUuX/rltTagtGOv5dc29jlDF+gBwpLsvSFgmta+2yLCvNoieNybUKDRJdJxuBnzl7hOyLLs8oRr5c3f/IGGR56PnrRLmTfS6W2HpplL3XnNmZvsCJxBqf7rTsEdXd3Lv9ZPr8Vuo9WwZPb/i7ksTXjOOcFu30Ir5PX3L3X9MmJ76fX4+YV5q+k6EY2osdeeZl9z9mwzLX0DyMZhNan3dsy1YToFAqpdApyzLfZlheiqQ+FWW13eNnleMnjPdD820nabIdN8zdVJvnzYtY7ncfYmZfZ3jNrtFz59nW9DMOhIOyG0I967vJdxzXRQtcjGhwVWuUp/F1tEjk65pf/8e+IRwn/rc6LHYzJ4AzswQSDVV0ueR9Fnke0zlItvxtmIj8+qJ7j2+QbiKe53QWOkbwnvpRgjs0j+3nI+JmF2i58cyBAFQt6+Oy7KufPZVknzeQ2pfZjrBpqZ3S5jX2Hc2r95YZnYaoR3EbMLV8xRCo0Wnrp1EPt+vXI/fvNfj7otDE4zcfpOyTG+uYn5PM/3+53tMNecYzKZz9Jy1h105BQKp4TpXaXSp8OVJkrpyPNDdH8lhe6nle2WYv2oO62gJ6eX6NH2GmbUn7J9cfgRTX6JcruQPJAQBd7j7wNg2VyMEAvlIvYfr3P0Pubwguvr6G/A3M+tJiKp/Q/hybxo1LkyK0FtSvsdULrIdb0k1IJmO+WMJQcAl7j4kfYaZbU8IBNLlc0ykO4jQqHC4mXV091sSlkmVewt3/1+e689HPu8hVaZM3+XVYssVnJl1IDSM/BLo67EugdHnVOrmRM+Zjt1M01tLS3xPs51ncj2mWvIYTJ0rsw51XU55BKYTrkI3bOLrUz0Nds5x+Q8IUfmWZpZ0FdavieVorlSr6l0T5u1M7sFdan/sk8Oy60XPDyTMSyoHhHvNkHwF8jrhvmyun0U97j7D3R90918TairWJVQHZ9NYmZoi32MqF32jKuu4ftFzo9XdMfl+bqnPZZc803dPJdQKTAJuNrOTE5Zpyr7K+/OKbmm8A/Qys0arU6NbMJ8Aa5jZ+gmL7BY9v5kwLx+NvY/uhKu9lxOCgK7UVTWXstQxub2ZJZ1TdspzfeXwPc0ktS/6ZZifmp46piZRd55Jum2TdAzmun9SKYYnNroUZRQIeGj9MBbobmbrZVs+wcOEL/3JZvbLpAXMbHsz6xJtbxHhvvTyhMYr6ctVE7qCFMOI6Hmwma2cmmhmnQgtenP1KKGh0AFmdnh8ppmlX1HVRs/9YsusA/wpw/pTtyh6x2e4+wzCvq02swujq6L49te1MJ47Zrasme0edQtLX6YjoTsQhC9TNhnL1ER5HVM5WpHQpzx9Hanj7TtCd6xc1UbP/WLr24rQWr4ed59J6Be/GnBt/EfdzLpmCIqJTmK7Ehrq3WBmZ8YWuZ1wtX6xmW0Tf72ZtUtIKdvUz+vv0fPN8fJG21ktbdJthLYu10Q1aqnlugMXpi3THI29jxmEY/fn0Yk/tf2OhBqwrPd3i83dpxLasKxH6Nf/k6jLaL7tA8rhe5rJS4ST+05mdmhsG4cSAuYPCe0mcPeFhN/CrsClseXXJXQTX0ToopqS6/5JtQMbna3Q5XRrAMKVzSGEbjp53RN290VmdjChD+njZvYyIVKaR2j1ujWhMdVq1J1Uziek4Dwj+jFO5RE4jNCo5IDmvqF8uftLZvYPQsO5d8zsfuryCMwmxwZF7r7QzH5F6Ft7t5kdT4icOxEazexO3fGR6vP9BzPbnBD19ib0M36c5ANyNOHq8ioz2ywqG+5+eTT/FEIjnUuBAWY2jnAvcfVo+1sDhxN6CnQm9MutNbPXCC2zOwF7Rss+4u7v5/C2s5UpL008prIZCxxrZtsSflRSx1s7QrfYOY29OOZOQkPB681sN0Kr7vUJn9uD0XrjTiHUrpxAyA/xNKFL4NqE790B1O+7/BN3nxlt52lCINHJ3a+I5n0d/RA+BLxqZs8REmItJRw/2xOqMtPbAD0Xlf+W6Dj/AfjW3W/I8r5vJVyFHgV8ZGYPE2oTVyf0L7+NuuD+WkKt2IHAW1Gbky6EW049CT1XxmXZXjavED7/M6LgPXXP/B/u/p2Z/Z3Q5uXtqKzLEK4EVyYcs7slrLPUnEw4Xm+MTrapPAKHEE7EBxI+61yUw/c007bczI4mtPW4N/o8PyDUZB9E6B1wVKxR5bmE2opTohwEo6nLI7A8oQfEZ2nLTyLc/v2NmS0ktClx4C53nwwh4CX8hk9y93dyKXjZPAhfkC+B1xLm9SO3LhU9CUlW3iF86D8QfiDvJ/SD7hBbflXCD8dM6jILDsx1e2nrGUjm7oO1GV4zhIRuNNRlFnyf0Mf0C0IilKZkFuxNGKTiM8IP/teERBmDY8utRYhcP4/2w7vA2YRgwYExCes+Mtpf86NlPOHzPIWQpCjVX3YK4QRwBrCK13WxOZvQw2EKoTvVTELgcgKwTB7HUMYy0Uh3tcY+73yPqQzlqorWP4IQ3DxM+AGcR/iB3SvXYyq2zCaELmgzCN2UxhPaDvy0vYTXLAcMpi7T3veElvzXAz2z7S9ghajMDlyW8D5voC6T5RzCD+VdwEEJZfkDdce5k19mwf6E3iXfRdv6jHAM940t14kQ9L8THRffE4L+pAyXGfdbNH8MseM8mr43ISD4IXXcUZd/oUP0Pt+Ltv9ltD/6JO3jTGXI9Hk0dvyS5+9FND/T930jQoD5bXSsvULILHhW9JoGn285f0+zLLdh9BlOJ1yoTSdk+Nwww/LdCLWrH0XH+reEYKJBRtFo+a0Jv5PfEYKmeucJ4BfRtDNy2d8WvahsmNl5hIGH+nqW7kEiIlJcZjYKOIKQq2FSsctTCczsAcKtunU9hy7W5RgIdCJUjfzP3fcvdnlERCpdVBXd092/jE3fnVAlP8ndNy1K4SqMmW1JaFx4mme/jQaUXxsB3H2BmQ0AdjOz5Ty/xCciIlJ4ywBTzWw04VbPYkLq6D0JtxyTepJIy1iN0NB1aK4vKLsaARERKS1Rj4vrCY0x1yQ0uJxFaPx6tW7jljYFAiIiIhWsbPIIiIiISOGVXRuBxnTv3t2rqqqKXQwREZFWMX78+Fnu3qM56yh6IGBmtxESnMxw982iaSsTBrapIvRz/bW7Zxt1kKqqKmpqkgakEhERaXvMbHJz11EKtwZGEJJtpDsXeM7d1yckTTi3tQslIiJSCYoeCLj7WOrGTU45kLox1+8gpGYUERGRAit6IJBBL49G4oqeexa5PCIiIm1SqQYCOTOzQWZWY2Y1M2fOLHZxREREykqpBgJfpYYKjZ5nZFrQ3Ye5e7W7V/fo0ayGkyIiIhWnVAOBR4Cjo7+PJozEJiIiIgVW9EDAzP5FGK5yQzObZma/IwwVuaeZfUTIVX11McsoIiLSVhU9EHD3w919NXfv6O5ruvtwd//a3Xd39/Wj53ivAhERkbIxahRUVUG7duF51Khil6hO0RMKiYiItGWjRsGgQTBvXvh/8uTwP0D//sUrV0rRawRERETKXWNX/IMH1wUBKfPmhemlQDUCIiIizZB0xX/ssTB2LLiH/5NMmdJ6ZWyMAgEREZFmSLriX7AAhg2Dbt2gU6fwf1zv3q1Tvmx0a0BERCSLpKr/GTPgppsyX/Gbwddfw623Qpcu9ed16QJXXNHSpc6NagREREQakVT1f9RRodrfHTp0gMWLG76ud+8QOKQaBA4eHG4H9O4dgoBSaCgIqhEQERFp1PnnN6z6X7oUll8e3noLRozIfsXfvz/U1obX1daWThAAqhEQERFJtHAh3Htv5kZ9338PP/tZeEDpXvFno0BAREQq3qhRdSfyNdaAHXaAcePgiy8ar/pP6d+/fE78cbo1ICIiFS3VBmDy5HDPf9o0uO8+WGkleOIJuP320m7s11yqERARkYp27rkN2wAA/PAD7LNP+NusfKv+s1EgICIiFWnuXPjLX0INQJL0tgHlXPWfjW4NiIhIm5eeB6BPn5D5b/314eKLoXPn5NeUSsKflqZAQERE2rR4G4ApU2D4cFhuudAg8JZb2nYbgGwUCIiISJuWlAIYQvfAHXcMVf7DhoWaArPwPGxY270VEGfuXuwyFEx1dbXX1NQUuxgiIlIiFi2CZZZJnmcWEvyUMzMb7+7VzVmHagRERKRNGj8ett468/xKaQOQjQIBEREpa/EBgW6/Hc45B7bdFr76Cs44o7LbAGSj7oMiIlK2kgYE+t3vQqPA3/0OrrkmJAaqrm67eQCaS20ERESkbFVVJQ8D3LNnqA1o69RGQEREKlqmAYFmzmzdcpQzBQIiIlKW3KFbt+R5agiYOwUCIiJSdmbPhkMOCc/tYmcyNQTMjwIBEREpafFeAUOGwFZbwaOPwrXXwh13VG4yoEJQY0ERESlZ8V4BKd27w+OPwzbbFKdcpaLNNxY0s9+b2btm9o6Z/cvMOhW7TCIi0noypQfu3FlBQKGUbCBgZmsApwHV7r4Z0B74TXFLJSIirSlTr4BMQwdL/ko2EIh0ADqbWQegC/BFkcsjIiKtqHv35OnqFVA4TcosaGbdCVfrWwMOvArc4O7fFKpg7v65mV0LTAHmA8+4+zOFWr+IiJQud7j88pAPoF27+oMDqVdAYWWtETCzb8ysb9r/awETgPOA3kAVcBEw3sx6FapgZrYScCCwNrA6sJyZHZmw3CAzqzGzmpnKICEiUvbmzoXDDoOLLoIBA2D4cPUKaEm53BroRv2ag6uBZYBt3H1Td98E2A7oCgwpYNn2AD5z95nuvgh4ENghvpC7D3P3anev7tGjRwE3LyIirSG9e+Caa8Kmm8L998Of/xy6Bg4cCLW1oVagtlZBQKE15dbAXsCl7j4hNcHda8zsasLtgkKZAmxnZl0ItwZ2B9Q3UESkDYl3D/z88/B81lnwxz8Wr1yVpCmNBQMs4h8AACAASURBVLsRbg3EvQms2rzi1HH314D7o/W+TSjrsEKtX0REii9T98B//7v1y1Kpcq0RqDazrtHfM4EVEpbpBiR8nE3n7hcDFxdynSIiUjoydQ/MNF0KL9dA4B/Rs0XPuwKPx5bpCyQMBikiItLQ0qXQtSt8/33Deeoe2HpyCQR2S5j2XcK0tYF7mlccERGpBD/+CEcdFYKADh1g8eK6eeoe2LqyBgLu/kIuK3L3Bl37RERE4r77Dg46CMaMgWuugdVWC20FpkwJNQFXXKGeAa2pSQmFREREcjVqVN2JfvXVQzfB6dNh5Mi6E75O/MWjQEBERFpMpu6B55yjk3+pKPWxBkREpIxl6h54j1qUlQwFAiIi0mLUPbD0KRAQEZEWkynzu7oHlo4mBwJm9ryZrVnIwoiISNvx6KPw9ddhsKB06h5YWppTI9AP6FKgcoiISBvy73/DwQdD374wdKhGDyxl6jUgIiIFNXIkHH007LADPP44rLBC6DkgpUltBEREpFnShxFeZRUYMAD69YOnngpBgJQ21QiIiEiTxfMEfPNNCAiOPBKWW664ZZPcqEZARESaLClPwNKlcMklxSmP5E+BgIiINJnyBJQ/BQIiItJk3bolT1eegPKhQEBERJrk5pth9mxo377+dOUJKC8KBEREJG+33QYnnAD77QfDhytPQDlrTq+BPQHdBRIRqTB33gnHHgt77w333w/LLhvyBkh5anKNgLs/5+4LClkYEREpPel5Anr0CCf93XeHBx8MQYCUN+UREBGRjOJ5AmbNCgHBEUdA587FLZsUhtoIiIhIRsoT0PYpEBARkYyUJ6DtUyAgIiIZ9eyZPF15AtqOvNoImNl2wN7AdsDqQGdgFjAJeAH4j7vPLnQhRUSk9b31FsyZE7oFutdNV56AtiWnGgEzO9rM3gZeBs4AugAfAa8Bs4FtgVuBz81shJmtXYjCmVk3M7vfzD4ws/fNbPtCrFdERBr30Ufwi1+E0QSvv155AtqyrDUCZvYW0BO4EzgKmOieHhv+tNyKwH5Af+BdMzvG3e9tZvn+Bjzl7oea2TKEAERERFrQ1Kmwxx6hUeB//wsbbQSnnVbsUklLyeXWwO3A0Gw5A9z9O2AUMMrMtgBWbU7BzGwFYBdgYLT+hcDC5qxTREQaN3Mm7LknfPstjB4dggBp27LeGnD36/NNHOTub7n7000vFgDrADOB281sgpndamYa3VpEpMDSEwatsQZ88gk89hj07VvskklrKOVeAx2AvsBN7r4VMBc4N76QmQ0ysxozq5k5c2Zrl1FEpKylEgZNnhwaBC5aFAYRUvfAymEJt/sbLmS2FnAosAi4x91nmVlvwol5PeBj4K/u/nHBCma2KvCqu1dF/+8MnOvu+2Z6TXV1tdfU1BSqCCIibV5VVQgC4vr0gdra1i6N5MvMxrt7dXPWkUtjwY2BV4AVoknnmtnuwLNAV0IQMAA4zMy2cveCxJHu/qWZTTWzDd19ErA78F4h1i0iIoESBkkutwaGANOAjQi9B14DHgG+BKrcfWtCrcAMEqrum+lUQuPD/wFbAlcWeP0iIhVt+eWTpythUOXIpdfADoQq+Q8BzOxcQgKhw6OeArj7V2Z2PSHHQMG4+0SgWVUeIiKS7LrrQsKgDh1g8eK66UoYVFlyqRHoAaRXEtVGz5/GlpsErFWAMomISAu7+274wx/gkEPgttuUMKiS5VIjMJsQDKQsAcYDc2LLrYD6+YuIlLz//hcGDoRdd4WRI6FTJxgwoNilkmLJpUbgPUIKYQDcfam7bx014Ev3M+CTQhZOREQKq6YGDj4YNt4Y/vOfEARIZcslEPgTMCGH5foC9zWvOCIiUmjpCYO23RaWXRaefBK6dSt2yaQUZL014O7P5LIidz+4+cUREZFCSiUMmjcv/O8Oc+eG9MFqByBQ2pkFRUSkmQYPrgsCUhYsCNNFQIGAiEibpoRBko0CARGRNso95ARIooRBkqJAQESkjRoyJLQH6Nix/nQlDJJ0CgRERNqg4cPh0kvhmGPg9tuVMEgyyyWhUCIzex44yt2nFbA8IiLSTE8+CccfD3vtBTffHGoEdOKXTJpTI9APyHD3SUREiuHNN+FXv4LNN4d//7vhbQGRON0aEBEpc+kJg7beOiQMeuKJzCMLiqRTICAiUsZSCYMmTw69BJYuDXkDnn++2CWTcqFAQESkjClhkDSXAgERkTKmhEHSXAoERETKWKZ2AEoYJLlSICAiUqaGDYM5c6BDrCO4EgZJPhQIiIiUoaeegpNOgn32gdtuU8IgabomJxQSEZHieOutulwB994bbg8MGFDsUkm5ak6NwJ6AmqOIiLSiadNg332hWzd47DHlCpDma3KNgLs/V8iCiIhI4+bMCUHAnDkwbhyssUaxSyRtgdoIiIiUsPSsgb16wdtvw/33w89+VuySSVuhNgIiIiUqlTUwlTBowQJYZhmYObO45ZK2Ja8aATPr21IFERGR+pKyBi5cqKyBUlj53hoYbWa7tUhJMjCz9mY2wcwea83tiogUm7IGSmvINxC4G3jCzA6JzzCzncxsXGGKVc/pwPstsF4RkZLWs2fydGUNlELKKxBw9xOBq4B7zOwEADPb3MweBcYCKxWycGa2JrAvcGsh1ysiUuo+/TTcFjCrP11ZA6XQ8u414O6XAicAfzezF4AJwGbAb4HNC1s8rgfOBpYWeL0iIiVr9mz45S+hY0e45hplDZSWlXevATNbGdgAWALsDLwM9HP3xYUsmJntB8xw9/Fm1q+R5QYBgwB6q75MRMrcwoVw8MHw2Wfw7LOw885w5pnFLpW0Zfn2GrgY+BQ4GfgLoRagGvhr4YvGjsABZlYL3AP8n5mNjC/k7sPcvdrdq3v06NECxRARaR3ucNxxMGZMGD9g552LXSKpBPneGhhMaDC4rrtf4O4jgF8CR5vZvWbWsVAFc/fz3H1Nd68CfgM87+5HFmr9IiKlID1h0EorwZ13wqWXqvpfWk++twY2dvdP0ie4+/NRl8IngKeA3QtVOBGRtiyeMOi776B9e1hnneKWSyqLuXthVmS2HvC0u69bkBU2QXV1tdfU1BRr8yIieamqgsmTG07v0wdqa1u7NFKOzGy8u1c3Zx0FG2vA3T8GdijU+kRE2jolDJJSUNBBh9z9q0KuT0SkLcs0eqA6QElryhoImNnDZrZVris0s05m9odUwiEREWlowQLo1KnhdCUMktaWS43AFOBVM3vNzE4zs75mVq+RoZmtbmYHmdlwYDqhW+GbLVBeEZGyt3QpHH00fPwxnHaaEgZJcWXtNeDup5rZ9cAZwBBgRcDNbA7wIyGtcEfAgNej5e5yd2UDFBFJMHgw3Hcf/PnP8Mc/wt/+VuwSSSXLqftg1GXwVDM7E9ge2BZYHegEfA18AIx194T2ryIiknLLLXD11XD88XDWWcUujUieeQTcfSHwQvQQEZEsRo0KNQBTpkCPHjBzJuy1F9xwQ8MBhUSKIe+xBgDMbFl3/7HQhRERaUviCYNmzAgn/0MPhQ5N+vUVKbx8xxroZ2aTgXlmNtvMXjCz68zsKDPbzMwK2h1RRKScDR5cFwSkuMPllxenPCJJ8o1J/wnMA04BugNbAQcBp0fzFwBdClY6EZEypoRBUg7yDQTWBn7l7o+nTzSzbkBfYMtCFUxEpNyttVbySV8Jg6SU5FuV/wGhq2A97v6tuz/v7i0xHLGISNlxh/XXbzhdCYOk1OQbCPwVOLYlCiIi0pZcey089xzsu68SBklpy/fWwC7AxmZ2L3CRu09qgTKJiJS1++6Ds8+GX/8a/vUvaKdm1FLC8g0EdgR6E9oKHGpm04DxhHTC44E3NfCQiFSyl16Co46CHXeEO+5QECClL9+EQpua2bLA5oQeA1tGz+cAywEOtC90IUVEysGHH8IBB4TGgA8/nDyokEipyRoImNm27v5a6v8okVBN9EgtY8CGwBYtUUgRkVKVnjmwfftw8n/ySVhllWKXTCQ3uVRavWJm2wCY2XAzO8XMdjSzrqkFPPjA3e9tsZKKiJSYVObAyZNDL4HFi8Pj1VeLXTKR3Jm7N76A2dbAp+7+tZm9BWxE6EK4FPgEmJD2mOjuM1q2yJlVV1d7TU1N9gVFRAqgqioEAXF9+kBtbWuXRiqRmY139+rmrCOXYYjfSPt7CzPrCGxGaBuQevwS6IraCIhIBVHmQGkL8h72wt0XUVcDAPzURmADlFlQRCrIiivCt982nK7MgVJOCjL+lYf7C5Oih4hIm3f33SEIaN8eliypm67MgVJu8h198HYz+1ts2nZmtlNUKyAi0uaNHg0DB8Kuu8Lw4cocKOUt3xqB3YCfYl0zOxn4e/Tva2a2p7vPLVThRERKzdtvw0EHwQYbwEMPwUorwdFHF7tUIk2Xb86rXsBHaf+fDjxISD3cEzirQOUSESk5U6fCPvtA164hV8BKKxW7RCLNl28gMBvoDGBmmwHrAhe6+0vAlcBhhSqYma1lZqPN7H0ze9fMTi/UukVEcjVqVOgm2K4drLsufP11CALWWqvYJRMpjHwDgTeAQ6K/jwY+c/cPov8/BvoUqmDAYuBMd98Y2A442cw2KeD6RUQaFU8YtGhReH777WKXTKRw8g0ErgQGmNkHwO+Bu9LmrQEUrH2Au0939zejv78H3o+2ISLSKgYPhnnz6k/78ccwXaStyHfQodfMbCegPzADuCZt9h60UPdBM6siJC56rfElRUQKRwmDpBLkFQiYWd8o0+AbCbPnAwUfayAa0+AB4Ax3n5MwfxAwCKC3sniISAEpYZBUgnxvDYw2s92SZrj7Ke5+QwHK9JMonfEDwCh3fzDDdoe5e7W7V/fo0aOQmxeRCjZ8eF3CoHRKGCRtTb6BwN3AE2Z2SHxGlFRoXGGK9VPa4uHA++7+10KtV0Qkm0ceCY0E99pLCYOk7cu3jcCJZjYduMfMTnX3oWa2OaER4b6EBn2FsiMwAHjbzCZG08539ycKuA0RkXpefBEOOwyqq+H++0POACUMkrasKYMOXWpmnwM3mdnhhBP2VOC3wJ2FKpi7jwOUtlhEWs3bb8MBB4Qr/8cfD0GASFuX760BzGxlwkiDS4CdgVeB9d19hLsvLXD5RERaVHrCoK22CtOefhq6dy9qsURaTb6DDl0MfAqcDPyFUAtQDegevoiUnXjCoCVLYMECGFew1k4ipc/CCMI5Lmy2ELgVuMTdv4qm/R/wEPAUcKS7L2qJguaiurraa2pqirV5ESkzVVUhCIjr0wdqa1u7NCL5M7Px7l7dnHXk20ZgY3f/JH2Cuz8fdSl8ghAM7N6cAomItBYlDBLJ89ZAPAhIm/4msBNQVYAyiYi0uEWLoFOn5HlKGCSVJO/Ggpm4+8fADoVan4hIS1m6FI45BubPh44d689TwiCpNAULBABS7QZEREqVO5x6amgoeOWVcPvtShgklS3vPAIiIuVs8GC48UY4+2w499wQAOjEL5WsoDUCIiKl7E9/gquuCl0Gr746BAEilU6BgIi0aamEQWahBmC77UKNgIIAkaBJgYCZtTOz581s/UIXSESkUNITBqX8739wzz3FK5NIqWlqjYAB/YDlC1cUEZHCGjwY5s2rP23evDBdRALdGhCRNksJg0SyUyAgIm3S00+HroJJlDBIpE7O3QfN7KK0f1MBxCAz+yJ9OXe/tBAFExFpquefh4MOCif8mTND4qAUJQwSqS+fPALHJEzbH1iY9r8DCgREpGjGjoX994f11oPRo0PNwODB4XZA794hCFDeAJE6eY0++NOLzDoQAoDqaJyBkqDRB0Uq28svwy9+EU74Y8ZAz57FLpFIyyrE6INNbSOQf/QgItICUnkC2rWDnXaCrl3huecUBIjkSo0FRaRspecJcA+POXNCGwERyY0CAREpW0l5AubPV54AkXw0KRBw9yXAbsCkwhZHRCR3yhMg0nxNHn3Q3V8oZEFERPLxyiuZ5ylPgEjudGtARMrOuHGhd0CPHtC5c/15yhMgkh8FAiJSVkaPhr32gjXWgAkT4JZboE+fMJpgnz4wbJjyBIjkI69bA2bWt5TyBohIZXnmGTjwQFh33dBFsFevcNLXiV+k6fKtERhtZru1SEkSmNneZjbJzD42s3Nba7siUhrGnTSKaR2qWGrtmNq+ijv3HsUGG4RagV69il06kbYh30DgbuAJMzskPsPMdjKzcYUpFphZe+CfwD7AJsDhZrZJodYvIqVt3Emj2OqmQay5ZDLtcNZaOpmbfRB/+fkoevQodulE2o68AgF3PxG4CrjHzE4AMLPNzexRYCywUgHLtg3wsbt/6u4LgXuAAwu4fhEpYVXDBrMc9ZMELMc8NrpTSQJECinv7oPufqmZfQ7cZGaHAzsCU4HfAncWsGxrROtNmQZs29gLJk2aRL9+/epN+/Wvf81JJ53EvHnz+OUvf9ngNQMHDmTgwIHMmjWLQw89tMH8E088kcMOO4ypU6cyYMCABvPPPPNM9t9/fyZNmsTxxx/fYP4FF1zAHnvswcSJEznjjDMazL/yyivZYYcdePnllzn//PMbzL/++uvZcsstefbZZ7n88ssbzL/55pvZcMMNefTRR/nLX/7SYP5dd93FWmutxb333stNN93UYP79999P9+7dGTFiBCNGjGgw/4knnqBLly7ceOON3HfffQ3mjxkzBoBrr72Wxx57rN68zp078+STTwJw2WWX8dxzz9Wbv8oqq/DAAw8AcN555/FKrD/YmmuuyciRIwE444wzmDhxYr35G2ywAcOGDQNg0KBBfPjhh/Xmb7nlllx//fUAHHnkkUybNq3e/O23356rrroKgEMOOYSvv/663vzdd9+dCy+8EIB99tmH+elD2AH77bcfZ511FkCD4w507DX32HtsyWQAbgTqHXlLJkO/fjr2dOzpd2/QoAZla4q8AwEzWxnYAFgC7Ay8DPRz98UFKVHaphKmNRjjwMwGAYMAll122QIXQUSKZQ4r0JU5Dab/yLLomy5SOHmNPmhmFwO/JwQQ1wMfA0OBYe5+WkELZrY9MMTd94r+Pw/A3a/K9BqNPihS/pYuXsqL25zJrhOuZzHt6cCSn+bNpQsTThzGTjeqm4AIFGb0wXxrBAYDtwKXuPtXUSGmAA+ZWS/gSHdf1JwCpXkDWN/M1gY+B34DHFGgdYtICVo0bxGvbX4su356Jy9seTrtt62m6tYLWH3JFL5o35vaQVcoCBApsHwDgY3d/ZP0Ce7+fNSl8AngKWD3QhTM3Reb2SnA00B74DZ3f7cQ6xaR0jDupFFUDRvM6kumML39mnzXsTs7LZjAmP+7lF3/ewHWzmDokQCsGT1EpLDyCgTiQUDa9DfNbCfCSbtg3P0JQoAhIm1MqntgqmfAGkumsvqSqYxddyD9nruwuIUTqSAFSzHs7h8DOxRqfSLStiV1DzRgndrRxSmQSIXKORAws35m1t/M+maYvwbQsB+JiEiC1ZckjxWcabqItIysgYCZdTWzl4HngLuAN8zsKTNbPbbomsDFLVBGEWmDPm+XPFbwF+01hrBIa8qlRuB8YGNgICHV78nAVsBrSvkrIk0xciScu/QK5tKl3vS5dKF2kMYQFmlNuQQCBwMXu/td7v6Buw8F+gJfAWPNbOsWLaGItBnucNllMGAAfNGvP6/+dhjT2vdhKca09n2UI0CkCHLpNdAbmJA+wd0/N7NdgUeB58zsAGB+0otFpHKNGgWDB8OUKbDWWlBVBWPHhkDg1lthmWX6w/Bw4lf3QJHiyCUQmEHC99Pd55rZPsCDhC5+1xa4bCJSxkaNgkGDYF7UMWDKlPA4+GC44w6wpCTiItLqcrk1UEOGUf/c/cdo3mPABQUsl4iUucGD64KAdOPHKwgQKSW5BAL/AvqY2SpJM6PBhg4DbgbU70dEgHD1n890ESmOrLcG3P0B4IEsyzhwYqEKJSLlzR26dYPZsxvO663egSIlpWCZBUVEAH78EY47LgQB7dvXn9elC1yh3oEiJUWBgIgUzPTpsNtuMHw4XHABjBgBffqENgF9+sCwYdBfvQNFSkq+ow+KiPwkvXtgr14wfz4sXgz33w+HHBKWOfLI4pZRRBqnQEBEmiTePfDLL8OV/5VX1gUBIlL6dGtARJokqXugOwwdWpzyiEjT5B0ImFk7M3vezNZP/7slCicipUvdA0XahqbUCBjQD1g+9reIVIh77sk8T90DRcqLbg2ISM7mz4cTToDDD4f11oPOnevPV/dAkfKjQEBEMho1KgwU1K4drLEGbLgh3HwznHMOvPsu3HKLugeKlDv1GhCRRPFeAV98EZ7/+Ee4+urwd//+OvGLlDvVCIhIokyDBt13X+uXRURajgIBEUmkXgEilUGBgIjUM28enHZayAmQRL0CRNoWBQIiFSy9MWBVFVxyCWy1FfzjH7DXXuoVIFIJFAiIVKhUY8DJk8PV/+TJMGQIzJoFzz0HTz2lXgEilcA8U/1fYy8y2xWocfe56X8XrFBm1wD7AwuBT4Bj3P3bbK+rrq72mpqaQhVDpE2rqgon/7i11lI7AJFyYWbj3b26OetoUo2Au7+QOvGn/11A/wU2c/efAR8C5xV4/SIVL9PJftq01i2HiBRXSd4acPdn3H1x9O+rwJrFLI9IW/PEE6FdQBI1BhSpLCUZCMT8Fngy00wzG2RmNWZWM3PmzFYslkj5+eIL+NWvYN99oVcvWHbZ+vPVGFCk8hQtEDCzZ83snYTHgWnLDAYWA6Myrcfdh7l7tbtX9+jRozWKLlI20nsFrLwyrLMOPPYYXH45fPYZDB+uxoAila5JjQVbg5kdDZwA7O7uCfnNGlJjQZE68RTBEAKCa6+F3/++eOUSkcIpRGPBvAIBM9sO2BvYDlgd6AzMAiYBLwD/cffZzSlQtJ29gb8Cu7p7zvX9CgRE6qy5Jnz+ecPpffpAbW2rF0dEWkCr9Rows6PN7G3gZeAMoAvwEfAaMBvYFrgV+NzMRpjZ2s0pFHADsDzwXzObaGZDm7k+kYrxww9w0UXJQQCoa6CI1Jd19EEzewvoCdwJHAVM9IRqBDNbEdgP6A+8a2bHuPu9TSmUu6/XlNeJVJpRo8LgQFOmhP7/e+4Jjz8OX34ZGv4lDRqkXgEiki6XYYhvB4a6+4LGFnL37wiN+kaZ2RbAqgUon4hkEG8DMGVKaPy33nrwyivwyScN2wioV4CIxGUNBNz9+nxX6u5vAW81qUQikpNMwwQvXAjbbRceqeWmTAk1AVdcoV4BIlJfLjUCPzGzvu7+ZksVRkRy8/rryemBAaZOrfu7f3+d+EWkcfnmERhtZru1SElEpIGk0QH32Qe23VaZAUWkMPINBO4GnjCzQ+IzzGwnMxtXmGKJSKbRAV98Ea6+OiT/6dKl/mvUBkBE8pXXrQF3P9HMpgP3mNmp7j7UzDYHrgT2Bd5viUKKVKLzz09uA7DyynDOOeHvTp3UBkBEmievQADA3S81s8+Bm8zscGBHYCphTIA7C1w+kYqzYEGoDchldEC1ARCR5so7EDCzlYENgCXAzoQkQ/3SRgsUkRyk5wDo3RvOPjv0/x86FGbOhI4dYdGihq9TGwARKaR8ew1cDPw+et1fgI+BoYR0wKcVvHQibVQ8B8DkyXDyyeHv/feHM86A6dOVB0BEWl6+NQKDCamEL3H3rwDMbArwkJn1Ao5094RrGBFJd955yff/V18dHnmk/jS1ARCRlpRvr4GN3f2kVBAA4O7PA7sBuwJPFbJwIuUs3vVv1KjQ4v+YY+r39U83fXr9//v3DwMELV0anhUEiEih5dtr4JMM0980s52ApwtSKpEyl1T1P2BA6Aa4/PLQtWsYHChO9/9FpLXlWyOQkbt/DOxQqPWJlLNzz21Y9e8Oq6wSrvqHDlUOABEpDVkDATN72My2ymVl7v6VmXUysz+Y2QnNL55I6YpX/Q8bBrffHjL/pXfxS/fNN7DccqGKf9gw6NMHzMLzsGGq+heR1pfLrYHJwKtmNpGQWfBF4H/p3QXNbHVgG2B/4GDgc0JeAZE2Kanq//jjw99rrw0rrABz5jR8XXrVv3IAiEgpyOXWwEJCY8DXgYuBN4AFZvaNmU03swWEhEIPApsCZwA/c/fXW6jMIq0iqbGfO0ycCKecktzqf9VVw/C/N96oqn8RKQ/m7o0vYLYY2MHdXzezKwk9A7YHVgM6AV8DHwBj3T3DeGito7q62mtqaopZBGkj4lf8AO3bh4Z+336b+XVmoYV/ah3q+iciLcnMxrt7dXPWkUuNwDfAStHf5wAL3P1P7n6Gu5/g7oPd/a5iBwEi+Uq64oeQ4vfMMxte8S9ZAgsXhnYAa66ZvM541b+6/olIqcslEBgHXGtmRwIGNF6FIFIGkkb2O+YY2Hhj6NYNvvoq+XXz58PAgWH0P1X9i0hbkEsgcArwJXAHIQh41sxeNLO/m9kxZralmXVs0VKKNEGme/xTpoQUvvEr/kWL4NNPw/3/Hj2S15m64lerfxFpK7K2EfhpQbNVgS8IKYa7AVsC60azFwHvARPc/XctUM6cqI1AZWnsHnyme/zLLZfcmj8ldY8/6fVduuhkLyKlpRBtBHLOLOjuX5rZQ8B17v5+VICuhIBgq+jRtzmFEclVUve9444LLfq7dYMrr0y+x79kCfzjHyFo+PLLhutNv+IHNfYTkbYv5xqBcqAagbYl0xX/Dz/A+usnn8iz0RW/iLQlrdVrQKRFZGq1n5p33HH1G/MddRR07x668GUKAsxC974+fZLn6x6/iEh9CgSkRTR2kk/NT2q1v9tusMsu4aQ/f3791yxdGqZdfnnjjflWXDHUHmRr1a/ufSIiJR4ImNlZZuZm1r3YZZH6sl3Nx0/yxx4LZ58Nf/87nHZa+D+p1f7YseE1qaQ8cfPnh9sF113X+IleV/wiIrkp2TYCZrYWoYfCRsDPhQZOAQAACBlJREFU3X1WtteojUDh5Nsiv0sXuPTS0A//yCNh9uzM6840BC/U3cOvqgoBRFyfPuHqPVsZRUQqQSHaCJRyIHA/cBnwMFCtQKD1ZDrRDx0K228PO+wAM2fmv16zMARvz55hYJ7GTvRqzCcikl2bbSxoZgcAn7v7WzksO8jMasysZmZTzk4VKlPVvjuce27Davt588J9+/XXbzwIeOmlxtPv9uoVAoJs9/BVtS8i0krcvSgP4FngnYTHgcBrwIrRcrVA91zW+fOf/9wlGDnSvU8fd7PwPHJk/XlduriH0354tG/vvt567iutVH96/HHbbe69eiXP69Mn8/q7dKlfhmxlFBGR7IAab+b5uORuDZjZ5sBzQOqadE1CRsNt3L3RnuO6NRAkVasvuywcemh4HjkyDJ4Tt+yyoeX+vfcm3+PPp9pe9+9FRFpem7w14O5vu3tPd69y9ypgGtA3WxBQaZKq9pcuhY8+Ss6j/+OPYZnHHksOAiBMv+mmkHmvudX26ponIlIeSi4QkCDf7nlHHRVO1htsALMyNKs0C6PqFSLZjk70IiJtQ85jDRRLVCtQUTLl0X/3XejcOeTRX7Cg/muWLoWOHeGf/4QLLwyt8+NSJ/orrkiu2o8n29HJXUSk7VONQJE0dsU/eHDDqv358+Gqq+CiixoGASlz58LvfgfXXKMW+SIikpuSrxFoizJd8b/ySrhPn9S/HsJJ+5tvYMstk5fJZ+Q8XfGLiAiUcEKhpiiXXgOZsuYBrLBCSLUbz7MPSrYjIiL1tcleA21BUrW/O7z3Xqi2b+yK/+uv4ZZbVLUvIiKtQzUCBZZ0td6hA3TrVteav2PHcNUfpzz6IiKSD9UIFEljDf3OO69hQ7/Fi0NDvptvDif222/XELkiIlIa1FgwT5ka+r34InzxBUydmvy6BQvC6yC3xnwiIiKtQbcG8tRYQ7/evUNq3u+/bzgvvdpfRESkEHRroIVkqvpftChcwScxCyf6m27KXu0vIiJSKhQIxCSl7z32WNhnnzC8bqYKlN69QzCgFv0iIlJOdGsgprGq/4MPDvOHDlUffhERKb5C3BpQY8GYxqr+H3gg/N23rxr6iYhI26BAIKZ378bT94LS84qISNuhNgIxV1yhxn4iIlI5FAjEqLGfiIhUEt0aSKCqfxERqRSqERAREalgCgREREQqmAIBERGRCqZAQEREpIIpEBAREalgCgREREQqmAIBERGRCqZAQEREpIK1qdEHzWwmkGHswCbpDswq4PoqkfZhYWg/Np/2YfNpHzZfofdhH3fv0ZwVtKlAoNDMrKa5wztWOu3DwtB+bD7tw+bTPmy+UtyHujUgIiJSwRQIiIiIVDAFAo0bVuwCtAHah4Wh/dh82ofNp33YfCW3D9VGQEREpIKpRkBERKSCKRDIwMz2NrNJZvaxmZ1b7PKUGzO7zcxmmNk7xS5LuTKztcxstJm9b2bvmtnpxS5TuTGzTmb2upm9Fe3DS4pdpnJlZu3NbIKZPVbsspQrM6s1s7fNbKKZ1RS7PCm6NZDAzNoDHwJ7AtOAN4DD3f29ohasjJjZLsAPwJ3uvlmxy1OOzGw1YDV3f9PMlgfGAwfpOMydmRmwnLv/YGYdgXHA6e7+apGLVnbM/n979xf61xzHcfz5yr80yQWh/dTU5BJr7cKvpCV/ZpE7ihvKjT/5U2o3ajcutRtE3JG1GlLCVkiy/JlGNGUk1mgXyNyYP28X56il78XOj/rs+Dwfdfqec65efS9Or/P5nM85eQBYD5xZVZtb55mjJF8D66vqhHoXgyMCi20ADlTVV1V1FNgO3Ng406xU1dvAD61zzFlVfVdVH437R4D9wOq2qealBr+Mh6eMm3c/EyVZAq4Hnm6dRf89i8Biq4Fvjzk+iBdgNZRkDXAZ8F7bJPMzDmnvAw4Du6vK/3C6bcBDwJ+tg8xcAbuS7E1yZ+swf7MILJYF57yLUBNJzgB2AvdV1c+t88xNVf1RVZcCS8CGJE5VTZBkM3C4qva2zvI/sFxV64DrgLvGKdTmLAKLHQQuOOZ4CTjUKIs6Ns5r7wSeq6oXWueZs6r6CXgLuLZxlLlZBm4Y57e3AxuTPNs20jxV1aHx9zDwIsM0dHMWgcU+AC5KcmGSU4GbgZcbZ1JnxgfdngH2V9WjrfPMUZJzkpw17p8OXAV83jbVvFTVlqpaqqo1DNfCN6rq1saxZifJqvGhX5KsAq4GTohVVRaBBarqd+Bu4HWGB7R2VNVnbVPNS5LngT3AxUkOJrmjdaYZWgZuY7gD2zdum1qHmpnzgTeTfMJQ8HdXlcvf1MK5wDtJPgbeB16pqtcaZwJcPihJUtccEZAkqWMWAUmSOmYRkCSpYxYBSZI6ZhGQJKljFgFJkjpmEZAkqWMWAUmSOmYRkDRZkrVJfkuy9R/nn0hyJMn6VtkkTWMRkDRZVR1g+Db9/UnOBkjyMHA7cFNVfdgyn6Tj5yuGJa1IkvOAL4HHGT7k8xRwS1XtaBpM0iQntw4gaZ6q6vsk24AHGa4l91oCpPlxakDSv/EFcBqwp6oeax1G0nQWAUkrkmQj8CTD56aXk1zSOJKkFbAISJosyTrgJYYHBq8EvgEeaZlJ0spYBCRNkmQt8CqwC7inqo4CW4FNSa5oGk7SZK4akHTcxpUC7zKMAFxTVb+O508CPgV+rKrLG0aUNJFFQJKkjjk1IElSxywCkiR1zCIgSVLHLAKSJHXMIiBJUscsApIkdcwiIElSxywCkiR1zCIgSVLH/gJJR4ozznTt7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_root_bracketing(f, a, b, dx, ax, xbounds=(-0.1, 1.4), ybounds=(-5, 6), flabel=''):\n",
    "    x = np.linspace(a, b, int((b-a)/dx)+1)\n",
    "    y = f(x)\n",
    "    # plot the sub-intervals in blue\n",
    "    ax.plot(x, y, 'bo-')\n",
    "    for i in range(1, len(x)):\n",
    "        if np.sign(y[i]) != np.sign(y[i-1]):\n",
    "            # plot the sub-interval where the sign changes in red\n",
    "            ax.plot([x[i], x[i-1]], [y[i], y[i-1]], 'ro-')\n",
    "    ax.set_xlabel('$x$', fontsize=16)\n",
    "    if not flabel:\n",
    "        fl = '$f(x)$'\n",
    "    else:\n",
    "        fl = flabel\n",
    "    ax.set_ylabel(fl, fontsize=16)\n",
    "    xlim = ax.get_xlim()\n",
    "    ax.plot([xlim[0], xlim[1]], [0., 0.], 'k--')\n",
    "    ax.set_xlim(xlim)\n",
    "    ax.set_title('Root bracketing\\n' + '(red indicates the bracket containing the root)', fontsize=20)\n",
    "\n",
    " # let's also use our plotting function from above.\n",
    "fig, ax1 = plt.subplots(figsize=(8,5))\n",
    "plot_root_bracketing(f, a, b, 0.1, ax1, flabel=r'$f(x) = 2x + xsin(x-3) - 5$')   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 4.4: </span>\n",
    "\n",
    "For $2x + x \\mathrm{sin}(x-3) = 5$, use the subinterval $x \\in(a,b)$ you found in Exercise 4.2 and complete the code below to implement a bisection algorithm. Derive the concept from a <a href=\"https://en.wikipedia.org/wiki/Bisection_method#Algorithm\">pseudo-code description</a> and compare the result to `scipy.optimize.bisect`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.790355086326599\n",
      "2.7903546180675676\n"
     ]
    }
   ],
   "source": [
    "def bisection(fct, a, b, atol=1.0E-6, nmax=100):\n",
    "    n = 0\n",
    "    while n <= nmax:\n",
    "        c = (a+b)/2.\n",
    "        if fct(c) == 0. or (b-a)/2. < atol:\n",
    "            return c\n",
    "        n += 1\n",
    "        if np.sign(fct(c)) == np.sign(fct(a)):\n",
    "            a = c\n",
    "        else:\n",
    "            b = c\n",
    "    raise RuntimeError('no root found within [a,b]')\n",
    "\n",
    "def f(x):\n",
    "    return 2*x + x*np.sin(x-3) - 5\n",
    "\n",
    "a, b = 0., 5.\n",
    "print(bisection(f, a, b))\n",
    "print(sop.bisect(f, a, b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 4.5: </span>\n",
    "\n",
    "For $2x + x \\mathrm{sin}(x-3) = 5$, use $a$ from the subinterval $x \\in(a,b)$ you found in Exercise 4.2 as initial guess $x_0$ and complete the code below to implement a Newton algorithm. Compare the result to [scipy.optimize.newton](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Newton (analytical derivative) used 10 function evaluations\n",
      "0.5671432904097811\n",
      "0.567143290409784\n",
      "Newton (analytical derivative) used 10 function evaluations\n",
      "2.7903546180673837\n",
      "2.7903546180673837\n"
     ]
    }
   ],
   "source": [
    "def newton(f, x0, dfdx, atol=1.0e-6):\n",
    "    \"\"\" Function to implement the Newton-Raphson method\n",
    "    \n",
    "    f is the function we are trying to find a root of\n",
    "    \n",
    "    and dfdx is another function which return the derivative of f\n",
    "    \"\"\"\n",
    "    x = [x0]\n",
    "    fevals = 0\n",
    "    while True:\n",
    "        x.append(x[-1] - f(x[-1])/dfdx(x[-1]))  # two function evaluations (f and dfdx)\n",
    "        fevals += 2\n",
    "        if abs(x[-1]-x[-2]) < atol:\n",
    "            print('Newton (analytical derivative) used', fevals, 'function evaluations')\n",
    "            return x[-1]\n",
    "        \n",
    "        \n",
    "###### case 1\n",
    "def f(x):\n",
    "    return x - np.exp(-x)\n",
    "\n",
    "def dfdx(x):\n",
    "    return 1 + np.exp(-x)\n",
    "\n",
    "x0 = -1.  # initial guess\n",
    "print(newton(f, x0, dfdx))\n",
    "print(sop.newton(f, x0, dfdx))\n",
    "\n",
    "\n",
    "###### case 2\n",
    "def f(x):\n",
    "    return 2*x + x*np.sin(x-3) - 5\n",
    "\n",
    "def dfdx(x):\n",
    "    return 2 - np.sin(3-x) + x*np.cos(3-x)\n",
    "\n",
    "x0 = 0.  # initial guess\n",
    "print(newton(f, x0, dfdx))\n",
    "print(sop.newton(f, x0, dfdx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 4.6: </span>\n",
    "\n",
    "Extend the Newton algorithm to compute $f^\\prime(x)$ using a finite difference approximation. Compare the result to `scipy.optimize.newton`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7903546180673837\n",
      "2.7903546180673837\n"
     ]
    }
   ],
   "source": [
    "def quasi_newton(f, x0, dx=1.0E-7, atol=1.0E-6):\n",
    "    \"\"\" Function to implement quasi-newton\n",
    "    \n",
    "    f is the function we are trying to find a root of\n",
    "    \"\"\"\n",
    "    x = [x0]\n",
    "    while True:\n",
    "        dfdx = (f(x[-1] + dx) - f(x[-1]))/(dx)\n",
    "        x.append(x[-1] - f(x[-1])/dfdx)\n",
    "        if abs(x[-1]-x[-2]) < atol:\n",
    "            return x[-1]\n",
    "\n",
    "def f(x):\n",
    "    return 2*x + x*np.sin(x-3) - 5\n",
    "\n",
    "x0 = 0.\n",
    "print(quasi_newton(f, x0))\n",
    "print(sop.newton(f, x0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <span style=\"color:blue\">Exercise 4.7:</span> \n",
    "\n",
    "For $2x + x \\mathrm{sin}(x-3) = 5$, use $a$ from the subinterval $x \\in(a,b)$ you found in Exercise 4.2 to find $x_0 = a$ and $x_1 = a+0.1$ and complete the code below to implement a Secant algorithm. Compare the result to `scipy.optimize.newton`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7903546180673446\n",
      "2.7903546180673837\n"
     ]
    }
   ],
   "source": [
    "def secant(f, x0, x1, atol=1.0E-6):\n",
    "    \"\"\" Function to implement the secant method\n",
    "    \n",
    "    x0 and x1 are the two required guesses\n",
    "    \n",
    "    f is the function we are trying to find a root of\n",
    "    \"\"\"\n",
    "    x = [x0, x1]\n",
    "    while True:\n",
    "        dfdx = (f(x[-1])-f(x[-2])) / (x[-1]-x[-2])\n",
    "        x.append(x[-1] - f(x[-1])/dfdx)\n",
    "        if abs(x[-1]-x[-2]) < atol:\n",
    "            return x[-1]\n",
    "\n",
    "def f(x):\n",
    "    return 2*x + x*np.sin(x-3) - 5\n",
    "\n",
    "x0 = 0.\n",
    "x1 = x0+0.1\n",
    "print(secant(f, x0, x1))\n",
    "print(sop.newton(f, x0))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "427.513px"
   },
   "toc_section_display": false,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}