{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical Methods \n",
"\n",
"\n",
"# Lecture 2: Numerical Differentiation\n",
"\n",
"\n",
"## Exercise solutions"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# some imports we will make at the start of every notebook\n",
"# later notebooks may add to this with specific SciPy modules\n",
"\n",
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 2.1: Compute first derivative using forward differencing\n",
"\n",
"Use the forward difference scheme to compute an approximation to $f'(2.36)$ from the following data:\n",
"\n",
"$f(2.36) = 0.85866$\n",
"\n",
"$f(2.37) = 0.86289$\n",
"\n",
"You should get an answer of 0.423."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.42299999999999693\n"
]
}
],
"source": [
"dx = 2.37-2.36\n",
"df = (0.86289-0.85866)/dx\n",
"print(df)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 2.2: Compute first derivative using central differencing\n",
"\n",
"Use the data below to compute $f'(0.2)$ using central differencing:\n",
"\n",
"$$f(0.1) = 0.078348$$\n",
"$$f(0.2) = 0.138910$$\n",
"$$f(0.3) = 0.192916$$\n",
"\n",
"You should get 0.57284"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.57284\n"
]
}
],
"source": [
"dx=0.1\n",
"df = (0.192916-0.078348)/(2*dx)\n",
"print(df)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: Write a function to perform numerical differentiation\n",
"\n",
"As covered above, the expression\n",
"\n",
"$$\\frac{f(x+\\Delta x) - f(x-\\Delta x)}{2\\Delta x},$$\n",
"\n",
"can be used to find an approximate derivative of the function $f(x)$ provided that $\\Delta x$ is appropriately small. \n",
"\n",
"Let's write a function `diff(f, x, dx = 1.0e-6)` that returns the approximation of the derivative of a mathematical function represented by a Python function `f(x)`.\n",
"\n",
"Let's apply the above formula to differentiate $\\,f(x) = e^x\\,$ at $\\,x = 0$, $\\,f(x) = e^{−2x}\\,$ at $\\,x = 0$, $\\,f(x) = \\cos(x)\\,$ at $\\,x = 2\\pi$, and $\\,f(x) = \\ln(x)\\,$ at $\\,x = 1\\,$, i.e. functions we know the exact derivative of.\n",
"\n",
"In each case, using $\\,\\Delta x = 0.01$, let's write out the error, i.e. the difference between the exact derivative and the result of the formula above."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The approximate derivative of exp(x) at x = 0 is: 1.000017. The error is 0.000017.\n",
"The approximate derivative of exp(-2*x) at x = 0 is: -2.00013. The error is 0.00013.\n",
"The approximate derivative of cos(x) at x = 2*pi is: 0.00000. The error is 0.00000.\n",
"The approximate derivative of ln(x) at x = 0 is: 1.00003. The error is 0.00003.\n"
]
}
],
"source": [
"def diff(f, x, dx=1e-6):\n",
" numerator = f(x + dx) - f(x - dx)\n",
" derivative = numerator / ( 2.0 * dx )\n",
" return derivative\n",
"\n",
"dx = 0.01\n",
"x = 0\n",
"f = np.exp\n",
"derivative = diff(f, x, dx)\n",
"print(\"The approximate derivative of exp(x) at x = 0 is: %f. The error is %f.\"\n",
" % (derivative, abs(derivative - 1)))\n",
"x = 0\n",
"\n",
"def g(x):\n",
" return np.exp(-2*x)\n",
"\n",
"derivative = diff(g, x, dx)\n",
"print('The approximate derivative of exp(-2*x) at x = 0 is: {0:.5f}. The error is {1:.5f}.'\n",
" .format(derivative, abs(derivative - (-2.0))))\n",
"\n",
"x = 2*np.pi\n",
"f = np.cos\n",
"derivative = diff(f, x, dx)\n",
"print('The approximate derivative of cos(x) at x = 2*pi is: {0:.5f}. The error is {1:.5f}.'\n",
" .format(derivative, abs(derivative - 0)))\n",
"\n",
"x = 1\n",
"f = np.log\n",
"derivative = diff(f, x, dx)\n",
"print('The approximate derivative of ln(x) at x = 0 is: {0:.5f}. The error is {1:.5f}.'\n",
" .format(derivative, abs(derivative - 1)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 2.3: Compute the derivative of $\\sin(x)$\n",
"\n",
"Compute \n",
"\n",
"$$\\frac{d(\\sin x)}{dx}\\qquad\\textrm{at}\\qquad x = 0.8$$\n",
"\n",
"using (a) forward differencing and (b) central differencing. \n",
"\n",
"Write some code that evaluates these derivatives for decreasing values of $h$ (start with $h=1.0$ and keep halving) and compare the values against the exact solution.\n",
"\n",
"Plot the convergence of your two methods."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact derivative at sin(0.8) = 0.6967067093471654\n",
"Forward differencing Central differencing\n",
" 0.256492 (error= 0.44) 0.586258 (error= 0.11)\n",
" 0.492404 (error= 0.2) 0.668038 (error= 0.029)\n",
" 0.600269 (error= 0.096) 0.689472 (error= 0.0072)\n",
" 0.650117 (error= 0.047) 0.694894 (error= 0.0018)\n",
" 0.673843 (error= 0.023) 0.696253 (error= 0.00045)\n",
" 0.685386 (error= 0.011) 0.696593 (error= 0.00011)\n",
" 0.691074 (error= 0.0056) 0.696678 (error= 2.8e-05)\n",
" 0.693897 (error= 0.0028) 0.6967 (error= 7.1e-06)\n",
" 0.695304 (error= 0.0014) 0.696705 (error= 1.8e-06)\n",
" 0.696006 (error= 0.0007) 0.696706 (error= 4.4e-07)\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAFYCAYAAAAmzDckAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1zV1f/A8dcBxL0HGlZqWu5JKo5CrdTU0rShtqwEXGhqpllqWj+zcmeiado3zS3mnok5KFc5Sc1BriTci31+fxxAQFDAC5/L5f18PO6DPud+xvtyP/H2nM8ZSmuNEEIIIe5wsjoAIYQQwt5IchRCCCGSkeQohBBCJCPJUQghhEhGkqMQQgiRjCRHIYQQIhlJjsLhKaU8lVILlVLnlFKRSqmLSqkNSqm3lFLOVscnHoxS6pRSanYGjntbKfVOJoQkHIAkR+HQlFL9gO1AMeBD4BngHeAoMBVoa110wmJvY+4FIe7iYnUAQmQWpdRTwDjgG621X7K3f1ZKjQPyZ31kDy6uxqu01tFWxyKEI5Kao3Bkg4FLwKCU3tRaH9da74/fVkrVV0ptVErdUErdVEptUkrVT3yMUmq2UuqMUqqOUmqrUuqWUuqYUso32Xm0Uqpd8msqpaYqpf5TSuVKVNZdKbVPKRWulApTSs1UShVLdpxWSn2ulBqslDoJRAI14t6rGxdLuFLqtFLqI6XUp0opnewcLkqpIUqpv5RSEXHNzGOVUnkS7VMu7lo+SqmRSqnzSqkrSqkVSqmyKXye7kqpvUqp20qpy0qpLUqpRonez6eUGqOUOhnXpH1SKTVUKXXPvz2J4uiplBqnlAqN+12vVEqVu9exib6DVL9LpVQg8DTQOO46Oq5MCENrLS95OdwLcAZuAT+lcf+awG1gD9AJ6AjsiiurlWi/2cA1IBjwAZ4FfgI00CzRfn8BC5NdwxW4CExOVPYFEAWMBZ4DugFngd8B50T76bjyrXGxtQLcgBLAZeAQ8ArQHtgChJj/vZNcfz5wExiGaV7uA1wBliTap1zctU7Ffa7WwFtAGLAl2fm+jtt3BtAOaAOMAl6Le98lLt6LQD+gBTAUCAfG3uf7iI/jNLAi7tzdgPOYJvFcifY9BcxOz3cJVAX2AvuAhnGvqlbft/Kyn5flAchLXpnxikscGhidxv0XxyWKIonKCmFqnksTlc1OIRHmjkse0xOVDY37Y1w4UVn7uGPrx22XA2KAYcliaRy3X/tEZRo4B+RNtu//YWqRZROV5QUuJE6OQNO4c7yZ7PiuceW1E8WkU0iEA+PKH4rbrhgX+7h7/E7fiDvmqWTlQ+NiLnWPY+PjOAw4pfC7eTdRWfLkmNbvMhDYZvW9Ki/7fEmzqhDGU8BKrfWV+AKt9TVgOab5LbFbWuvNifaLAI4BjyTaZw4mab6cqOwN4IjWemfc9rOYRxtz45o8XZRSLpha47W4mBJbq7W+naysIRCktT6TKJ7bwKpk+7XCJKQlya61PtHnTyz58QfifsZ/xmfiYp9O6lpharA7UrhmrrjY72ex1jo2fkNrvR04A3je45j0fJdCpEiSo3BUFzE1t0fTuH8xTJNdcv8CRZOVXU5hvwgg4dmd1joE+BV4HUApVQTTNPhjomNKxf38G9O0mvhVCCie7BopxVcGCE2h/EKy7VKYZt0bya4Tf2zya11Kth0R9zP+M8bvf4bUlcL8/pN/tvh/HCS/ZkqSf474Mvd7HJOe71KIFElvVeGQtNbRcR0snlVK5Y6r3d3LJaB0CuWluTtRpNWPwHdKqUeBlpjkNDfR+xfjfj5Hygn3YrLtlNaXO8+dJJuYWwrnCsc0r6bkXCrlqQmL++kOHElln4vAScyz0JScSsN1kn+O+LI/73FMZnyXIoeRmqNwZF9gaidfpfSmUqq8Uqpm3OYWoI1SqmCi9wtiOppsyeD1F2ESUldMk+qvWutTid7fAMQCj2itd6fwOpmGa/wGeCbuSaqUyouppSa2FlPrK5zKtdKbHDfGxe59j33WAg8DN1K5Ztg9jo3XKXHPVqVUY6AsEHSPY9L6XUZgns8KcRepOQqHpbX+VSnVHxinlKqC6UzzD6ZprQXwHtAF2I/pZdkW2KSUGoOppX0I5ANGZvD615RSy4FemObP7snePx53rW+UUk9g/nCHYxLKs8CMxM82UzEO6AGsU0p9ivmD3z/uZ0JNU2sdqJSaByyOG9+5E5PcygHPAx9qrY+m47MdV0qNB/rHJZ7lmA469YG/tNYLMLXkbpjf6VhMz1BX4DHgBUyHo1v3uVRBYJlSahpQEhiNeb77v3sck9bv8jDQUyn1KnAcuK61Tq0WLHIaq3sEyUtemf0CGmFqcecxz7wuYTqFvE7SnpANMDWiG5ghD5uI61maaJ/ZwJkUrhEIBKZQ3gbzxzlJz9Vk+7yBqQHejLt2MPANSXugauCzVI6vC2zDJNazwCfAROBysv2cgL6YJBUOXI377y/jY+NOL9H3kh3rFVfulazcF/OPi4i432sg4Jno/TzACMzQlvh9dsWVudzjO4uPoyfmHwD/YYbmrALKJ9v3FIl6q6bjuywNrAaux13rru9PXjn3pbRO6TGGECK7ips9Zy8QprVuYXU8GRE30P8k0F1rPcPaaEROJM2qQmRzSqlRmB6vIZhnrO9hBsI/b2VcQmRnkhyFyP40Ztabh+L+ez/med4aS6MSIhuTZlUhhBAiGRnKIYQQQiSTY5pVS5QoocuVK2d1GCIDbt68Sf782XJlKWFjci+IeLa6F/bs2ROmtS6ZvDzHJMdy5cqxe/duq8MQGRAYGIiXl5fVYQg7IPeCiGere0EpFZJSucM3qyql2imlpl+9etXqUIQQQmQTDp8ctdYrtNbehQsXtjoUIYQQ2YTDJ0epOQohhEgvh0+OUnMUQgiRXg6fHIUQQoj0cvjeqkqpdkC7ihUr3nO/a9euERoaSlRUVNYEJtKscOHCBAcHWxpDrly5KFWqFIUKFbI0DiFE1nD45Ki1XgGs8PDw6J7aPteuXePChQu4u7uTN29elFJZGKG4n+vXr1OwYMH775hJtNbcvn2bs2fPAkiCFCIHkGZVIDQ0FHd3d/LlyyeJUdxFKUW+fPlwd3cnNDTU6nCEEFnA4ZNjWnqrRkVFkTevLAgu7i1v3rzS7C5EDuHwyTGtvVWlxijuR+4RIezD/PkwfXp5goIy7xoO/8xRCCFE9nf6tEmKM2bA0aMAj7BsGWzaBJ6etr+eJEchhBB26b//YPFi+Okn2LbNlJUtC0qB1orISAgMzJzk6PDNqjJDTvbSu3fvdE8mvHv3bpRSnDp1CjATEiulCAsLS9jn559/plKlSri4uPD222+nWiaEsNa1a/C//0Hr1lCmDPTsCZcvw2efwd9/w8KFkCcPODnF4uoKmTUPvcMnR0eeIeftt99GKXXX688//7Q6NEs1atSI8+fPU7x48YSy9957j44dOxISEsLEiRNTLRNCZL3wcFiyBDp1Ajc3eOst+OsvGDQI9u+Hgwdh6FB47DFTS9y0Cd5551SmNamCNKtme8888ww//vhjkrISJUpk+HxRUVHkypXrQcO6p9jYWLTWODs7Z8r5XV1dKV26dML2lStXCAsLo2XLlri7u6daJoTIOtHRJsnNmwdLl8L16yYxdu8OnTtDw4am+TQlnp4QEfEPnp4VMi0+h685ZrWgIBg9mkztRZVY7ty5KV26dJKXi4v5N09ERAT9+vXDzc2NPHny0LBhQ7bFN9xzp/lx9erV1K9fH1dXV9atW0euXLn4/fffE/YrW7YsVapUSdjesGED+fPnTxjWMG7cOGrWrEn+/Plxd3fnvffe48qVKwn7z549mwIFCrB69WqqV6+Oq6srwcHBxMTEMHDgQIoWLUrRokXp168fMTEx9/3Ma9eupXLlyuTJk4emTZty1Dydv+tzhYWFERgYSNGiRQFo3rw5SqlUy4QQmSs21jw77NULHnoIWrWCZcvg5ZdhwwY4cwYmTTLJz+rO4VJzTEW/fpDe1smrV00TQGwsODlBzZqQntbc2rVhwoT0XfNeBg0axMKFC/n++++pUKEC48aNo1WrVhw7dowyZcok7Pfhhx8yduxYKlasSMGCBalbty6bN2+mQYMGHDt2jKtXr3Lx4kXOnz9PmTJlCAwMpFGjRgk1TCcnJyZMmECFChUICQmhT58+9OnTJ0mNNjw8nM8++4xp06ZRsmRJypQpw9ixY/nuu+/47rvvqFmzJlOmTGHu3LnUrVs31c90+vRp2rdvT/fu3enVqxf79++nf//+qe7fqFEjDh06RLVq1ViyZAmNGjWiWLFiKZYJIWxPa/O3dN4809v09GnImxfatTM1xNatIXduq6O8m8PXHLOyQ87VqyYxgvmZFX2A1q5dS4ECBRJerVu3BuDmzZtMnTqVMWPG0KZNG6pUqYK/vz9ubm5MmTIlyTlGjBjBc889R4UKFShZsiReXl5s3rwZMLWwJk2aUL9+/YTaVfIVuPv160fz5s0pV64cTz/9NF9++SULFy4kNv6XAcTExDB58mQaN27M448/TsGCBZkwYQKDBg3ilVdeoXLlykycODFJc2hKpk6dyiOPPMKkSZOoXLkyr7zyCr6+vqnu7+rqSqlSpQAoVqwYpUuXTrVMCGE7R4/Cp59C1apQty6MH28qDHPmQGgoLFgA7dvbZ2KEHFBzTMvcqinJSA0uKAhatIDISHB1hblzM+9hcbynnnqK6dOnJ2zHz/Rz/PhxoqKiaNy4ccJ7zs7OeHp6cvjw4STn8PDwSLLt5eXFlClTiIqKIjAwkGbNmnHz5k0CAwN58cUX2bVrF19++WXC/r/88gujR48mODiYq1evEhMTQ2RkJP/++y8PPfQQAC4uLtSuXTvhmKtXr3L+/Hk8E/2CnJycaNCgAadPn0718wYHB9OwYcMkA/I9M/uXLIRIkzNnTNKbNw/27DFNo08/De+/Dx07QqI+cnbP4ZNjVorvRRUYaLoXZ8Xf7Hz58pHSiiNaayDlWV2Sl+XPnz/JdtOmTYmIiGDXrl1s2bKFfv36cePGDXx8fNi+fTu5cuWifv36AISEhNCmTRu6d+/OyJEjKV68OHv37qVz585ERkYmnDN37tw26YAT/7mEEPYhLMyMRZw3D7ZuNc2oHh4wdiy8+ipk1/5ukhxtzNMza5Li/VSsWBFXV1e2bdtGhQqmR1dMTAxBQUF06dLlnscWKFCAunXrMn36dK5fv07dunWJiorin3/+Ye7cuUmeN+7evZvIyEjGjx+fkPxWrlx53/gKFy5MmTJl+O2332jevDlgEt/OnTuTPA9NrmrVqixZsgStdUKS/+233+7/CxFC2Mz16/DzzyYhrl9vep5WqWKaUV97DSpVsjrCB+fwzxxzqvz589OjRw8GDx7M6tWrCQ4OpkePHly4cIGePXve93gvLy/mzJlD06ZNcXZ2Jk+ePDRo0IA5c+Yked5YqVIlYmNjmTBhAidPnmTevHlMSGObdN++ffnyyy9ZvHgxR44coV+/fpw/f/6ex/j6+nLq1Cn69evHkSNHWLx4Mf7+/mm63v38888/VK5cmRUrVtjkfEI4kvBwCAiAV16BUqXgjTfM+MMBA0yHm0OH4JNPHCMxgiRHhzZmzBheeeUVunXrRu3atdm/fz9r1669Z80sXrNmzYiJiUmSCFMqq1mzJhMnTmTcuHFUrVqVGTNm8PXXX6cpvgEDBtCtWzfee+89GjRoQGxsLF27dr3nMY888ghLly5l7dq11KpVi/Hjx/PFF1+k6Xr3ExkZyZEjR5DZlEROFz8kbetWUzPs1s2MQXzpJfPY6N13zZCMkyfhiy+gVi3rh17Ymsopz3A8PDz07t27U3wvODg4yTg+YV+sXuw4MblXrJW8p7Swve3b73QsjE8PhQqZxNi5MzRvDi528EDOVveCUmqP1tojebkdfEQhhBBW0hp27TI9TWfMgIiIO+917WrK8uSxLj4rOHxyVEq1A9ql1KNTCCFyqvjB+QsWmMm8T56EXLmgfn2TKGNizJC0Xr1yXmKEHJAcMzrOUQghHNGhQyYhLlhgBuo7O8Mzz5jONO3bQ9Gi5pljVg5Js0cOnxyFECKnO3r0TkI8dMhMb+nlZXqavvQSJF+rwF6GpFlJkqMQQjigkydNc+n8+XfmiW7SBCZPNktD3WemxhxPkqMQQjiIM2dMQlywAHbuNGUNGsC4cWbli7JlrY0vO5HkKIQQ2di//5rp2+bPN8MwAOrUgTFjzID9cuUsDS/bkuQohBDZTFgYLFliaohbtphVgKpXh1GjzHymjjJLjZUkOQohRDZw+bKZvm3BArPAQUwMPPEEfPyxSYhVq1odoWOR6eOEJUaMGEH16tXTdUxYWBhKqYR1JU+dOoVSisQzH23fvp2aNWvi6uqaMHtGSmVCZAfXrpn1D9u1M9O3vfsuHDsGH3wAf/wBwcF31kzMSYKCgpg7dy5BQUGZdo1smRyVUhWUUjOVUoutjsVqFy5coG/fvjz22GPkzp0bd3d3WrduzerVq216nYwks8z28MMPc/78+STrRPbt25datWpx/Phxli5dmmqZEPbq5k1TO3zppTsTfO/bB35+ppPN8eNm3tPatR1vPtO02Lp1K08//TQzZ86kRYsWmZYgs7xZVSn1PdAWCNVaV09U3gqYCDgDM7TWqc4mrbU+Abyb05PjqVOnaNy4MQULFmT06NHUqlWL2NhYNm3ahK+vL//880+WxxQZGYmrq2uWXMvZ2ZnSyfqj//333/Tq1YuHH374nmVC2JPwcFizxiTFFSvg1i0z1MLb2ywB1bChGZuY0/3++++8+uqrREVFAebvTWBgYOYseK61ztIX8BRQFziYqMwZOA5UAFyBfUBVoAawMtmrVKLjFqf1uvXq1dOpOXz4cKrvpdeOHTv0//3f/+kdO3bY7Jypad26tS5Tpoy+fv36Xe9dunQp4b+vXLmiu3fvrkuWLKkLFCign3rqKb1r166E92fNmqXz58+vN27cqKtVq6bz5cunvby89IkTJxLeB5K8Zs2apbXWGtDffPON7tChg86XL58eMGCAjo6O1u+8844uV66czpMnj65YsaIeM2aMjomJSbjm8OHDdbVq1e75+Xbu3Knr1q2rc+fOrWvXrq1XrlypAb1582attdYnT57UgN61a1fCfyePMbW4M8qW94pIv/jv3hFs2aL1m29q3bKl1gULag1alyihta+v1ps3ax0dbXWE9uPixYvax8dHK6V0iRIltKurq3ZyctJ58+Z94L+1wG6dQs7I8pqj1vpXpVS5ZMX1gb+1qRGilJoPvKi1Ho2pZWaIUsob8AZwc3NLeFaVXOHChbl+/XqSsg8//JADBw6k63rXrl3j4MGDxMbG4uTkRPXq1SlUqFCaj69RowZjxoxJ076XLl1i7dq1fPLJJ2it74rfxcWF69evo7WmVatWFCpUiAULFlC0aFF++uknmjdvzp49eyhdujTh4eFERETw2WefMXnyZPLkyYOvry/vvfcey5Yt4/nnn6dPnz6sXbs2obm2UKFCCdccMWIEw4cPZ8SIESiluHr1KiVKlGDWrFmUKFGCPXv20LdvX/Lnz8+bb74JQEREBLGxsXfFHe/mzZs8//zzNGnShMmTJ3PhwgX69+8PwK1bt7h+/To3btxI2LdixYocO3aMWrVqMWzYMDp27EiBAgXuKkscd0aEh4eneh+JzHfjxo1s/fuPjFTs2VOMn39+iN9/LwYoQNOw4UU6djxLnTpXcHY2S2Fs3WppqHZBa826deuYNm0a165do1OnTrz99tucPHmSnTt3Ur9+fSIiIjLlnrCX3qruwOlE22eABqntrJQqDnwO1FFKDYlLonfRWk8HpoNZsiq1zhjBwcF3LYnk6uqasLJ9Wl2/fp3Y2FiAhD/8RYsWTfPxrq6uaV6aKTg4GK01tWvXvucxv/zyCwcOHOC///4jb968ANSuXZv169cTEBDAoEGDyJMnD9HR0fj7+/PEE08AMGjQILp160b+/PkpWLAgxYoVw9XVlZQmcH/ttdfo3bt3krLESb569eoEBwcTEBBAr169AMidOzdOTk6pxj5v3jyioqL48ccf0VrTqFEjbt68yRtvvEG+fPkoWLAgBQoUAMzCzkWKFKFIkSIopXBzc0sSZ0plGZUnTx7q1KnzwOcRGZMdl6wKDzdrIi5aBMuXm042iSfydnZWvPBCCQYOLJH6SXKgQ4cO0aNHD7Zu3Yqnpyf+/v7UrFkz4f1q1apl6r1gL8kxpcfKqS40qbW+CPim6cQZXJUjravZJxYUFESLFi0SnrvNnTs3c9rCIb5Z+b727NnDrVu3KFmyZJLy8PBwjh8/nrCdO3fuhMQI8NBDDxEVFcWVK1coVqzYPa/h4XHXUmj4+/szY8YMQkJCuH37NlFRUTz66KNpihlM8q9ZsyYFChRIqOll1u9SCFu7fRvWrjWD81esgOvXzYTenTqZV/780KqVWTPR1dXMcyqMGzduMHLkSMaPH0+hQoWYMWMG3bp1wymLH7raS3I8AyTuLVEWOGeLE+ssXJXD09OTTZs2JfzrNjP/mFeqVAmlFMHBwXTo0CHV/WJjY3Fzc2NrCm00iZt8XZKtXqriusHF14TvJX/+/Em2FyxYQL9+/fj6669p1KgRhQoVYsqUKQQEBNz3XPHSmvyFsBe3bsHq1SYhrlxpep0WL27GIL78MjRrZpaEirdpk6x8kZjWmmXLltG3b19Onz7Nu+++yxdffEGJ5LOiZxF7SY67gEpKqfLAWeA1oIstTpzV6zl6enpmSQ2nWLFitGzZkm+++QY/P7+EJsZ4V65coUiRItStW5cLFy7g5OREhQoVMnw9V1dXYmJi0rTvtm3baNCgQZKm1sS11LSoWrUqP/zwAzdv3kwo++2339J1DiEy240bJiEuWmR+3roFJUvC66+bGqKXF7ik8ldWVr644+TJk/Tp04dVq1ZRo0YN5s2bR+PGjS2NKcs7Byul5gFBwBNKqTNKqXe11tFAb2AdEAws1FofssX1tNYrtNbehQsXtsXp7Mq3336L1hoPDw8WLVrEkSNH+Ouvv5g6dWpC2/wzzzxD48aNefHFF1mzZg0nT54kKCiI4cOHp1ibTE25cuUICQlh7969hIWFEZF4qfBkHn/8cfbu3cuaNWs4duwYo0aNYsuWLen6bF26dMHFxYV33nmH4OBgNmzYwOeff56uc6QmJiaGypUr4+/vb5PziZzl+nWYNw86djTjEF991XSeeest+OUXOHcO/P3NGompJUZhRERE8Pnnn1O1alW2bNnC2LFj2bt3r+WJESxIjlrrzlrrMlrrXFrrslrrmXHlq7XWj2utH9Na2+avIKbmqJSafvXqVVud0m6UL1+evXv38uyzz/Lhhx9Ss2ZNmjdvzvLly5k2bRpgmkdXr15N8+bN6d69O0888QSvvPIKR44c4aGHHkrztTp27Mjzzz9PixYtKFmyJPPmzUt1Xx8fH1555RW6dOnCk08+yalTpxgwYEC6PluBAgVYuXIlx44do2nTpgwcODDNPXnvR2vNkSNHCAsLs8n5hOO7ehXmzjWLAZcsCV26mAWB333XzG169ix8+61pOpWEmDa//PILtWrV4uOPP6Zt27YEBwfTv3//ux7xWCal8R2O+MqqcY7C9q5du2Z1CAnkXrFWVo5zvHxZ6x9+0LptW61dXc04RHd3rfv21XrrVq0TDdsV6XD+/HndpUsXDegKFSroNWvWZOg8troXsJdxjkIIYa8uXYKffzadajZsgKgoePhh6NXLdKpp0EBmqsmomJgYpk6dytChQwkPD2fYsGEMHjw4YYiZvXH45JjVHXKEENlLWJhJiIsWmR6k0dHw6KPQt6/pVFO/fs6cw9SWdu3aha+vb8JjoClTplDJztfVcvjkqLNwKIcQInv47z+z/NPixaYTTUwMVKgA/fubGmK9epIQbeHy5csMHToUf39/SpcuzYIFC3j55ZcThorZM4dPjkIIAXDhgkmIixaZ8YWxsVCxIgwaZGqIdepIQrQVrTVz5sxh4MCBhIWF4efnx8iRI9M1nabVHD45prVZNX4+VCFSk5YJEYT9CAoys9NERMDevfDrryYhPvEEfPSRSYg1a0pCtLXDhw/Ts2dPtmzZQsOGDVm3bl2SZeWyC4dPjmlpVs2fPz9nz57Fzc2NXLlyZYsqv8g6WmuioqK4cOHCXbMBCfsTEgLjxsE335hkCFCuHHz8sWkyrVZNEmJmuHnzJqNGjWLs2LEULFiQadOm8d5772XbSofDJ8e0KFu2LGFhYYSEhBAdHW11OCKZ8PBw8iSeqdkCLi4uFC5c2LKprMS9HT0KS5aY1549Sd9zdjbrIg4ZYk1sOcHy5cvx8/MjJCSEt99+my+//PKu+ZyzG0mOgJOTE6VKlaJUqVJWhyJSEBgYKCthiCS0hgMH7iTEQ3HzadWvD2PGQPnyZsYamdg7c506dQo/Pz9WrFhBtWrV+PXXX2natKnVYdmEwydHGcohhGPQGoKDC7JmDSxdCn//bZpHmzaFiROhQwczJjFe2bIysXdmiYyMZNy4cYwcORInJye++uor+vbtS67EM6tncw6fHGUohxDZV0wMbN9uaodLl8KZM/VwcYEWLeCDD+DFF8HNLeVjZWLvzBEYGEjPnj0TVgSaOHEiDyf+V4mDcPjkKITIXiIjYfNmkwyXLYPQUMid26x/+PrrwQwaVIV0rCEubOTChQsMHDiQOXPmUL58eVauXEmbNm2sDivTSHIUQlju9m1Yv97UEFesgCtXoEABaNMGXnoJnn/ebAcGXqBo0SpWh5ujxMTEMG3aND766CNu3brF0KFD+eijj8iXL5/VoWUqSY5CCEtcv27WQFyyxPy8eROKFjVNpS+9BM89BxZ3Us6xgoKCCAwMxM3NjalTp7J7926aN2/Ot99+yxNPPGF1eFlCkqMQIstcvgzLl5uEuH69GaBfqpRZHPill8ySTw7UpyNbCgoKokWLFoSHh6O1plixYvz000+89tprOWoMuMMnR+mtKoS1Llwwzw6XLDHPEqOjTa9SX1+zYHCjRmYsorCe1ppJkyZx+/ZtwKwH2xQSELcAACAASURBVLt3bzp37mxxZFnP4ZOj9FYVIuv984+Zx3TJEti2zQzDqFgRBgwwCdHDQ2apsTcnT56kV69erFmzBqUUTk5OuLq60qpVK6tDs4TDJ0chROYJCrozlrBECdPDdMkS2LXLvF+9OgwbZhJi9eqSEO1RVFQU48aN49NPP8XZ2ZkJEyZQt25dtm3bhpeXF545dDyMJEchRIbs2GHGG0ZEmG2tzU8PDxg92jxDfPxx6+IT9xcUFISPjw8HDhygffv2TJo0KWHMoqPMdJNRkhyFEGkWGws7d5oa4syZEB5+5702beDbb+GRR6yLT6TNlStXGDJkCNOmTcPd3Z2AgADat29vdVh2RZKjEOKeoqLMck/xg/LPnQMXF6hbF/74wyRMV1cYOlQSo73TWrNw4UL69u3Lf//9R9++fRk5ciQFCxa0OjS7I8lRCHGX27dhwwaTEJcvN0Mw8uaF1q3NHKZt2pgxiYmfOebQR1PZxsmTJ+nZsydr166lXr16rF69mrp161odlt1y+OQoQzmESJurV81g/KVLYc0aMyi/SBFo184kxJYtIfmkKDJ/qf1L3uFm4sSJ9OrVC2cZP3NPDp8cZSiHEKkLDTU1w6VLYeNG04RaujS88YZJiF5epslUZE9BQUF4e3tz8OBBOnTowKRJkyhbtqzVYWULDp8chRBJhYSYMYgBAWYMYmysWf/Qz8/0MG3YELLp4u0izuXLlxM63Dz88MP8/PPPvPDCC1aHla1IchQiBwgONrXDgADYs8eU1agBH39sEmLNmjIG0RForVmwYAH9+vXjv//+4/3332fkyJEUKFDA6tCyHUmOQjggrU0SXLrUvI4cMeUNG8KYMabJtFIla2MUtnXixAl69uzJunXrpMONDUhyFMJBxMSYZtL4GuLp02bOUi8v02T64ovg7m51lMLWoqKiGDt2LJ9++ikuLi5MmjSJnj17SoebByTJUYhsLCLCdKQJCICff4awMLMwcMuWMGoUtG0LxYtbHaXILNu3b8fHx4dDhw7x0ksvMXHiROlwYyOSHIXIZq5fN0MtAgJg1SqzXajQnYWBW7UyCwMLx3X58mUGDx7M9OnTpcNNJpHkKIQdix9kX7s2/PuvaTLdsMHUGEuWhFdfNQmxeXNTYxSOTWvN/Pnz6devH2FhYfTv359PP/1UOtxkgmybHJVS7YE2QClgitZ6vcUhCWFTy5bBK6+YsYfxHnkEevQwHWoaN5Z1EHOS48eP07NnT9avX8+TTz7J2rVrqVOnjtVhOSxLkqNS6nugLRCqta6eqLwVMBFwBmZorb9I7Rxa62XAMqVUUeBrQJKjyPaCg++MQdy9+065UtCzJ0yeLEMucprIyEi+/vprRo0aRa5cuZg8eTI9evSQDjeZzKqa42zgG+B/8QVKKWdgCvAscAbYpZRajkmUo5Md/47WOjTuvz+OO06IbEdrs/ZhfEKMH3JRvz74+sLs2abm6OoKXbtKYsxptm3bho+PD4cPH6Zjx45MnDgRd+lynCWUjl+ELasvrFQ5YGV8zVEp5QmM0Fq3jNseAqC1Tp4Y449XwBfABq31xlT28Qa8Adzc3OrNnz/fxp9CZIUbN2441DOVmBjFvn2F2bq1JNu2lSAsLDdOTprata/QpMl/NGlykZIlzSKJhw4V4s8/i1C79hWqVbtmceTWc7R7ITXXrl1j+vTprFq1Cjc3N/z8/GjUqJHVYdkVW90LzZo126O19khebk/PHN2B04m2zwAN7rF/H+AZoLBSqqLW2j/5Dlrr6cB0AA8PD+3l5WW7aEWWCQwMJLt/d7dvw/r1pna4YgVcumRWuWjZ0jw/bNtWUaxYUaBokuOy+ce2OUe4F+5Fa828efN4//33uXjxIgMHDmT48OE54h8E6ZXZ94I9JceUGoxSrdZqrScBk+57UlmVQ1jk8mUz1CIgANauhVu3kq5y8dxzkD+/1VEKe3H8+HF69OjBhg0bqF+/PuvWraN27dpWh5Vj2VNyPAM8nGi7LHDuQU8qq3KIrHTunBmMHxAAmzdDdDSUKQNvvXVnlYtcuayOUtiT5B1uvvnmG3x9faXDjcXsKTnuAioppcoDZ4HXgC4PelKpOYrMduzYnQ41v/1myipVgv79TUKsX19WuRB3CwoKYvbs2axfv55Tp07RqVMnJkyYIB1u7IRVQznmAV5ACaXUGWC41nqmUqo3sA7TQ/V7rfWhB72W1ByFrWkNf/xxJyEeirtL69Y1U7Z16ABVq0rPUpG6devW0aZNG2JiYlBK8dVXXzFw4ECrwxKJWJIctdadUylfDazO4nCEuK/4Sb0DAszg/JAQUxts2hQmTID27eHRR62OUti7+A433t7exMTEAODk5ERU4pkehF2wp2bVTCHNqiKjwsPvTOq9fPmdSb2ffRaGDTMda0qWtDpKkV0knuGmSpUqnDhxgujoaFxdXR26B2525fDJUZpVRXpcvQqrV5uEuGYN3LhxZ1LvDh3MpN4FC1odpchOUltSaufOnQnDETw9Pa0OUyTj8MlRao7iflatghkz4Px52LvXzEhTqhR07mwSokzqLTIqKCgIb29vDh48SIcOHZg0aVLCklKenp6SFO2YwydHqTmKlJw4YWqHP/wABw6YMqXMKhe9eoGnp0zqLTLu6tWrDBkyBH9/f9zd3Vm2bBkvvvii1WGJdHD45CgEmB6m+/bd6WEanxBLlzZJUWvTwaZmTWjSxNpYRfaltWbJkiX4+flx4cIF/Pz8GDVqFAWlLT7bkeQoHFZMDGzffqeH6alTJhE2bgxjx5oephcuQIsWEBlpJveWfhEio0JCQujVqxerVq2iTp06LF++HA+Pu6bsFNmEwydHeeaYs4SHw6ZNd3qY/vefSXrPPANDh8ILL5jnifEqVDD7BwaaxCiPgER6RUdHM2nSJD755BMAxo4di5+fHy4uDv/n1aE5/LcnzxwdX0o9TAsWTNrDtFCh1I/39JSkKDJm9+7deHt788cff9CmTRumTJnCozLg1SE4fHIUjunff80cpsuWmZpfVBS4uUGXLqa5VHqYisx0/fp1PvnkEyZPnoybmxuLFi2iY8eOKJkWyWFIchTZxt9/m2QYEABBQaYTTYUK4OdnaogNG0oPU5H5fv75Z3r37s3Zs2fx9fVl9OjRFC5c2OqwhI05fHKUZ47Zl9bw55/w/ffl6NMHDh405bVrw4gRJiFWry5zmIqscebMGfr06cOyZcuoXr06CxculHGKDszhk6M8c8xeUp7D9FGaNIHx402TablyVkcpcpKYmBi+/fZbhg4dSlRUFKNHj2bAgAHkkrXHHJrDJ0dh/8LDYcMGkwxTmsO0WLEdtG/f2OowRQ70559/4uPjw86dO3nuueeYOnUqFSpUsDoskQUkOQpLXL1qpm2L72F686bpUdq2rakdJp7DNDBQViwQWevmzZuMGDGC8ePHU7x4cebOnUvnzp2lw00OIslRZImgIFixAqKjzUw1mzebHqalS8Prr5vnh82amTGJQlhp9erV9OzZk5CQEN577z3GjBlDsWLFrA5LZDGHT47SIcdax47BxIkwdSrExpoyd3fo188kxAYNzLRtQljt/Pnz9OvXj4ULF1KlShV+/fVXmjZtanVYwiIOnxylQ07W0tqsbBHfoebQoaTvOztDz57w0UfWxCdEcrGxsUyfPp3BgwcTHh7OyJEjGTRoELlloGyO5vDJUWS+6GjYuvVOQjx92tQGn3rK1Boffhi6dr0zf2mzZlZHLIRx8OBBfHx82LFjB82aNcPf35/HH3/c6rCEHZDkKDLk1i3TwzQgwDxLvHQJ8uSB556DkSNNx5oSJe7sL/OXCnty+/ZtRo0axVdffUXhwoWZPXs2b775pnS4EQkkOYo0u3TpTg/TdetMgixSxCTCDh2gZUvInz/lY2X+UmEvNm7ciK+vL8ePH+ett97i66+/pkTif8kJgSRHcR9nzpim0mXLTM0vJsZ0qOnWzQy5ePppkLHQIjsIDQ1lwIABzJkzh0qVKrFp0yaaN29udVjCTklyFHcJDr7z/HDXLlNWuTIMGmQSooeH9DAV2YfWmlmzZvHBBx8kTBj+0UcfkSdPHqtDE3bM4ZOjDOW4v9hYkwTjE+KRI6a8fn0YPdokxMqVrY1RiPQKCgpi0aJFbN68mT///JMmTZowbdo0qlatanVoIhtw+OQoQzlSFhkJW7aYhPjzz3DuHLi4mA4zfn7w4oum+VSI7GjLli0888wzREdHAzB48GA+//xznKTJQ6SRwydHcceNG6YjTUAArFxppnDLlw9atza1wzZtoGhRq6MU4sFs2bKFTp06JSRGZ2dnChUqJIlRpIskRwcXFmaGWgQEmKEX4eFQvDi89JJJiM8+C3nzWh2lEA/u4sWLfPDBB8yaNYsyZcqQO3duoqOjcXV1xcvLy+rwRDYjydFBBAXdGUdYpsydHqZbt5pnio88Aj4+ZshF48amCVUIR6C1Zu7cufTv359Lly7x4YcfMmzYMPbt20dgYCBeXl6y7qJIN/kT6QB27IAWLSAiwmxrbX7WqAFDh5qEWLu2LAosHM/x48fp0aMHGzZsoEGDBmzcuJGaNWsC4OnpKUlRZJgkx2wqJsbUFpctg1mzTHNpvNatYdIkkA66wlFFRUXx9ddfM3LkSHLlysU333yDr68vzs7OVocmHIQkx2wkIsJMwxYQYBYFDg01c5XWrQt79pjmU1dX+OQTSYzCcQUFBeHt7c3Bgwfp2LEjEydOxF26Vgsbk+Ro565eNYsBBwTA6tWmx2nBgvD886a5tHVrs0hw4meO0pIkHNGVK1cYP348K1asoGzZsixfvpx27dpZHZZwUNkyOSqlqgB9gRLAJq31VItDsql//zVjD5ctMzXFqCgoVQo6dzYJsXlzSL6ajsxdKhyV1prFixfj5+dHaGgoffv2ZeTIkRQsWNDq0IQDS1dyVEoVBsK11hEZvaBS6nugLRCqta6eqLwVMBFwBmZorb9I7Rxa62DAVynlBHyX0VjsybFjJhkGBMBvv5lONY89Bn37miEXDRuatRCFyElCQkLo1asXq1atok6dOowYMQIfHx+rwxI5QJpHxSqlXICLwHMPeM3ZQKtk53YGpgCtgapAZ6VUVaVUDaXUymSvUnHHvABsAzY9YDyW0No8J/z4Y6heHR5/3MxdGhEBn34KBw6YhPnVV2bohSRGkZNER0czbtw4qlatyubNmxk7diw7d+7kiSeesDo0kUMoHd/vPy07K3UW6K61Xv1AF1WqHLAyvuaolPIERmitW8ZtDwHQWo9Ow7lWaa3bpPKeN+AN4ObmVm/+/PkPEvYDi4lR7NtXmG3bSrB9ewlCQ/Pg5KSpWfMKjRuH0aRJGKVLZ7hS7rBu3LhBgQIFrA5DZJEjR44wduxYjh07RsOGDenbty+lS5cG5F4Qd9jqXmjWrNkerbVH8vL0PnOcA7wHPFByTIE7cDrR9hmgQWo7K6W8gJeA3PeKRWs9HZgO4OHhoa2YJePWLVi//s6UbYkXBe7QAdq2VZQoURQoClTK8viyg/iB3MKxXb9+nWHDhjFp0iTc3NxYtGgRHTt2TLIAsdwLIl5m3wvpTY6ngC5KqV3Az8B5IEnVU2v9fQbiSGl4eqpVWq11IBCYphNbsCrHpUsmEcYvCnz7dtoXBRYiJ1q+fDm9e/fmzJkz+Pr6Mnr0aAoXLmx1WCIHS29ynBL30x2ol8L7GshIcjwDPJxouyxwLgPnuTugLFqV4/TpO1O2bdlyZ1Hgd96RRYGFSM3Zs2fx8/Nj6dKlVK9enQULFsisNsIupDc5ls+UKGAXUEkpVR44C7wGdLHFiTOr5qg1HD58p4fpnj2mPH5R4A4doF49WRRYiJTExMTg7+/PkCFDiIqK4v/+7/8YOHAgueRfkMJOpCs5aq1DHvSCSql5gBdQQil1BhiutZ6plOoNrMMM5fhea33oQa8Ftq05bt8Oc+bAtWtmceBjx0x5gwayKLAQabV//368vb35/fffefbZZ5k6dSqPPfaY1WEJkUSGJgFQSlUHngaKYYZ3/Kq1PpiWY7XWnVMpX43tO/rYzJIl0KnTne0nn4Rvv4UXXpBFgYVIi1u3bjFy5EjGjh1L0aJFmTNnDl26dEnS4UYIe5HeSQBcMOMUO5O0E41WSv0EvK21jrFdeA/OVs2qR46YVS20NmMOO3SAHj1sE6MQjm7dunX06NGDkydP8s477/Dll19SvHhxq8MSIlXpfSI2HHgFGIZ5/pg37ucw4NW4n3ZFa71Ca+39oD3fmjUzQzCcnc3k3tKbXIj7u3DhAl26dKFVq1a4uroSGBjIzJkzJTEKu5feZtXXgVFa688TlYUAn8fNctMNk0AdjqenmedUJvcW4v5iY2P5/vvvGTRoEDdv3mT48OEMGTKE3MknBRbCTqU3OT4EBKXy3g5g6IOFY3u27K0qk3sLcX/BwcH4+PiwdetWnn76afz9/aksPdVENpPeZtVzQONU3muEjcYm2pKtmlWFEPcWHh7O8OHDqVWrFgcPHmTmzJls3rxZEqPIltJbc5wLDFVKxcb993mgNGZc4lBgjG3DE0JkB5s3b8bX15ejR4/StWtXxo0bR6lSpawOS4gMS2/NcQSwGPgUOAbcAP4GPk9UbleUUu2UUtOvXr1qdShCOJy1a9dSr149mjdvTnR0NOvXr2fOnDmSGEW2l95JAKIxc6t+DjyFGed4CdiitT6cCfE9sKyaPk6InERrzfDhwxk1ahQALi4uzJgxg2bNmlkcmRC2kebkqJRyxTSb/qS13gXYZAYbIUT28vfff9OjRw82btyYUKa15rfffpPkKBxGmptVtdaRgA9mbKMQIoeJjIzk//7v/6hRowY7d+5k4MCB5M2bF2dnZ1xdXWUpKeFQ0tsh5w+gBvBrJsSSKaxYskoIR7Njxw68vb05dOgQnTp1YuLEiTz00EO89NJLCevqyWoawpGkNzkOAOYppUKAVVrrVNdctBfyzFGIjLty5QpDhgzB39+fhx9+mBUrVtC2bduE9z09PSUpCoeU3uS4CCiMWeg4WikVStJFibXW+lFbBSeEsIbWmsWLF+Pn50doaCjvv/8+I0eOpECBAlaHJkSWSG9y3ETSZCiEcDAhISH06tWLVatWUbduXVauXEm9eimtbS6E40rvUI63MykOIYTFoqOjmThxIsOGDUMpxbhx4+jTpw8uLhla2U6IbC3NvVWVUq5Kqb1KqecyMyBbk0kAhLi/3bt3U79+fQYOHEjz5s05fPgw77//viRGkWOldyhHeSA688KxPZlbVYjUXb9+nX79+tGgQQP+/fdfFi1axPLly3nkkUesDk0IS6V3+rgNQLaqOQohUrZ8+XKqVavGpEmT8PHxITg4mE6dOqGUuv/BQji49LaZTAbmKKVcgGWYiceTdNDRWp+wUWxCiExw9uxZ/Pz8WLp0KdWrV2fBggUyHEOIZNKbHLfE/ewPvJ/KPs4ZD0cIkVliYmLw9/dnyJAhREVFMXr0aAYMGECuXLmsDk0Iu5Pe5NgtU6IQQmSq/fv34+3tze+//86zzz7L1KlTeeyxx6wOSwi7ld6hHD+k9p5SyhkzQYAQwk7cunWLkSNHMnbsWIoWLcqcOXPo0qWLPFcU4j7u2yFHKXVJKVU30bZSSi1XSlVItqsH8J+tA3xQMpRD5FTr16+nevXqjBkzhjfffJPg4GC6du0qiVGINEhLb9UiJK1hOgFt48rtngzlEDlNaGgoXbt2pWXLlri6uhIYGMjMmTMpXry41aEJkW2kdyiHEMJOxcbGMnPmTCpXrszixYsZPnw4+/bt4+mnn7Y6NCGyHZn+QggHEBwcjI+PD1u3buWpp55i2rRpVK5c2eqwhMi2pOYoRDYWHh7O8OHDqVWrFgcPHmTGjBls3rxZEqMQDyitNUf3RB1wnBOVXUm0T1nbhSWEuJ/AwEB8fHw4evQoXbp0Yfz48ZQqVcrqsIRwCGlNjotTKFuWbFshy1kJkekuXrzIBx98wKxZsyhfvjxr166lZcuWVoclhENJS3KUgf9C2AGtNXPnzqV///5cvnyZwYMH88knn5AvXz6rQxPC4dw3Od5r4L8QInMFBQURGBhIxYoV+e6779iwYQMNGzZk+vTp1KhRw+rwhHBY0ltVCDsVFBREixYtCA8PR2tNvnz5mDJlCj4+Pjg7yxTGQmSmbNtbVSmVXym1RynV1upYhMgM//vf/7h9+zZam0f5fn5+9OzZUxKjEFkgy5OjUup7pVSoUupgsvJWSqkjSqm/lVKD03CqD4GFmROlENa5evUqPXv2xN/fH6UUTk5O5M2blxdeeMHq0ITIMaxoVp0NfAP8L74gbtLyKcCzwBlgl1JqOWbYyOhkx78D1AQOA3myIF4hsoTWmiVLluDn58eFCxfo168fbdq0YdeuXXh5ecmai0JkIRXfZJOlF1WqHLBSa109btsTGKG1bhm3PQRAa508McYf/zmQH6gK3AY6aK1jU9jPG/AGcHNzqzd//nybfxaR+W7cuEGBAgWsDiNT/fvvv0ycOJHffvuNSpUqMWDAAJ544gmrw7I7OeFeEGljq3uhWbNme7TWHsnL7aVDjjtwOtH2GaBBajtrrYcCKKXeBsJSSoxx+00HpgN4eHhoLy8vG4UrslJgYCCO+t1FR0czadIkhg0bhtaasWPH4ufnh4uLvfyvaV8c+V4Q6ZPZ94K9/B+Y0ho6963Saq1n3/fESrUD2lWsWDEDYQmRefbs2YO3tzd79+6lTZs2TJkyhUcffdTqsIQQ2E9v1TPAw4m2ywLnbHFiWbJK2JsbN27w/vvvU79+fc6dO8fChQtZsWKFJEYh7Ii9JMddQCWlVHmllCvwGrDcFieWxY6FPVmxYgVVq1ZlwoQJeHt7ExwczMsvvywLEAthZ6wYyjEPCAKeUEqdUUq9q7WOBnoD64BgYKHW+pAtric1R2EPzp07R6dOnXjhhRcoVKgQ27dvZ+rUqRQpki3WDBcix8nyZ45a686plK8GVtv6evLMUVgpJiaGadOmMWTIECIjI/n8888ZOHAgrq6uVocmhLgHe2lWzTRScxRW2b9/P40bN6ZXr17Ur1+fAwcO8NFHH0liFCIbcPjkKERWu3XrFoMHD6ZevXocP36cH3/8kfXr1yOtF0JkH/YylCPTSLOqyErr16+nR48enDhxgm7duvHVV19RvHhxq8MSQqSTw9ccpVlVZIXQ0FC6du1Ky5YtcXZ25pdffuH777+XxChENuXwyVGIzKS1ZubMmVSuXJlFixYxbNgw9u/fT7NmzawOTQjxAKRZVYgM+uuvv/Dx8eHXX3+ladOmTJs2jSpVqlgdlhDCBhy+5ijNqsLWwsPDGTFiBLVq1WL//v189913BAYGSmIUwoE4fM1RCFsKDAzEx8eHo0eP0rlzZ8aPH4+bm5vVYQkhbMzha45C2MLFixd55513aNasGVFRUaxZs4affvpJEqMQDsrhk6PMrSoehNaaOXPmUKVKFf73v/8xaNAgDh48SKtWrawOTQiRiRw+OcozR5FRx48fp2XLlrzxxhuUL1+ePXv2MGbMGPLly2d1aEKITObwyVGI9IqKimL06NFUr16d3377jcmTJ7Njxw5q1apldWhCiCwiHXKESCQoKAhvb28OHjxIhw4dmDx5Mu7u7laHJYTIYpIcRY4XFBTEmjVrOHjwIMuWLcPd3Z1ly5bx4osvWh2aEMIiDp8cZRIAcS87duygWbNmREZGAvDyyy8zc+ZMChYsaHFkQggrOfwzR+mQI1ITEhLCO++8k5AYnZ2dqVOnjiRGIYTjJ0chkouOjmbcuHFUq1aNU6dOkStXLpydnXF1dcXLy8vq8IQQdsDhm1WFSGzPnj14e3uzd+9enn/+eaZMmcL58+cJDAzEy8sLT09Pq0MUQtgBSY4iR7hx4waffPIJkyZNomTJkixYsICXX34ZpRTlypWTpCiESEKSo3B4K1eupFevXvzzzz/4+PjwxRdfUKRIEavDEkLYMXnmKBzWuXPnePnll2nXrh0FChRg27Zt+Pv7S2IUQtyXwydHmVs154mNjeXbb7+lSpUqrFixgs8++4w//viDxo0bWx2aECKbcPjkKEM5cpYDBw7QuHFjevXqhYeHBwcOHGDo0KG4urpaHZoQIhtx+OQocobbt28zZMgQ6taty7Fjx/jhhx/YuHEjlSpVsjo0IUQ2JB1yRLa3YcMGfH19OXHiBG+//TZfffUVJUqUsDosIUQ2JjVHkW2Fhoby+uuv89xzz+Hs7Mwvv/zCrFmzJDEKIR6YJEeR7Wit+f7776lSpQoLFy7kk08+Yf/+/TRr1szq0IQQDkKaVUW2cuTIEXx8fNiyZQtNmjRh2rRpVK1a1eqwhBAORmqOIluIiIjg008/pWbNmuzbt4/p06ezZcsWSYxCiEwhNUdh9/bt20ePHj3466+/eO211xg/fjylS5e2OiwhhAOT5Cjs1qVLlxg0aBAzZ86kXLlyrFmzhlatWlkdlhAiB8iWzapKKS+l1FallL9SysvqeIRtaa356aefqFKlCrNnz+bVV1/l4MGDkhiFEFkmy5OjUup7pVSoUupgsvJWSqkjSqm/lVKD73MaDdwA8gBnMitWkfVOnDhBq1at6Nq1K48++ii7d+/G19eX/PnzWx2aECIHsaLmOBtIUgVQSjkDU4DWQFWgs1KqqlKqhlJqZbJXKWCr1ro18CHwaRbHLzJBVFQUY8aMoXr16gQFBTF58mSCgoKoXbu21aEJIXKgLH/mqLX+VSlVLllxfeBvrfUJAKXUfOBFrfVooO09TncZyJ3am0opb8AbwM3NjcDAwIwHLjLN4cOHGTt2LCdOnKBp06b06dOHkiVLsnXrVsCsxSjfnQC5F8QdmX0v2EuHHHfgdKLtM0CD1HZWSr0EtASKAN+ktp/WejowHcDDw0N7AoqiWgAADpJJREFUeXnZIlZhI9euXeOjjz7i22+/5aGHHmLZsmW8+OKLd+0XGBiIfHcC5F4Qd2T2vWAvyVGlUKZT21lrvRRYmqYTK9UOaFexYsUMhiZsTWtNQEAAffr04fz58/Tp04fPPvuMggULWh2aEEIA9tNb9QzwcKLtssA5W5xYlqyyL6dPn6Z9+/Z07NiRkiVL8vvvvzNx4kRJjEIIu2IvyXEXUEkpVV4p5Qq8Biy3xYllsWP7sG3bNtq2bcvjjz/Oxo0b+eqrr9i9ezdPPvmk1aEJIcRdrBjKMQ8IAp5QSp1RSr2rtY4GegPrgGBgodb6kC2uJzVH682ePZunnnqKVatWERkZyY8//sjAgQNxcbGXVn0hhEjKit6qnVMpXw2stvX15JmjdW7cuMHw4cMZP348WptHyEopjhw5YnFkQghxb/bSrJpppOZojVWrVlGtWjXGjRvHCy+8QN68eXF2dsbV1VV6Gwoh7J60awmbOn/+PH379mXRokVUrVqVrVu30qRJE4KCghK6Xnt6elodphBC3JPDJ0dpVs0asbGxTJ8+ncGDBxMeHs5nn33GBx98gKurKwCenp6SFIUQ2YY0q4oHdvDgQZo2bUqPHj2oV68eBw4cYOjQoQmJUQghshuHT44i89y+fZuhQ4dSp04djhw5wg8//MDGjRupVKmS1aEJIcQDkWZVkSGbNm3C19eXv//+m7feeouvv/6aEiVKWB2WEELYhMPXHKVZ1bb+++8/3nzzTZ555hnAJMnZs2dLYhRCOBSHT47CNrTWzJo1i8qVKzN//nw+/vhjDhw4QPPmza0OTQghbM7hm1XFgzt69Cg+Pj4EBgbSuHFjpk+fTtWqVa0OSwghMo3D1xxlbtWMi4iIYOTIkdSoUYM//viDadOm8euvv0piFEI4PIdPjvLMMWO2bt1KnTp1GD58OB06dOCvv/7C29sbJyeHv2WEEMLxk6NIn8uXL9O9e3eeeuopbt26xerVq5k/fz6lS5e2OjQhhMgykhwFYDrczJs3j8qVKzNr1iwGDhzIoUOHaN26tdWhCSFElpMOOYKTJ0/So0cP1q1bx5NPPsm6deuoXbu21WEJIYRlHL7mKB1yUhcVFcWXX35JtWrV2L59O5MmTSIoKEgSoxAix3P45CgdclL2+++/4+HhwYcffkjLli0JDg6mT58+ODs7Wx2aEEJYzuGTo0jq2rVr9O7dG09PTy5evMjSpUsJCAigbNmyVocmxP+3d/+xVdVnHMffDzWlZBkEFY1DMmnmkKKMJhekLF3IggOdkx8jCiEMGKOtAhI3dS4yCGHEQXRBixWKU0BTwChUQRBcBasLBhokKlRZA04aQ4QFUcEMvDz7o1WasxZoufeecw+fV3L/uN9zzvc+KQ/99Ht6eo5IZCgcLyHr16+nb9++VFRUMH36dPbt28fo0aPDLktEJHIUjpeAQ4cOMWrUKMaMGUOPHj145513KC8vp2vXrmGXJiISSQrHGEsmkzz++OMUFBSwdetWFi1axK5duxg0aFDYpYmIRJr+lCOm3n33XUpKSqirq2PEiBFUVFTQu3fvsMsSEckKWjnGzIkTJ7j//vsZOHAgn3zyCatXr2bTpk0KRhGRdoj9yvFSedjxjh07WLZsGVu2bOHw4cNMmzaNhQsX0r1797BLExHJOrEPR3ffAGxIJBLTwq4lXTZu3MioUaNIJpOYGRUVFdx9991hlyUikrV0WjWLnTlzhmXLljF27FiSySQAnTp14vPPPw+5MhGR7KZwzFJ79+6luLiYsrIy+vXrR15eHjk5OeTm5jJ06NCwyxMRyWoKxyzz9ddfM3v2bAoLC/nwww959tlnqaur44033mD+/PnU1NRQVFQUdpkiIlkt9r9zjJOamhrKyspoaGhg4sSJPPbYY/To0QOAoqIihaKISIpo5ZgFjh49yqRJkxg2bBjuzuuvv86qVau+C0YREUkthWOEuTsrV67khhtuoKqqiocffpj333+fYcOGhV2aiEis6bRqRO3fv5+ysjK2bdvGkCFDqKyspF+/fmGXJSJyScjKlaOZdTKzBWZWbmaTwq4nlU6dOsX8+fPp378/u3fvZunSpbz11lsKRhGRDMp4OJrZM2b2mZl9EBgfYWYfmVmDmT10nmlGAj2B00BjumrNtLfffpsBAwYwZ84cRo4cSX19PaWlpXTqlJU/w4iIZK0wvuuuAEa0HDCzHOBJ4FagABhvZgVmdpOZbQy8rgL6ADvc/fdA1t8K5tixY5SUlFBcXMzJkyd59dVXWbt2Lddcc03YpYmIXJIy/jtHd681s+sCw4OABnc/AGBma4CR7v4IcHtwDjNrBE41v0229VlmVgKUAFx99dVs3779YstPKXdn27ZtLFmyhOPHj3PnnXcyefJkunTpErlaw/TVV1/p6yGAekHOSncvROWCnJ7AoRbvG4Gbz7H/OqDczIqB2rZ2cvdKoBIgkUh4lO4cc/DgQe655x5ee+01EokElZWVFBYWhl1WJG3fvl13/RFAvSBnpbsXohKO1sqYt7Wzu58Epl7QxBF7Ksfp06dZvHgxc+fOJScnh8WLFzNjxgxycnLCLk1ERJpF5UqPRqBXi/fXAp+mYmJ33+DuJd26dUvFdBdl586dDBw4kAcffJBbbrmFffv2MWvWLAWjiEjERCUcdwHXm1lvM8sFxgGvpGJiM/uVmVUeP348FdN1yBdffMG9997L4MGDOXLkCOvWraO6uppevXqd/2AREcm4MP6UYzWwA+hjZo1mNtXdvwFmAFuAeuAFd9+bis8Le+VYXV1NQUEBS5YsYfr06dTX1zN69GjMWjuTLCIiURDG1arj2xjfBGxK9eeF9TvHxsZGZs6cSXV1Nf379+ell17i5pvPdY2RiIhERVROq6ZNpleOyWSS8vJyCgoK2LJlCwsXLqSurk7BKCKSRaJytWos7Nmzh9LSUnbu3Mnw4cOpqKggPz8/7LJERKSdYr9yzMQFOSdOnOCBBx4gkUjw8ccfU1VVxebNmxWMIiJZKvbhmO7Tqps3b+bGG2/k0UcfZcqUKdTX1zN+/HhdcCMiksViH47pcvjwYcaNG8dtt91GXl4etbW1LF++nMsvvzzs0kRE5CLFPhxTfVr1zJkzVFZW0rdvX9avX8+8efPYs2cPxcXFKZlfRETCF/twTOVp1aqqKvLz8yktLWXAgAG89957zJkzh86dO6egUhERiQpdrXqB1qxZw4QJEwDIzc1lwYIF9OnTJ+SqREQkHWK/ckyVgwcPfneRTTKZ5M033wy5IhERSZfYh2Oqfuc4dOhQ8vLyyMnJITc3V4/NERGJsdifVnX3DcCGRCIx7WLmKSoqoqam5rtniBUVFaWoQhERiZrYh2MqFRUVKRRFRC4BsT+tKiIi0l4KRxERkQCFo4iISEDswzETNx4XEZF4iX04Zvp5jiIikv1iH44iIiLtpXAUEREJUDiKiIgEmLuHXUNGmNkR4N8pmq4bkK4rfC527o4c355jzrfvxWxva9uVwNELqi7zotwLHZlDvdBx6oX0bE93L/zQ3Xv836i769XOF1AZ1bk7cnx7jjnfvhezva1tQF3Y/+bZ2AsdmUO9oF5QLzS9dFq1YzZEeO6OHN+eY86378VsT+fXNV2i3AsdmUO90HHqhfRsD6UXLpnTqpK9zKzO3RNh1yHhUy/It9LdC1o5SjaoDLsAiQz1gnwrrb2glaOIiEiAVo4iIiIBCkcREZEAhaOIiEiAwlFERCRA4ShZy8xGmdlyM3vZzH4Rdj0SHjPLN7O/m9mLYdcimWdm3zOzlc3fDyakYk6Fo4TCzJ4xs8/M7IPA+Agz+8jMGszsoXPN4e7V7j4NmAzclcZyJY1S1AsH3H1qeiuVTGpnX4wBXmz+fnBHKj5f4ShhWQGMaDlgZjnAk8CtQAEw3swKzOwmM9sYeF3V4tDZzcdJdlpB6npB4mMFF9gXwLXAoebdkqn48MtSMYlIe7l7rZldFxgeBDS4+wEAM1sDjHT3R4Dbg3OYmQF/BTa7++70VizpkopekPhpT18AjTQF5B5StOjTylGipCdnf/qDpobveY79ZwLDgLFmVpbOwiTj2tULZnaFmS0FCs3sT+kuTkLTVl+sA35tZk+RonuxauUoUWKtjLV5Cyd3fwJ4In3lSIja2wv/AfQDUvy12hfufgKYksoP0spRoqQR6NXi/bXApyHVIuFSL0hrMtYXCkeJkl3A9WbW28xygXHAKyHXJOFQL0hrMtYXCkcJhZmtBnYAfcys0cymuvs3wAxgC1APvODue8OsU9JPvSCtCbsv9FQOERGRAK0cRUREAhSOIiIiAQpHERGRAIWjiIhIgMJRREQkQOEoIiISoHAUEREJUDiKxJSZPW1mbmZ/C7sWkWyjmwCIxJCZdQEOA98HjgA9m+8uIiIXQCtHkXgaDXQFFgFXEXhorIicm8JRJJ4mAQeBP9O0cvxNy41m9iMzO21m8wLjT5nZl2aWyFypItGjcBSJGTP7AU0PgX7e3U8Da4A7zKz7t/u4ewPwNHCfmV3ZfNwc4LfAaHevy3zlItGhcBSJn4k0/d9+vvn9KqAzcFdgv3lADvBHM5sKzAUmuvs/MlWoSFTpghyRmDGzvcCX7j64xVg9cMzdhwT2XQD8AbgMmOXuT2a0WJGI0spRJEbMbCBQADwX2PQcUGRmPw6M/4umVeUOBaPIWQpHkXiZBJwG1gbGnwecFhfmmNnPgWU0PVD2p2b2k0wVKRJ1Oq0qEhNmlgt8CvzT3Ue2sn0bkA9cBxQC22laUd4H7Af2uvsvM1WvSJRdFnYBIpIytwNXAIfMbFQr2w8AQ4HfAX8BtgIz3f1M8590PGNmP3P32kwVLBJVWjmKxISZvQzccQG7OlALDHf3/zYfmwN8QCsX7YhcihSOIiIiAbogR0REJEDhKCIiEqBwFBERCVA4ioiIBCgcRUREAhSOIiIiAQpHERGRAIWjiIhIwP8ALfQUjTGxAH4AAAAASUVORK5CYII=\n",
"text/plain": [
"