Lecture4conv.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Programming in Python, Mathematics, Statistics, and Modeling\n",
    "## Lecture 4:\n",
    "This lecture contains the contents: \n",
    "* Nonlinear equations\n",
    "* Iterative methods to solve nonlinear equations\n",
    "    * Bisection method\n",
    "    * Newtons Method\n",
    "* Examples\n",
    "\n",
    "Second part of the lecture covers:\n",
    "* An intro to statistics\n",
    "    * Random Variable\n",
    "    * Probability Density Function\n",
    "    * Variance and Covariance\n",
    "    * Linear Regression (covered in Linear Algebra lecture)\n",
    "    * Simple Kriging\n",
    "    * Variogram\n",
    "* Some examples\n",
    "\n",
    "This lecture starts with a brief recap of the previous lecture. Then we will introduce the concept of nonlinear equations and suggest two iterative methods that we can use to solve the nonlinear equations (using Python of course!). We will then code together both methods and apply them to a example (our own square root evaluator, exciting!). We will then deal with a 2x2 nonlinear system, which we will solve using Newton's Method. In the second half of the lecture we will introduce some concepts from statistics and do some examples. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recap of lecture on numerical differentiation and integration:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_05/Slide29.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+sAAAF6CAYAAABhpcmGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXhM59vA8e9kkURIgiAhEksspdRO7NRWKlV0UUvttVNU+9ZapdVWlVqKokpRLVpb+dVOtGqrfSeEIAmRENkz8/4xmenMZNZsk8T9ua5ckznnOec8Z+LyzH3uZ1GoVCoVQgghhBBCCCGEyDMc7F0BIYQQQgghhBBC6JNgXQghhBBCCCGEyGMkWBdCCCGEEEIIIfIYCdaFEEIIIYQQQog8RoJ1IYQQQgghhBAij5FgXQghhBBCCCGEyGMkWBdCCCGEEEIIIfIYJ3tXQAhTUpOSWBYQQGpCAslPnlDIwwMnN7cM5VLi4qjSowevrFqV+5W0A2VaGkvKliUlLo6UZ88YHBqKZ/nyJsuH7d/P9p49SY6NJTUxETdvbxSOjhnKxUdEWDyXMRdWr+bgxIkkREVRtnlz3j5wwLYbymaX1q1j/7hxeaY+ppxauJB/FywgMToa1+LFqTNqFHVHjrT6+ORnzwiZNImrGzeSlpxMscBAms6YQUDbtnrlEh494tyKFVz++Wee3L6NMiWFYlWq8NJ771Fz0CAUCkV235oQQugx1Z6rlEpSnj6liJ8fFTp2JGjqVAqXLJmjdbC27czL7h4+zNY33iAxOpoiZcow5NYts+WPzprFqQULSIiKAsDNyGesTE6mkIeHxXMZs3voUK5u2kTCw4cETZtG0+nTbT5HdtozciRXNmzIM/Uxxto23JyYmzc5MH489/76C4AyTZrQ6uuv8apYUa9c9NWrnF68mJs7dpAYHY1KqcS3YUMafvQR/q1bZ+t9iRygEiKPO/fDD6qvQHXuhx+M7g+ZNk31x7vv5mqd8oKQadNUX4EqJjTUqvJ/vPuu2fK2nMuYpQEBqvUtW2b6+OxmrD6af0u39++3S500Dn38sWp+0aKqO4cPq1QqlSrs4EHVPHd31aGPP7bq+LSUFNX6li1VP7z4ourJ3bsqpVKpOj53rmqOg4Pq5q5demXXBgWp5nt4qG7u2qVSpqWpUhITVUc++UT1Faj2jBqV7fcmhBCmGGvPUxITVce++kr1Fah+atxYlZaamqN1sLXtzMvWt2ypWhoQYHX5pQEBJsvf3r/fpnMZigkNVX0FqpBp0zJ9juxkqj6a70L2ZEsbbsqTO3dUi0qVUm1/5x1Vcny8KvnZM9XWt95SLSpVSvXk7l1tuZTERNVXCoVqVa1aqocXL6pUKpXqWVSUanNwsOorhUJ1+ddfc+QeRfaRbvBCCJGLoq9e5dgXX1B37Fj8mjUDoFyLFtQdPZpjX3xB9NWrFs9x4ccfuXvwIG3mz6do2bIoFArqv/8+vo0asWfYMJSpqdqyKqWSRh99RIUOHVA4OODk4kKTqVMpExTEvwsXEnPjRo7dqxBCWOLk4kKDCRPwadCA+0ePEnHqlL2rJESOsqUNN+Xw//0fKc+e0fa773B2c8O5cGHaLVlC8tOnHP7oo/8KqlSgUtFuyRJKvPACAIW9vXll1Sqc3Nw4OGFCTt2myCYSrIt8r8nUqXRcudLe1cj3xqWk5NsugfnJ+R9+QJWWRuWuXfW2B77+Oqq0NC5YMZzj7PLluHh6Uq5VqwzniA0NJWzfPu226n37UvXNNzOcw7dxY1CpeHDyZKbuQwghspNnhQoAPMlEN2yRdeVatmSQPLzNFba04cYkP33KlV9+wb9NG1w8PLTbXb28KNe6NVd+/ZXkp08BcHByotH//R8+DRvqncO1WDGKV63Kk9u3eRYZmT03JnKEBOsi3wo7cIBl5cujcHBA4aD+p6xSKjnxzTesrF6dxaVL812ZMuzs35+4e/e0x11YvZrFPj587ejIz61acTckhHVNm7KwRAnmKBT8FhzMopIlmaNQML9oUZbpBLBhBw4wR6Eg6uxZ7ba1QUHML1KEbz08ODlvHgB3Dh5k29tvs6xCBRZ5e7O0XDn2jRlDclyc9rjH16+z2MeHeW5uzFEoeBoezpZu3VhStixzFAp+1vlP/PLPP/NDjRosKFaMldWr8++iRdn2Oe7s148j06fj4PTfFBa/tG3Lt56ezFEoCNMZ861bX2vdO3qUX9u3Z2Hx4iwsXpy1jRtzddMm7X5bPgdjkp89Y++oUSwqWZLvypRhQ+vWen8fjU2dO7NvzBgAtnbrxmIfHxb7+BB55ozV95Id7h46BAoFJWrU0NtesmZNQP1vx5zUxEQeHDtGiRo1tP/uteeoVSvDOeoMH45XpUoZzqN5cu/q5WX7TQghRDZ7fO0aAF6BgdptzyIjOfR//8eqWrX4zteXBV5ebOzQgQcnThg9R+iuXfzcooW6PfD15aeGDfnn88/NBiOhu3ax2MeHOQoFC0uUYG1QEKcWLtS2gQtLlGD7O+9oyx+ZPp1FJUuiUioBiLt/n8U+PswtVIhF3t7cP3YMgPOrVvFL27YsLVeOBV5erKxenVMLFqBSqbTnsvR95Ej6WGtlaiohU6eyxM+PxaVL81OjRoTt35+5D9qIZeXLc+fgQRzS57NJevLEaHsf/tdfevW11oXVq1ldp476e0CJEmzu0oXI06f19lvzOZgSc/Mmm7t04VsPD5b6+/NH374kRkfrlUlNSmKxjw9XNmwA0H4H+N5I+5iTbG3Djbl39Chpycl4p39vMDxHWlIS944eBdTBevPPPtP+bXUpU1NRODhQqGjRzN6OyAUSrIsC5c8hQzgyZQptFy5keEQEfU+fJvbmTdYGBWkb6xp9+zL8wQOKlivH07Aw/l24kG7btzPswQP8WrakVN26vHfnDk6FCxMYHKw32cqNrVsBuJ7+CtDr778p17o1wZs2UW/sWABOLVhAfEQEvY8fZ8TDh7y+bRvXfvuNXQMGaI8rFhjI8AcPqPrWWwDsHTmSZp99xtDwcBr93/9py13++We29+xJpeBghkdG0vfUKeLu3dM2ODnhzT17aDN/fobtuvW1Rti+fWxo2RLvGjUYeu8ewyIiqPb222zt0YMzy5YB1n8Opmx74w0urF5N199/Z9i9e3RctYqDEydqnyprdN+xQ3tPwZs3M/zBA4Y/eECpl14ye/4tOoG9pZ/LVvxNHl+7hmuxYji5uOhtd3J1xbVYMe0XVlNibt5EpVTi7uOTYV8RX1/tNSx5cPw4biVKULZ5c4tlhRAipyTGxPDXJ58Q+e+/1BoyhNJ16mj3PTh2jFPffsvLCxYw7P59hoSF4VayJL+8/DLPHjzQO8+FNWvY1KkTFbt0Ydj9+wwND6fWkCEcnjSJy+vWmbx+hY4deaFnTyq//jrDIyLo9fff1B05kte3bQOgxezZvKpz/I2tW0l4+JDw9Em9ivj6MvTuXdxKlGDwrVv4pmcwD06YQMlatRh04wYjo6NpOmMGB8aN4/R332nPZen7iMbeUaM49vnnvLxgAcMjIui2Ywen5s/PsV4ILh4eRtv7sk2aaOtrrb9nzmTXgAE0/PBDRjx8yODQUAoVKcK6pk21Abu1n4MxCY8e8XOLFjy9c4f+Fy8y5PZtXhwwgL2jRumVc3Jx0bsnzXeAwVb0JrD2O8BiHx+e3Llj9lzZ0YZr9mflHEmxsTy+epXyHTrgbGTyZpF3yGzwIt/YN2YMh3TG4WhmLtW4c+gQ51asoN7Ysfi3aQOAe6lSvLxgAT++9BIhH39Mh+XL9c759M4deoaE4FqsGKBulBUODji5uhLQti2hu3ahTEvTPpG8uWMHhUuV4sa2bQRNngxASnw8D44do5xOg+JVqRJVP/qIwt7eAJSqXZvGkyaxe+hQHl+7RrHKlTPc3wvvvEOJatUAeGnYMHwaNECZmsqB8eNx9/Gh2cyZ6no4O9Ns5kyubtyYqc/xpwYN9GaDT46NpcGHH2bqXOaolEr+fO89ivj50XLOHO1nWG/sWK5u3EjIxx9T4913MwStxj4HU27u3Enozp3UnzCBsk2bAuAZEEDQlCmsTx8PnlWvbd6cLefRSIqJwb10aaP7nAoX1s7Wa+54AOfChY0eD+q/qTkPTp7k/tGjtPn2W2mkhRC5TtOeq5RKEh4+VI9bnziRpp98olfOxcuLBhMmaNtXFw8P2i1dypUNGzi1YAHNZ80C1N2C940aRclatWj4wQfa42sNGsS1334zW5eDEyfy+Pp1Xtu8Wa+HWdmmTXEtXpwb27dTa/BgAJ7evat9SHBj2zbtvCN3Dx+m5EsvUahIEe3xPg0b0mL2bBwLFQKgao8eXNu0iX9mzaLO8OEZ6mHq+8jDixc5s2QJVXr0oPLrrwPqMcdtvv2W7ytWpKifn6WPO8N1FhsEeZbancyKCQ3lr+nTeaFXL6q9/TaQ/jdcsoTrW7cSMnky3bZvz1A/Y5+DKce+/JK48HDaLV2q/Sz8W7WiYqdO2lnSs2q4wYOhrMiONjw7znFuxQpUaWk0mznTcqWFXUlmXeQbbebP1z4JHf7gAcEGQdTVX38FoGLnznrbS9aqRVF/f65u3KjttqbhFRhIkTJltO/LNG6sfSpeqUsXEqOjCT9yBICHFy/i4uVF1bfe4sHx49oG+/bu3fi1aKFtkAFafvEFPvXr611LE6A/unTJ6P35tWih/d2jXDkqv/46ESdPEnfvHgFt2+p1YVIoFHrlbdH7+HG9z9GWTLktIv79l5jr1/Fv0yZD96vS9eqR8OgREUa6Mhr7HEzR9HSo0LGj3vYyTZro/T3yHBPDCGxaRs1IWWuOT01KYvfQoVR45RXq2LBUnBBCZBdNez4iMpKx8fF0++MPrm3axJp69XgSFqYt59esWYYAvpC7O+6+vkTrtKW3/vyTpNhYKr36asZrzZtHoIl2JGTyZO4eOkTwxo04Ojvr7XNwdKTCK69we88eUhISAHWvuup9+uBds6a2/QF1W1SpSxe947v/8UeGdsircmXi7t0jyUggZer7yI30DL9hO+fh768d52+LouXK6X0HsDVTbotrmzejSkujfPv2ettdPD0pFhhI2N69GSZTM/e9zJgbW7eicHDIcA3D8eB5Tibb8Ow4R0xoKH998glNZ86kdN26tl1T5DrJrIsC4/H16wC4p3cB0lXE15f7YWHEP3yIe6lS2u3m1nOt9OqroFBwY9s2yrVooW2MfRs14t8FC7ixbRu1Bg/mupFG+ml4OMdmzyZs/371E2uFAmVyMkCGcVTm6hIbGmrynox1f8pLYtL/HpfXr9d+2dBQJifj7O6eoRsjmP+bGDL1+SgUCgrr/J3zEhdPT1Lj443uS3n2jEKenhaPB3WPDmPHA2bP8b+BA1GlptJlwwZZY10IYXdOrq74t25NmwUL2NypEwfGjyc4/eG7Mi2NcytWcGHVKp7cvo0yLQ1QZ4ITdNaSNtf+G+vJBvD3p59yYdUqChUtSuKjR3oBokalLl24tHYtYfv2UalzZ25s3UrjyZNxcHLin88+0/aUu7FtG2/s3at37P3jxzn+5ZdE/vuveliWQkFK+rw1iY8fa/8v1zDV9j2x8D3gqYVu1/ak+R6wb9QoDowfr7cvNSEBhaMjCdHRVn8vMyY2NBTXEiUyPGzJq9+RstqGZ/UcSbGx/Pbqq1R9800a5UCvSpH9JFgX+ZZ/q1Z648nRmbTFama6Vrn7+OBTvz43t22j1VdfcWPrVtp+9x0lXniBQh4e3Ni2jZqDBnFr1y5afPGF9rjkp09ZFxQECgWvbd5M6bp1UaRP1PZL69Ymr2esm5dmIhpVZu7NSq9YMfu4LsPeCaYLquv80rBhtPrqK6vPb667m6lrZOpvb6Ut3bppxyZa0mb+fKpZ6KlQrHJl7eQwulmX1KQkEh8/pkxQkNnjPStWROHgwLP79zPs0zz8KKYzQZOu/ePGEXHiBG8dOiQTyggh8pQyjRsDcEdnUtMD48dzav582syfT60hQ3BydQXQm/gVyFQbEHP9Oj3+/JON7duze9gwXt+yJUOZCh074uDszM1t2yjXsiUPz52jTFCQNli/sW2besyvuzueAQHa4+4fP876pk3xb92atw4c0HbPPjJ9On8b9BTQMtH25cb3gCE2jn239XvAK2vWUMmg16NJtnwHMLhOTjEcNmBO7+PH8TDTUyErbbiG5gGUredIiY9nc+fOlKhRg3ZLlpi9hsg7JFgXBUaxypW59b//8ez+fbwNZtqOu38fF09P7Rhya1Xq0oUjU6cSfuQIcffvaycjK9+hAze3b+fO/v14Vqigd97be/fy9M4dWnzxBT716mXpnrzSMwfxRjLQxrLS2Ukzdk9l0EUt3sqxbcWqVAEgLjw8w760lBTuHjyIb+PGemP8bOWZ/vk8e/BA72+uUqmIj4zEMxtmec3uMet+LVpw7++/eXj+vF73M80M9pYm03F2c8OnQQMeXbyISqnUe7ihmdne2DlCpk7l+u+/8/bhw9oshmaVBGMZJSGEyFXpPX10s4UXVq2i+AsvUHf0aLOHmgtekmJjSUtJydD+d1ixgmKBgdQdPZqT8+Zxaf16XujZU6+Mi6cnfs2bc2P7dvzbtsU/fUiab6NGFC5dmutbt5KWlJShd92ltWtRpqTQ9NNPbR5TbsgzD3wPUKam6o3nT3j4UFsvc8x9D3gWGUn0pUt68/1khmeFCjy+ejXDA/Ds/Gyyc8x6ZttwXb6NG+NYqBAPz53LsC/yzBkcCxVSL8+qIzUpid+7dqWQpyed167VDk+MvnKFImXKyAP8PEzGrIsCQ7OW9M0//tDbHnX+PE/DwqjcvbttWVvQNsB7RozQGwtfqUsXUhMS2D9uHBUNGmlHzYRpBl2MM9NVrXS9erj7+nJr925t9z9QB6N3Dx2y+Xzm/PLyy1xYvVr7XhPAPdVpZNNSUog8dcqq85WqXVv7ACU1MVFv383t29n+9ttZHleu+fuE7typt/3eX3+Rlj7sQJezuzvw3wOIGzt2cCWTE/VlVo1+/VA4OHDt99/1tl//7TcUDg7UePddve0Jjx5px0tq1Bw4kKSYGL0MlOYcHgEB2gkWNY7PmcP5lSt5Y88eipYtq91+Ztky7az8QghhTw/SlzzTfYjp6OKSYbhOWnIy8RERetvKt2+Pi6dnhrYA4LfgYE5+802G7Zrgs9msWXhVqsS+0aONPoyu1KULceHhHJ05U9vmKBQKKnbuzL0jR7i4dm2GYD07vwdoxuHfNLi3J2Fh2qFg2eVrJydib9/Wvnc38j0g8swZUg3aJFMqd+uGwtGR6wbtHcA/s2Zx4uuvs1hj9d9HpVRy688/9bYbto8amu8BmrHyJ+fP11tGLjfY2obH3bun15vBxcODKj16ELZvH0lPnmi3J6afs0qPHnrrrytTU9n+9tsoU1N5bdMmvSEDf773Hg9OnszmOxTZSYJ1UWD4NW9OzUGDOPv999q1weOjotg3ahRF/fxo/tlnNp+zVO3aFC1XjqgzZ/Qa44qdOqFwdMywHdRLm7j7+nJ60SIeXb4MqJfQOJqJGTcdnJxo+dVXxEdEEDJpEmkpKaQmJhIyebLRCWqyQpmWptcY+AYF4VaiBOeWLyfpyRPSkpM5MmUKriVKWHU+hYMD7b//npRnz9gzYoR2KbUHJ06wb/Roms6cmeVgvWKnTpRv356zy5ZpJwKMvX2bQx9+iLORjH3x9FnmNU+0T82bl+vj/UpUq0aDDz7g1Lx52jrfPXyYU99+S4OJE7Uz4YN6Epilfn6sMZgA5sX+/SnbrBn7xozhaXg4KpWKk/Pmcf+ff2i7eLFeQ3xm6VIOfvABZZs148Lq1RyZPl37Y+rLjBBC5BaVSsX9Y8fYN3o0Dk5OehPKVe7enUcXL3Jm2TJUKhWpiYnsGzs2wwPgQkWL0mbBAiL+/ZcT33yDMi0NZVoaJ+bOJfryZWobmX1dw7lwYdp//z0Jjx5lWO4L/nsoHH3pEhU6dNDbrkxNJSEqCh+DCdAqv/46CgcH/po2jcTHjwH1JHjmlpAzxbtGDWoNHsyNLVu4mt7TK/7hQ3YPHaqdMT27qNLS9LqUax4UnJo/H2VaGgnR0Rz74gvcrPwe4Fm+PM1mziR01y5OL1mi/Z5xad06Lvz4I00srJ9ujQYTJ+Lu60vIpEk8vXsXUK9TfnHtWqPldb8HJMbEcPTTT3N0iIExtrThl9avZ0nZsuwZMULvHM1nz8apcGH2DB9OSkICKfHx7B46FGd3d5rPnq0tp1Iq2dmvHze2baN03br8M3u23veAnFr+T2QfhSq3/4UKYaXUpCSWBQSQmpBA8pMnFPLwwMnNjeaff07N/v2NHqNSKjk5fz5nly0jMToahYMDAe3b03zWLG1XtNBdu9jZrx8JUVEoHB1xLV6cmgMGmAzmdw8fzsU1axiRvryMxvrmzYkLD2fwzZsZjok6e5aDH3zA/WPHcC1enGKBgVR45RX2v/8+hTw8KBMURJdffmFFlSokx8aSmphI4dKlKV6tGm8bCaAurVvH3zNn8uzePdxKlqTGu++SlpzM0U8/xc3bmyrdu5scf3Rr925+69IFZUoKKqVSHSAbmVgsLTmZjitX8mK/ftpt948dY++oUURfuoRnhQo0njKFm9u3c+HHHylcujT13n+fIr6+HJw4Ue/zDP71V/zS1+++f/w4f02bxv1//sHJ1ZUiZctSf9w47TIuSU+eWP05GJMcF8ehDz/kyoYNKJyc8KxQgZZffskfffoQd+8ersWL8+r69finzxcQMnkyZ5cvR6FQUKZJE15ZvZpC6U/ac4tKpeLfBQv4d+FCEmNicPXyos7IkdQZNUovi/QsIoKf6tenRI0a9Ni1K8N9H/74Y65t2kRacjJelSrRdMaMDDPiLitfnic6mRJDQdOm0TQbvjAJIYQpptpzUE805uLpSZkmTag/fjy+Ost1psTHc2TqVK788gup8fEU9fen1pAh/PP55zy7fx/X4sV5++BBiletCqjb96MzZxJ95QoOzs741K9P89mz8a5eHYAVVasSFx5OyrNnuHl789KwYVTv1YufW7YkISoKlVJJ4dKlM8w/srJ6dYr6+fGGTvY2+dkzFnt7U61nTzquXJnhnq/99ht/f/opj69do2i5cpRp3BhHV1fOfPcdbt7e1Bk1Ct+GDa36PqJMTeWvTz7h3IoVKFNTKVKmDE2mTePk/PmEHz6MW8mStPzyS2r07Wv08w+ZMoXjX31FWlISoJP5N5CWlMTg0FA8deYEOL9qFf989hkJDx/iXbMmLy9YwG/Bwdr2NfjXX7m0di1XN20i4eFDnN3dcS5ShKHh4dqu1pfWr+fE118TGxqKs7s73jVr0mT6dO3f2tbvZYZibtxg35gx3D10iEIeHvg2bEi999/n5xYttPUZcvs2Ti4uJD99ys533+XuoUM4ubnx4oABGVYcyA3WtuG3du9ma/fuNPzwQxpPmqS37/H16xycMEE7r07ZJk1o9fXXeOkMAYy9dYvvLawa8Ob+/fjn9dnzn2MSrAshhBBCCCGEEHmMdIMXQgghhBBCCCHyGAnWhRBCCCGEEEKIPEaCdSGEEEIIIYQQIo+RYF0IIYQQQgghhMhjJFgXQgghhBBCCCHyGCd7VyAvcHFxoWTJkvauhhBCCKEVFRVFUvpSSyLrpK0XQgiR11hq6yVYB0qWLMndu3ftXQ0hhBBCy8/Pz95VKFCkrRdCCJHXWGrrpRu8EEIIIYQQQgiRx0iwLoQQQgghhBBC5DHSDV4IIUSmKZVKVCqVvauRLykUChwc5Jm5EEII+5A2POdlta2XYF0IIYTNHj9+TFRUFGlpafauSr7m6upKQECABO1CCCFyjVKp5Pbt2yQmJtq7Ks8FR0dHSpYsSbFixWw+VoJ1IYQQNnn8+DGRkZGULVsWV1dXFAqFvauUL6lUKsLDw4mMjMTHx8fe1RFCCPGciIyMxMHBgcqVK0sbnsNUKhWJiYmEh4cD2BywS7AuhBDCJlFRUZQtW5YiRYrYuyr5XunSpbl16xalS5eWL0xCCCFynEqlIiYmhvLly+PkJKFgbihSpAhly5bl3r17Ngfr0u9OCCGE1ZRKJWlpabi6utq7KgWCs7MzKpVKxgwKIYTIFZo2x9nZ2d5Vea64urqSlpaGUqm06TgJ1oUQQlhNE1RKFjh7SbAuhBAiN0h7Yx+a7022fv75IlhPTEyka9euVKlShdq1a9OxY0du3bpltOyKFSuoXLkylSpVYsiQIaSmpuZuZYUQQghh1OjRoylfvjwKhYLz58+bLGeuLd++fTvVqlUjMDCQ7t27ExcXlxtVF0IIIXJdvgjWAYYMGcKVK1c4ffo0r776KkOGDMlQJjQ0lClTphASEsL169d58OABK1assENthRBC5KYtW7bwwgsvULt2bc6dO2fv6mitWrWKHj162LsaeUaPHj0ICQkhICDAZBlzbXlcXBwDBw7k999/5/r16/j6+jJr1qzcqr4QQogcolAosu3h661bt/D29s6Wc9kqu9v9fBGsu7q60qlTJ233gcaNG3Pz5s0M5TZu3Mjrr7+unahn6NChrF+/PrerCwkJcPTof++Tk+H+/dyvhxBCPCeWLFnCjBkzOH36NDVr1rTqmOzueSU9uSxr0aIFfn5+ZsuYa8t37txJ/fr1qVatGgDDhw+3TztvhEql4tT9UzyKf6Td9jTpqR1rJIQQIiflRrufL6cA/Pbbb+nSpUuG7WFhYXpP68uXL09YWFhuVk3t++9hzBi4ehUqV4YPP4QVK+DRI3B2hmHDIDISNm1Sl583T12uc2fNjYCvr7qsEELkdcHBcONGzpy7UiXYutVskdGjR3P48GGuXLnCN998w9SpU/n4449JTU2lWLFifPfdd1SvXp0DBw4wduxYWrRowfHjxxk+fDhjxowhIiICZ2dn6tatywsvvMDatWu5efMm7dq148aNG+zdu5fJkyeTmJhISkoK48ePp3///gC0atWKpk2bcjT9Ae3OnTsZNWoU+/fvp2zZstqgUljPXFtubF94eDhKpTLDWvVz585l7ty52vc52V1epVIxZtcYFhxbwG9v/VCMTRwAACAASURBVEZXn1akzP8GL9UMPJMdqBzjSJ1nRannUJZ6ZRtQq1Y7nGrXhcBAcMgXeRMhREGVk204WNWOa8yZM4fdu3cTFRXFJ598Qs+ePQHo3bs3ly9fJjk5GX9/f1auXEmpUqUA+OGHH5g/f7520ryNGzfqnTM5OZkBAwZQpEgRvvjiCypUqJCv2v18F6x/9tlnXLt2jSVLlhjdrzvpkakB/DnegMfEqF9jY9WvERHw9CkkJakD8GPH1MG6xpQp0KSJOliPjVU33pMnw9Sp6uM+/FD9o/mColKBTO4khBCA+gHu2bNnmTBhAg0bNqR69ers37+fmjVrsnbtWt58803t+OizZ8+ycOFCvv32WwCWLVvG33//TY0aNUhLS+P48eMA7N69m7Zt2wJQt25dQkJCcHR0JDo6mrp169KxY0d8fX0BOH36NLt27cLZ2ZkFCxYQGhrKhQsXSElJoUWLFpQvXz73P5R8zlxbbu3khuPGjWPcuHHa95Yy+lkxZf8UFhxbQD3fetS4GQeNAohPekLfrq7cLOXEJa9EjpWOBqKBc9wbuBLfOIjy9SC0VW3qNnwNp5fbQY0aErwLIZ5bCoWCI0eOcPPmTRo2bEizZs0oV64c8+bN03Zrnz17NjNmzGDhwoUcOHCAWbNmcfjwYXx9fYmPjwfU68gDPH78mO7du9O+fXs++ugjAGrUqJGv2v18FazPmTOHzZs3s2fPHgoXLpxhv7+/v97Ec7dv38bf3z9DuRxvwDVdIsy96nabSElR/4A6WE9Jgago9ft//oHvvoOaNdUZ+d274a234ORJqFBB3eU+LAyqVs3eexBCCGtZ+cQ8N/zzzz/Url1b2xW+V69ejBgxgvvpQ5GqVKlCs2bNtOXbtm3Lnj17iIiIoEOHDly6dInz58+zZ88e3nzzTQAePXrEwIEDuXr1Kk5OTjx8+JALFy5oG+0+ffpol8DZv38/7777Ls7Ozjg7O9O7d29CQkJy8yPI98y15f7+/uzbt0+779atW5QtWzZDVj03Hb17lFmHZ9GwbEP2vLmDovWCIC0Nz+9X80PPnuDkhEqlIvxpOCfC/ubSpRB8F9eH06f5/fZvDKl6iOL3D9FlKrweWZx2L3WjcNc3oFUrKFTIbvclhHhO5KE2fNCgQQBUrFiRZs2acfjwYd555x3Wrl3LmjVrSEpKIiEhAR8fHwB27NhB3759te2xbnyYmJhI06ZNmTx5Mu+88452e35r9/PN49u5c+eyfv16du/ejZeXl9Ey3bt357fffiMiIgKVSsWSJUt4++23c7mmQFqa5VfN74bvDQN7TRCveX/pEjx+DKGh6vezZqkD+SdP1O+vX9cfLy+EEM8RlUplNPOq2VakSBG97ZpGe+/evbRt25a2bduye/duDh48SJs2bQAYOnQoLVu25Ny5c5w+fZoqVaqQmJioPYfuOWVJnKwz15Z37NiR48ePc/nyZQAWL15sn3Zex/lIda+NqS2mUnTR9+p2eMoU6NMHnNQ5EYVCgZ+HH11ffIP/e2O+et/XXxO0eCuTGk6gXIkK/FgburaPxrv4cqbO7gAlS8KAAXDgANi4Lq8QQhQECoWCkJAQFi5cyM6dOzl37hxz587Va4NNcXFxoWnTpmzbtk1vbHl+a/fzRbB+9+5dxo8fT0xMDK1bt6Z27do0atQIUD+B2Zr+RKhixYp88sknNG3alEqVKlGqVCkGDhyY+xW2NbOu+95SoG94rkeP1AH90/RJbN5/H7p2/e/cf/0FZ89m/Z6EECIfCAoK4vTp01y6dAmAn3/+GT8/P+1TeEONGjXi8uXL7Nmzh+bNm9O2bVvmz59PuXLlKFGiBKDuRhcQEIBCoeDQoUOcOXPG5PVffvll1qxZQ2pqKgkJCaxbty77bzIfGzFiBH5+fty9e5e2bdsSGBgIWN+WFy1alOXLl9O1a1cCAwMJDw/n448/ttv9ANyOuQ1AQJIrzJypnoNm7Firjn2x1IvMfOUrTk+8yc3RN/mm/Vzq+9ajdJsu6vP88AOLJrbmYh0/9fC427dz8laEEMKuVq5cCah7TYWEhNCsWTMeP36Mh4cHxYsXJzk5maVLl2rLd+nShdWrV/PgwQMA4uPjtV3hFQoFy5Yto3Tp0nTr1o2kpCQg/7X7+aIbvJ+fn8mnFsuXL9d7P3jwYAYPHpwb1TLNMKA2l1nXPC03FYxbOpfh/qdPQXcMfp8+UK6c+sk8qIP3wEBIn5RBCCEKkpIlS7JmzRp69epFWloaXl5e/PLLLybLOzk50axZM+Li4nBzc6NGjRqkpKRox62Benzc8OHDmT17NtWrV9c+LDZmyJAhnD17lurVq+Pn50fz5s25LQGW1qJFi1i0aFGG7ba05cHBwQQHB+dI/TLjw2Yf8maNN6myYQ/Ex8Nnn4GLi83nqVCsAmOD3mds0Pvq7zyTFdw9sY9R219GpbjPyzdnMfKVz+hSpQuOI0fByy/L/DVCiAJFkw2PiopiwYIFlCtXDl9fX3766SeqVauGn58fTZo04X//+x+gXmFk8uTJtG/fHoVCQaFChfQmmFMoFMybN48pU6bQuXNntmzZgru7e75q9/NFsJ7vWAqodTPp1gT25s5lqYv9s2fqH4DERPUYuMGDQfNl6d499czz0uALIfKxA5oHkqi7Snfs2DFDmVatWnHixIkM27cajNcLDw/Xe9+uXTuuXbtm8boAhQoVYtmyZVbWWhQERQoVoWbpmnD6G/WGli2zfE7NsI2y9VpzqNRhFv6zgE2KTeytmEZAzFbGTd7KqNHVUEz8EHr1ktVjhBD5niYxO3HiRL3tTk5ObNiwQW/brFmztL/369ePfv36ZTjfw4cPtb9/+umnevvyU7ufL7rB5zvWZMOtzZRn5lW3i71u8J6UpN9lPjQUypaF9C4nQgghhLDNkbAjhMWGwb//qtvUkiWz7dwKhYJm/s34+Y0N3H4/jGktpxFfujjb2wWgCL8H/fujCqwECxeqJ5wVQghRoEiwnhOsyYZbCs4N91vKtJvKrJvL4kdE6L+eOQMvvAAmniQJIYQQ4j9pyjRa/diK0TtGwoULUKdOjl2rTNEyTG81nVvj77BiwmEIC0M1axbt2z3g8w2jeFa5vLrXXHJyjtVBCCFE7pJgPSdYk/1OS1Ovl24pw25rUJ+aqj6vZiy85lrmjtW8nj4Nly9D+nrE3LwJc+aozyeEEEIIPfee3iNVmYp/qru651rt2jl+zcLOhSnnWQ48PYkYPYCbL/rxcVsIfOchi38cSXL1qrB2rcwgL4QQBYAE6znB2nHoSqXt3d+t7TavWy6zx65aBR98oF6GBtRBuxVLJQghhBDPg7DYMAACHqcHxjmYWTfGp4gPl0ZdZlGnRahKejOiM9R4NYztM3pDvXpw+HCu1kcIIUT2kmA9J9gyDj0nJpizdA1rs/TpSxxou9SNHw/Vqpm+byGEEOI5cjs2fdm2O+lzweRCZt1QIcdCDG8wnBtjbjKz9UwelHJj8xsvwsWL0KIFvP023LmT6/USQgiRdRKs5wRbsuGZnWDOlqXeLD0AsPZcN2+q13jVdIvfuRP++MP4ZyCEEEIUcNo11i/dB09PqFDBbnVxL+TOpBaTuDLqKl/83z64cIGk1zrz2b0NPHuxinpJORnPLoQQ+YoE6znBlmx4VpdwM1dOqVT/ZLaLvaXXDz8Eg+UVhBCiIFi1ahU9evSwWO7333/n2LFj2vcnTpygV69eOVk1kYckpSXh7uyO/8lr6qx6HlgGtUzRMpR0LwmBgfw09XUmvQw1h6Tx58pJ6q7xR4/au4pCCJGjClIbLsF6TrAlC56TS7jZ2oXe1mskJek/pe/UCebPRwghnheGDX39+vVZu3atHWskctP0VtN5+sa/lIp4Zpcu8Jb0q92PeR3mEVmsEB36QJ+qF4lqGwQjR8KTJ/aunhBC2FV+aMOd7F2BAslccKxSGV9mLbuWbjO3PFxmg3dz19Z0iVep1N3iCxWCMWPU2x48gFKlwEGeCQlR0NVYXMPo9i1vbyGweCDXo6/z2s+vGS1zYfgFAP53/X+M+3Oc0X3W6N27N5cvXyY5ORl/f39WrlzJxYsXGTt2LE2aNOHIkSOkpqby448/Ur9+fVJTU+ncuTOPHj0iISGB2rVr8/3331O4cGG983bu3JnevXvTs2dPdT3/9z+mTp3KtGnT2Lp1K3v27GH58uWMHDmSwMBAJkyYwIkTJwDYsWMH06dPJzk5GYVCwdKlS2nUqJHV9yTyPoVmudMXX7RvRYxwdHBkTOMxdK3WleF/DOcn/uDPKs7sXLWIutu3w48/QsuW9q6mEMLOcqoN191vibThxkkUlRPMZdR1l1KxJrOe2fXYszKxnC3X1p3ZXveYqCjw85NMuxAi18ybN48TJ05w9uxZmjVrxowZMwC4cOECAwYM4MyZM4waNYpJkyYB4OjoyLp16zhx4gTnz5/Hw8ODxYsXZzjv2LFjWbRokfb9woULGTlyJJ06dSI4OJiPPvqI06dPM2jQIL3jrl69ysCBA1m7di1nzpzh+PHjVJNJOgsMlUrF54c/Z094+ozr3t72rZAZAV4BbO+5nbXd1hJQoTZVJn2tbqdbt4Zx4yAhwd5VFEI856QNN04y6znBXMZa87vhe1u7oNsSSGc2w27NtTWZdcMyDx+qyz148N/2P/+El18GFxeEEAWLpSfngcUDLZbpENiBC4HWZ9INrV27ljVr1pCUlERCQgI+Pj706NGDqlWrUr9+fQCCgoKYM2cOoA62vvnmG3bs2EFqaiqxsbG0aNEiw3nbtWvH2LFjOXPmDB4eHpw4cYKNGzdarM/u3bvp1KkTVapUAcDZ2RlPT89M35/IWx4nPubjfR8z0CWItgBeXvauklkKhYJ3ar5Dzxd7olAooH0XvpjelqDN39Bi1y712uy5vPScECJvkDY8o7zShktmPSdYs4ya4fvsnmDO2Lkzm6W3ZjZ7S/U/cAA6d4atWxFCiOwWEhLCwoUL2blzJ+fOnWPu3LkkJiYC4Orqqi3n6OhIavr/U+vWrePgwYMcOnSIc+fOMWHCBO0xhkaPHs2iRYtYvHgxAwYMwEUeOj737sSql0Mrl5z+7yuPB+saivRJ8KLKeDHzxWha91cww+cKaUGNYMGC/x7CCyFELpE23DQJ1nOCNbO/a/Zndhk1W9Zwz+r4d3N1sDZbr5nI5mn6WrQPH8Ly5frDAoQQIpMeP36Mh4cHxYsXJzk5maVLl1p1TIkSJShatChPnz5l1apVJsv26dOHnTt38uOPPzJ06FDtdg8PD2JjY40e06FDB3bu3MnVq1cBSElJMVlW5D8Jqequ40US0tuxfBKsa5R0L8mJwSd4sXRNprVU0ra/I/cmjYZu3SA62t7VE0I8R6QNN02C9ZxgTSZasz+zWW5rymd1/HtmsvfW9hRYswYGD4azZxFCiKx65ZVXCAwMpFq1anTo0IHaVszM3bdvX+Li4qhevTrdunWjefPmJssWLlyYrl270rx5c8qVK6fd3qdPH9atW0ft2rVZvny53jGBgYGsWLGCnj17UqtWLRo2bMiVK1cyf5MiT0lVqts3p4T0TE4+C9YBqnpX5Z9B/zCiwQgO+CTy0hgX9p35XT2zvc4MyUIIkZOkDTdNxqznBHNBsKlu8CqVOsucExPMKZXq8+fEBHOa7nLWdtvXvGoms9G8Xrmi/mLQpw9CCGErJycnNmzYoLdt1qxZANpZXQFefPFFbt26BYCnpyd79uwxer5+/frRr18/7fu0tDQOHTrEwoUL9co1aNCACxf0x+jpXq9Tp0506tTJ5vsReZ82WI9PUq+v7uFh5xpljquTKws7LaRNhTYM2zEM9wkDYfw30Lw5LF4MAwfau4pCiAJO2nDTJLOeE8wFv6YmmNO8z66u69ac29Zl48w9ELA1s264/euvoW9fePYMIYTIS7Zu3UrFihVp0qSJ2Sf34vni5epFx8COBMSo1IF6Pl+mtNsL3QgdE0qj4Z/BkSNcqFacuOGDYNgwSE62d/WEECJT8nsbLpn1nGAuM20qs64pn9WZ2TVjwI1dy9b11a2ti7lzWftAQDMhRHIyuLvDyZMQGwtt2iCEEPYUHBxMcHCwvash8pjaPrXZ2WsnfFM/X3aBN6aws3p94ofV/GnbW4l3lDu/rVhC4NmzsHEj+PrauYZCCGGb/N6G5+/HwHmVuWDYXGbdXEBtyxh1zXtTDwKsDaStmYE+M5PSWVNu3DgYMAAhhBAiT4uJKTDBukYJtxKMCXqfC+7xNBjlws7Iv6BePfjrL3tXTQghnisSrOeEzCzdZvje1gnmzGXSDd9nNpDW3a9U/jcWXqm0vd6WxrjHx/83nh3gzBkIC0MIYV+aZZ9UsrxTttB8jprPVeQfJ+6doP+W/vxTKAqKFbN3dbKVQqHgo2YfsbPXThRuhencW8G3lR5C69bq9diFEPmWtOP2kdn2XrrB5wTdIFelsm6CObBtBndjQbBhcG4ps57ZBwPmeghkddI63Vfda3TuDHXryjrtQtiZg4MDrq6uhIeHU7p0aZydne1dpXxLpVLx6NEjnJ2dccjn452fR6GPQ1l1ehWvOCpoVMAy6xodAjtwfPBxOq/rzJi2V6igLEGX3r0hNBQmTVJPrCeEyFccHBxwdnbm0aNHlChRQh4W54KUlBQiIiJwdXW1ub2XYD0nGGagdbdb6gaflQDaVOBveO7s6Lpuqsu9rePhreliD+p12jVrtYN6nfZixcDRESFE7goICCAyMpJbt27Jk/kscnZ2xt/f397VEJmgnQ0+TVXgusHrqlS8En8N/IvFxxfTeeC70CUYpkyBmzdh6VKQB3ZC5Dv+/v6EhYURHR1t76o8FxQKBV5eXpQqVcrmYyVYzwmmAuicnmDOUmY9OyeYs9R939Yu97YMGXj2DPz94cMPYdo0hBC5y8HBAR8fH0qXLo1KpZKAPZMUCoVk1PMxbbCupEAH6wDF3YozucVkAJQHDzBtQj1GrP8Bn7Aw9cRzBfz+hShoChUqRGBgIEqlUtrwHKZQKLQ/mSHBek4wFUDn9ARzpgJ/TbnsnGAuO+tt6tqmHj7ExqrHs9+/jxDCfrLS+AiR3z1PwbquPVH/MLPsDda8X4Q/luylerNmsGMHBATYu2pCCBvJA+O8T/5COSE7MuuWJmDLzgnmrM2Gm7qWNePhrQ3SzV3D1Gdx/z7UqAFHjiCEEELkhuc1WG9fqT0rg1cS7pxIsxFuHI29AEFBcOGCvasmhBAFjgTrOcGasd3G9mcmQ20qE20u057ZzLqp8eTG6m3rGHVz3fo1M8+b+iwuX4aLF9VrswshhMjTrl27RpMmTahSpQoNGzbk4sWLGcqsXr2a2rVra3+8vb3p1q0bALdu3cLJyUlv/40bN3L7NqhZuibjS3WlfAzPVbAO0L9Of3a8s4MkZwVtB7uwp0gktGgBx4/bu2pCCFGg5ItgffTo0ZQvXx6FQsH58+eNljlw4ACFCxfWa7wTdJf+yk3WBLXG9mcmQ21LZt2ajLm1r6buIzP1NvaalvbfsnC2Zufv34ePP1Yv/yaEECJPee+99xgyZAhXr15l4sSJDBw4MEOZvn37cvr0ae2Pr68vvXr10u738vLS21+pUqXcvAUAGvs1Zo7761R5xHMXrIM6w76nzx6cXdwY0L84Scnx0KYN7N9v76oJIUSBkS+C9R49ehASEkKAhfFQ1atX12u83dzccqmGBkwFtZa6weu+NwxUTQWz5sbHm+tGrlLZNkbd8NVUvW3tEWBNJt7WLvTbtsHnn8PffyOEECLviIyM5NSpU/Tu3RuA7t27Exoayq1bt0wec+zYMSIiIggODs6lWtogJkb9WsDWWbdWULkgDvY7yNb+f+Kyc7d6hZZXXoEtW+xdNSGEKBDyRbDeokUL/Pz87F0N65maNd2WCeY0+80Fs4bLwlmbWQfzDwKMXVM3uDdWz8xm1K05zlJAb/iakqL/mpQEdugiKYQQQt+dO3coU6YMTk7q+W0VCoV2CSFTVqxYQZ8+fXDWWSLsyZMnNGjQgLp16zJjxgzSdNukdHPnzsXPz0/7ExcXl633svrMaho9/pJL3jyXmXWNWqVrUdunNjRrxtHN3/JtMyfo3h3WrLF31YQQIt/LF8G6ta5cuULdunVp0KABixcvNlkupxtwveDS2ky6qfLmssi2BOfmZle3ZhI73QcD5u4jszPMmwvmrcn0mzvnF19AtWr/ZUCEEELYjeEKAuaWDYqPj2fDhg16XeV9fX25e/cux48fZ8+ePRw+fJivv/46w7Hjxo3j7t272p8iRYpk300Ad5/c5RjhJDrxXAfrGiqVig9uf8+Y5s/4pHMR6NsXzHwXE0IIYVmBCdbr1q3L3bt3OXXqFL/99htLlizhl19+MVo2pxtwqzLQxt7bMs7clgnlDM9leLw1XdatnSgvs5l1c+Vtzawblo+KUm+LjUUIIYT9lCtXjrt375Ka/v+0SqXizp07+Pv7Gy2/ceNGXnjhBapXr67d5uLiQqlSpQAoXrw4AwYM4PDhwzlfeQPP62zwpigUCn5/63fq+dZjet1YpnUvjmrECFi0yN5VE0KIfKvABOseHh54enoC4OfnR8+ePe3SeOuNI89qZt3SDO6Wzm3teHlrgmBz57bUXd/aaxgrr3vuzHShN/a6bBkMG4YQQojcVapUKerUqcNPP/0EwKZNmyhfvjzly5c3Wn7lypUZJqCLjIwkJX2YU1JSEps3b6ZOnTo5Wm9jUtLUdXBCAdn90D+fKlG4BHv67qF+mfrMqBmtDthHjoSFC+1dNSGEyJcKTLB+//59lOldtZ8+fcr27dvt0nhn6C6e2QnmDI+3JrNu7QRzhvstBbvWZNY17zUT49maSTcXaNs6wZyla27cCD/+iBBCiNy3dOlSli5dSpUqVZg9ezYrVqwAYNCgQWzdulVb7saNG5w8eZK33npL7/iQkBDq1KnDSy+9RN26dfHx8WHSpEm5eg+gk1l3LwoOBebrVJZ5uXqxu89uGpRpwKc1o/mjdVkYNQoWLLB31YQQIt9xsncFrDFixAi2bNnCgwcPaNu2LUWKFOH69esMGjSI4OBggoOD2bRpE9999x1OTk6kpqbyxhtv0L9//9yvrC3jxm0JuM11Dzf13tKDAlsmf7OUWTf13taMurFrm3uIYOy9pWsb1vfiRXB1hYoVEUIIkbOqVq3K30ZW61i+fLne+0qVKvH06dMM5bp166Zdc92eNMG6s7uHnWuS93i5evFnnz9Z+e9KOg1+G15+GUaPVu8cNcq+lRNCiHwkXwTrixYtYpGRMU+6DfvIkSMZOXJkblbLOHPBeXZk1k2NIzfXNd1YeXNd7I1d01x3fXP3YW1GXbecUvlfDwVbxtZbm1k3rO8bb4C3Nxw8iBBCCGGNN2u8yQtf/YC3a3F7VyVP8nL1YlzQOABS9+5mc79GvDF6NAqV6r/AXQghhFn5IljPV7KSWbd1LHhmA2jDa9mS3bb1vqzNrJu7li1d5k2d03C7Uqnusq9QwJMnUKgQQgghhLUalG1AgxNpUNfb3lXJ8z67tpJpTe9xwdObT8aMAScnGD7c3tUSQog8TwZZZbesBNDmylvTNd3WbvCZzazbcm5z2W1Tr5bG1tuSQTe239xnqVLBW2/B+vUIIYQQJimV6lVGZCZ4i8Y0GqOedO7Fh8wJ9oYRI2D1antXSwgh8jwJ1rObpbHd2TWuPCcmmMtMAG2q3tZkt029mvoMbM2w2xK8a7YlJMAvv8D27QghhBCmjNjyHkU/UhFTzM3eVcnzPF092dVrFzVK1uCDug9Z2q4Y9O8PmzbZu2pCCJGnSbCe3awNvjVlbRlXrum6bercmX0QYCzjrJnV3ZprmXtwYOv48uzIrNuaeTfWbV+3DufOQXIyQgghhEZ8fAxxLuDkWczeVckXShQuwe4+u6lUrBLDmsbwc1BR6NkTdu2yd9WEECLPkmA9u5kLYq0JsM0F3LoBY05k1k2NG7d0Llu6wVsTOJv6DDQPK7IatJvLrBu+hoZCrVrqtdmFEEKIdCmJ8QA4e0iwbi3for7s7buXhmUbUuvLH6FoUejWDQ4dsnfVhBAiT5JgPbtlZYI5SwF3UpLpc2flQYCxtdHNTXxn7L5M3WdWAmdrz20u8De339j9Gr5GR+u/CiGEEEBqqvoBulPhInauSf4S4BXA3wP/pnqT1+B//yPRxRFefRWOH7d31YQQIs+RYD27mQs0s5pZ1w3WLZ3bmky7uSDY2rLG6mkus25ukjdT5zZVl+zqDm9LV/uICBg2DKKiEEII8fxKTUsBwKGQi51rkv8oFAoAQkolUnGiC/+UTIaOHeHSJTvXTAgh8hYJ1rObuSDWUibdUlbZMLNuy7lsyZbb+lDBXHCfmbHr1tbb1sy6uUDc2qz8vn2wZImsyS6EEM+51LQUnNJAIUt/ZpqTgxOxqgQ6D3DhilOsOmAPD7d3tYQQIs+QYD27ZSWzbingtrUbvC3j5Q0DbFu655sLqC0F6Wlp+pPZWcrqmxtfbm3GXbcOhkMALAX2KSn673WPEUII8dz4yq8/R1YCEqxnWmO/xvzS4xdi0p7RYXQx7j0OUwfsMTH2rpoQQuQJEqxnt6wEzJYC7sTE/37PzgnmDK9la73MBdS2jBc39T6zWXtrusPb0j3f2Ou6deDtDY8eIYQQ4vlR2bEUDcMBZ2d7VyVf61ylM8uDl3M79SGdJvgSe+08vPaa/nceIYR4Tkmwnt3y4wRzhteyZoy6tT0ETGWudYNhWz4za7qs2xK8W5uNN7X/0iV1BuDePYQQQjw/7sdHEl4UCdazQb/a/fiszWecUd1n9ZiW6tnhe/eWnmtCiOeeBOvZLStBrqWA2twEc7Zm0i1l1jMbQJsqby4YtuUzsyazbkvm3dK4d0uZdcNyd+/C4cMIIYQo2HqGXiVo2QAAIABJREFUzaXOUKQbfDb5qNlH7Oy1k5Gf7YW33oJNm2DMGPVDfyGEeE5JsJ7dsppZNxdQm5tgzpou99aWt2XcuLXlzQW71mbpDfdbG0hnJrNu6UGAqfKTJ0O7dpINEEKIAi5VmYqTEsmsZxOFQkHHwI4oHB1JWbmc9W9Vh0WL4PPP7V01IYSwGwnWs5stQa6lgDolxfrMujUPBizVTTPJm6XA3lwwbuq9uSDX2vHvhuWtDaTNBfeWAnpbu9Y/far+O+nWWQghRIGTIsF6jvnk6GzeeeEi33YrC5MmqeeHEUKI55AE69ktq0Gu7v7kZP1zZ/fSbeay9rY8VLClS76xANqWcxv2ANA9l7UZdt1ra7YpleZnhbfU/d1UuRUrYPFihBBCFCypqvRgXbrBZ7vxQeOpWqIqY2vdY2tQcejfH0JC7F0tIYTIdRKsZ7esBrmmMumQcTZ4zTJixq5lTVbflvHwtmS/LXWxN3y1tgu9YV2s7bJuLjtu6tq2ZtZNlZs/H775BiGEEAVLqipNMus5pJhbMf7o9Qfehb3p2Smek+UcoWtXuHbN3lUTQohcJcF6drN2vLU17w2DdcP3upl3a4Jcc0GxpfHw1gbQlq5l7NXWsfeWst/WBtrm7sPWBwLmMvC69ZeJcoQQokCQYD1nVSxWkS1vbyENFV0GuHEnLRo6d5alUoUQzxUJ1rObbvCnVOpnv40FteaCXEvBelYCbFsz67Y8ZDB37bQ0dcBqbWbd3AMPazPputdSKv8LmM3dV2ay9Ka269a/SROYOhUhhBD520nH4RxfhnSDz0FB5YL4seuPuBXx4tknk9WZ9ddfz/h9SAghCigne1egwNEN/kA/+21LUAu2Beu2dl03fK/bxd5YvbMywZxhZl0zkZ2pY819Rob3pVTaFkBbOre13dwtdX83lVk/cQJKlEAIIUT+5poKpCCZ9Rz21otv8Vq113B1coVbceqhZQMHwpo1oFDYu3pCCJGjJLOe3TRBmkP6R6sJghWKvJdZt/ZaWe1ib5hZt+WBhS1ZfWsCaEtd7rMy/t3UfnOf+4ULkJCAEEKI/GVv6lX+9kMy67nA1ckVgD+HtuPDYYGwdi3MmGHnWgkhRM6TYD27aQIxFxf1qybodXHRD9R032ueymsCO82TYsNst6Vg3VJ38sxm8TMTQJvLrGc20NfstzThnLXjx20ZW5/ZzLruOTU9CjT7oqKgZk2YOxchhBD5Sx/F70xoj2TWc4lKpWLusfl8Wfo633UrB9Onw6+/2rtaQgiRoyRYz26mgnVXV9PBuqasJrBzddU/VvNFwNLs8LZk0m3J4lsKzi31EDAcC57ZLvSmjreU9c5MZt1SMG5txt1ctj4mRj1+PjoaIYQQ+UsqSplgLhcpFArWdV9HpWKVGPXSPfbXLQbvvgv//mvvqgkhRI6RYD27aQIzw4DbxUU/WCxU6L/3mmDdMHjXPRb+C84N9zs4WBfkWjvBnO61dOule25bxpkbBsG2ZtZtzYabyoJnJbOemYy6tfXTvCYkwJw5ErwLIUQ+oA3WpRt8rinuVpytPbdS2LkwPbqnccMjDV57DSIi7F01IYTIERKsZzdrusE7OoKTU8bgXBPUmgrWzb3XDQodHS0H7ykp+hO9ZWU8vC1j8VUq2ybdy0w23FzQbmuW3tbg3FiQb21X+kOH4IMPYNs2hBBC5G2SWbeP6iWr83OPn3mc8pTg0SWJf3BHZogXQhRYMht8dtMEXsay4TEx6v1OTuofY5l1pdJycO7qCk+e6O9PTta/tu6DAU2Qqhvk6gbMuuc29t7WLvaWurobnttcoG8pG27Yxd7Sa2bHv2c2s26sG7ypYzWfi3zhEEKIPC8VJc4SrNtFp8qdmNN+DgoUuMXcg6/mwNChsHKlzBAvhChQ8kVmffTo0ZQvXx6FQsH58+dNlluxYgWVK1emUqVKDBkyhFTdoC63WJtZ12S/zXV7z2xmXfe94Xh4zSz15rq9G+63lKW3ZlI43f3mJsazFOiby7wbC5jNreluTWbdUoBtab8tXesNt//7L3TuDLGxCCFEQXHt2jWaNGlClSpVaNiwIRcvXsxQ5sCBAxQuXJjatWtrfxJ0Vs7Yvn071apVIzAwkO7duxMXF5ebtwCAX2phSsch3eDtZFzQON4Peh/F57NRdu4Eq1bJhK1CiAInXwTrPXr0ICQkhICAAJNlQkNDmTJlCiEhIVy/fp0HDx6wYsWKXKxlOsPMum42XBNYGmbWjY1vt/W9bqCpeW9sPLzmWrpLyumeS8OWCeYsBeeG5S1l1i2NKzcVvBsLeg3XdLc1S5/ZIN3w1dh68Ja6xe/bB3/8AefOIYQQBcV7773HkCFDuHr1KhMnTmTgwIFGy1WvXp3Tp09rf9zc3ACIi4tj4MCB/P7771y/fh1fX19mzZqVm7cAwNXrr7BiK5JZt7NHSTG06xbH2o5lYeJE2LXL3lUSQohsky+C9RYtWuDn52e2zMaNG3n99dcpXbo0CoWCoUOHsn79+lyqoQ4zmfVlNZP5tlQoq2ukss0/gRDPWK66J5DkZjDbu2FAbSmY1zwIMJdZNwzeTZ1Lw9KycLZ0g09Jsf7ctnaxNzeBmy3j3Q3LWwqkMzOGPatZe8PfhRAin4mMjOTUqVP07t0bgO7duxMaGsqtW7esPsfOnTupX78+1apVA2D48OH2ae9TUtQPvB0dc//aQitFmcKVmBsMbvqI05Xc4Z134OZNe1dLCCGyRb4I1q0RFhaml3kvX748YWFhRsvOnTsXPz8/7U+2dp/TCZgfu8IFp2jt+xktVIwJvMa7HeIJbnmP5q1DqfpmJHtLPQOFgoSUeDr0hgF1wpjVHH51v82Z0hDvlj61gKXZ4TXj0HUz65rueYbBu+GDAc17TXlz2W9Lwbvhe0td7C0F45mdYM5SJt1cZl2pNJ4NtyU4VyrV3fAN78OarvPGtn/9Nfj7Z5xvQAgh8ok7d+5QpkwZnJzU7ZpCocDf399oe33lyhXq1q1LgwYNWLx4sXa7sfY+PDwcpW5PKnK2rU9VpvJl8UvsrCaBur35FPFh05ubSPt/9s47TKrq/OOf2Z0tFOlF+oo0lW5BRbGAXbFFjb1i1xgSS2I30RCTaPJTUNSoGAz2FmOJxohgRxQVbEhdqgJSt035/XH37Ny53HvPOTPDzi77fp5nn7vv3DPnnilw93u+73kPCY4/rzlrqtY5Bec2b8730ARBELJmuyowF3EVFUkqkeTD+PHjGT9+fF2sc+2tUEKrtJSbDoL3BnzFh69DtKSEF56AzUN25ad1K1jXoQXrIpWsZjO7Ve0AhYX8EN/IzJ6wpXg1dAOYA/0BZvHjZ9C+spKlrWB62QaGLYX+VRXOB+jnlldUOCIvGk2vDm9SvK66OiWw3VXr3a/Rxv3W7Q+fTYE53ZpwnWuvKzhn66y729uk6/s97r3WV1/B8uWwcSO0b48gCEJjJOIpAOZ3vx4+fDjl5eW0bt2a8vJyjjzySDp06MDJJ5/s24cf2/JeXxWr4toeX3P2blGOyFmvQqaM6D6Ce4+4lwtfvpDTru3DK7d8TuG4cfD441JwThCERs12I9Z79uyZlka3ePFievbsWf8DcTnrb3WBtSXVFCaceI/lQOcS+KEFVHV0qsOvi8HuLSEapeeGCJvugJ+OHsP82W/y7d59+HbtfJbs0pV2Fcuhqoo3dobzB8+HwVAaf4NB3WF4ciljtsDPqqqcm1JRkSPolFhX6+P9nHXdHu4lJc7stGkquordQtUmxT6biQE/QW2zlt4kxd7GYQ+6tqlI111zxQrYtAn69kUQBKEx0KNHD8rLy4nFYkSjUZLJJEuXLt3qft2qVau637t3786pp57KjBkzOPnkk+nZsydvvfVW3flFixbRrVs3CgrqL1kwlnD+H45uPwmKjZ5xu4/jo2Uf8dCnD3HnuF35zQPTYM894Ze/zPfQBEEQMma7ucuceOKJPP/886xatYpkMsn999/Pz3/+8/ofSK3Q+rE5zOsEB69sTgTSRbC3wJxyv6uqiABtoy3ZczmcvqgVt74Njywd7vRRVcUR38ET3w3h2pkwal0rFrWByd1X8ejQVN9Teq3j7gHrmL3DJhLRHFeeLy21X2duU7xOJ6jDnPh43Ek7D3PWbdPiM1mbrutLty7eVMyPGwdHiKcjCELjoVOnTgwbNoypU6cC8Oyzz1JWVkZZWVlauxUrVtSltW/cuJGXX36ZYcOGAXD44Yfz8ccf8/XXXwMwadKker/f14n1iKTBNyTuPfJerhpxFeNuexkGDICrr4b//S/fwxIEQciYRuGsX3bZZbz44ousXLmSMWPG0LJlS+bPn88FF1zA2LFjGTt2LL179+bWW29l5MiRJBIJDj744MAKs9uUWkE1s9VPAIxaUVs8zu1ge7duKyx0BLuBYO6yCU7ZVMYpb86Bn/qQnDWLJcceyOb33oa+VVBYyMSy1XzcejOwng5VhYxZGOGQwiUcXRKjU662iYvHnW3g1NZoYaJX56zrUtNt3HD3OnOTFHvTa2dSWC5onLbp70Ht1q1zsjMEQRAaEZMnT+acc87hjjvuoFWrVkyZMgUg7Z7+7LPPct999xGNRonFYpx00kmce+65AOywww489NBDHHfcccRiMQYNGlTXR30hznrDpCRawt2H3+0Ezz9P9d57UnzKKfDJJ9CjR34HJwiCkAGNQqxPnDiRiRMnbvX4Qw89lBaPGzeOcePG1dew/KkVVO+0+BGA/ZfW3sjdhdtatkxfC66c9gwEdAToFWkDPwA9KiEa5c1PB/HBhnm81TvCf7pX8cSACp5gHs93h+NKS0kC70VXsHcECrPZJq6wMOVkh4nebJ11W8GdS2fdJu3de14n/G1T7P3Ou8f7xBPQuTMcdBCCIAgNlf79+/P+++9v9bj7nn755Zdz+eWXB/ahRH2+qEk4u5wUibPeYPm/9f9h8nVtee+WpbQ+4QSYMSNVVFcQBKGRIFPCuaZWPH1Wso4Om2GXVfGUcw6OUFXOuhJ0Ks5m6zYVR6O0ipRy6KJCJnzcmtkzd2PVw+3555y+HLzQee7sLrDfwA/p8ms4b/hSXuoPFdVb0vv2Xtu9pl0JTzXJoESjmpCwSYPXCX1bwe1NNw9Lodel83tT7FUbv6PfeV1KfVBavKmL7x3vL34Bt92GIAiCsG0pLizmqOUt2XVTs3wPRQggmUwyr2opZ13bl8Qns+Cyy1I7tAiCIDQSRKznmlpB9WbVKcx6ACJV1SlRC+lr1mtqnLTtTJ11v6Jw3hT7aJRONcWcungHWlU5bTtugRu/786Om+CRbqs59lToMPwNxh1jeG0lEtVaeyWKlbBXIldVYA3buk3nbtsWiXNva6YT/jZ9m7riYRMFpk66qYvvvUZ1tWzrJgiCUA90atGJl9/uyrglHfM9FCGAK0dcyemDTuelgu+4/dw+8PDDMHlyvoclCIJghYj1XFMrqApKm9FrPekCGtLFuhJW6nxYRXbQn9cUr1Nte66H2+Z15vP74PvPDuCu12D4hpbUFFInuJ9uuYT/7gTxEs++66WlKZGoMgaUKFbjUJMQ3nEqvFu32axJ16XFh/WdSTV43fpxb3G7INfb6/j7HTNJiw+ajIjHYexYeO45BEEQhG1AdXUqo0xocEQiER445gGG7jiUm3t+z2v77QhXXgkffpjvoQmCIBgjYj3XxGJ8tiM8XvQ1P5XiCFe3s15dvbU4z2LN+laxZfG63rTllx/AjHf78dBLzvlEBH5Z9jVjzobuoz/j2jHwTWRN6vlBznpYej6kT1i4CUuLtxXYtsXrbNa/q+e4Y/dR5+LbOug6UR+WhbB2LfzrX/DqqwiCIAi5Zen6pZy97ype6rYp30MRQmhe1JxnT36WViWtOOvIKjaXRODkk2HNmnwPTRAEwQgR67kmHufZXeCM5DOsbFn7mHvNOqSLd3eshJbJ9mlB593ryAOc9aDnRmv3gy9IwpsfDeCmt6FZMsqd+8GAoTPZ/1xY2gq9s246yaAIS4vXCWxvbLOHu8n69yCR7iekdRMBpqn0Ns56MulkMQRdwz2GjRsRBEEQsmdNxRoe61fBF62r9I2FvNK7bW8eP+Fxpp78BC3unghLlsBZZ6XunYIgCA0YEeu5JhajslaHN6upfUwJZoU7Ld4vVqLWGyshpkuDV253TY3VtnDueMDaAm59G+Z/cRBvPQqnL+/A0tawY6QlxGIsKNrEnE7J1CSDn7Oui/2c9lw66ybr4U3T4E2EtG1fNg562LX90uK9bb7/Htq0gaeeQhAEQciOuq3bChrFpjpNnqP6HcWhOx8K559P1VmnwyuvwB//mO9hCYIgaBGxnmvicSpqt1YvVVotyEkPiouKnD3MFUHiXc0KqzVzag931Zc7Bd+73t1wfXxBSSkHLYKps3oy//+gqLgZxGL8ofcyhp64mv3H/siTHVc529hks9ZeYSKgw2LbSvOmbriu2JvtJEMiYeeg+13bfT6RSG2j59dm2TKnzZIlCIIgCNmhxHqRiPVGxbwfv2Lg8A947tAecMMN8Pbb+R6SIAhCKCLWc43LWa8T6zon3S92i3fvvqBu8R6WYu922pWIVRMBJin1nvPRBM7EQDzOhQvaceb3LfmoYzU/H/Y9vc5dz627/sCqFpgXyvNe270NnMJkXfm2LDCnc7eDqr+biPdMnXS/MZi67+q4ebOzL3uNSv8QBEEQTBFnvXHSuqQ166s3cO4B65jftRR+/nNYuTLfwxIEQQhExHqu8RPrts66N20+yFn3e26Q8HennmdazM61Rn3PNSU89mFXlj5Xxh1zdySagFv6ljOvo2FfQXGYuw25ddbDCrR526u14WHOum5iIOjapk66SQq+qfB//nk49VSYORNBEATBjljcmeiMFopYb0x0a9WNf574TzbWbOZnl3WkYs0q517ovl8LgiA0IESs55p4nB7rYUhpGcVKm5kUmAtz1sPEum3fKg5yv3XF6woLHeFa7ewf3yleym/mtWfBpCJen7c7By5y2v5nZzh67wVM7wXJKost6dzCU629DxPvaps4hc2adT8xbzoxEOZu+/UVFpsK7LA17rZ9bN6cfhQEQRCM6d6iC1e/C8NqOuR7KIIlY3qP4eYDbmZO1WKu/NWuTir8zTfne1iCIAi+iFjPNbEYf/gvfDZoIhH1mEmBuWycdZu+gyYCbPZwh/Qt6OJxojVxDq3s5rzmqio+6gavdFjHgefCXv1n8NwukPDu2a5z1r3iHcKrv3tjmzR3vzhMzPuJ4rBxhrn2pgI7G2dd9/jMmXDmmeIuCIIgGNCnZU/ufAP2iXfN91CEDLhh1A0c0vsQHmo2j5eO6gN33OEUnRMEQWhgiFjPNfG4sybcNu09LPaKdfca9kyddW/fyp02SYNXsXebuOKUGL/hHfj2s/255GP4svkmTjwFhh7wFR92C+hbxW7hWVqailXBPa8oDhPryWT6mmzb9e9hafLxeKqom1/bbNbDmwrsMGddt3bd+/gLL8DUqVKAThAEwQR1bykqyu84hIwoLCjk8RMe5+YDbuaI/3sF2rZ1JqzlHigIQgNDxHquicW4Z+8If1zxdOqxbVlgzk+Mhznp3vNBfas/RPzS4CFVeb6w0HHZ3deqFcx9Iu2Z9G9Y+L/B/PpdWNysinYVTl+xAoj7Oe1uUVtbzI54PDVOr6AOE+veOBNnPSg2Efqm1eEzFdhhznpY5fiwvlWcSEBFBYIgCMLWvLXwLfYcB/9pIcXJGisdW3TklgNvoah3X3jsMWI/rYWTT5bCq4IgNChErOeaWIxHByd4YKUrnSqfBeZ058P6dseVlf6V5t3F6ryV52ufu+PGJH96A5Z/Npq+a52+Hh0Kux44l38MhlhVRaq9ErkFBY5joYSoGocSvZFIahxuwmJb9ztMgGci9MP68jvaiHdb4a8T71dcAUOGIAiCIGzNmk2rmdUN1hcn9I2FBk1VrIpz488y7jcD4cMPnS3dBEEQGggi1nNNPE5FUYTSwuLUY7YC2rbAXFjfOuddJ9a9KfeqLz+xrqk036KoRV1cGYWVJdWcdQIMGPg20wbWrmlXwlONW4la1ZdXvNs666bryNX5MDdcJ8ZNxbypsHYL7GTSPw3epoJ8WPsFC2DhQgRBEIStidU495ZooaTBN3aKCotYtWkVjxZ9ydTjd4Y774TXX8/3sARBEAAR67kn5mzdVlqYhfu9LbZuC4rDxHnYtWKx1LVUcTpdpXlXfPlHsPjz0dzyP/ihqJrTfgbDR33FvHbxVN/ufdfdzrpbvHvFuXufde95m3XkuvYmae66SvOmBeb8BLa7An6Y8M9GxLvH/7//yTZvgiAItdTUifViTUuhoVMQKeDR4x6lS8suXLL7Sr7r1dJZv75iRb6HJgiCIGI958TjW4v1TIrAbeut20z69ruW3zjdafDuOGwbOKBNtCU3T4cFbzlr2leV1LDjJqC6mkQ0VWk+TawHOevFnvXvCm8avK2zHhSbpNAHiXPveZ1oD0t79/adaaG5IPdeTQpceilcdRWCIAiCy1mPiljfHujUohNTT5jK5tgWTrmkI1Vrf4Azzki/bwuCIOQBEeu5JhajsjBJs0KXQ53J9mpBglqt5Q7qW7dnu81EQCSy9bX8xqkqomvS4IPEe/tNcf70Biz8/OC6AnS3jKzhuFHLmdt8c7qT7i04F9S3wu20b0tn3SaF3ns+E0EdJPxN16br3Hzv8ysqts5aEARBaKLEYkqsSxr89sLBOx3MDaNu4NPKhVw/fgi89RZMmJDvYQmC0MQRsZ5r4nFOXtiMQ3fcN/VYtlu3RaOprct04lu3Z7tNgTkT197vWqo6fNg2cJAS2LUisLS4ed35pa3hxW4bGXz4Qs7bexXLW9UWlKupcdxe3f7w7kJ4iky2VzN11k3S4INiXaq6ibOuYuWGZ5oOHybu1blkEg45BCZNQhAEoSlyYOuhPPo87FbULd9DEXLITQfcxMW7X8yFv3gMBg2Cm26SJWCCIOQVEeu5Jhbj/g86cs0u41KPZbt1m1sU68S37URAcXGqsjqEbwvndy2TwngGa9i98SPT2/LRO/05cFUpj/TbTL+DvuR3B0SorN5i3VcdNmnuqn2QALd16cNiW2ddNxFgkxZvKt7d462pgTffhLffRhAEoSnSr7gLZ8+BriUd8j0UIYdEC6Lcd/R99Os2GJ580jEVTj0V1q7N99AEQWiiiFjPNbFYbpz0IFFs4qzrhL+p064bp65v9zpzd2ySJl9YyJ6bWvPmKx341ytt6FZVwpQhSSJVGtc+aGKgoCAzZ90dZ7P+PUz4Q7obbiKwTdP3M02LD3PzvW0B3ngD1q9HEAShSaD24i6SNPjtlacTX3LMTX2JLS+Hc89NLfkTBEGoR0Ss55hNVHPk6BVM+v6J1IPZFphzPz8TZz3TNey6tfS6cQdVlg8S6+4icbV9R+IJjv6+kC8+3oPXnimlpNL5A+nO/j/yQXeCxbnf+njdunJdKns2leVtXXv3UZcGH3Zt0zXpJuI96DkLFsChh8LkyQiCIDQF/rLonzS/Hj6JSMXw7ZV3l77Ly5VzuPWygfDSS3DPPfkekiAITRAR6zlmU6SGV7tu4asNC1IPZltgzv18E7fbxln3xqUhhfF0zrpuPbzq289pdxfOq6xM37otHqe4sJg+G4ugspLlO8ANfZawzwVwWs+PKW+FmdOeybryIGcdUs4KmIlx07R43dp1Exffdk26yZr2IEG/YYNz3LgRQRCEpkBlTQUVRVAgBea2W/445o8M7zKc29vP5a2RXeHqq2H27HwPSxCEJoaI9RxTgSNgSotcojeTQm1BothmfbuKM9133Vbo21Sad8e1ae91fbmc9ToR66o033UjzPliP478FqZ1XMmAy+HOAWuoLiR8DbtbcBYX260rh62roYdVmvdLkzedKPAT0olEKgUvF8LfNOXer8Cc7vjdd7DnnvDNNwiCIGyPxOLOZG1RUammpdBYKYmW8MSJT9CiuAWnH13F6uZJOOUUmZgWBKFeEbGeYypxbuDNipqnHjRx1k0Ludmsbw9qb7LOPGjcNmvYvWv5gsR6PJ7elxLryllXdQAKC+vE9y7J9vz7n/DvWf3ZcRNc238J5x1LuFh3C2i/tHhdrNvDXbceXifeTdPesxX+ftewKTCnO376KcyaJQ6EIAjbLbGYUz8lWlSiaSk0Zvq278vkoyezsmoNZ/16ZxLfz4dLLpH164Ig1BuNQqx/99137LvvvvTr14+99tqLefPmbdXm7bffpnnz5gwdOrTup6Kiot7HWokjbkqLmqUezEXBObezbiOYTQrMBbnh2RaY0znt3pR7t1hXExZKYKvzHjF+5LLmfDkJfrekD798H6isJAn82LIg1Zdqn0ymbyunBKjaFs+mwJw3znY9vJ8YNk17N3HWTYV2mBMflpYfdly7Fp55Rv64EQQBMLunv/XWW4wYMYJdd92VgQMHcv3115Os/T9k0aJFRKPRtPv9999/X6+vQTnrIta3f04bdBqX7XkZJx/xayLHHAOPPw6PPZbvYQmC0ERoFGL9oosu4sILL+Tbb7/lmmuu4fzzz/dtt+uuu/LZZ5/V/TRr1sy33bakotZZTxPruiJwJue3ZYG5IGfdxvE3ubZuT3fV1i3OlYBW532c89IY3PDjruy+won/1R/KDprDn/aF6uqKrdrXxapvNWlgU2DOG9uuhw9rn4mzHhSbCusw8e7du912AuChh+CkkyQtXhAEwOye3rZtW6ZNm8a8efOYNWsW06dPZ9q0aXXn27Rpk3a/33nnnevzJYhYb2Lce+S9nDf8fCIPPwLdusHll8P8+fkeliAITYAGL9ZXr17N7NmzOeOMMwA48cQTWbhwIYsWLcrvwALouR7u/K6MUT33Tz2YiTh3ny8o2LYF5tR5d5E3k2vZXjuXLJmgAAAgAElEQVRMrAf1VVPjiMUAZ90vblENnWtKuOZQGLrrO8zoGdBeOdbuLebCBLVXrOvWrJsKam/7TJx10zT4bBx1Xep8ULx5c/oRHPEvCEKTw/SePmzYMHr37g1AaWkpQ4cOZcGCBd7u8sYvmh3ERw9Alxad8z0UoR5ZEt3MuTcMZHP1Jjj99PRCs4IgCNuABi/Wly5dSteuXYnWCrlIJELPnj1ZsmTJVm2/+eYbhg8fzp577smkSZMC+7zrrrvo3r173c+mTZtyNt7uG+Dq8l7s0X2v1IPZbN1WWAiRSHiBuWycddvidbbXNhXrQX0pwed11r1bs7ni0Qth7rwDufV/sKBkC6POg4uGlfNTKf7OutoyTglmFauJAoXOWbdNVbcR2JkKf52Tnsla9UxFuzrecQf07y+CXRCaIDb3dMXKlSt55plnOPLII+se27BhA3vuuSfDhw/ntttuI+7+/7SWbXqvT7Rkz+VQUtIiZ30KDZ/nvnqOR1e9zvirh8BHH8Gtt+Z7SIIgbOc0eLEOzs3cTdJn7evw4cMpLy9n9uzZPP/889x///089dRTvv2NHz+e8vLyup+WLVvmbrDuYmgKkwJzQQ61W0j79WUj/IPOmxavy6TyfNC2cKBfH+/3OsK2ZoO686XFzblpOnzxzkAOWgj/6PIDa5rhL9ajUSerQAlmrwuv0K1ZtykwZyKw3XGYa+8n/E2ddJPHM3XWg45ffOGkD3rfT0EQmgQm93TFhg0bOOaYY7jmmmsYPnw4AF26dKG8vJyPP/6YN998kxkzZvCXv/xlq+duy3v9D1XrWNQG4tFG8WeUkCOuHHElY3qP4YGiOTx/VG9n8vmdd/I9LEEQtmMa/F2mR48elJeXE6v9Qz+ZTLJ06VJ69uyZ1q5Vq1a0bt0agO7du3PqqacyY8aMeh/vcztVUrbnu7yx5H+pB7MpKOc9KqddFUUzKTBnuh4+WzGue74uDd7E1VcpZyFp8O6477oI/50CH83dl53XAZWVvNEblrRm6+J1SgR7JwIKPMXqFNkWmDNdZ65b/+4n/G3WpPu183PWvWvXM+3T2/6xxyAkE0YQhO0H03s6wMaNGzn88MMZO3Ys48ePr3u8pKSETp06AdCuXTvOO++8er/f/3bzS+x0FWwq2NrRF7ZfCiIFTDluCu2bteeCkWtY1rEEzjgD1q3L99AEQdhOafBivVOnTgwbNoypU6cC8Oyzz1JWVkZZWVlauxUrVpCoTavduHEjL7/8MsOGDavv4bKuOM7i0kriuJyCbLZq8x79xLtNKnpYanu2ae66jAK3OI9EUunmQeMKG3eQWPdx2iPAwEQHADbENnPKSbDb/p/zt13WE4/VpMapRK9uIkBhW2AuzA33CuxEIpUmbuKs2xaYs3Hcg15XNn26H//rX+HOOxEEYfvH9J6+adMmDj/8cA477DBuvPHGtHOrV6+mpnbitqqqiueee67e7/c1idoCc8Wyz3pTo+sOXfn72L+ztno9Z/9yJxLlS+Hii2XHE0EQtgkNXqwDTJ48mcmTJ9OvXz8mTJjA3//+dwAuuOACXnrpJcC54Q8aNIghQ4aw9957c8ghh3DuuefW+1grCxyBVRotTRfWtoXawpx17+O5KjAXjdq59jbr4WHrNeo2Qt87lqIiZ6zurdggeE177flWW+JMeyFKh5pirhq5kX1Gfcuc9rGUs+6XBu+dCFCEpcHbpKp7Y53Q16XYu68V5IbbCO6gjIFM17LrJgRmzoR63oZJEIT6w+Se/re//Y2PPvqI559/vm57tttvvx2AmTNnMmzYMIYMGcLw4cPZcccduf766+v1NcSSzv9Z0ZL633VGyD/HDjiWi3a/iDa9d6XihLHw1FMwZUq+hyUIwnZIVN8k//Tv35/3339/q8cfeuihut8vv/xyLr/88vocli8VbrGu3NpsC8ypx9xHt+OebYE5Pxc/kdBv1eb+3eTaNmnvpmJeibwAce4n3g9bWsyXn+zJraUfcNfwCnY/fDG3tS7lt0o06voqKnLS8XUF5nRF4ILaZ1L9vaAgPTZdD28i4oNeR7Zr2YOq3R97LBxwADz3HIIgbH+Y3NOvv/76QAF+wgkncMIJJ2yz8ZkQi8egUJz1psw9R9xDtCBKZMw6+GCws53bfvtBnz75HpogCNsRjcJZbzQkk1QWOmlQWznruSww5z3arFHXOeu6vm3XtAc567Yp9aaufUDBubQ4GqVFQQl3vlPCx6/1YOj6ZuxY6bNm3WCbuDps0tx17ZPJ9O1gdOvd/a4VdN62GFzYRIGts26aar9pk/OjWLkyfQyCIAh5JpZ0/k8qKJZ91psqRYVFRCIR4m1ac++EE9hcsxlOO022cxMEIaeIWM8l8TiVtfqxWbRZ5oXbwpx1vzXsugJyNgXmdH3bFpgLqgZv4tpnuh7eHfsJbHXtWIxhPxbxwYz+nPtdC4jFWLBDjN8MWkVFFDPhr7ApIGfSXle8Tufih7n26jF37H7cvaWayXp406NNqr2K16+HXr1kTbsgCA2KWCJGNA4R771HaHI88eUTXPH9PfzymiHw8cdwyy35HpIgCNsRItZzSSzG2Z/By8sPpHur7uaF23RF31Qbdxw2EWAj3v1cfJtJBptrhznrtuvhw/p2x37rztVza4VotLCISGEUYjH+tmeCCb2XM+xieL9opfOcoLR4t5j3c8PD3G6dAHevjzdx1k2ddxNBbdqX7TEsDV49lkymx+vWOXUJVq5MjUH2ZxcEIc88/dMh/DQBZ1mU0KQ5ddCpjOk9hgejtdu5/eEPMH16voclCMJ2goj1XBKP03ctHFVdRoviFsHp5LrYpMCcW2CHid4wZ72gwPmxSYO3ccO9sbf6u01BOZsUe3esRF+As048nuqrupq//Af+uKgvi9rAyF3e51eHwpZmtdexcdpNisDZOOuZFpjztjcR0qZr6zN11v2uaSr4v/0W2raFt95CEAQhXxTVJGhRg4h1wbOd21rZzk0QhJwiYj2XxGJsKYKKaJJkMmleBM59XlVjDxPn7qPJGvWga5kcTVPo1XlTMW8zzqD2Jmnw3rh2zXqds17j2rqtspJoAq5ZM4DP7ocR63fgrn1hyIFfsbYZRsXr6sZukhYfJt7DiteZuPZBzruJkNatrdcJ/0wKzJmupV+wADZscES7IAhCnngvUs4bvUmfhBaaLF136MrDxz7M2uqfOOuXZSSWlcNFF8l2boIgZI2I9VwSizHuGGjefYqzB6uNQ20qzjNxv4NEcS76zlTM65YGmKTvmzjrCvUHVXV1+rWVWI9G08T3gB9h5rv9+fPrsO/GNrSrwLzgXEmJvZPubW/rrOvS4m2EtM5Zt0mp94vDnHVT0a6O06c7lXeXLUMQ8kIyCeefD88+m++RCPXIb9p8wuknIs66UMfY/mO5ePeL+aZ4I0tOORyefhoefTTfwxIEoZEjYj2XxONUFEEkCUUFRenC23Std1Dae5iYz7TAnM5ZtxXjOqfdmzFguhTA77y371LP9jnuOKxvd1q8R3wXVlbxq/dhyqJhAMSqKzjjBPhgx5rUc13t02KTNeru82q/eEXYmvWwau9+1/JLXQ8T0rbOei4KzsXjqfXqQePye/zTT5092b/7DkHIC9XV8PDD8Mwz+R6JUI/EiBNNIGJdSOMvh/2Fzy/5nLKJj0P37nDFFXJ/EgQhK0Ss55JYjMoolBIlEomkC2/TAnM6kW6SBm8q3nXOuq1gDrt2YaGT4h+21t62Orypsx42TnfBOU2a+6ziNTy9K4zcex43HgTVVVv82ytnXQlLlXJv6qR7Yz/xbSqovbFOtNu4+Kbp8JkUmrMZr/v42GOw++5bT34IQi753e/g7rud373fTaFJEEvWivUC+TNKSNG8qDntmrWDdu2YNekG2c5NEISskbtMLok7W7eVJn2Etmlquk5I+7nfap27etzUWTeZGLBNVTedGLBda68bS5hYDxP6icTWzrpKmfeI971XF/PJAzBk8w78/gDYe8BM5nbEvzq8EtTq2jZr1L2xzlnXifewAm5+YjjMtferBp+rQnNhrr2J0Ad47z2YPRvWrk2NWdYMCrnm4YdhyhTnd+93UmgSxEgQTUbyPQyhgfJB+Qfs/ell/PqaoTBrljPBJwiCkAEi1nNJrbPejNq0uDBBHSQ8dSnqOqFt46zr3PtMUtVt1uLr+rYR8zbOetDrUgJZtfdx2geuhg++HcUN0+Hz5hvZ/SJ4tduWuvN17ZXwVH3p1qy7097dfam2tk560MSAn9hNJlPboZkI/0zXqOsEt8nadduJgl/+EkaPRhCyYuJEOOusVBw2ASY0CWIkKEqIWBf82avbXhxQdgD3Rz/j34ftBLff7kwmC4IgWCJiPZfE41REoTSiSS/Pxlk3SYs3FdAmR3eKXzYF5vyO2RaYM3XWTdL3o9FU6rQmLb64uBm/+x+8+/6u7L8Y9qlon3Y+0FmPxVJ9eQW1rbNuKqhVbFLcze+5uvXwiUTmxeHCXP1MnXXv8dNPnR/F5s1bv9eC4Ic7I+Nf/4KnnkrFYZNWQpOgS2URPTcV6hsKTZKCSAGPHvsobUrbcN6BG1jdpgjOPBM2bsz30ARBaGSIWM8lsRivTYW342c6cZgrHiQedWnvJg67qfttUmAubJ25zbpzk0mFXDrr7gJzYTUCgs4HOOvueMTqIt74B7QpbAnAhE7f8PAwSJYU+zvr8Xj6/u+mafA6MR6Wqq7iMDdbl2Jvej4TgR3k6pv2YTIR4B7vqFHO+kFBCOORR6BXL2dyB8yWg4iz3qR47ZNdeOOV9vkehtCA6dG6B5OOnMTqqjWM+1VfkgsWwFVX5XtYgiA0MkSs55JYjM6boUe09gYeJrxtC8yZOuw268ozce9t3HBdGnw2zrr3fFFRehZAJs66Nzbdqq20lC1FcP+Oyzj/WDh2xEJWllQ7IlHnrKsx2xSYy8RZV7ESx27RG+ac68R7tqnrQdeyTX8Pc9jd1ygvT9/m7Z13YNUqhCZOIgFff51y0+fMgaVL4YcfnNhmOYjQNKipkT3WBS2nDjqVUweeyr9q5vLxCSOcehcvvJDvYQmC0IgQsZ5L4nE+6gbfFNQWt8rEWc+kwJw3DhK5BQWOU27r3rsft3HDwwrM6ZxykxT7MLHv/iMqk4kA5YZDsFh3FZRrXgOfTh/AGXPgXx3XMvCcCp5rtyp1La9Y98ZesR62dZtJqrrNGnadc27q2pu63W4hbSqAdOLc1Fl3j3/LFjj4YCn8I8Bzz8EuuziFoCD4+6ayQMIKLQpNgondljO1X6W+odDkmXTUJN4+5232uuc5aNcOLrgAVqzI97AEQWgkiFjPJbEYh54JlydfceJcOOu5LDBn2mdY35mmqvtdO2z9u83ae29cUJC+963tuL3X1jnrtXHbTXH+8Tw8/Z2zL/uJg+Zxwz4V/mnwNTWOixck1t1xMpm+7YttWnzY+Vw466aC2k9YB43Lpg+/x/2K1nmvWVHhnFNrCCsr4aKL4JtvELZzqqvhD3+AlSud+Mcf048mk0FB30GhSXB7/1Xc11/WHwt62pS2YVSvUdC1Kxsn/ZXYujVw/vmyU4kgCEaIWM8l8dqt2wpqhWJYmrn6XVWJz0X6uzoGiVzbvnSTDDp33MZZN3Htw8S7+3wmQj9s4kCtf9fsw67in20p44v74MjVrTl8WbN0Z13XlxqXruCcjaDWOdjZrFkPE9gm7niQkNb1YSPaEwn/dfHe42efwQMPwIsvpsbkrdIvbB/MmgW//S0884wTm2ZyuOOgc0KToCaSICp/QgkWzF87n6GrbuH2K4bAq6/Cfffle0iCIDQC5E6TQ5I1NVRFoVmkNgU7rICbaVq7rRseJnJtUumDrqUTzDYp9+415jrBbCPedS59WGaD37WCKs37Ce5aV7/LRvj3+73Zb00LKCxkZpv13DZsA7FS/z3c0/Zod59X6ArO6da0hxXHMl2THhSbuNvuOKi4XZhbmQvR7n29uuPbb8MOO6RXkxcaJ9XVMHAgPPigE6t/T2oHiEwyOILaCE2CWAQR64IVPVr1oGVxS37X7ks+3KML/PrXTq0MQRCEEOROk0MqayoAKC2oFWRhgtrUQde54TYF5mxS6XXjVs8LErm6cUaj6ZXm/QS1jfvtfZ1h4ttk0iFTse5+/6uq6q71p51Xc/OILRyw33wWtA14rl+s0G3llo2zbromXRdnIqjDJhH8+jA9bzJBEFaUDmDBAqft4sVO/N13Tuq0pC42Dt54w8mWAFi3DubODV6TbvPdTSb9v/eSBt+kiEWSRCPyJ5RgTkm0hKnHT6WwoJAzTipkU6ISzjgjNWkoCILgg9xpckhl9RbAJdZNCsxl6rCHCW5dKrquj0xS7L3XymRiwEZAh7nlYe+BSd/e80FiXf2R7hbY7mu7xPuT73Zh/Kwi3mu7kaEXw2MtvycJ+rR4VSjPnY4d5m6Dvkic+w+DXDjrOtEbJnx0znqmEwEmqfemfavjo486qdMLFzpxTY0jAoWGgRLRirPOgmuucX7XfbY2jrp7SYXfc4QmQawgSRGF+oaC4GJQ50FMGD2B+RXl/OqaIfDJJ3DbbfkeliAIDRgR6zmkuqaSVpXQsrBWgIWJ4WxEra6dqbOeSYE5W2ddd23T1HUbZz2TFPqw82F7uLvPV1amf7aVlXV9l1Yn+MtbRfxn3u60rIazd/6c00+s3ZMd7NLidWJcJ+ZtU+qzddbDhJK3L1sRZZOCn0g4gs5UuHnbqUkOdbziCthtN4QGwnnnwUEHpeItW5wigmD+PTLJ1NB954QmwRELo+y1uXW+hyE0Qn6x9y8YvdNoHij8jDcP7eNkbL37br6HJQhCA0XEeg7pXNiK9RPgrnanOQ+ECXLvmm7TNekm7rcqWqfO265ND4tN08lt1tirdjbOuu7aOqFvMxEQJM69cYizrv7IP6SiK1/cByesak/fNRAp0TjrfmI9G0EN4S69idMeJO5Njkow+13b3bepsLYR7aZ9mLZbvjx9+53f/hZuuAGhnvjPf5zq/er79PXX6ZX8w7IoTLI+gs7rnis0CZ57vpibVg7I9zCERkhBpIBHj3uUWw64hVF/ecbJoDvzTNiwId9DEwShASJiPZeoP9ZM08z9UtSzSX8Per6t2x0WmwrqTK5pI6i9qe3ea9u49LqJgCBx7o2Tya3FuoqVWCwspH11lGc+LuOm6c5zfyqF37efS1UherFeWpqZs+6OdS59kNvt15ets26Tvh/mYpo4oN4418JNHVVa9DPPOPt1K+6+G955ByFHrFzpvKfqfX/ySad6/08/OXHYxJLN98TvcRNnXcR606K6On2LUEGwoHur7tx84M0UDxwCf/4zyYUL4aqr8j0sQRAaICLWc8jqyjU8OBy+jNe6bTrhGpY2nqno9XvcVvCbCuowEWyaCRA0EeC3ztw0dT0TcR42EeD9g8zrtLtj93O9Yj2RqIsjVdUUJoGSEv66N9zYcz57XwBftU+mngv+4l0ncnXndS697Xp4U0FtOxEQ5mLarlXPRLhl4sSrNskk/OpX8H//l4offBCWLUMwpKYGnn7aSWcHmDYNxo9PVefXfb/ccaZHv8/e/RknElu3EbZ7amLVnHlMjAc6Lc33UIRGzoaqDZzd9SMePGs3eOSR9AlfQRAERKznlPkV5Vw4Ft6onOc8YJIKrnPSM00rj0a33sPd1I0PmyBQ6fvq8SARnMnEgG6dedC11JiycenDJgKKitLbhznt7udWV6fG7d5H3R2XlnL9O3DT/O583hmG7zeXiXtCsipgzboS60oYFBTYi/dsnfUgVz8WSy/yFbbeN+hamU4EmBwzFWqZTgi40/i//x4uvDC1dVhlpfP7pk0ILv79b1izxvl9+nQ4+eTUvvfqe6uOfp+xW0j71UAwddTDPnvdhICw3VNdtYWpQ+C9HX7K91CERk4ymeTtxW/zy34L+K53Gxg3zlliJQiCUIuI9RxSGXP+iGwWNSgwp4424jusL5P0d1PX20ZwB4lcW2c9zCkPupbJJIjqKxtn3dveK9bdzntQRoASGKovVwG5ogTcOq8TMx6BLjUlXH4UHN37A9aX4F9wzi08VaxLVQ9Kg7cV597zXoGdSJinwevS+W1S102OtsItkwmBoHGpz08d//tfR7y/8ooTl5fDH//Y9NzZF16AL75wfi8vh6OPhnvucWI1kbF5s3O0mSzxVmw3/axNJoOCvqNN7bNrwtRUOt/JaKRQ01IQwmld2prHjnuMilglZ1zUkZqf1jrFMmWLUEEQarES6zF3YSphKypitfusK7FuIlRNhXMmrrhpmrttWrxfir1y8bOZAPCK3DABrZsYKCjwb+vXty7F3ts+yEkPe26Qs+5yzvddCp/N3Z+zPoN1BdW0qEGfBp9JWrzXWc807T2ovYm7rZ6rE++2Ij3MIbUVbtk467qjqlKujtOmwXXXweefO/Grr8Kll6b+YFMp142dRx5JZRcAnHoq3Hqr87sS5+qomxzx+zxt2trE26Gzns/7+Xfffce+++5Lv3792GuvvZg3b55vu7///e/07duXnXfemQsvvJCY6z1++eWXGTBgAH369OHEE09kUz1mqcSqnfeuqCCqaSkIeg4oO4Cr972ajyq+4/Yrh8Lrr8PEifkeliAIDQQrsf5g797bahyh5OLGXh9UxpwbeGlhrcDyE5PuSu1hzrpp+ruJ+M1UlJs48blIrff2ZRLr1uK70+J1699Nzgc566br49V3UcUBa9JbFe/AlBfgzQ/6EU1AsrSER4bClpotqfZuUeKO1Ti8gqKmJn2W3i3WEwnnvMLG7faLdYLVVJx7z2eTtuwdZ7bCzWRCwFRoqlh9JmpbuGnT4L77UsL1/PPhiCOoY/p0p/p5Q+Onn5yJBsWECbDnnqn4rrvgz392fk8mnUwDJRpNP8ug91i9/za7CGTitAd9hxuZWM/X/Rzgoosu4sILL+Tbb7/lmmuu4fzzz9+qzcKFC7nxxhuZOXMm8+fPZ+XKlfz9738HYNOmTZx//vm88MILzJ8/ny5dunD77bfX2/hjVc4kW1TEupAjbjvoNoZ0HsLv237BB3t1hauvhq++yvewBEFoAFiJ9c0rV/L4Pvsw54EHqN64cVuNaSuyvbHXF5Xx2jT4ombOA2Hi1htnmgZvI6RNU9IzcdZNRXuYwLYR0GHOum4iIFvx7nXWbVLsg5x1z77qzSucP/xf7byR846DPXaZwWc7Eu6sq77C0t51sVdQ65x0E4Ed5qyHTQy442zcbu+1TIVbJq6+rdA0bffFF/Dll9Rx0knw6187v1dXw/Dh8PjjTlxVBXfeCatXp/qaMyfdpVeOvord73tNTXq8YEEqFT2RcNLU1ZZ1c+fCkCHw3XdOPGkSHHkkLFzoxB9/DLNmpbICampykyUR1t5d9M30/ddN5PhNyKi4kabB5+t+vnr1ambPns0ZZ5wBwIknnsjChQtZtGhRWrtnnnmG448/ns6dOxOJRLj44ouZNm0aAK+++ip77LEHAwY4W6ddeumldefqg1i1iHUht5RES3j8hMcZ2XMk7X9/l/P/+umnpyZwBUFosliJ9eadOrH/7bezbMYMHuzdm1fOOoslb7+9jYbmkIsbe33RJlnC3kuhY0lb5wG1llkd/dZCB4lY3dFE1Nuuhzd1wU2cddNrZ+p+hznruomATOKgNHhTZ93bt+G+6odVdGPCGzC/dAt7jYM/9f+RBElH9EQizv6sXrHuFRRece5Nf7VZw+4V3H7OulvYeF1O275txZT76N7T3d2Xdw93nXNuI/RsnXWb9u73ZvPmlID+4QenQvr77zvxhx/Ctdemqgo/9RQMHQqzZzvx3XdDly6p5191FeyzT6rv0aPhtNOc32tqYNAgZw95cCYNrrwS/vEPJ54920ndV+n7SvSpjICwCQ3b98jm8zFtm8lkStCEUyNz1vNxPwdYunQpXbt2JVr7/2IkEqFnz54sWbIkrd2SJUvo1atXXVxWVlbXxu/csmXLSNTTUpEWyWKumwGjE2X1cj2habBbp92Yfs50+h5yivN/7qefwi235HtYgiDkGatp4VETJtDz4IPpefDBVG3YwNfTpvHOtddS8eOPDDznHHY75xxa9eiR0wGG3djLysrq2oXd2L3cdddd3HXXXXVxrta6HTXkJI6aUwG7HOY8cPbZzrZe3bs78UkngTv1sE8faNfO+V1VHG9W68rbiNygY5jItnk8THgrMWTj+AcdTdeoe2ObDADTa4UJ7mycdXVe/VGpEeuFJaVc+y4cUroLp+0yj2v6L+HVs+Af69bTTfUdjzvC3eusFxc7s/I2znom7neYwHb/8RyWMq+7VibOdNAkgnecpsLNJIVa7ZSQaZ8mrr039rb1ptSr6upr1zrH77+H9euduEULxx13Ly36/vvUdmmVlc7vqg810RNWkT3s6Cduc5V94G5vO9FiM5kS9Dk0Mmc9H/dzRUT9O6klGVBMy93O28bbhx/b6l7ftmtv/nD2Y7DLLjnpTxC8TDuxH61m9+GoCROcJVD775/vIQmCkCesnPWB55xT93tJq1YMuegizvjwQ4578UW+/9e/eHCnnXj60EP5+okniOcwdScXN3Y348ePp7y8vO6nZcuWuRnoqFHw0EOw665OPHQo/OEPqT/gL7jAWYeqeP11ePJJ5/fiYmdPYeVglZU5wn7YMCdWQqy42DmauOWZOug2Tnwma+mDjt5t4UzSyU2u4efShznpOsHt3lc9Vy69xmkf/kOU2ZPh4lU9+KQLVMcq09fD+znr3hR7hU0afJjb7Rd7xZhO2AcJH2/7TJzpsNfhJ+hykXJv48KHvQ4bRzcT4e9tH/S52Lx202tvi6wDb3vbiRabyZTtxFnP1/28R48elJeXo2rKJJNJli5dSs+ePdPa9ezZMy2DbvHixXVtvOcWLVpEt27dKChI/746obIAACAASURBVJNmm93rd9gBzjwT9tgjN/0Jgou1FWu57PVfcM7+a1jZvsT5rm3YkO9hCYKQJ6zE+vIPPkiL1y9axLu33MLzxxzDylmziBQUEG3WjG+efpoHevXineuuo1qlWmZILm7sDRbv/t3HHw99+zq/t23rOFzHHOPEhx4Kt98OY8Y4cfv2qXYQLqQzSW8Pi03S4HMxMRDmWPtlDphMKnjHHXTedM16pmnwCiX8XfuuB8XNa+C+5cP5+l7YaX0BRKPM7FTJxkjtGuMgsa5z1t1iPkzUuvsOOu8VeqbCXtc+E2daN05TQW3jxtqKQ117v/fSu4d4tuI27DPeFsLadDLBZgLA+7hpW5PPKSw7xB27l100AvJxPwfo1KkTw4YNY+rUqQA8++yzlJWVpWXKgbPk7fnnn2fVqlUkk0nuv/9+fv7znwNw+OGH8/HHH/N1bZHFSZMm1Z0ThMZOu2btmHz0ZH6sWsf543cmuXixs/xIEIQmiZVYf/GEE4hVVTFv6lSeGj2ah/r04f3bbqOoRQsO+NOfuKi8nONffJFjn32Ws7/4gkQsxr9OPjmrAebixr5d0Lq147q3auXEl1zirB9VExJ9+zq/+znryrE2TW+vzwJzYRMCQSJYvR6bcevEeZiz7m2vshu849KN268v757tAQXn3OK9yybnfHmbAo4YuZBhRyzmg9YbU+NSAsLnuWmxwtZZ1znWfi6qX99hjq63valI9HNXg8aZS9fbuw7eVOhn4kh79xDP1llX7417fX82Lr3umKm7bfP52H6WYZ+T7QRTIyEf93PF5MmTmTx5Mv369WPChAl1xWAvuOACXnrpJQB69+7NrbfeysiRI9l5553p1KlTXXHZHXbYgYceeojjjjuOPn36sGzZMn6rstIEYTvgpN1O4szBZ/JK9Vwmn70bTJkCzzyT72EJgpAHovomKSrXrOG+zp2p2rCB0jZtGHzhhQw891y6uLfmqaV5hw4c+Oc/M2nHHbMe5OTJkznnnHO44447aNWqFVOmTAGcG/vYsWMZO3Zs2o09kUhw8MEH+1aN324oKYGBA1Px738PN9+c2hZujz1SRZ7ALC0+Ewde/YGfzVr1oHZBAjoTN1+3/t3UWVdF3XLRFwSLdYMCdN0qotzxTXeu7r+E/Q4t56ZvYvz260KiQc56SUn6VlkK3Zp1nbMeJoJ1wkYn3m2Fs59g0/Vt45wHxV53Ndu060wc6WyFdCLhfFf9BHWmkyZhkw7ZivKwMWQ7OeI3uaBi3QSU+993AyZf93OA/v37874qiOjioYceSovHjRvHuHHjfPtQ935B2F6554h7eGfxO4zv+z0H9WlD/4sugn33ha5d8z00QRDqEau/KhKxGF1GjGDguefS5/jjiXqFhovqzZt564orcpIWmIsb+3aPV0Q+/HD6ez90qLPFEpg76ybtlMtn66ybiPYg0ZtJobxsUtW9GQQ2hfC85/1c+kgk9Vmpf1Nq7/MQsR6JFnFFeTcOmlfB6aN/4uZdVvL62RGmrtnATvH41k66ty9VnE5XDT5MrIQ57TpX3uRa7jiRsBPtppMOtuJXNzEQNsmQqYh0v3emR9sJgXg89Z3I1lk3mXTI1v0OG0M274H3eUGfrV/s/r2Bk6/7uSAIZrQubc1jxz/GMdOO4ZsbL6X/2XfAOefAa6+ljBlBELZ7rMR6806d+Nnrrxu1jRQU0P2AA+h30kkZDUzIAe7CfB9+mIo7dHCqzqtKv7bryd3F4Gyccr/zYUI7SEBn66xnUmAuyPHPdv27aq/+yLdx2mufO/DHAj58vgO/HQ13l63g/c0/sJPuuSresmVrZz0bQe0VwbZ9hU0EuM/rxKFOZPn1ZSrg/OKg9yBTh90vTib9JyxyeXRvBWgyvkwmPEw+O5P3JGwsmU6K+D3PdNLH/bxGgNzPBaHhM6rXKJZctYTWpa1hxmqniPE998AvfpHvoQmCUE9YifVLVqwwblvUrBkDzz7bekDCNsIt3Nu0gZUrnYq2kFoHr44m68vd6+GD2vjFNscggW3rrOtS1U0Etvsa7hntXBWvCxLrXne8qCgVl5TUPbc0Wchd8/txzjsbGNzd2Q5waWtoWQptbcS6n/gOEydhAtvEhdQJ6rCJAPUcd+wn2PxiE/fb5ujtOxtXWNc+WwfdZAym4zN11P36NhHKYe3c7dVEhordfbhrCdi+Ht131G9CqZEg93NBaBy0Lm0NwKJbrmLtnNcYfu21MHp0+lJIQRC2WySPpqnSqlVKwI8cCW+84VSjh5QgNBXSYevhdU56JqnrthMBfqnoYe64yXr5sLX12RSvC3LWY7H0vqqqUu+JEhTRKIN/KoGqKmIFcNKw+Qy5BKYXL3eeE5QWr9awFxXZixOdY20rzsOc9TCB7SfYdKLKVMyaOLxBExrZush+j+dC3Aa1D1pzb9qnSTvbdf0mExo6tzvT99tmAsn9PEEQhBxSUVPBPk+M4cSfxdlQUAOnn751sVhBELZLRKwLjuAbMyYl3o49Fi67DHbf3Yn9xG+QkDZd755N6nomqfVhzrrNOnSTtfU21/K2DxLr6ppusa5iJSjUuKuqKEzABT/0YE0zOKj/B/xmNFSXulx58BfvtuLEG4dVmrdNKfaLTR1QEwFnKmZNROW2ctbDrmUryk0EdSZOuek1vQ63bkLDZmIg00kek6Pue9VInXVBEBoPzYqaceOoG1lUsYIrrh0En38O11+f72EJglAPiFgXtqZ3b7j3XmddO0CfPtCxYypt3tZ5dh917riJwLadEHCvsdeN3/TaJhMXtn1DuFh3tw1w1olGneJzwAXrd+bTybD7hpZM2B/2GfUt37QneA17aWm48DSJddvCBQlmv9hEYAc5oCbOejYCzn3tsIkBW2fdRHjrBLSp2PU75tKl97Zz71fuNzGgE+VhEwSmE0qm751fJoPftfyuLQiCsA24ZI9LOLLvkTzGHJ46ri/85S/w3//me1iCIGxjRKwLeq69FsrLU+KuuDh4+zITl9s0tnXWbZx3G2fd79pBkwom4j1MzOvEumqbSKRiJRjUtVziu98aeO/d/lz/Dny6w2buHEl4wblkMlWJ3k+8Z+us65z0sImAXDvrtsItzAUPeh3ZTAAEvc5sJxnCjrkS5zqH2s9ZN01V1wnqsM86k8/Dxll3/y4IgpBDIpEID499mI7NO3LxXqsp37E5nH02rF2b76EJgrANEbEumOEW5zfc4MzoKtwC1LQKvI2zrlxxnbNu47jbOOsFBeGTDro0eBtnXaWmK4LS4N2xEicqVmvQa59bVFnD79+CGXNHcNfrQGUlNQXwY8vaf/5e8V5dnYrDUopVnI2zHtZ3Lp31MEGn1lDnQuTqxmkqfsPa656TzdFWnGfSt3r+tnhP/OIwZ133OnTOut/3ShAEYRvRuWVnHj72YdZVr+cf1x0By5bBxRfLtoqCsB0T1TcRBA+jR6fHXbrAjjs6v+scdhsXPJtK8roxmK4rLyz036IuzFm3cdJt1qwHjVMJTa+z7nHOR9bsCFVOfOuB8PcRc5myCg71iPu05ysx4p4UCHOsVV8KW2c9TMz7ia4gwWzSV5DAMxFySuDbjjMTF9l9LT/3O1PXOMz9Nh1fps66iYA2nejQffYmEwNhznpQX+q8OOuCINQjR/c7mvfPf58RXfeCGSfD00/DUUc5LrsgCNsd4qwL2TNzJkyd6vyuE+cmVeDV47pt4WxT703dcF3au/eYTYE502rw3nH5xa4162nPVQLaVVBu1x9gS2GCw86Eq7rMoSLq095dcE6d0znYNs66TviECW6dGDTpK0jgmQg523Fm4yLr0sfDhGbYNcKEaTYud9Dzgt7/TN1808/eJNMh089anHVBEPLA3t33JlJQQOy+iSzr0xkuvxwWLMj3sARB2AaIWBeyp02bVDG6tm3h8MPh0EOd2NT19juvHO1MRXk268xt0vh1a9Rt0uLD0uBNhX/YmvTa+LQvYM7cUYxcAn/ruoTdL4JZHWu2bq/EitrOTye4sykwZ3PeJM0602vphFwmade5EqJ+LnEmojyovfpdbYOnc8wzmQjQvQ7bCQOT1HTTSQeVyRDUV9hEgLsfQRCEbUw8EWfMv0/hsHGlVFRugjPPlP+DBGE7RMS6kFsKC+HVV+Hkk504aM9207Xruag0b7LOPCwNXjdenTjPtMBcQYH9+vdoNL1AHAQ67WWRtkx/BP7wTQ/mt4PThy0gHiE4Dd5deV7VMLBx1m2cc3VeJ1Tdx7BJBBPn3UTIBV1L56zrxK6NyNX1aSukg97DsAwB08mGXGQIZDOh4b2WSZ+672TY90zS4AVBqCcKCwo5dOdDmVuxmOuuHgrvvQcTJuR7WIIg5BgR68K2ZfhwuPpqOOEEJ7atCu8nnHWuvMnjYc667QSBLu1dJ7CD0uCzXQ+vBHWI016YhOu+7sCsB+Af8wdTmHTO/9gcR+wrsaLGooSPmgiork4vbOMV6+417JkIoSD3OxZLub9+fXlFron7aioWTSYG3OcTCb37bTox4H6d8fjWe5f79aVz1N1H0wwBU9fbxlm37TOo77DPWn1nwvqyTXsXZ10QhDxx7chr2a/nfvxfyWe8PqYMbrkFPvoo38MSBCGHiFgXti0lJXDnnVBW5sQqXV4dTdaf27ryOmddl07ujk1S73Xr33XXCnLWbVPog16HLi2+spLBq2Cvmk4AvFO6il5XwV391pCI1fg7696+1HWz3brNVMy7t5gz6UvnrHtFsHqOOw4StbmcCPC7tulEgM1kQ9AxbHLEb3w2x7CJgGTS3K3XufY2Trtf32HfQZNYEAShnigsKOQfx/+DViWtOGf0Jn7coRBOPx02bcr30ARByBEi1oX65dhj4YEHYP/9ndhk/XkmqfPePrznw0Sv+7zJenidONe5+O72atmAeq7ppILftdR5JSCC0uI9cWlljG4b4Vf9F3HQsetYWLQ5NW4ldELWw6dhs2bdNk1ebTHndy4Xgtom5d5U4GXiIuveM53wD7uGNzshaILDNB0+7Jph75F7fbztNW3eI+/rMhH+4qwLgtCAKWtTxsQjJ7Ky6keuv2Z3mD8fxo/P97AEQcgRItaF+qVVKxg3LiXuWrZ0ji1aOEc/YZ1J6rzuqEsnN3XWdeLcVmAXFWWWrh/0OvxceyWgA8T7XquL+PTBAi5b3p13usUYfMBcHhywhWRhQUqcZCLWw0STik3Tmb1924qqMEFnktIdNjEQdm3TVPUwZz0o1vVt60jbOOsm19R9PjZ9e7fOM+3b294mXV+d1332giAI9czpg07n3iPuZcKV/4IxY+DBB+GFF/I9LEEQcoCIdSG/XHEFvPQSdO3qxIWFTmG1SMSJM3HWTZz2IGe9oCC9sJvJRECB65+RiaAOE/NeV992IsBkPXzYvuq1cYtkEfcu3o03niqlTXUhl4z4ga87RlLiJCSlHth6vTyEiyYV68SjOw6bCDCZGAiKTUS7jbPu59qbilybSQdbZ91kciSb1HRvarvus7e5hulnafq6TDMZTD4PQRCEeiYSiXDZXpfRtkV7mDKFWPu2cMEFsGJFvocmCEKWiFgX8kvHjnDMMam4ZcuUyw7+LrdJ0bewx8PWketS6f2uobaY8/alYlNx7m1vUmDOpu+gPdz93PHaa49ZAF+83pun3+/OLhtLIRZjZUkNyRJ98TogvcBcIrH1OvNcOuu2gi0bZ91UnHvP5zp9PBtn3UT0mopzv9fj3gLNZvLEZCLAVPjr2ptMBNhO+giCIOSR2757iAOv7UzNujVw3nnpRWAFQWh0iFgXGhYTJsCLL6ZiP2fd1EE3LTBnUv3d5ho2a9jDUvJzIfyDKs27Y7+0eDWOeJw2FUmOX9UOCgtZH6lijzMrOH7vxaxoiV6sewvQedeZ24hz73ldpXmbdGad2+rdezvTSQaTo9uRNhGPNn17j7aufaaC2uY9M7mWrYAOep22kyMm31lBEIQ8srFqI+9u+ZrbrxwKr70GEyfme0iCIGSBiHWhYbHzznDQQam4tDQlIk0ddBPBrdLd1Xnb9fBhj9tUh3cL7Ehk6xT8bNLeveeDxLoS0D7Oep14qR1LPB5jv6URXuy8jl0vgykt55P0Phf0xexAL7p0zq+u0nymzrqfaLRJhQ57HTpH3W9iwFQ8BjnsYUfTdPFMnHXd+52Ni28zqWMyMeM+2jrr3msJgiDkkd8f/HuGdB7C79p+zvt7dXW2z503L9/DEgQhQ0SsCw2biRPh7rud303dbdtU9myK2OmcdRuB7XetMOfcNg4S695YifXCwtQWabVxuwp44ukkz30zlJI4nLPzFxx1OixtE0k919uXO/a64TaC2hvnssCc+7yqkB7k5PqJXFv3O0yQmopaFevS3m2cdV26uM21cvU6/B63mdTJdmJG56xLGrwgCA2IkmgJj5/wOMWFxZxxImyMVMNpp21dAFYQhEaBiHWhYXPwwTBqlPO7Ep5qezNTZ91EvGfizgedN12T7o39RLtO+NvEXnGu3G/veZUGr67tdtprXfjjN/dk3kQ4a3lHXu0LT/Tc4LQ1KF5XF5uIc68QcrvObuGvS082EXRBAtDEbbV1v8OEqK1rHCY8/Y5hEwNB7793AsMkU0AnmDOdGLCd1NFlDNgUwjO5liAIQp7ZrdNu3DnmThZULOeqawbDnDlw4435HpYgCBkQ1TcRhAZCNAp33AH77JOK/Y4mgtp7XufO2xSc8wpsVdletQtbo67ryyb2iv0gce6NY7H098S9hr2ioq5tuwqYMrsX5/7nB/bbeSdgET/EN7KpDexkUnle55B6hZJuD/dM3VVdnMlEgM41NhW5NpMMupR1v2vrJiyCXqeJsNYJ5jABbTNu2/fMOxbvJIRuwkgXC4IgNAAu3+tyvlj9BWfudhq8dh38+c9w+OGOCSIIQqNBnHWhcfGb38CBBzq/ZyOovY/bOuthferS3oNiv2OQsDeJTQvMQfpaedWXn7Pu45QfuAiiJc0AuLzXXAZeCvf0XUciQrCzXlqaLoyKivSupc0e7iausakoy3YiwH3eVuTauMY6Z92vnakD7W1vK6hNJgJsxm/zHuk+65qaVLVkm3H6xSLWBUFoIEQiER445gH2730gTJ3q7LRz5pnw44/5HpogCBaIWBcaL717w8iRW6fJZ+Kwm66D14l2v3XjOkEdds1snHXTAnNBz1XjqKxMjVsJZLWPusc5H7ekA502w5W7LmLkefBFQe0fBX4F59wiqrhY75hm66zbuPi2znqYGHRfy289vGrjNy4/QR0moL0p3UHXMLmWycRALlLs/eJs3iOb9wzsih6afK8EQRAaGO+VrOaMmwcRX7Eczj9ftnMThEaEiHWh8dKiBcycCWPGOHEm6e8FBen7pGdaYd5v/btJ7DdRkKmzHtY3hIt1XV9eZ90b14rxMctK+OI+uGp5Lz7qBsMHvstvRkNNaW2dAe+a9Xg8NSmhBJ0aRy6d9WzcWFunVnftsHRyW6FpInLdR6+Yt30dNin2pi69rr1Jir0uDT6sfdj3SJd9IM66IAiNgKfmPsXjm9/n9ssGwUsvwX335XtIgiAY0uDF+pYtWzj11FPp06cP/fr147nnngtsG4lEGDx4MEOHDmXo0KHMmDGjHkcq5J1u3dKPNunvNs65rk9dgTmdsx5WYM5U+PudV264wr2GXbfe3SvOvU67K829ZTXcvXIIHz0Igzc158PuEC0OcdbVtVSs+rJx1m0Esjpvuu5Z91xb0et17d1HXd82Lr/OSc/E/dY56u7zNoI5rL3fNW0mXmydddMJi6D2QiCm9/Ply5dz2GGH0b9/fwYPHszJJ5/M2rVr686XlZUxYMCAunv9k08+WV8vQRAaJRPGTGBw58Hc2nEu7+7VBX71K/jyy3wPSxAEA6L6Jvnlz3/+MyUlJcyfP5+FCxeyzz77cNBBB9G2bVvf9u+99x4tW7as51EKDYKRI2HtWlDfDZ3D7rde3HTLtmwKzLnPFxSYb90Wlqrud94r3ouKnFgJCm8afNh6d9VXTU36tUKqve++Aj6cMYCf5n5CZJ9SksAfOnzF+S2gs7savLqWElUlJbBlS0pURSKOIxwm1r2xiVgME3g1NcEusZ/ItRHvYQ61Sd+msZ/7bZvSbZoRYOJ+2wpsd99hFdt1Ey+6segmfWwzAoRATO/nhYWF3Hjjjey3334AXH311Vx33XU88MADdW2eeeYZBg4cWK/jF4TGSmm0lGknTmOPB/bgtBOSzPmymjY//zl8/DE0a5bv4QmCEEKDd9affPJJLrvsMgB22mknRo0axYsvvpjnUQkNFvcffeoGpBxkP1Fsm/Zush7eVGAXFqan4PsdC1z/RHWufTZp8Zm6+Jp91aOV1XTYApSW8t/ecH2v+Qy4HB7c+ScSsZr0NHivs+6NbcR6JmLRLcpqt6gLfK6ta6xzv4NEbiZp2GHOum1fOmfdPW61Hj/oubpMhqBr22YbmAjqoIwNmwmLoGsLgZjezzt37lwn1AFGjBjBggUL6m2cgrA9smvHXfnr4X9lSeVKLrx2V5Jz58Kvf53vYQmCoKHBi/UlS5bQq1evurisrIwlS5YEtj/wwAMZMmQI48ePZ/Pmzb5t7rrrLrp37173s2nTppyPW2gAnHMOPPgg7LabE/uJYdu0d5P18KYC2uRakUi4a59p5XkILzBnuj5eI9bd8egF8MzsvjSLwYW7LeCA06qYV7w+1ZfbWYdU7C1Op3Dvs+49n61YtBVwNuLQm+KdSKSK/diIWN3rytZZD5uE8JsQyOV75G5v4pxn47S7v0c2n13QtYRAbO/nAPF4nIkTJ3LMMcekPX766aczaNAgLrjgAn744Qff58q9XhDSGTd8HCftehIDRh5L8sADYNIkEANMEBo0eRfr+++/Px06dPD9Wbp0KeCsRVckQypYLl68mFmzZvHee+/xww8/cPXVV/u2Gz9+POXl5XU/kja/ndK5M1xwQSoN3U9wZ5v27tfe7YaHieBMXH2btHid4M6kwJw3Vk6qV1Cr2FUtPgKcuKQFX90Ll63sybs9YMjwj3hiQMy5tp+z7hbvQXu2K3JZcE5XIdxUaAZd29Q1zsZV9jvavAcmrn2uJgLCXrfus7P9rL3iPZfLJ5q4s57L+7k6f+mll9KmTRuuuOKKusffeecd5syZw+zZs2nfvj1nn3227/PlXi8I6UQiEZ742RPcNvr3FEx9HNq1g/POg2XL8j00QRACyLtYnzFjBj/++KPvT48ePejZsyeLFi2qa7948WJ69uzp25d6vEWLFlx66aVSYE5I55BD4PDDU+nx2RSY0611NxHYOmfdr53O/TZ11iMRZw17UFs/8W4i/HVOe2Ulravg3hXDeP8h2O+nHThgVSlEo8TjMZKxGv9t3sJihY3A9nOCTQVcJu53kODT9W0rFv1S7LNxv8PcbnffmYzT9HXbOuW27d2fdSKRqssQ9Fzdd6EJk8v7OcCVV17J0qVLefLJJylwTYKq5xQVFXHVVVfJvV4QLCiIOP+W1rdvyQ23j6Zq/Vo444wm//+XIDRU8i7WdZx00klMnDgRgIULFzJ9+nTGjh27Vbt169axZcsWABKJBE8++STDhg2r17EKDZxzz4VXX0057WFp8Nk+bpIWn4l7b1twLqzyvG2xOpOUerXW20C8j1gG/5vZhy7VJVBYyF97r2bU0T8wp0M8/Lk6MQ/2InhbOuthLrPtpIKNg23irOuEaFB7737xJpMMNu+/u302KfQmY/FO+rjrFZhmI6idFpq4s67D9H4OjlCfP38+zz//PMWunSw2b97MTz/9VBdPmzZN7vWCkAEPzn6Q21c9zW/HD4G334Y//jHfQxIEwYcGL9avvvpqKioq6NOnD4cddhgTJ06kXbt2ANx///3cdNNNAHz99dfsvffeDBkyhEGDBrFmzRr++te/5nPoQkOndWvnBxwB767Mbpoeb5LSrtuqzeZatkXgwiYKdGJcJ/yDUurdcZCgVn98V1XVXXt1cQ3vda5h+P5zueIIWBfbmP7coDR473nILOU7SMAlk+luq0lKt6mDbZueb5Py7ees2/ZlKsBNJgJy5axnkppu6qxD+vfI7z3yG7d76YYQiOn9/N133+Wee+5h0aJFjBgxgqFDh3L88ccDsGrVKg466CAGDx7MoEGDmD59Oo899ljeXpMgNFau2vsq9u2xL3e1mMNrB/eEm26CDz7I97AEQfAQ1TfJLy1atAjcQ/Xiiy+u+32fffbh888/r69hCdsDU6ak/6Hul6Juk+7ud14ViTO5hknqvc2adF0Kvi7NXSf8g4rVuWM/sV5YmErBr6qqu/YfP2nLzz+t5vKjC7l3xHqeqPkft5fDeZFi5z8qi2J21unLXgEXVrzOdo10mFg0Sau2dcPDnHVbl97UHc9kIsB0UiKZ3Lo6v83r0F1bt5zC7zUXFztjUnFJCWzcKGmkGkzv5yNHjgxcz967d28+/fTTbTI+QWhKRAuiPH7C4wy9fyhnH7KJzz8tpPNpp8Gnn6aMDEEQ8k6Dd9YFYZvRqxf065eK/Yq/2brffueDBHUmfWfjrLvHYuLS26TcB4l1vwJ07vdAifVoFGIxhq2MMHPuXjz2QoTCBNx6AFQ1cwl7d99hsa376j1vux7eZi14LlO8bZx1W5feJLU9m4kAm+wEm235bNPgdbsK+D3XnfYuzrogCI2UsjZlTD56Mqur1nLOr/uQWLQQLr00tUOJIAh5R8S6ICjOPBNqUy21DrtNGrz6vaAgfF/1bLeFMykKF+Ssm1SDD3Peg8S6N3aLcxWrvmpFViRaxJnzivju7cG8NA1aRJsD8GDr75nbEbO0eLdgCxLvYeJRV2ne1v0OE9S6NGvb9HIVJxLOj1u026wbt3nPMhHMNhMDuortNlkUNTXpf4iaftbFxalxqX8P6nWIWBcEoZFyysBT+P6ovAAAIABJREFUOHfouazrtAPrjzsc/vlPePTRfA9LEIRaGnwavCDUG5Mnp34PKxjnPpoUmHMXtAt7rol4z2R7NXcctGbdz5XXVZo3cda9cU1ttXevs662bovH68ayw5YYu68AhpaysiVc1vsr4pfAuGXLuPUb6BzmrG/YkBJZQWnxbmFlI+Aycb9NJwZsUuiDzrtjXbG0bNLH3bGtYLadGLCt9J9p2nvYtdSkTzye+vfhToNX1xYEQWhkTDxyItGCKEXHb4BZw+Cyy2CvvWC33fI9NEFo8oizLgh+tGjhHFu1co7ZOOu6dHddGnyYwNZVcA8S3CauvM21iovT1+artHeFW7xrnPW62CW+d9wE736wGyOXwuQeq+h7JdzR4Su2FOHvtPuJrFgsNS6btHd33+q5uXK/vdfKxJG26ds2pd50DXu2Lr2Ns26bJeGNw9Le3efd4lzF6j1S32F17drvYDxWjSAIQmOjWVEzigqLoH17pt19HmsLquDkk6F2lyVBEPKHiHVB8KNDB2crk1/96v/bu+/wKKr1D+DfTQ9JaKETUiihpgBSvfTOtYAURVBBFJGLiijX8hOv5SrqRURFBRVURJEiXvSC7SIiCheEUER6SaN3EkjP/P4YJjtZppzZXTKb5Pt5njybk52dOTu7Orznfc8Zue1uZl1rQTmz1xq1zRaBEw3e3Znv7o0F55R9Kf3IzXW2Cwrksm2l7ZI573AqAOs+AlbsS0bdbOD/og9g5Aho38rNtQzeJagqlSFV+qFmZR6zSEbatQxbmcuvtW+rAwHezNqLBueuz1sdRPAks+5Oll5kYEa5j7fIdArl+361fbi6hFHDgNERP4KIqLxae2Qt7tz1PO59ogWk3buBhx6yu0tElR6DdSI9PXoAV28r5FbJuvJ30QD/emTW1ftQz5f393cGJ3r7srrSvGhZvN5AgFKyrbQ1ytwdAIZeicGf7wJv747FYxucz38fchRFDtX2yv5CQrTLldVtd0ujRdquwWK+S/bVm4vXGWWkrQbUVqYOWB1EECnf15uWIFJt4M56BKJrHyhTNVRTN9KrAV8kANnIR0FRAYiIyqOesT0xqs0orCzejbfvaQEsWAAsWmR3t4gqNQbrRCJatACaNgWSk+W2lQDcndJ5vX1aLVU3OpbDYW11eCuZd9cyeHVbbyBACZJc2xpBVFARMPlADfRKldu/xAADW21FwiTgy+hsSOrt1QuDKf1Qgi7XfQe6rDyv8OYCc9dz30ZZZKuZdKvZbiuDCFbK9wHjaQjunhP1VAxA/y4D6kGfq/s6ESZhao9c/Fj1DFBYiJ45dbH9PeA/J3rLpaREROWQw+HA3JvmomnNpni8ySFsSYgEJk4E9u+3u2tElRaDdSIRMTHAgQNAly5yWzQAdy09F3mtUYCtzoaLLAInWgEgMt/dSubdNbMeqApg9AYClIDLNbOuVebu0k48CTxzqCHSqwHD2x5AhwnAd6GZkPwc8rGVAE99yy2tzLpIpt2d7LdIZtcb+7ZaBm82j1w0ay/SbysBtieVDWYDHK6ftd5dBZTvner5E4F5mNoqA3HDMvFGcg4W1zkpT2kIDkbSSXCBOSIq96oGV8XS4UvhcDhw++hgXCy+Is9fd50eRkRlgsE6kTuUf9grwZ9RVtwsKLdSYq/Ohptl1o1K8LX+7ulK83pz1q0OBChtsyBKFXRVzwVe/LMuDr8JTDkWgz/qAINapuCZPo7Sq3Yr71OvDN7sWIDn2W+j+fHemOutdyyrwbhZJt41oC5QlX+7k1m3MqDhzkCAyC3/dJ4/UwV4LPEEGo86hTeanUGLrCD8+4dIzE+JKr2t+j0REZVTbeu3xaz+s3C6OAs7nxgH7NgBTJ1qd7eIKqUA802I6BrDhgFZWUCnTnJbKwh2OK6dF671KJpZV7eVxdG8lVnXeq3Vleb1Musi5ftabfVcYcA8oM7NRZ3LwBsnkjD1izS8Orwu7th7AagfgMLiQqysdx5DTtaAv2qusW6Apm7n5pYunb7emXVvzv329gJzVrLfZgG1aNbetS1J5rekKyx0DsqoB2by8szL3jXaF4OBN2OOo81Zfzy3py5uPVUDjuzLQHh+6W2ZWSeiCmJSh0kY0mIIGobVA35JBd57D+jZU86yE1GZYWadyB116wJPPOHMrGsFwXoBqjtz1fXaopl1kX17stK81cy6WbCuFfhbCLIaXQLmbG+IhAtBgL8/lsRkYXi3Y2g58AA+SgYKigtKz2E3K7lXB+/qoNXhsJ79Npuz7slcb6OA2hsLzBkNOqiz+GbnwGygwOz2amaVDuqBGNeBGYFM+sGawMS2R/FhO7nd5IIDm3d0RMqnVTAkLRSOgMDS1R/qBeeIiCoAh8OBhlUbAv7+OPTuP/Fn85rA/fcDhw7Z3TWiSoXBOpE31K0r/2NfuS+7Uaba04Xm9OaZG7VFsvaeZNb1Amy9fpmVwWsF/vkuWUyBTLty7IGpgZi+vRpOBRfh3sH5aJq4Dm91KEZ2lQDjfWm1JclZ8u16T3fAvNzcakbak2y4p1l70RJ712O5MzAgeo603pfyWvXK/8r3QG/Kg8ZAzZYGwIiYzYh/CJjX6CR+bHz1eX9/tMurCb+iYmdVi3pdBeU7zGCdiCqYrLwsdFp5M24bVwWX8i5x/jpRGWOwTuQN06cDhw8D4eFy2yib7G45vNHidWbt651Z1wuwRV5rJUuvbosE2Ff7FXm5GC9srYq0tUl4+bcQ5PgV4ZFBwMz40+L7ArQz8UqQq/TZLItsNbNuVmLvur3ePdy9kVl3d9DBk/J9o3ZgYOnXKpUPWpl1g4UKj1QH+naRFyZcXuskBh4Efv6zA75YfnV79dQMZe0DdbCutFkGT0QVTERwBKZ3n479uZkY/2RLSCkpwJQpdneLqNJgsE7kDQEBQIMGzrZR5tmorF3k0agMXmkr8+XNBgSMsvTeCrCtznd3bQcGyu/Hdd8iAbbyPq4GotUQgqe2hSPtpyTM+wZ44FQMAOB8UTbuGQL8Xqeg9L4EVqIvCXLVpdBGwbvRAnNm2WyrJfauxxKZ623U9mTQwZPMurrcXP283r3Q1Zl0ncz6pbAAnAuV25E5QErVyxi9E9jxawJWfwb0yKsPh7K9y/eo5DvKzDoRVQIPd3oYI1qNwHLHHrx5dzwwbx7w6ad2d4uoUmCwTnQ9DB4MDBnibGtl1kXL4o0CbL2A2uxYWvs2C86NAmqzMvjrtdI8YBxQK+fdJSMamlOACVuB+n7ytIXvwk9iYTLQseuf6Hwf8EnEYeQEQH8gQD3PWR2cK9lVdRm2ku12DTQVrm3XgNpooTYri9e5tkX2fT0XmBPJ2ivfG6NKB/Vr1Z+HzvSJg9JZPDIQiOq9Ha/dKO+rah6Qvr0XFq0AEs/4O/etHEvne8TMOhFVBg6HAx/e8iHiI+Mxrelh/NauFvDAA8Aff9jdNaIKj8E60fUwbRrw8cfOttb8cSsZdK2/G5Wqmw0ImM1/19q3UQBtVgZvZb67WVsJolyPpQROyqJ/6oyoEsC5zjW++tpR6VWx/T3g3pMNsLMuMLbpH4iaCrzV7Jy8ndHtvvQy62YVACLBuzcXr3Nn3+6WqrubtVc+Z7Pbq7mef6VfQUHO96VUZOTmIt8fWBqdhX53Ac3arsdbnYHYvFAknXDuKzwwTPtYytoHSt8KCpzfYXXGXwnmiYgqoKrBVfHlyC8R5B+Erx8eIP+/cNgw4NIlu7tGVKExWCcqC0FB+vdk96QsXnTFeZFsvZWA2SjgVpffm/VT2ZeVtt4cduXYrsG6Otut9EUj+Es6CczPaIejrwOz98Si9hUgyC+w5PkfGwM5oc52yeuVwFPpp96CZu4E795eRd215P563bpNNEuvvGer58y1sqGwUP7slTns6s8jLw/fNQVuTz6An2OB4Sci8dPHwI7d3TFqF8QGBtTfQfX3SsEyeCKqBNrUaYOdE3filbs/BWbOBA4cAMaPlxdeJaLrgsE6UVmYNQt47jn5d9G56yILzlnNrBtl8a1kzo2e9/eXAyflHvNar7VScm9UYu/aVp8DJcgSyKyrA7YaucAjB2pizxzgvnOxAIBUx0X0vxuo13cbHrgJ2Bh4ApLyeiUwVWdX3QnWXduui6cB3l1F3eoCc+ryfK19m926TdmXco7U1QcFBfI/9tw5Z+ppCP7+OO/IwweNL6B790PY2MgB5OVh8AFg9pEWODoLWPZ7LHqlAo6Q0NL7MppOof6Oqu4yUIJl8ERUSTSp2QQOhwMFf3sQ343rBixfDrz5pt3dIqqwGKwTlYURI4C+feXf7cisi+zbypx0o7bWsVz3bVYW725mXS+IkiQ52HTNrBvcZ90BICBYDuiqXy7C7G+BxrmheP8GoGubTWj+EPBq89Ol5zEri5BZvK+3bltr8TR1ubgnZfBWS9WNFsZT79vfX56jr7Xom5L9Fs2kCy7wl1eUh69aAsO6ZqLe0IOYcONZ/F79MnbWk4P1gGLgkfPxqHPZZF96xxZdFJGZdSKqJB75fgoGxazHtz0aylP/fvvN7i4RVUgM1onKWmioHBjUqCG3rc5dV7f1MtJG5e5afxdZYE4006716ElZvLuZddd2QYHzWBaz3dUvF+GRTcC2vT2xbS4wJa0+LoQAa2tllwTQ/6tXiMM14H5m3ei2cMr2rkGu3sCAa0m9UuVgVqqulcV3fa3rvtXTDvTep7Jv10oHD6oPSiobJAmPRe3GbUPz8VWDi/jLyRAs+LkaTqzrgAf+DHG/skEvs66cA602g3UiqiT+fuPfUSOkBsYMyEZarQDg9tuBU6fs7hZRhcNgnaishYQAKSnAiy/KbbNyd6OgXS9jLRr4qx8dDv3g3SyAVm/vzkCAWdbe08y6az+tBJou7eQTwBs76uPo68CCI4klQe2krufR5M7T6NBxB2a0u4JdkUVyQGmWSXd9Xt02uy2Z2cCA3rFdV1FXt9Wr2IsE1CLZb9dV1N2co36lSiC+agGMid6CZg8DhSHy5zc2sxZe/TUUaWuTsObb2hj3ZxCqIcTawIze+wDEMussgyeiSiS2eiw+HfopzuVfxPBH6iPvxFHgzjv5/0EiL2OwTmSHVq2AiAj5d08CbHXA4OdnPfB3J2uv97xoZt3qgnKugZFZsC7Sz4KC0tu7EYgGwg8NEAEUFkIqLMCLO2vh7oNh2B+Wg6e75SGh+y40fgTYFnDadF+6bSXIdTic5ePqANqs5F5krrfSdt2Xt+be683r1+q3co5UgwwFfsC71Q/gllFA7X7bcNsdwGe1j6NmDnA8wgEAuOFMEP6+PQyNCsP01yfQ2Ldmv11u9QZALLPOMngiqmT+Gv9XPP2Xp7El7wgefqwlsGYNMH263d0iqlAYrBPZLTYWqFULaN1abpuVv6v/7nA4V193ONwP/LX+bjbPXC940eqn1Uy6lcy6+lZuZoG+XoWAO/PKVXPUHUXF+OupavhkYz2cWt0a3y8OwEPHGiGwGIi7JB9/R2QBbh4FzK6xDzvrAsUhbga5rtlvrYBacK63ZtZeqTbwZCDAtS/qFdqVoFbdb9X5zgoCvgvNxLZ6cjugGHip4SGsbgbckF0Vb3wHpP3cFps/ABo5qjlfrw6Y1WsIGFUXaPVbYRSs6/33wIwSEVUyL/R6Af0a98OaBnk437UdMGMG8OWXdneLqMIIMN+EiK6revXkeV4OOUsovOCca9Cu/pvoInVGAwJWAmx1UKx1DHcDf6223n3W9fopUiGgZFP1AjitgE8doCmBZ3Exggsk9D8ooX+zJOA/GUBnOYD7o3o+vmsK/Md/H/AgUDt/O3pFA0ODUnGH676VtlaQ67qKenGx+CCD6y3P1Jl1pR0YKA8AKceuUkV7X1ba589rl8GHhOBCCLC+SgZ+6Qesa7MHKYOAIr8UjMsFFmQHwwHgm63N0WTdH6jW90bgf6uAZMn481Bn1i9fNu6nesBC4XA4By0A7akcWt8jZtaJqJLx9/PHF8O/gAMO1Bh+BWjfHrjnHqB5c6BNG7u7R1TuMbNO5AuUYBtwLwtuFoyblcWbldibZRbV91Z3p6TebD6wlQXmjDLrWu9Da9+ii44pr83Pdx5bp8R7zKVYnH8VWL2lOR7/DYjKD8Gy1sC/qx8v2ffrXYD/q7sL38QDJyP8St+WzHVgQDmuSL+1BiGUFfKVtuu96JWg12DF/FLHUgYCtMrNr+4rN8gPm+sWYkntkyXn6M1OwC0ttmHmjcChKrm4aT/w+t4YTPvNue92ZwJRLU/w83DNrJt9lgqjOepm33+WwRNRJVYztCZqhNYAGjbEj/OewLHAXGDIEHmglog84vPB+oIFC5CQkICAgADMmTPHcNtNmzYhOTkZ8fHx6NOnD44fP15GvSTyoh495Fu9dekit0XmnVvNqJsF82ZBrTpQ9fe/dnE612OIzncXabtmQAMDxV9r1he9AE703tv+/s4AOCBAHsRQvTY8Hxh0tAr+9SOQsqc7Tr8GvLw3quT5xQnAy1FHcMudQL0+Kaj7aBH6tdiMdTEoGQjIkfKvXRjPSkWAent1sO66CJxBqbrhvlXts6HAv5qfwb19stGu4zZE3HcKnW45ibuSjyBPKgCCg3HbHmDO7jjsfBc4va0f/v0FMHV/JFqesXasawY0lLZ6FXs/P3mQQv1ahVnZu9kAE8vgTV25cgWjRo1C06ZNER8fjxUrVuhu63A4kJiYiOTkZCQnJ2P9+vUlzx04cABdu3ZFfHw8OnbsiN27d5dF94nIxL4z+zBwx+MY9ng08lIPccE5Ii/w+TL49u3bY+nSpZgxY4bhdpIkYfTo0fjwww/Rs2dPzJw5E1OnTsXixYvLqKdEXhIbCyxd6myLZNaVzLxIEC7yaCXAFimxNwp8ROaZq9vq8mSzLL3VCgG9+fDqVdMBOQAMD3fuW12GrQSWrpldjYA6MgeIvOAoef5/HwK7+iTgf+f/wJYODfFHwVFsqH8R+YFhQKE/rhTnIeLuo6hTcApNEoEmIVvQpAfQOPIShu0EQvPyUOwAHMFBcLgcy/XY17Szs6+9vZrJAnPb6gFHoi8jPRhIq78D6SOBw41248NMoH1eHgr9gL83TwcA1M91YFBmKG64UAU3ZEXAUXgCCApCwikg4WAN4NQRINgki280Fz8iQj7nBQWly+CVxQSVttYCcq5tdwaYmFk3NXPmTAQHB+PgwYM4cuQIunTpgl69eqGGchtLFxs2bEB4ePg1f3/ggQcwYcIEjB07FsuXL8f48eOxcePG6919IjLRvFZzTO08FTM3zsTDU5pj3uvfyQvOvfyy3V0jKrd8PlhPSkoCAPj5GRcBbNmyBcHBwejZsycA+WJep04dFBQUIFCdeSMqb0JD5UclUNEKipX/Pqxm0o3+blS6rjVH3WhAwJ155Gbbq8vD1f0yaxsF80YBnLqdmwtUr+7ct5LZLS4unclVZ3YF5pEHFAPJp/2RvB1A/XbAN0dR3LE9pBNpQL0AXCrIwu0Hg3GobhD21irAb1WOA3UA4CBuWSUH62vigJsHbEG9rkBd//WoWhOoGn4EMZHArKvH/q1uHn75C+AXuQd+NwJ+jc7A7/xFNCssxk3HwoGAAHxUKwObBxQgp0E6LjQFztT/FWfigL+dPYOHNsvvY9jtwJEa+6+enHT41QOi8otwtor8fJ3LwK+7OqLFDymIbNwauHRJft8OB5BfJF6qLlLer+zLdaV/9WdtFKyr10Jw53vFzLqpJUuW4OOPPwYAxMXFoXv37li5ciXGjh0rvI9Tp04hJSUFP/zwAwBg2LBhmDx5MlJTUxEbG+v9ThORJTP6zsD2k9vx/uH/ov2wWEyYMQNo1w4YPtzurhGVSz4frItKT09HTExMSTsiIgIRERE4fvw4oqOjS207a9YszJo1q6SdnZ1dZv0ksuz++4HoaDnjDmgH1lYXpzN73soicGaBv978XpF9a22vLncWWfjLykCA0eJ1rm31sbWCdeV5kRJug7ZfXj7gJ7+PeueBz78PBxISgE2bcLFdKxw+tBWHeyUhIm8HkJeH8Hygz8WaOJF7Fkdr5mJPFJAddBbN/YBZf8h9WVv7MqZ3BIADQAwAyNnvWzOLcNPJaoC/P36segaL4wHgFALqApGFV1ArD/D3Dyzp5/NrASQkIObXPxAT1xYNftmGwIG9gUOrgRp5cAC4Mb8ekOVSmp6fb3x7NZFzpJ4Coa5sUMrczb53VtY+EP1eMbNuyPU6HRsbi/T0dN3te/bsiYKCAvTp0wcvvvgiwsLCkJGRgQYNGiDg6ufjcDgQHR2N9PT0a4J1XuuJyl6AXwC+GPYFOnzQAZOTMtHmQA10HTsWaNGCC84RucH2OevdunVDrVq1NH8yMjIs7cuhXqQLcmm8lqlTpyIzM7PkR6vMjshnREUB48cbl7qbBeGiZfHqv5uVj7ubWRcJ3kXL5K1mQM2OJZpZdz22euBANLNuJTBVD1Cogt5ql4vQ9gQw7HJMSdl7l0xg1eEu2Po+kPlTMi7NAAo3D0TKPOe+7zsXi63zgN83tMGmT4KwcWtb/La8GmZsrVGy73d2xeD0a0DWlv7IfxE4saYtdr0LTDrftKRfd+0E7sqKQ/c0IOaChMBinfclSXLGWzlHyir2yufhzoCG+rNVLyin/nyMvkdmn63Zd1DrO1nJg3WR67n6Oq13jQaAtLQ0bNmyBRs2bMDp06cxbdq0kud4rSfybZFVIvHV7V8hNDAUu5++X/5/NhecI3KL7Zl19aIxnoiOjkZqampJOysrC1lZWahfv75X9k/kM2rVAiIjSwfvShm8u5l1re3MysfNBgbUbatl7qJl8iKZddEKAcA8WHedL++6r+Lia+dIX7lSel9WS7617iHuOjCgtQK7al9+wSEILXS26/lXQ73jAOoFAqcCgcs1gEwHUM0PqCX3u0Z2IXAFQGCYfr/M3of6LgH5+fI8f/UAht6Ahtm+8/Pl92wloHb9bhhVUYhk1rW+R5W8DN7seq5cp2vXrg1ADsgHDx6suy0AhIWFYdKkSZgwYQIAoFGjRsjMzERhYSECAgIgSRIyMjKuqaAjInsl1UvCkUeOoGZoTeBcLDBpEjBqFLBqVen/nxKRIdsz697Svn175Obm4ueffwYAzJs3D0OGDOF8dap43nwTSElxtq1kua0E76IBttm+AgLkgQX1gIJRdttKJt6TQF/r2EbzmM3mv+u9D3dvr6Ze8V09J1o9V9+dud2uz6sDzSKBeeSiFQHKIINyztTvQz1VwEr1gdHnYfZd8CSzbva9YWZdyIgRI/DOO+8AAI4cOYJ169bhlltuuWa78+fP48rVQa7i4mIsWbIEbdu2BQDUqVMHbdu2xaJFiwAAX375JWJjYzlfncgH1QytCQA4c/dw/GtKB0jffw888YTNvSIqX3w+WF+0aBGioqKwbNkyTJ8+HVFRUdi2bRsAYO7cuXj22WcByAvQLVq0CI888gji4+OxatUqvP7663Z2nej6CA+X57ArRLLcon/Xy6wblY+LDgQY7dtKRtS17N8siLISvFu5nZfVtsACcyXPq/utDqit3kNcJPutHMv11m1WbwNnVqqurhBQZ9ZF3ofrOVRea/bZeuuzFh0YqOSZdTPTpk1DTk4OmjZtigEDBuCdd95BzZryP+bV1/O9e/eic+fOSEpKQkJCAs6ePYvZs2eX7GfevHmYN28e4uPj8corr2D+/Pm2vB8iEvPM2un4e/XfMXN0LPD66wD/myUSFmC+ib3GjBmDMWPGaD43ceLEUu0uXbpgx44dZdEtIt/RqZPzdyvZbrO/q+/AoA5WlCy56ACB+vn8fOPAX9lONKspUgFgJXhXl7kDxqXRWscuLtZ/Xq/sWivzrpWR9vd3Zr/NMtJW2so5UO/b399Zvu96T3d39u36vswy667nSKttNjBj9nlZqaIQ+c4qn39xcen/dqhEWFgYlixZovmc+nrepUsX7Ny5U3c/zZs3563aiMqR1/q9hl/Tf8UT2I3m3erglgcfBJo2BXr0sLtrRD6P/6IgKu8WLpR/APNMulmQrn69w1H6eaMg3OgYWvu2klk3GijwZpZe63lPSqPNMrlGwbtekJufLy/WprQ9XGm+VKm6XmY9MFA+5+4E664VAso5UgYoXM+pldX4Pc2sGw3MiHxP9D57ZteJiEqpGlwV34z6BpFVIjF6QDZ2NgwAhg0DDh2yu2tEPo/BOlFF4mmQLhJgi2bvjYJ3s2y4XhCmdSwrWXqr85w9WXTMLJMreuswrbJ3Zd9Wy971Su79/Z2L4ynHVrZ1Pba7mXX1QIDW+7R6jsw+W0+mQFj9nqjbDNaJiK4RVyMOK0auQF5xAW6eEI6TeWeBm28GLl60u2tEPo3BOlFF4nAATz0FjBsnt82CdCuBtpVyd7PtzDLSekGY1vuwGlQZBfdGpesiwaGVzLo7gahyiypl38rK8964LRxQOquvt2K73srzSrZceV65NZuVANqs+sBKqbrZsdWfl/qe7Vr7MvueqNtcZI6ISFO3mG54/+b30bRhAgL/8QKwZw9wxx38/yaRgQDzTYioXHn5ZefvZiXrosG7OvAxe43IdlaCKvX2RoMIyvMic4tF21bLrtVz1l23Nyq7dm2LBKLqtrJvJZhX9qX8A0hpK8G9EvTm55c+lrqtvgWdvz+Qk1N6X67Hcud9aD3v5+c8j54sAif62buW6+v1S7TNf3QSEekamzwWdyfdDT84gD2pwIIFwOOPA6pFJInIiZl1ooos7Oo9sqtUkR9Fs99az4tkzLX2ofV3K5l19fZa+7QSVIkcSy/TqzUQYCVrbyWLLJLF9yRDbZRFNhtM8aRCQOTzET2W1YEYo4EcT17r2m+WwRMRGfJz+EECMGtcc7xwd4x8S9q337a7W0Q+iZl1ooosNhZYutS54qpo9lsrwBbNxovs291yZbPMuqfBoVGwaHVfes+W14sFAAAgAElEQVQrwZynt4nTep/u7Fu031r7AowXhRPJdmu9D9cV8kXehyffBSsDRmbHZmadiMhUflE+Ptu9BCmN0xAzoC7ueeQR+ba0t95qd9eIfAoz60QV3YgRQJ068u8iGXS9v5sF51aCd3eDLK3BBGWFeNd+Km0r2XB123Ues9WsvNmghDfLx12PZTUjbeUceTOzbvV9qKcSmH22VsriPR3AYGadiMiS4IBg/GfUfxBdLRr3dT2LNclVgVGjgM2b7e4akU9hsE5UmTRsWPpRNFvuWjKs9VqzMnh3AyN1wCcyMGCWERUN8NwJNM0CPPX+jDLSIgG0Udts3+5k7bX2pdVWB9Qin7XRsdTvw8/v2vJ9K/sy+jxEBjBEB5SYWSciElI/oj5W37kaYUFhuG1YIXZFFskrxB85YnfXiHwGg3WiyqRTJ+DcOaBvX7ltlLFWHpXstVn23UqW3iwI1gtUzQYErGRTlbbe81YGEdx5Xh3kqu8dr3dsK8Gjlay9SIbak31bGSwxeh+eZumNPk9PM+vqvjFYJyIS1rpOa6y4fQVyivNx66QayD97Chg0SP63ChFxzjpRpVOjhvN3kcy6aMbcSom9aEDt51d6oEBkQMBKdtssi28UIGsFbOoSaKuBqNWg1ih49HTRN9Fsd1m+D3c+S9G+iAwEmGXele1ZBk9EZEnvuN74eMjHqFWlFoJCdwOPPgoMHQr88MO11zOiSobBOlFl1qsX0Ls30KKF3DbKVLuTQdd63t/fOcdc+btekGWlTF9pi5Yru74/17bVzLrVLLKnAbSVUnXXleZFz4FW2+gWdCIDGmYVAnrvw52KDNFzqlSPKLeN82QggJl1IiLL7ky4U/5lSn+cSv0TNd7+EIH33AN8/nnpfzMQVTL89hNVZp07A2vWAFWrym2jTLVWibx6WysBteujXpAlOgCgF2CLlMXrtT3Jymsd2yiz7k5putH2niwwZ/a+AgNL/8PJ6grtZpl10eoDqwG10TnVenR3QIOZdSIitx3LOobOcWtw70ONULx0CfDII4Ak2d0tItswWCcip9BQOcOoZE+NAme9INwsqDfK3ivPe5K9txqIGmXWjSoARLLIRgGeJ/cjt5pZt3rrNrOBAvXz6qy9NwZLjDLrVs6Jle21vsOigzqu+2ZmnYjIbXXC6iCxbiIW1cjAY/dGQZozB3jpJbu7RWQbButE5PT3v8slZ1WqyG11EKK+PZr60Sg7qX40ep1eZt1sIMBqZl0kOytaASASQLsbiFoNoI2y9mbHsrpvo/dh9RyIBPPK987qvkQGHYwGhawMvKj7zWCdiMhtAX4BWDxsMbpFd8PsRpl4dWR9YPp0YN48u7tGZAsG60Tk1KwZcMcdzrY6CHE4jEuHPSmLv56ZdbNgUPTYIplaK8Gj6zxyoyDXk6y9a9udQQejY1vJfntyDkUCeyuZd/XzWt8v0UEd1+1ZBk9E5JHQwFB8PeprJNZNxFOtjmP+gNrAgw8Cy5fb3TWiMsdgnYj0xcQA0dHOtkggbbUs3igQMgvStQYCrASersf28zM+ttFrrWZ21Qu1Wc2sm72voCB5cEXrWCJZYivZb6v3h9erojDb3upnabZv9Wet9V11tyKAmXUiIo9VD6mO70Z/h7jqcfj38ARIdesAo0cDP/1kd9eoEku7kAapjNdQYLBORPoWLwZ+/tnZFilR9zSzblSeLLJPs9Jno+yq1m3iRAM4q4GqayBqZV9m78tqRtpqWy+z7mmW3iz77UmW3vV8qx9FBmZEj8XMOhGRV9SPqI9fxv2CFeO+g+O77+WKtCFDgK1b7e4aVUKbj25Gm/fa4LmPxwLr1pXZcRmsE5G+wEB50TmFVmZdNNNuFCiZBetmwbto0GsUZIkMNngzs65+vUg/rQTQZgMB3syse7LAnJUMtsiAhtnnIRKk6/VLtN/MrBMReU1U1SgE+gdCSkzEjDm3Y2PtPKB/f+CPP+zuGlUiu07twqDPBqGosAB9X1kqTxnNySmTYzNYJyJxjz8O3H+//LtRwGz0qLWd8jc/P+fceK1tPc2sGwVdIsewGlDrBYvKPb2tZL+9mVm32m8rmXWrgwqiwbzZ4IdIMC8SpOv1S7TfDNaJiLwu7WIa/pnxGQaNC0RK2CWgb19g3z67u0WVwKFzh9D/0/7Iyr2EFV8GoFumH/Dll6WTWdcRg3UiEjd9unMBOneCcr3n3d2XVlCvBPxK22q21egYVoJYo+ys1YEAd0q+jfbtSXCu7ovDIVdf6PXTk9J1swENd86ZaBm8v7/x98jo2CyDJyLyutjqsfj6jq+Ri0L0nxCKXX5ngD59gMOH7e4aVWCpF1LRe2FvnMw+ic9Wh2LgngJg5Uqga9cy6wODdSJyT0SEHKDUqCG3zbLhIiXH7gT6Zo+imXWRSgHRsnatY2tldvUGAkQCU9G2SFbeStZevW9Ps/RGFQFan7W7gyOu+xat3BD9fJhZJyK67vo07oOvbv8Kl6Rc9P1bBPblHwN69wYyMuzuGlVQRy8dxfnLZ7FgTRhGpOTKGfW+fcu0DwzWicg9VasC27YBTz8ttz3JhosES1b2qX60M7Mumtk1C77dCbC1BkNEjiUScKv77ednvK27peta3wV3s/Su+xbJrLtuZ7RvdV8YrBMRXTeDmg3CkuFLcKY4Gzc9UgcFGWlywH78uN1dowroxiuROPxBFdyzMQdYtgy46aYy7wODdSJyX5s2QFiY/LtokO7NzLpIht1qcGi0b0+y32ZZY9F+6rWNAlPRQFO0rT5HrmsMWDlHIpl1KwMaogMcopUZell99bGUcnnlNSyDJyK6roa2HIrPbvsM8+78HIHvzgUOHpSznSdP2t01qgCOZR1Dr096Yf/WH4DevVEr8xywZAlw66229CfAfBMiIgGJicDAgUDPnnLbrNRYKxASzZiLZkZFgkV3AzaR4N2TrLGVYxkFqlrPuRv4uz6vPkZRkWdZetfttT5LZR57cbG1823Ub61Hre30juW6L2bWiYiuu9vb3C7/8kBvHMrOgOOll9C4Vy9gzRqgfn17O0fl1vGs43KgfnY/fvtkO+JPZcm3Mb7tNtv6xMw6EXlHZCTw7bdAQoLcFs2sBwTIAZifn3jG3GwgQC8baxTkigTtVgJRkayxlYGAssysu3vO3MnS6/VN5HMRzaS7bi/6/RI5luu2zKwTEZWZ3MJc9A34DD2mVMfBk3uAHj2AzEy7u0Xl0Mnsk+izsA/2n92PtzdUx7h1l4BFi4ARI2ztF4N1Iro+atYs/SiS1bQSPIk+6gVsWquNaz26DiooK4RbzYYbBaKeZtbNssjuvtbseaMqCdF9WRnQcD2mt/dt9B0VrQBgZp2IqMyEBITgn73+iWO4hO5TqmHv+QNywJ6WZnfXqBzJuJiB7h93x54ze/DGxmqY/ONFYOFC5x2QbOTzwfqCBQuQkJCAgIAAzJkzx3Bbh8OBxMREJCcnIzk5GevXry+jXhLRNUaOlBega99ebotkw80y5aKl9epHq2Xvoll8T7LhZllikQDb6HnXQQjXReDcDc5djyWSWTfKfhsdW6Rqwso5MTr/IgM1esdyfWSwruvKlSsYNWoUmjZtivj4eKxYsUJzu927d5dcx5OTkxEbG4uayqAfgNjYWLRo0aLk+SVLlpTVWyAiHzQ6cTQWD1uMU1I2ej5cFX9mH5YD9iNH7O4alRNv/O8N7D+7H3PWhWPKmsvA558Do0fb3S0A8P056+3bt8fSpUsxY8YMoe03bNiA8PDw69wrIjIVEAAkJ5duqx+VIFIrUy0SqBn9XSSzLhLo6x2zoEA7GFQHakbBo7cz0mYDHcoicEVF3ilVdyeo1eq3p4MpRvsSOf9WBmz0jsUyeGEzZ85EcHAwDh48iCNHjqBLly7o1asXaii3f7yqVatW2L59e0l78uTJcCj/n7hq+fLlaNOmTZn0m4h838jWIxHoF4jbl9+OnpMjsHF2Gpp27w789BPQrJnd3SMf92rVYRi04n3025sv357tllvs7lIJn8+sJyUloWXLlvDz8/muEpGRO+8EJk8GGjSQ2yLZStEg3SgANyvhFs2siwSLeu9HKbk3ej+uq6pbybSLDHQYDQy4m2nXOqbVfZsNBBg9iu7LdcqD6ECN+n3p3aLOdV/MrOtasmQJ/va3vwEA4uLi0L17d6xcudLwNXl5efj8888xfvz4sugiEZVjQ1sOxYrbV6BXi4GIfusT4NgxoHt3YOdOu7tGPmjb8W14du2zkNauRWD/geh3UAJWrfKpQB0oB8G6VT179kRSUhKmTp2Ky5cva24za9YsREVFlfxkZ2eXcS+JKqHOnYG333YGPVoBqEhWW+vvRs/rHcMbWXzRrLE7gw56/dY6lpWBDpHA32gQwuh9eZq19zSzbjaAYfZ5GH2v1LdnMxocYWZdV3p6OmJiYkrasbGxSE9PN3zNihUrEBcXh2R1lQ6A0aNHIyEhAffddx9Onz6t+Vpe64kqn5vib8LSEUsRNPpuSEuWYJf/GTlg59RYUtmYsRG9PumF19a/gr13DZT/bfr99/ItAH2M7cF6t27dUKtWLc2fjIwMS/tKS0vDli1bsGHDBpw+fRrTpk3T3G7q1KnIzMws+WHZPJEN4uPlH4VWwCkaQIsGvUaZdZGMtOvfRTPrVvqrtW+zAFsduJqdQ6PXah3LKCgWyaybHcus8sFKZl1v3yLnxOh5o/fp+ppKnFkXuZ6ry9klSTLd54IFC67Jqv/yyy/YsWMHUlJSEBkZiXvuuUfztbzWE1Vur9U7hLYTJHwWnwv07w98/bXdXSIf8N3B79Dv034oyM/FqkXFaJkTLk+X+Mtf7O6apgDzTa4vby4CFx0dDQAICwvDpEmTMGHCBK/tm4i8bM6c0oGNVsZWNIAWDXrNAk2jR5HMut77Ed2nu5l1q1lks0EGvQBaWWPA6P1YzdIrbbOpAu5UCLgbnBv1IS/P+HxX4mDd7HoeHR2N1NRU1K5dG4A8wD548GDd7dPS0rBhwwYsW7bsmv0AQGBgIKZMmYJ49aAfEdFVtzS/BXN+n4Mxf83E6aphmHLbbcAHHwDjxtndNbLJx9s/xv3f3I+I4kD8Z34euhY1ANZ9D/jwGii2Z9a95fz587hy5QoAoLi4GEuWLEHbtm1t7hUR6XI4gMBAZzsiQv4BvJtZVy9iZ5QRFS2VNsrsupvF13o0y1BbzazrBdTeDnr15nbrvS/R0nStcyi6yr3VfYu8X73XsAxe14gRI/DOO+8AAI4cOYJ169bhFoO5gR999BGGDh2K6tWrl/zt8uXLuHDhQkl78eLFvNYTkaaWtVtiw70b0KJWCzza7TKevjUc0r33Aq++CghU9lDFcubKGUz5bgoaFITitzk56BrRCti40acDdaAcBOuLFi1CVFQUli1bhunTpyMqKgrbtm0DAMydOxfPPvssAGDv3r3o3LkzkpKSkJCQgLNnz2L27Nl2dp2IrPjiC+D99+Xf3clyG22nlbH2JKOuPIpmqK0OOhgFzGZtkWoE0VJ0vfdhZdDE06kCRudQvSif1qCDu4M8IgM2zKxbNm3aNOTk5KBp06YYMGAA3nnnnZJbsqmv54BcIv/xxx9fUwJ/8uRJ9OrVC4mJiUhISMC6deuwcOHCMn0fRFR+NKrWCL+O+xWdGnbCjMSLeGx0JPDkk8Cjj3JwtZKp5ReO1bvbYuPsLLRs3QP49VfgaqWWLwsw38ReY8aMwZgxYzSfmzhxYsnvXbp0wU6u9khUfpnd5k396E5gXVioXWbtbkCtFdS6lsGL7lOrv6IZade2yDHNAn9PSuzV58KsIkDkHIkMpii3pDPLrHtyDNe23rb8x5+usLAw3Xuiq6/ngDy3PTU19ZrtGjduXDJgT0QkIrJKJNbcvQZjvhqDocPvBzJeBd58Ezh0SL6ftlLVRxXO5fzLmLhqIp5s8yBaj38SXdevB0aOBBYuBIKD7e6eEJ8P1omoEgoNBWbMkFdwBTzLviptkQDM6Bhar3c3sy5ybKNst0ipuh2Zda221UEGs89Y6/3k57vXb6ufvVJ6bzS3npl1IiKfExYUhq9u/0pu/NAH2x8agSorvkH8X/4CfPNNuciwkjXHs45jyJIh2Hx0M+ovWY3X1p+TKypmziw9jc7HlZ+eElHl8uSTQNeu8u+ezGtW2lYCZZF9m2W/jQJokWMbZbuNAmyRc+WtzLqVDLRZ2xvZb9F9i34XjL5HevtmsE5E5NPy/IFbW+1A54dCsfbSTqBjR2DzZru7RV70+9HfccMHN2Dz0c34v80heHXZeeCNN4BZs8pVoA4wWCei8qBFC/m2K717y21vZNY9ya4qj2YZadFKAK1H0X2LrqLuTr9d2yJZe9e20YCGSGbdyucjum+rx9L6HikXe9dtWAZPROTTggOCMe+meSgKDkT/sf6YH3sO6NEDWLrU7q6RFyzauQjdPuqGc1mn8fkKB/75WzAcq1YDU6bY3TW3MFgnIt8XHg58/70z024166oVcFrNrmo9KoGy8jrltVZvQ6YVzJtluz2pFFCvkC+StfekssFoqoDWIIonlQ9mmXWrWXyjAQFlcTtm1omIyp2BTQdiw70bEFW9Ee4bVIApf/VDwajbgWef5aBrObcx7VfUyvPH+vcLMCovHti0CRg40O5uuY3BOhGVP4mJQLNmQOvWclukzNndMmsrWXFPssKuj6LZb3eCWr1+K8+LBtQi59Issy46SOJOZl3kvXv6fphZJyIql1rXaY1N921Ct+hueDPhChYMiQZefBEYPBg4c8bu7pEFF3MvIq8wDzhxArNf2Y6UN67ghsSBwP/+BzRvbnf3PBJgvgkRkY9JTAT273e2RQJuq0G41ayrOsvqaQDt7+/MzhcXGwfU7mSHlXZBgdjAgDvnwvWcWN23O4MretUIokG7yDnTq3CoUwc4cgSoVg1ERFQ+1AmrgzV3r8H8bfNx37S7gFpTIb3/Phzt2gHLl8vz2cmn7TixA8OXDUeP4Ob44PkUBB47jjrTpskLFauv2eUUM+tEVP4lJwONGgFRUXJbK9hyN4NrJXj0RpCut53ofGyrQa5Iyb07VQbKo2jW3t1Se633oZT5iwTfovvWej9RUaW/c7GxQI0aICKi8iPQPxATb5gI/yphkObOxejXOuFfccch/eVG4L33AEmyu4ukQZIkfLD1A3T6sBPSzh1B0kergQsXgM8+A157rUIE6gAz60RUEQwfDgwb5pyHbRTselJmbfY6b+9TL0g0qhSwmi02Ku1WsvsiAb9eW2+QwbVtdSBD6/24Lvrm6b7NBk22bAFCQkBERBXDuZxz2BGWhcW9C7GhaRV8NHUSqm/YIAft4eF2d4+uupx/GQ+uehCf7vwUMXmhWLqwCB2rtQJ+Xwa0amV397yKmXUiqhiUQB0AunSRf4wy61bL30VeZ6VkW2SfIpl1syBW5P3olXZ7mvVWl/MrbfXAgHogwBvZb9dF36wG7WafT1ycnD1XREQAgYEgIqKKIbJKJDbdtwmj2ozCv6OvIOnRUPyyfhHQrh3w++92d4+uGrdyHD7d+SluSgtByuwcdOx9l3z7vQoWqAMM1omoInrsMWDDBu0Vz0UDaasBnZUyeCsl6q77dm1fz8y6t0rt1Y/uvHejttbfrZ4D0QGOL74AfvoJRERUcYUHheOz2z7DBzd/gDPhDvQc58BH1Q7Ld6SZMYMLidpEUqYj5OfjuS3hmPkDsHKpH2rOmQ988gkQFmZvB68TButEVPENHSr/AKWzvCKPVoJJT+dEGx1bL/vtyRxv130pz3s6EGD06G4FgJXg3upnqffYti3QsiUQGur8e3AwiIioYnM4HLiv3X1ImZCC3nG90eftVUB8PPD000CfPvKColRmUi+kotcnvfDVD28CnTuj1Wsf4bHCjvDbvgO4997S1ZUVDIN1Iqr4XnxRXmwEcC485q1SdXcy6+4EwaKZdXeCX/W+jQYCvBEEu1vJYOXYVrP1eo+PPQb8+WfpwQwiIqo0mtdqjv/e/V9Edx4AbNmC+Y/1xquF61CY2AZ45x35ji103RRLxXhn8zto824brEtbh60zpwI7dgDTpwO//go0bWp3F6+7APNNiIgqmH/+E0hKkn8XnSstEmS6W7rtTkDtSebZtW2WWXd3QMOoQsDTgQyjY7ubre/bF/jtN6BFC+f5qMCj9UREJK44JBhzm1/Clgjgy+RCfPT8ZLRetgyYPx9o0sTu7lU4B88dxH1f34d1aesQcyUIXy0H+oW0ADYsADp1srt7ZYaZdSKqfJ54Ahg4UP5dWd1VeXQnsyu6aron2XDlOeW2ZN7IQCvPm2XW3Z0qYFQhYGUOutVjigbnSkm78tipE7B6NVf8JSKia/g5/PDL2F8wres0bK1TiLaT/PB//utwJbk18NJLQF6e3V2sMHILc3Hj/K5Yl7YOk7Y48Mfbheh3x9NASkqlCtQBButEVNk1bCgHaA8/LLfdKVW3mml2Jxvubvm4SFZfGQAQPZYnfbkeUwRc/y46kNGnD/D668CgQSAiIjITGhiK1/q9hg33bkCregl4uTvQZ6wD0jPPyBV7a9bY3cVy7XjWcUCSEPLNt3jt2yKs/Rh452gyItZvlgdEKuG6MSyDJyJSB2tNmsi35GrWTG6LBInqrLfWa7xZPu5O5lz02IWFYoMPVgJn10eRsnyRv4sMbLgev3Zt+fcGDeR2WBgwdSqIiIis6BTVCVsmbMGczXPQILQOHA3SgBdfxLmb+6LmzSOBV18tfatPMpR+MR1P/vdJrNr3DfZu6oD6X6/FPdWqAf98G3jwwUq9dgyDdSIitVatgAsXnIF3UJD8qNxP26yM3Ciw9mQ+ttVMszvz5QsLtftvNXi3UqLuSeBv9Rh16wLHjwO1aoGIiMgTAX4BmNJ5itxIAvYP7ojkFQNw/+aleLbtvxE5YYq8eny1avZ21IddyruEV359BW9sfAO5Rbm4aT9Q9PNaYMwY4F//AurVs7uLtmOwTkTkyk81Q2j4cDl4b9dObhsFh4B2ZtqTLLHyqMxTdzfwNzuGUcm9uwMF3qwQsHJsdaXDU08BN96IErVrg4iIyNvy6tVBYqMb8JbfJixsV4Sn1r2Gv7WYj7AnpwMTJjhvA0oAgC92fYGHVz+E0zlnkHjKD//6Huhftyvw/etA5852d89ncM46EZGRRo2A5593ZthbtJDLqNUL0qnLs6zM+baaDbeyQJtoVt+17efn3cXeXLdX5se7G/jrHVsZYGnQAKhfHyVefhn4619BRER0PSXUTcDG8RuxeNhiVKsdhSf6AY3vPo+lH0yRp9jNmQPk5trdTVsVFRfJv2Rno2jlVwg4fRYL/g2k/BCL/q8ul2/HxkC9FAbrRERWTJsGpKU5g/fwcHnus0IdvHs7i2xlIMCdrL6/v3EgbTVLb1RG72llgPIYEiIvOKMskLd6NfD11yAiIiprDocDd7S5A/sm78O7g99FYJ16qDl2IlBcDOmhh3CpVRNg9mwgK8vurpapy/mXMWfzHDSeHYtdLz0CxMXhjulLcWBFQ4x7cB78/9wDDBvG26VqYBk8EZFV6rL3N9+Uy+TVz3maWTcKfkUDf6vPG82L99ZiduoF+KwuMKe3/f/9nzxVQcHbrhERkc2CA4LxYIcHcW/bexHkHwRMeh1r356KIZfex/3fPopHXn8W0aMnAQ89JN+VpoI6lnUMczbPwdzN7+J8/kXUvgwcWfkW2oQ2gv97cxE2bpwz+UGaGKwTEXmidevS7aeeAqKj5d89zaxrbSca+LsTcHtzX3r79vR9DB4MHDgAxMTI7bZt5R8iIiIfExxw9VZjVaqgeMRwxHz7K2Z1/RNvds7CLftexf19Z6J/wlD4T3gA6N279Jo55dw7/3sLj/7wGAqkQjQ7C7y0EbhHSkSV6dOAkSMZpAtisE5E5E2PP+78PTIS6NcP6NtXblvNtGs9XxYLs1ktoRcp+1fP67fyPhIS5NvftGght3v0kH+IiIjKkb6N+2LnpD/ww6EfMPt/b+Dffj/gq5ZFeGr9crzcbznQuDFw333AqFHl8rZvGRcz8NXerzC5en/4LfoMrb59Dze2KcSjmxy4qelg+L32GNCzJ0vdLWKwTkR0vfj7Az/84Gw3ayb/KCvLi2aw1eXj3l6YzUpm3ZOBAtHMemwsUKMG0LKl3O7QAThyBEREROWdw+HAgKYDMKDpAKReSMWCbQsw8q+dgeW/IGfhAgza/zRuvvNpjAhui+hb75aneUVF2d1tXceyjmH57uVYunUhfju9FQDQciHQ7zDQKzoavRo/DLw4Tl6sl9zCYJ2IqKzUrw/s3+9st28vl9G3aiW3tTLZrkGupwuyiQTU3sjSG70PdWa9alU5MFcGMBo0AM6e5cg7ERFVaLHVY/FCrxfkxg2D8eeDt2Lbwr5YF3sFj2MbEvZvw6Bxj2KQ1AQ9OgyHY9BgoEsXIDDQ3o4DQG4uhn7QB/8+twEAUCUfuH0/cPuRKujRcwTw3p1Anz6lq+rILT4/MeLpp59Gy5YtkZSUhI4dO+Knn37S3XbTpk1ITk5GfHw8+vTpg+PHj5dhT4mILGrfHti1yzn/+uabgREjnCPQwcHyj0JrFfXrscCcksn31qr2QUGl/3ExciQwZIhzmz//BCZPdj7PQL1CWrBgARISEhAQEIA5c+YYbmt0PT9w4AC6du2K+Ph4dOzYEbt3777eXSciuu5uiOmCk0+dxYqRK3B369E42aAaXvsLMLFNKhyvvAr06IENravi2+HJOPPs48B338mD29fRhdwL+PHQj3h+1TQMmNUWk6e3l6ei1ayJ+K834PZdwNIfq+P0ibvwxbhVGLrhPILmfwz0789A3Ut8PrPerVs3TJ8+HaGhodixYwd69uyJ48ePIyQkpNR2kiRh9OjR+PDDD9GzZ0/MnDkTU6dOxeLFi23qORGRRT17yocya+gAAA2iSURBVD+KZ54B7rjD2a5bV85EA9ZL062Uwbu2rWbYa9WSA3RlZfYnnwSOHnXu++23S79vBueVQvv27bF06VLMmDHDcDuz6/kDDzyACRMmYOzYsVi+fDnGjx+PjRs3lsVbICK6rkICQjC05VAMbTkUxVIxUo6n4NzlM8DoSGD1arxx+h0sr70DwA7U++/raPM5kJBdBcMcrXFjzF+AJk2Q3aAWqsQ2g190DFC9umHQLEkSLuReQMaFNBzN3INBfvFARgb+c+hbPHjxc2T6Zzv7VgBEHACQEg507YpX+/YDBg4EEhN5Hb+OfD5YHzRoUMnvCQkJKCoqwpkzZxDlMn9jy5YtCA4ORs+r/9B94IEHUKdOHRQUFCDQF8pFiIis6thR/lF8953zohsUJGej+/WT27VryxfLunXltjKgqWTmRUrW1Rd0dcm9sg9l5dY6dZzHBOQVbG+6CUhKktvTpwP33usM1pU+UqWWdPX74Wey2rHR9fz8+fNISUnBD1fXghg2bBgmT56M1NRUxJbDBZmIiPT4OfxwQ4MbnH/o0AHTTw5Ft0NrsH3PWvxxfDt+CzuG//pdQcya7bhxxe8AgFaPApk7gKp5QLVcoGqBH/z8/DF/a0PccDEMJ4IL0LnPEVwKKMKlwGIUqf6XfGEGUC0PqNkIqD0Y6Hk+AO0dDdG1ZhKSE/ohaGR34PPWzJqXIZ8P1tU++ugjNGnS5JpAHQDS09MRo5SSAoiIiEBERASOHz+OaOU2SlfNmjULs2bNKmlnZ2eDiMjn1ahRuq3OUDdsCGRmAvXqye1+/YBZs+RSNEAOpvftcy7cNnAgkJPjDMRvuQW4eNG5v8mTgeRk+fcGDYC33gKUwdP27YGDB4G4OLmdmAh8843ztaGhQNOmnr9fqpSMruenT59GgwYNEHB18MnhcCA6Ohrp6enXBOu81hNRRZNYNxGJdROBro8CAIqKi5B6IRVVp4YDp7IhpaZi8PZ/4nT2SVzARVwMzkFWcB6koiIUOwAUFaFKXjFq5/qjSWEgIhCMmghFI0c1NAqqDf9negDRTdG1USOkNGkiT8tj1txWtgfr3bp1w549ezSf27ZtGxpdnbu5Zs0aPP/88/jxxx919+Vw+TJJkqS53dSpUzF16tSStlbwT0RU7jRo4Pw9PBx49FFnu3174PPPne2xY+UfxQsvlN7Xv/7l/N3hAB56qHS7SRNv9JgqENHruQij6zmv9UREMn8/fzSpefV6HFEXjiZNMLdPH8PXVAXw+/XvGnmJ7cH6+vXrTbdZt24dxo0bh2+++QbNmzfX3CY6Ohqpqakl7aysLGRlZaF+/fre6ioRERHpELmeizC6noeEhCAzMxOFhYUICAiAJEnIyMi4poKOiIioIvD51eB/+eUX3HXXXVi5cmXJfDct7du3R25uLn7++WcAwLx58zBkyBDOVyciIipHjK7nderUQdu2bbFo0SIAwJdffonY2FjOVyciogrJIenVj/mIZs2a4dKlS6Uy5J9++ikSEhIwd+5cHDt2DC9cLd/cuHEjJk6ciJycHDRs2BCLFi1Cw4YNTY8RFRWFzMzM6/YeiIiIrKpo16ZFixbhySefxPnz5xEUFISwsDB88803aNu2raXr+b59+zB27FicPXsWVatWxSeffILWrVubHr+inU8iIir/zK5NPh+slwVewImIyNfw2uRdPJ9ERORrzK5NPl8GT0RERERERFTZMFgnIiIiIiIi8jEM1omIiIiIiIh8DIN1IiIiIiIiIh/DYJ2IiIiIiIjIxzBYJyIiIiIiIvIxDNaJiIiIiIiIfAyDdSIiIiIiIiIf45AkSbK7E3YLDg5G7dq1vbKv7OxshIeHe2Vfdijv/QfK/3tg/+3F/tuL/Xc6ffo08vLyvLIv8u61HuB31W7sv73Yf3ux//Yqy2s9g3Uvi4qKQmZmpt3dcFt57z9Q/t8D+28v9t9e7D+VF+X9s2b/7cX+24v9txf7L45l8EREREREREQ+hsE6ERERERERkY/xf+65556zuxMVTZcuXezugkfKe/+B8v8e2H97sf/2Yv+pvCjvnzX7by/2317sv73YfzGcs05ERERERETkY1gGT0RERERERORjGKwTERERERER+RgG6244cOAAunbtivj4eHTs2BG7d+/W3G7+/Plo1qwZmjRpggkTJqCwsLCMe6pNpP8///wzqlSpguTk5JKfnJwcG3p7rYcffhixsbFwOBzYtWuX7na+ev5F+u/L5z83NxdDhgxBfHw8kpOTMXDgQKSmpmpu64ufgWj/ffkz6N+/PxITE5GcnIxu3bph+/btmtv54vkHxPrvy+cfAJ5//nnD/4Z99dyTNbze24fXenvxWm8/Xut9g+3Xe4ks69Wrl/TRRx9JkiRJy5Ytkzp37nzNNocPH5bq168vnThxQiouLpZuvvlmae7cuWXcU20i/V+7dq3Uvn37Mu6ZmHXr1kkZGRlSTEyM9Mcff2hu48vnX6T/vnz+c3JypFWrVknFxcWSJEnS22+/LfXr1++a7Xz1MxDtvy9/BufPny/5/auvvpLatm17zTa+ev4lSaz/vnz+t27dKg0cOFCKjo7W/G/Yl889WcPrvX14rbcXr/X247Xefr5wvWdm3aJTp04hJSUFY8aMAQAMGzYMR44cuWa0bvny5Rg6dCjq1q0Lh8OBiRMnYvHixTb0uDTR/vuy7t27IyoqynAbXz3/gFj/fVlISAgGDx4Mh8MBAOjcuTMOHz58zXa++hmI9t+XVa9eveT3ixcvws/v2v+V++r5B8T676vy8vLwt7/9De+++27Jd8iVL597Esfrvb14rbcXr/X247XeXr5yvS9fZ80HZGRkoEGDBggICAAAOBwOREdHIz09vdR26enpiImJKWnHxsZes40dRPsPAPv27UO7du3QoUMHvPvuu2XdVY/46vm3oryc/7feegs333zzNX8vL5+BXv8B3/4M7r77bjRq1AjPPPMMPvnkk2ue9/Xzb9Z/wDfP/7PPPosxY8YgLi5OdxtfP/ckhtd73+er596K8nLuea23B6/19vGV632AV/dWSbiOrkg6d79Tb6e3jR1E+t+uXTtkZmaiWrVqyMzMxODBg1GrVi2MHDmyrLrpMV89/yLKy/l/+eWXceDAAcydO1fzeV//DIz67+ufwcKFCwEAn3zyCaZNm4bVq1dfs40vn3+z/vvi+d+4cSN+//13vPLKK6bb+vK5J3G83vs+Xz33IsrLuee13j681tvDl673zKxb1KhRI2RmZpYsHiBJEjIyMhAdHV1qu+jo6FKlZmlpaddsYwfR/letWhXVqlUDAERFRWHUqFFYv359mffXXb56/kWVh/M/c+ZMrFixAt9++y2qVKlyzfO+/hmY9b88fAYAcM8992Dt2rU4e/Zsqb/7+vlX6PXfF8//unXrsHfvXsTFxSE2NhaZmZkYMGAAvv3221LblZdzT8Z4vfd9vnruRZWHc89rvW/gtb5s+dT13qsz4CuJHj16lFqwpVOnTtdsc+jQoWsWHHjvvffKuKfaRPp/7NgxqaioSJIkSbp06ZLUtWtXaf78+WXZTVNGi7b48vlXGPXf18//66+/LrVr1046d+6c7ja+/BmI9N9XP4OLFy9KR48eLWmvWLFCatiwYckiOgpfPf+i/ffV86+m99+wr557so7Xe/vxWm8fXuvtw2u9b7Hzes9g3Q179+6VOnfuLDVr1kxq3769tGvXLkmSJGn8+PHSypUrS7Z7//33pSZNmkhxcXHS+PHjpfz8fLu6XIpI/99++22pVatWUmJiotSqVSvpH//4xzX/gdll0qRJUsOGDSV/f3+pbt26UpMmTSRJKj/nX6T/vnz+MzIyJABS48aNpaSkJCkpKUnq2LGjJEnl4zMQ7b+vfgbp6elShw4dpDZt2kiJiYlSnz59pG3btkmSVD7Ov2j/ffX8q6kv3uXh3JN1vN7bh9d6e/Faby9e632Lndd7hyT52OQGIiIiIiIiokqOc9aJiIiIiIiIfAyDdSIiIiIiIiIfw2CdiIiIiIiIyMcwWCciIiIiIiLyMQzWiYiIiIiIiHwMg3UiIiIiIiIiH8NgnYiIiIiIiMjHMFgnIq+4mJqKmQ4HfnvuObu7QkRERNcBr/VEZYvBOhEREREREZGPYbBORERERERE5GMckiRJdneCiMq3/730ErbOno2cM2cQGBaGwPBwAEDXf/wDyQ8+aHPviIiIyFO81hOVPQbrROQVF1NT8UFcHLr84x+4kXPZiIiIKhxe64nKFsvgiYiIiIiIiHwMg3UiIiIiIiIiH8NgnYiIiIiIiMjHMFgnIiIiIiIi8jEM1omIiIiIiIh8DIN1IvKKwLAwAIBUWAgAOLdvHza98oqdXSIiIiIv4rWeqGwF2N0BIqoYQmvVQmhkJM7u3g0A2Ld0KU5s2WJzr4iIiMhbeK0nKlu8zzoRec3+FSuwbto05GdloXrjxhj40UeIbNnS7m4RERGRl/BaT1R2GKwTERERERER+RjOWSciIiIiIiLyMQzWiYiIiIiIiHwMg3UiIiIiIiIiH8NgnYiIiIiIiMjHMFgnIiIiIiIi8jEM1omIiIiIiIh8DIN1IiIiIiIiIh/DYJ2IiIiIiIjIxzBYJyIiIiIiIvIxDNaJiIiIiIiIfMz/A7B5jqCeXzyZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Problem statement:\n",
    "y' = -100(y - cos(t)) - sin(t)\n",
    "y(0) = 0\n",
    "\n",
    "Solve with Forward and Backward Euler schemes, stepsize of 0.2 for a total time of 4. Plot both solutions on the\n",
    "interval [0, 4], on top of the exact solution which is given by:\n",
    "y(t) = cos(t) - e^(-100t)\n",
    "\n",
    "Forward Euler: y_n+1 = y_n + dt*f(t_n, y_n)\n",
    "Backward Euler: y_n+1 = y_n + dt*f(t_n+1, y_n+1)\n",
    "\n",
    "y_n+1 = (y_n + dt *(-100(cos(t)) - sin(t))) / (1 + dt*100)\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "from math import e, cos, sin, ceil\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Define parameters:\n",
    "dt = 0.02\n",
    "tot_time = 4  \n",
    "num_steps = ceil(tot_time / dt)\n",
    "time_vec = np.linspace(0, tot_time, num_steps+1)\n",
    "\n",
    "# Allocate mememory:\n",
    "sol_forward = np.zeros((num_steps+1,))\n",
    "sol_backard = np.zeros((num_steps+1,))\n",
    "\n",
    "sol_forward[0] = 0\n",
    "sol_backard[0] = 0\n",
    "\n",
    "# Perform time-integration\n",
    "for ii in range(num_steps):\n",
    "    # Forward Euler:\n",
    "    sol_forward[ii + 1] = sol_forward[ii] - 100 * dt * (sol_forward[ii] - cos(time_vec[ii])) - dt * sin(time_vec[ii])\n",
    "\n",
    "    # Backward Euler:\n",
    "    sol_backard[ii + 1] = (sol_backard[ii] + 100 * dt * cos(time_vec[ii]) - dt * sin(time_vec[ii])) / (1 + 100 * dt)\n",
    "\n",
    "# Calculate exact solution:\n",
    "sol_exact = np.zeros((num_steps+1,))\n",
    "for ii in range(num_steps+1):\n",
    "    sol_exact[ii] = -e ** (-100 * time_vec[ii]) + cos(time_vec[ii])\n",
    "\n",
    "# plot with various axes scales\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "# Forward Euler\n",
    "plt.subplot(121)\n",
    "plt.plot(time_vec, sol_forward, 'r', time_vec, sol_exact, 'g--')\n",
    "plt.title('Forward Euler dt = ' + str(dt), fontdict=font)\n",
    "plt.xlabel('t', fontdict=font)\n",
    "plt.ylabel('y', fontdict=font)\n",
    "plt.gca().legend(('forward','analytic'))\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(time_vec, sol_backard, 'r', time_vec, sol_exact, 'g--')\n",
    "plt.title('Backward Euler dt = ' + str(dt), fontdict=font)\n",
    "plt.xlabel('t', fontdict=font)\n",
    "plt.ylabel('y', fontdict=font)\n",
    "plt.gca().legend(('backward','analytic'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide25.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide30.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple example of ODE here (in-class exercise)\n",
    "The example we will go through is the following:\n",
    "* $ y' = 3t^2 + 2t + 1 - y = f(t, y(t)) $, with initial conditions: $ y(0) = 1 $.\n",
    "\n",
    "This equation has the analytical solution: \n",
    "* $ y(t) = -4e^{-t} + 3t^3 - 4t + 5 $\n",
    "\n",
    "**!!!** Convince yourself that this is true by filling the differential equation in (take the derivative and substitute back in the equation, also check for the initial condition) **!!!**\n",
    "\n",
    "### Forward Euler\n",
    "$ \\frac{dy}{dt} \\approx \\frac{y(t_{n+1}) - y(t_{n})}{\\Delta t} = f(t_{n}, y(t_{n})) $ \n",
    "\n",
    "$ y(t_{n+1}) = y(t_n) + \\Delta t f(t_{n}, y(t_{n})) $\n",
    "\n",
    "### Backward Euler\n",
    "$ \\frac{dy}{dt} \\approx \\frac{y(t_{n+1}) - y(t_{n})}{\\Delta t} = f(t_{n+1}, y(t_{n+1})) $\n",
    "\n",
    "$ y(t_{n+1}) = y(t_n) + \\Delta t f(t_{n+1}, y(t_{n+1})) $\n",
    "\n",
    "\n",
    "### Solving the actual problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple ODE example\n",
    "\"\"\"\n",
    "Problem statement:\n",
    "y' = 3*t^2 + 2*t + 1 - y\n",
    "y(0) = 1\n",
    "\n",
    "Solve with Forward and Backward Euler schemes, stepsize of 2.25 for a total time of 100. Plot both solutions on the\n",
    "interval [0, 100], on top of the exact solution which is given by:\n",
    "y(t) = -4e^(-t) + 3t^3 - 4t + 5\n",
    "\n",
    "Forward Euler: y_n+1 = y_n + dt*f(t_n, y_n)\n",
    "Backward Euler: y_n+1 = y_n + dt*f(t_n+1, y_n+1)\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "from math import e, cos, sin, ceil\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Define parameters:\n",
    "dt = 2.25\n",
    "tot_time = 100\n",
    "num_steps = ceil(tot_time / dt)\n",
    "exact_steps = min(50000, 100 * num_steps)\n",
    "time_vec = np.linspace(0, tot_time, num_steps+1)\n",
    "time_vec_exact = np.linspace(0, tot_time, exact_steps+1)\n",
    "\n",
    "# Allocate mememory:\n",
    "sol_forward = np.zeros((num_steps+1,))\n",
    "sol_backward = np.zeros((num_steps+1,))\n",
    "\n",
    "# Set initial conditions:\n",
    "sol_forward[0] = ???\n",
    "sol_backward[0] = ???\n",
    "\n",
    "# Perform time-integration\n",
    "for ii in range(1, num_steps+1):\n",
    "    # Forward Euler:\n",
    "    sol_forward[ii] = ???\n",
    "\n",
    "    # Backward Euler:\n",
    "    sol_backward[ii] = ???\n",
    "\n",
    "# Calculate exact solution:\n",
    "sol_exact = np.zeros((exact_steps+1,))\n",
    "for ii in range(exact_steps+1):\n",
    "    sol_exact[ii] = ???\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot with various axes scales\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "# Forward Euler\n",
    "plt.subplot(121)\n",
    "plt.plot(time_vec, sol_forward, 'r', time_vec_exact, sol_exact, 'g--')\n",
    "plt.title('Forward Euler dt = ' + str(dt), fontdict=font)\n",
    "plt.xlabel('t', fontdict=font)\n",
    "plt.ylabel('y', fontdict=font)\n",
    "plt.gca().legend(('forward','analytic'))\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(time_vec, sol_backward, 'r', time_vec_exact, sol_exact, 'g--')\n",
    "plt.title('Backward Euler dt = ' + str(dt), fontdict=font)\n",
    "plt.xlabel('t', fontdict=font)\n",
    "plt.ylabel('y', fontdict=font)\n",
    "plt.gca().legend(('backward','analytic'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other concepts that remain unclear after your revision?\n",
    "### Remember that the exam is on Friday and this will be one of the last moments at which you can ask questions :)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide31.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve the equation for pressure (difference):\n",
    "Plug in all the relationships defined above. Use the following script to write a new \"my_func\" and simply plug into what we have developed above. \n",
    "\n",
    "I.e.:\n",
    "$$ \\frac{V_0}{E_0} = \\frac{V_0 - \\Delta V}{E} + W_i $$\n",
    "\n",
    "Where the dependencies on pressure are:\n",
    "$$ E_0 = E(p_0) = \\frac{\\alpha p_0}{Z(p_0) T} , \\quad E(p) = \\frac{\\alpha p}{Z(p) T} $$\n",
    "\n",
    "And for water (expansion):\n",
    "$$ W_i = c_w W \\Delta p = c_w W (p - p_0) $$ \n",
    "\n",
    "Rewrite the top equation, while plugging in all the relationships, such that we obtain the following equation:\n",
    "$$ R(p) = 0 $$ \n",
    "\n",
    "**Question:** What is $\\Delta V$???\n",
    "\n",
    "This is what we call residual form. We can find the root of the residual equation, i.e. the solution p that satistifies this nonlinear equation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide10.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide11.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![NewtonMethod](lecture_04/NewtonIteration_Ani.gif \"NewtonMethod\")\n",
    "source: https://en.wikipedia.org/wiki/Newton%27s_method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide12.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bisection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Bisection method:\n",
    "# First we need to define the function for which we want to calculate the root\n",
    "def my_function(x, a):\n",
    "    return (???)\n",
    "\n",
    "# Set some convergence parameters:\n",
    "curr_iter = 0\n",
    "max_iter = 100\n",
    "threshold = 1e-12\n",
    "\n",
    "# Define number for which to calculate root:\n",
    "a = 2\n",
    "\n",
    "# Get initial solution:\n",
    "x_left = ???\n",
    "x_right = ???\n",
    "\n",
    "# Always test if there truly exists a root, i.e. is (f(a) * f(b) < 0) True?\n",
    "does_root_exist = ???\n",
    "\n",
    "if does_root_exist:\n",
    "    print('There is a root on the interval :)')  \n",
    "else:\n",
    "    print('Change left and/or right boundaries, no root on interval :(')\n",
    "\n",
    "# Find first guess:\n",
    "x_new = ???\n",
    "\n",
    "while ???:\n",
    "    ### your code here ###\n",
    "\n",
    "print('For bisection method:')\n",
    "print(\"The square root of {:4.2f} is {:8.6f} (found after {:d} iterations)\".format(a, x_new, curr_iter))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Newtons method:\n",
    "# First we need to define the function for which we want to calculate the root\n",
    "def my_function(x, a):\n",
    "    return (???)\n",
    "\n",
    "# For Newton's method, we also need it's derivative:\n",
    "def my_func_der(x):\n",
    "    return (???)\n",
    "\n",
    "# Set some convergence parameters:\n",
    "curr_iter = 0\n",
    "max_iter = 100\n",
    "threshold = 1e-12\n",
    "\n",
    "# Define number for which to calculate root:\n",
    "a = 2\n",
    "\n",
    "# Get initial solution:\n",
    "x_old = ??? \n",
    "x_new = ???\n",
    "\n",
    "while ???:\n",
    "    ### your code here ###\n",
    "    \n",
    "print('For Newtons method:')\n",
    "print(\"The square root of {:4.2f} is {:8.6f} (found after {:d} iterations)\".format(a, x_new, curr_iter))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve the equation for pressure - equivalent to the challenge assignment\n",
    "Plug in all the relationships defined above. Use the following script to write a new \"my_func\" and simply plug into what we have developed above. \n",
    "\n",
    "I.e.:\n",
    "$$ \\frac{V_0}{E_0} = \\frac{V_0 - \\Delta V}{E} + W_i $$\n",
    "\n",
    "Where the dependencies on pressure are:\n",
    "$$ E_0 = E(p_0) = \\frac{\\alpha p_0}{Z(p_0) T} , \\quad E(p) = \\frac{\\alpha p}{Z(p) T} $$\n",
    "\n",
    "And for water (expansion):\n",
    "$$ W_i = c_w W \\Delta p = c_w W (p - p_0) $$ \n",
    "\n",
    "Rewrite the top equation, while plugging in all the relationships, such that we obtain the following equation:\n",
    "$$ R(p) = 0 $$ \n",
    "\n",
    "**Question:** What is $\\Delta V$???\n",
    "\n",
    "This is what we call residual form. We can find the root of the residual equation, i.e. the solution p that satistifies this nonlinear equation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_E (p):\n",
    "    return (???)\n",
    "\n",
    "# input parameters\n",
    "water_fraction = 1  # assume the reservoir volume of aquifer = volume of gas reservoir\n",
    "depth = 2000\n",
    "p0 = 200  # initial pressure in bar\n",
    "poro = 0.3\n",
    "h = 120\n",
    "\n",
    "T = 10   # in years of extraction\n",
    "Q = 500e6  # gas flowrate [m3/year]\n",
    "gas_extracted_surf = ???  # delta_V\n",
    "\n",
    "# pore-volume occupied by gas\n",
    "pore_volume = 2000*3000*h*poro / 2 # here we assume only half of reservoir is available\n",
    "E0 = compute_E(p0)\n",
    "\n",
    "# surface volume of gas\n",
    "volume = ???\n",
    "\n",
    "# Compressibility factors:\n",
    "rock_compr = 1e-5\n",
    "water_compr = 1e-5 \n",
    "\n",
    "def compute_residual(p):    \n",
    "    return (???)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set some convergence parameters:\n",
    "curr_iter = 0\n",
    "max_iter = 100\n",
    "threshold = 1e-6\n",
    "\n",
    "# Get initial solution:\n",
    "x_left = ???\n",
    "x_right = ???\n",
    "\n",
    "# Always test if there truly exists a root, i.e. is (f(a) * f(b) < 0) True?\n",
    "does_root_exist = ???\n",
    "\n",
    "if does_root_exist:\n",
    "    print('There is a root on the interval :)')  \n",
    "else:\n",
    "    print('Change left and/or right boundaries, no root on interval :(')\n",
    "\n",
    "# Find first guess:\n",
    "x_new = ???\n",
    "\n",
    "while ???:\n",
    "    ### your code here ###\n",
    "\n",
    "print('For bisection method:')\n",
    "print(\"The solution is {:8.6f} (found after {:d} iterations)\".format(x_new, curr_iter))\n",
    "print('Size of residual is {:4.3e}'.format(compute_residual(x_new)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample space: Set of possible outcomes of an experiment\n",
    "<img style=\"float: middle;\" src=\"statistics_pngs/sample_space.PNG\" width=\"60%\">\n",
    "\n",
    "\n",
    "*-taken from Lecture of AESB2440-*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Event: Subset of the sample space\n",
    "I.e.: Event, E, is the sum of two dices equal to 4 --> $E = \\{(1,3),(2,2),(3,1)\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: middle;\" src=\"statistics_pngs/set_theory.PNG\" width=\"60%\">\n",
    "\n",
    "<img style=\"float: middle;\" src=\"statistics_pngs/set_theory_2.PNG\" width=\"60%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability: \n",
    "The extent to which an event is likely to occur, measured by the ratio of the favourable cases to the whole number of cases possible.\n",
    "\n",
    "Take previous event as example: probability that you throw a sum equal to 4, i.e.:\n",
    "$E = \\{(1,3),(2,2),(3,1)\\}$ --> $P(E) = P(\\{(1,3),(2,2),(3,1)\\}) = 3 / 36$, since favorable cases are 3, total number of cases are 36 (6 * 6). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional probability:\n",
    "Given that B is true, what is the probability of A, i.e. $P(A|B)$\n",
    "<img style=\"float: left;\" src=\"statistics_pngs/conditional_prob.PNG\" width=\"100%\">\n",
    "\n",
    "*-taken from Lecture of AESB2440-*\n",
    "\n",
    "Example:\n",
    "* E1 = {exactly one of the eyes of the dice is equal to 1}\n",
    "* E2 = {sum of the eyes of the dice equals 4}\n",
    "\n",
    "What is P(E1) given that E2 holds --> $P(E1|E2)$?\n",
    "\n",
    "E1 intersection with E2 (i.e. shared elements in two sets):\n",
    "$$E1 \\cap E2 = \\{(1,3),(2,2),(3,1)\\} \\cap \\{(1,2), (1,3), ..., (1,6), (6,1), (5,1), ..., (3,1), (2,1)\\} = \\{(1,3), (3,1)\\}$$\n",
    "\n",
    "Probability E1 intersection with E2:\n",
    "$$P(E2\\cap E1) = 2/36$$\n",
    "\n",
    "Now, probability of E2 we already know from previous example:\n",
    "$$P(E2) = P(\\{(1,3),(2,2),(3,1)\\}) = 3/36$$\n",
    "\n",
    "Therefore:\n",
    "$$P(E1|E2) = \\frac{P(E2\\cap E1)}{P(E2)} = 2/3$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random variable: \n",
    "**Informal**: Variable whose values depend on outcomes of a random phenomenon\n",
    "\n",
    "**Formal**: Measurable function defined on a probability space whose outcomes are typically real numbers.\n",
    "\n",
    "<img style=\"float: left;\" src=\"statistics_pngs/random_variable.PNG\" width=\"100%\"> \n",
    "\n",
    "*-taken from Lecture of AESB2440-*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"statistics_pngs/prob_mass_function_discrete.PNG\" width=\"100%\"> \n",
    "\n",
    "*-taken from Lecture of AESB2440-*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: middle;\" src=\"statistics_pngs/prob_mass_function_discrete_example.PNG\" width=\"60%\">\n",
    "\n",
    "*-taken from Lecture of AESB2440-*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a continuous random variable we cannot simply define the probability of $ p(a) = P(X = a)$ since there are infinite of values that a can attain, hence the probability $p(a) = 0$. Instead we have to define the probability over a certain interval. Probability is given in this case (a continuous random variable) by: \n",
    "$$ P(a \\leq X \\leq b) = \\int_a^b{f_X(x)dx} $$\n",
    "\n",
    "where $f_X(x)$ is the probability density (i.e. mass) function (**PDF**).\n",
    "\n",
    "Also note that the distribution function, most often refered to as cumulative probability function (**CDF**), for a continuous variable is calculated in a similar manner:\n",
    "$$ F(X) = \\int_{-\\infty}^{x}{f_X(\\tau)d\\tau} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of continuous random variable (PDF and CDF)\n",
    "\n",
    "<img style=\"float: left;\" src=\"statistics_pngs/pdf_cdf_example.PNG\" width=\"100%\"> \n",
    "\n",
    "*-taken from Lecture of AESB2440-*\n",
    "\n",
    "Above is a special PDF, particularly the standard normal distribution for a continuous random variable, defined as follows:\n",
    "\n",
    "$$ f_X(X|\\mu, \\sigma^2) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}}e^{-\\frac{(x - \\mu)^2}{2 \\sigma^2}} $$\n",
    "\n",
    "where $\\sigma$ and $\\sigma ^2$ is the standard deviation and variance respecitvely.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other important concepts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean and Expectation\n",
    "#### Population mean: mean value of the whole population (often unknown!)\n",
    "#### Sample mean: mean value of the samples  (taken from the population)\n",
    "$$ \\bar{x} = \\frac{x_1 + x_2 + ... + x_n}{n} $$\n",
    "#### Expectation: can be seen as the probability (relative frequency) weighted average/mean\n",
    "Continuous case: \n",
    "$$ E[X] = \\int_{z\\in R}{zf(z)dz} = \\mu $$\n",
    "Discrete case:\n",
    "$$ E[X] = \\sum_{i=1}^{n}{p_i z_i}$$\n",
    "\n",
    "#### Please note that mean and expectation are often used interchangably. Expectation is indeed simply the arithmetic mean. Take the following example:\n",
    "outcomes = {1, 1, 2, 2, 2, 3}\n",
    "\n",
    "mean = sum_outcomes / number_outcomes = (1 + 1 + 2 + 2 + 2 + 3)/6 = 1.83\n",
    "\n",
    "probability/relative_frequency = {2/6, 3/6, 1/6} (for 1, 2, 3 respectively)\n",
    "\n",
    "expectation = sum( p_i * x_i ) = 2/6 * 1 + 3/6 * 2 + 1/6 * 3 = 1.83\n",
    "\n",
    "If you look closely you can simply write this sum 2/6 * 1 + 3/6 * 2 + 1/6 * 3 as (2 * 1 + 3 * 2 + 1 * 3)/6 since they have a common denominator and this further can be expanded as (1 + 1 + 2 + 2 + 2 + 3)/6 **which is exactly how we caclulated the mean, showing that the mean and expectation in this case are exactly the same thing!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercises on: Mean and Expectation\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "# Load data set from porosity_data_shale.txt, assuming this is the whole population!!!:\n",
    "population = ???\n",
    "size_population = ???\n",
    "\n",
    "# Extract a few samples using random integer indexing (take 5 samples)\n",
    "num_samples = 5\n",
    "\n",
    "samples_id = []\n",
    "samples_added = 0\n",
    "while (samples_added < num_samples):\n",
    "    rand_id = np.random.randint(0, size_population)\n",
    "    \n",
    "    if rand_id not in samples_id:\n",
    "        samples_id.append(rand_id)\n",
    "        samples_added += 1\n",
    "\n",
    "samples = population[samples_id]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate and compare sample mean versus population mean:\n",
    "sample_mean = ???\n",
    "population_mean = ???\n",
    "\n",
    "print('Sample mean = ', sample_mean)\n",
    "print('Population mean = ', population_mean)\n",
    "\n",
    "# Expectation:\n",
    "# Calculate the expectation of a discrete variable, result of a single dice roll:\n",
    "result_die = ???\n",
    "probability = ???\n",
    "\n",
    "expectation = ???\n",
    "\n",
    "print('Expectation of fair dice roll = ', expectation)\n",
    "\n",
    "# Now what about an unfair die:\n",
    "# i.e. a die that is way(!) more likely to give a 6 as a result\n",
    "# i.e. p_6 = 0.5, p_1,2,...,5 = ???\n",
    "### your code here ###\n",
    "\n",
    "expectation = ???\n",
    "\n",
    "print('Expectation of new unfair dice roll = ', expectation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variance and Standard Deviation\n",
    "Variance is often refered to as a the spread, or in more technical terms the average squared deviation from the mean.\n",
    "\n",
    "#### Sample variance: the average squared deviation from the sample(!) mean\n",
    "$$ s_n^2 = \\frac{1}{n - 1}\\sum_{i = 1}^{n}{(x_i - \\hat{x})^2} $$\n",
    "\n",
    "Technical remark: the term $1/(n - 1)$ is used to make it an **unbiased sample variance**. In practical cases, the sample size (n) is often large enough to neglect this difference (as we see in an example below). See https://en.wikipedia.org/wiki/Bias_of_an_estimator and https://en.wikipedia.org/wiki/Bessel%27s_correction#Source_of_bias for more information.\n",
    "\n",
    "#### Sample standard deviation: the square root of the sample variance\n",
    "$$ s_n = \\sqrt{s_n^2} $$\n",
    "\n",
    "#### Variance is often also referred to as the second stochastic moment, i.e. second moment about the mean:\n",
    "Theoretical case:\n",
    "$$ \\text{var}(X) = E[(X - E[X])^2] = E[(X - \\mu)^2] = \\sigma^2 $$\n",
    "\n",
    "#### NOTE:\n",
    "Variance is not a unitless quantity, just like the expectation or mean! Therefore, having a high number for the variance doesn't necessarily mean that the spread is large and vice versa (e.g. mean = 1e6, var = 100 vs. mean = 0, var = 1)! Always consider the relative quantities!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercises on: Variance and Standard Deviation\n",
    "# Extract a few samples using random integer indexing (take 5 samples)     \n",
    "def mean_samples(sample_array):\n",
    "    return ???\n",
    "\n",
    "# Calculate sample variance and population variance:\n",
    "def sample_variance_biased(sample_array):\n",
    "    return ???\n",
    "\n",
    "def sample_variance_unbiased(sample_array):\n",
    "    return ???\n",
    "\n",
    "num_samples = 5\n",
    "samples_id = []\n",
    "samples_added = 0\n",
    "while (samples_added < num_samples):\n",
    "    rand_id = np.random.randint(0, size_population)\n",
    "    \n",
    "    if rand_id not in samples_id:\n",
    "        samples_id.append(rand_id)\n",
    "        samples_added += 1\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = np.array(population[samples_id], copy=False)\n",
    "\n",
    "# Biased:\n",
    "biased_sample_var = sample_variance_biased(samples)\n",
    "\n",
    "# Unbiased:\n",
    "unbiased_sample_var = sample_variance_unbiased(samples)\n",
    "\n",
    "# Numpy var:\n",
    "numpy_sample_var = np.var(samples)\n",
    "\n",
    "# Print results:\n",
    "print('Biased Sample Variance: ', biased_sample_var)\n",
    "print('Unbiased Sample Variance: ', unbiased_sample_var)\n",
    "print('Numpy Sample Variance: ', numpy_sample_var)\n",
    "\n",
    "# What about population:\n",
    "print('Population variance: ', np.var(population))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "e) Plot the histogram, fit the analytical continuous PDF, compute the following probability: \n",
    "* $ P(|X - \\mu| \\leq \\sigma ) $\n",
    "* $ P(0.0 \\leq X - \\mu \\leq 0.10 ) $\n",
    "\n",
    "Recall that the probability is given by (for a continuous random variable): \n",
    "\n",
    "$$ P(a \\leq X \\leq b) = \\int_a^b{f_X(x)dx} $$\n",
    "\n",
    "Also recall that the trapezoidal rule gives us:\n",
    "\n",
    "$$ \\int_a^b{f_X(x)dx} \\approx \\sum_{i=1}^n{\\frac{f_X(x_{i-1}) + f_X(x_i)}{2} \\Delta x} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now, fit gaussian distribution, plot, and calculate probabilities\n",
    "from math import pi, exp\n",
    "\n",
    "\n",
    "# Histogram:\n",
    "plt.figure(num=None, figsize=(10, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "plt.hist(population, 20)\n",
    "plt.title('Histogram Porosity Data', fontdict=font)\n",
    "plt.xlabel('porosity', fontdict=font)\n",
    "plt.ylabel('frequency', fontdict=font)\n",
    "plt.show()\n",
    "\n",
    "mu = ???\n",
    "sigma = ???\n",
    "\n",
    "def my_gaussian_distribution(x):\n",
    "    return ???\n",
    "\n",
    "my_distribution = np.zeros((size_population,))\n",
    "\n",
    "for ii in range(size_population):\n",
    "    my_distribution[???] = ???\n",
    "    \n",
    "# Plot distribution:\n",
    "plt.figure(num=None, figsize=(10, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "plt.plot(population, my_distribution)\n",
    "plt.title('Probability Distribution', fontdict=font)\n",
    "plt.xlabel('porosity', fontdict=font)\n",
    "plt.ylabel('probability', fontdict=font)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate probability that a new porosity measurement in the region is between 0 and 0.10:\n",
    "# Use integration\n",
    "def my_func(x):\n",
    "    return ( ??? )\n",
    "\n",
    "def my_integral_func(sigma_factor):\n",
    "    ### your code here ###\n",
    "    return ???\n",
    "\n",
    "print('Probability 1: ', my_integral_func(1))\n",
    "\n",
    "num_p = 1000\n",
    "int_low = 0\n",
    "int_high = 0.1\n",
    "int_len = int_high - int_low\n",
    "dx = int_len / num_p\n",
    "\n",
    "int_sum = 0\n",
    "a = int_low\n",
    "\n",
    "for ii in range(num_p):\n",
    "    ### your code here ###\n",
    "    \n",
    "print('Probability 2: ', int_sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional topics to read upon:\n",
    "1) Joint Probability Distributions\n",
    "\n",
    "2) Covariance\n",
    "\n",
    "3) Correlation\n",
    "\n",
    "4) Covariance Function and Variogram\n",
    "\n",
    "5) Kriging\n",
    "\n",
    "In the three resources given on the BrightSpace you can find all the necessary details on these subjects. \n",
    "\n",
    "You will use these concepts in further classes, so it is worth to investigate some time into learning these concepts. However, it will not be mandatory to know these additional topics on the exam (the basic statistics covered above will be possibily tested on the exam). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}