{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using DSolve \n", "\n", "In this lesson, we will use the sympy package and in particular, the dsolve function. \n", "The dsolve function will take a differential equation as input, and the output will be an equation involving the unknown function x(t) but none of it derivatives. In the best case scenario, the equation will explicitly give x(t) though sometimes it might only be implicit. \n", "\n", "To begin with, we do from sympy import *. This imports all functions from sympy." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from sympy import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The important thing about dsolve is that it does not find a solution *numerically* but finds a solution *symbolically* (i.e. it \"knows\" the methods that we will be learning about in class.) Therefore, we must make an equation that tells sympy that certain things are symbols and should be treated that way. \n", "\n", "There are a few things going on, and we will go step-by-step looking at some more basic examples (that have nothing to do with differential equation) first." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Eq(y**2, 2)\n", " 2 \n", "y = 2\n" ] } ], "source": [ "#eq is an equation\n", "eq = Eq(symbols('y')**2, 2)\n", "print(eq)\n", "pprint(eq)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first line encodes the equation $y^2 = 2$ in a format that sympy can understand. The Eq (equation) class is what sympy uses to represent equations. The part involving symbols tells sympy that y is a symbolic variable (instead of a normal programming variable that contains a literal value). So the \n", "expression Eq(symbols('y')**2, 2) stores the equation $y^2 = 2$ in a way that sympy understands. Next, \n", "we issu the print command so we can see what Eq is. The pprint command means \"pretty print\" and prints the equation in a slightly more readable format. \n", "\n", "Below, we see what happens if we don't tell sympy that x is a symbol. If x=1 (as we have below), then the line eq = Eq(x**2, 2) means \"$1^2 = 2$\" and this is of course not true, so the output given when we print (using either pprint or print) is False." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "False\n", "False\n" ] } ], "source": [ "x = 1\n", "eq = Eq(x**2,2)\n", "print(eq)\n", "pprint(eq)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The package sympy has a function called solve that can solve this equation. It is pretty straight forward. The solve function returns the list of solutions to the equation:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-sqrt(2), sqrt(2)]\n", "[-√2, √2]\n" ] } ], "source": [ "#eq is an equation\n", "eq = Eq(symbols('y')**2, 2)\n", "\n", "# soln will be the solution\n", "soln = solve(eq)\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The block directly below solves the equation $\\sin y + \\cos(x^2) = 1$ for $y$. The expression\n", "\n", "soln = solve(eq, symbols('y'))\n", "\n", "says to solve eq for the variable y. In other words, it says that the solutions to this \n", "equation are: \n", "$$\n", "y = \\arcsin(\\cos(x^2) - 1) + \\pi\\\\\n", "y = -\\arcsin(\\cos(x^2) - 1).\n", "$$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ⎛ 2⎞ \n", "sin(y) + cos⎝x ⎠ = 1\n", "[asin(cos(x**2) - 1) + pi, -asin(cos(x**2) - 1)]\n", "⎡ ⎛ ⎛ 2⎞ ⎞ ⎛ ⎛ 2⎞ ⎞⎤\n", "⎣asin⎝cos⎝x ⎠ - 1⎠ + π, -asin⎝cos⎝x ⎠ - 1⎠⎦\n" ] } ], "source": [ "eq = Eq(sin(symbols('y')) + cos(symbols('x')**2), 1)\n", "pprint(eq)\n", "soln = solve(eq, symbols('y'))\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The block below is similar, but instead of solving for y, it solves for sin y. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ⎛ 2⎞ \n", "sin(y) + cos⎝x ⎠ = 1\n", "[1 - cos(x**2)]\n", "⎡ ⎛ 2⎞⎤\n", "⎣1 - cos⎝x ⎠⎦\n" ] } ], "source": [ "eq = Eq(sin(symbols('y')) + cos(symbols('x')**2), 1)\n", "pprint(eq)\n", "soln = solve(eq, sin(symbols('y')))\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In other words, this says: \n", "$$\n", "\\sin y = 1 - \\cos(x^2).\n", "$$\n", "\n", "To add an assumption, like the variables are real, add the \"real = True\" expression. The block below asks solve to find a real solution to the equation $y^2 = -2$. Since there isn't a solution to this equation among the real numbers, an empty solution set ([]) is returned." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[]\n", "[]\n" ] } ], "source": [ "#eq is an equation\n", "eq = Eq(symbols('y', real = True)**2, -2)\n", "\n", "# soln will be the solution\n", "soln = solve(eq)\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But if we don't specify real solutions, it will give solutions over complex numbers: " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-sqrt(2)*I, sqrt(2)*I]\n", "[-√2⋅ⅈ, √2⋅ⅈ]\n" ] } ], "source": [ "#eq is an equation\n", "eq = Eq(symbols('y')**2, -2)\n", "\n", "# soln will be the solution\n", "soln = solve(eq)\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you might have noticed already, typing \"symbols('y')\" and such gets annoying. So, just like any data type, we can store these as python variables. In the code below, we first assign to the variable x the symolic variable symbols('x') and similarly to the variable y, we assign the symbolic variable symbols('y'). That way, when the code runs, sympy interprets the line: \n", "\n", "\n", "eq = Eq(sin(y) + cos(x**2), 1)\n", " \n", "\n", "as \n", "\n", "\n", "eq = Eq(sin(symbols('y')) + cos(symbols('x')**2), 1).\n", "" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ⎛ 2⎞ \n", "sin(y) + cos⎝x ⎠ = 1\n", "[asin(cos(x**2) - 1) + pi, -asin(cos(x**2) - 1)]\n", "⎡ ⎛ ⎛ 2⎞ ⎞ ⎛ ⎛ 2⎞ ⎞⎤\n", "⎣asin⎝cos⎝x ⎠ - 1⎠ + π, -asin⎝cos⎝x ⎠ - 1⎠⎦\n" ] } ], "source": [ "x = symbols('x')\n", "y = symbols('y')\n", "eq = Eq(sin(y) + cos(x**2), 1)\n", "pprint(eq)\n", "soln = solve(eq, y)\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In principle, the names of the variables don't have to have any connection with their symbolic representation. For example, the code below does the same thing as above (notice that the output is exactly the same), but note what names we've given the symbolic variables: " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ⎛ 2⎞ \n", "sin(y) + cos⎝x ⎠ = 1\n", "[asin(cos(x**2) - 1) + pi, -asin(cos(x**2) - 1)]\n", "⎡ ⎛ ⎛ 2⎞ ⎞ ⎛ ⎛ 2⎞ ⎞⎤\n", "⎣asin⎝cos⎝x ⎠ - 1⎠ + π, -asin⎝cos⎝x ⎠ - 1⎠⎦\n" ] } ], "source": [ "rob = symbols('x')\n", "cara = symbols('y')\n", "eq = Eq(sin(cara) + cos(rob**2), 1)\n", "pprint(eq)\n", "soln = solve(eq, cara)\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It would be incredibly unwise to do so, but you can also do this (note that we are setting x to be symbols('y')): " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ⎛ 2⎞ \n", "sin(y) + cos⎝x ⎠ = 1\n", "[asin(cos(x**2) - 1) + pi, -asin(cos(x**2) - 1)]\n", "⎡ ⎛ ⎛ 2⎞ ⎞ ⎛ ⎛ 2⎞ ⎞⎤\n", "⎣asin⎝cos⎝x ⎠ - 1⎠ + π, -asin⎝cos⎝x ⎠ - 1⎠⎦\n" ] } ], "source": [ "y = symbols('x')\n", "x = symbols('y')\n", "eq = Eq(sin(x) + cos(y**2), 1)\n", "pprint(eq)\n", "soln = solve(eq, x)\n", "print(soln)\n", "pprint(soln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we get the same exact output, but the naming convention is very very bad. \n", "\n", "Now that we have some basic understanding about sympy and the Eq class, we now turn to using sympy to solve some ODEs. \n", "\n", "Below, we solve the equation: $x'(t) = x(t)$. First, we set the symolic variable t. The next line, \n", "x = Function('x') tells sympy that x represents the function whose name is x. (If, \n", "for example, we had done x = Function('z') then x would represent the function whose name is z. It's a good idea to alter the code to see what I am talking about.) \n", "\n", "The expression diff(x(t), t) is the derivative of x with respect to t. So, the line: \n", "\n", "\n", "deq = Eq(diff(x(t),t), x(t))\n", "\n", "\n", "is the sympy Eq for $x'(t) = x(t)$. Then xsoln = dsolve(deq, x(t)) says to solve the differential equation for the function x(t). If you do dsolve(deq), this will work on most examples (bascially, dsolve guesses at what we want to do) but it's probably best to be in the habit of telling dsolve exactly what we want. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " t\n", "z(t) = C₁⋅ℯ \n" ] } ], "source": [ "#deq is the differential equation\n", "t = symbols('t')\n", "x = Function('z')\n", "\n", "deq = Eq(diff(x(t),t), x(t))\n", "xsoln = dsolve(deq, x(t))\n", "pprint(xsoln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now solve a slightly more complicated equation: $x'(t) = ax(t)$:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " a⋅t\n", "x(t) = C₁⋅ℯ \n" ] } ], "source": [ "#deq is the differential equation\n", "t = symbols('t')\n", "a = symbols('a')\n", "x = Function('x')\n", "\n", "deq = Eq(diff(x(t),t), a*x(t))\n", "xsoln = dsolve(deq, x(t))\n", "pprint(xsoln)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's graph some of these solution curves." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9e3RU93nv/dlzv2tGMxrdRkIgMBbIYBPAhhjHAY5p8iY4Pj21MUmcto7dNOe0WfVpVpI6ddPGPm7fnuZtHCc57mtnxU7Smtcrx0nrpLfE6bJz0hQDxg43gUBCGt01M5Lmft3vH3u0pUECBEgIxPNZa699md/svWeDvs+zn9/ze36KqqoIgiAISxPDYt+AIAiCsHCIyAuCICxhROQFQRCWMCLygiAISxgReUEQhCWMiLwgCMISRkReEARhCSMiLwjnQVGUbkVR0oqixBVFGVMU5ReKonxKURT5uxGuG+Q/qyBcmA+rquoGlgF/DnwOeGFxb0kQ5o6IvCDMAVVVx1VV/XvgAeATiqK0L/Y9CcJcEJEXhEtAVdX9QBjYttj3IghzQUReEC6dfqB6sW9CEOaCiLwgXDqNQHSxb0IQ5oKIvCBcAoqibEIT+Z8v9r0IwlwQkReEOaAoikdRlA8BLwPfVVX1V4t9T4IwFxSpJy8Is6MoSjdQCxSAEnAM+C7wv1RVLS7irQnCnBGRFwRBWMJIuEYQBGEJIyIvCIKwhBGRFwRBWMKIyAuCICxhTIt9A9MJBAJqS0vLYt+GIAjCdcXBgwdHVVWtme2za0rkW1paOHDgwGLfhiAIwnWFoihnz/eZhGsEQRCWMCLygiAISxgReUEQhCWMiLwgCMISRkReEARhCSMiLwiCsIQRkRcEQVjCiMgLgiAsIplMhrfeeoujR48uyPmvqcFQgiAINwKqqtLT08OhQ4c4evQohUKB9vZ21q5dO+/XEpEXBEG4SsTjcd555x3efvttIpEIFouF9evXs2HDBhoaGhbkmiLygiAIC0ixWKSzs5NDhw5x8uRJVFWlubmZbdu2sWbNGiwWy4JeX0ReEARhAYhGo7z99tscPnyYeDyO0+lky5Yt3HbbbdTUzFpLbEGQjldBEJY0AwMDfOhDH9L3n376aVauXMnq1av553/+Z/34zp07icViV3StfD7PV7/6VZqamlizZg0f//jHGR8f54EHHuCxxx7jnnvuuaoCD9fYHK8bN25UpQqlIAjzyWc/+1nuvPNO7r33Xo4dO8aDDz7I/v376e/vZ+fOnZw8eRKj0ciLL75IOBzm8ccfv+RrDAwMcOjQIX71q18xMTFBMBhkw4YNmEwmfvu3f5sTJ04swC+bQlGUg6qqbpztM/HkBUG47nnrrbdYt24dmUyGZDLJ2rVrOXLkCADf//73+bVf+zUAfvjDH7Jnzx6sVivLly9n5cqV7N+/H4Ddu3fzd3/3d3O+Zjqd5q233uK5557jueee49ChQ6xatYpPfvKT/P7v/z533XUXiqKgKMr8/+BLYF5i8oqifAv4EDCsqmp7+Vg1sA9oAbqB+1VVvbJ3IUEQhFnYtGkTu3fv5otf/CLpdJqPfexjtLe309XVhc/nw2q1AtDX18cdd9yhfy8UCtHX1weAz+cjm80SiUTw+/0V53/ggQfo6OgAIJfLkU6nyWQy3HHHHdxzzz184AMf4JZbbsHhcADw6quv8oUvfIHh4WF+9KMfXY1HcF7mq+P128CzwEvTjn0e+Kmqqn+uKMrny/ufm6frCYIgVPDEE0+wadMmbDYbzzzzDKCFUabHwGcLT0/3tIPBIP39/TNE/vnnn+fw4cO8/fbbRKNRrFYrt9xyCxs2bKC+vn6Gt37fffdx33338cYbb/DHf/zH/OQnP5nPn3pJzIvIq6r6hqIoLeccvhe4u7z9IvBviMgLgrBARKNREokE+XyeTCaD0+nEbreTyWT0NqFQiN7eXn0/HA5X5KdnMhnsdjtQmfr41FNPMTo6isViwW63Y7PZ2LdvH4899hgPPfTQee/prrvu4vTp04yOjhIIBBbgV1+ceet4LYv8a9PCNWOqqnqnfR5TVdU3y/ceBR4FaG5ufs/Zs+edxUoQBOG87N69mz179tDV1cXAwADPPvusHp/v7u4G4OjRo+zdu1fveN2xYwenTp3CaDSiqiqhUEjvQD18+DCJRAKn08mtt97KbbfdNieh7uzspLW1FUVROHToEB/+8IcJh8MLGpu/UMfroufJq6r6N8DfgJZds8i3IwjCdchLL72EyWRi7969FItFtm7dyuuvv8727dtpbW2ls7OTlStXsnbtWu6//37WrFmDyWTi61//OkajkXw+zyuvvEJ9fT3f/OY3URSFVatWsWHDBlatWoXRaJzzvXz/+9/npZdewmw2Y7fb2bdv36J2vi6kJ98B3K2q6oCiKPXAv6mquvpC55AUSkEQ5ptXX32VgwcP8uSTT874bDL18d133+UHP/gBGzduZM+ePaxfvx6Px7MId3t5LJYn//fAJ4A/L69/uIDXEgRBmJX77ruPSCSi76fTaX71q19x6NAhBgcHMZlMtLW18ZGPfITPfe5zGAxLK7N8Xjx5RVH+Dq2TNQAMAX8C/AD4/4BmoAf4DVVVoxc6j3jygiAsBKqq0t3dzaFDhzh+/DiFQoG6ujo2bNjALbfcone2Xq8suCevquqD5/lox3ycXxAE4XKYmJjgnXfe4dChQ8RiMaxWK7fddhu33XbbglV9vNZY9I5XQRCE+SSTyXDy5EmOHDnCqVOnUFWVZcuWcffdd9PW1rbgVR+vNUTkBUG47kkmk3R0dHD8+HHOnDlDsVjE7Xbz3ve+l9tuu23G4KYbCRF5QRCuSyYmJjhx4gTHjh3j7NmzqKqK1+tl8+bNrFmzhsbGxiXXiXo5iMgLgnDdEI1GOX78OMePHyccDgMQCAS48847WbNmDXV1dYteEOxaQ0ReEIRrmuHhYV3YBwcHAaivr2f79u20tbVd9frs1xsi8oIgXFOoqsrAwIAu7KOjowA0NTVxzz330NbWhs83o0KKcB5E5AVBWHRKpRLhcFgX9rGxMRRFoaWlhc2bN3PzzTdfVyNQryVE5AVBWBSKxSLd3d0cP36cEydOkEgkMBqNrFixgve9733cdNNNOJ3Oxb7N6x4ReUEQrhr5fJ4zZ85w/PhxOjo6SKfTmM1mVq1aRVtbG6tWrcJmsy32bS4pROQFQVhQstksnZ2dHD9+nJMnT5LL5bBaraxevZq2tjZaW1tvuAFKVxMReUEQ5p10Os3Jkyc5duwYp0+fplAo4HA4aG9vp62tjeXLl2MyifxcDeQpC4IwLyQSCU6cOMHx48fp6uqiVCrhdrvZsGEDa9asobm5WQYnLQIi8oIgXDbj4+N6RkxPTw+qquLz+diyZQttbW00NDSIsC8yIvKCIFwSkUiE48ePc+zYMfr7+wFtAuy77rqLtrY2amtrZdTpNYSIvCAIF6RQKBAOhzlz5gwnTpxgeHgYgIaGBnbs2EFbW9uiTVItXBwReUEQKlBVleHhYc6cOcPp06c5e/Ys+XweRVFoampi165dtLW14fV6F/tWhTkgIi8IAuPj45w5c0ZfkskkAH6/n1tvvZUVK1bQ0tJy3c+gdK1RLBQYPH2K3iPv4A7UsPZ98z/Pkoi8INyAZDIZuru7dVGfrA/jcDhYsWIFra2tLF++XLz1eUYtlRg+20XvkXfoOfou4eNHyWfSAKy9e6eIvCAIl0exWNTj6mfOnCEcDqOqKiaTiZaWFjZs2MCKFSsIBoNLLhtmYGCARx55hNdeew2Ap59+mhdeeAGj0cgzzzzDrl27ANi5cyevvPLKFRc/+973vsdf/MVfAOByufjzP/0SnlKenqPv0Hv0V2QScQB8DSHWbHs/zbesJ9TWjsNTdUXXPR8i8oKwBFFVlZGRkYq4ei6XQ1EUGhoauPPOO1mxYgVNTU1LflDSV77yFR555BEAjh07xssvv8zRo0fp7+9n586dnDx5EqPRyMc//nG+8Y1v8Pjjj1/R9YI+L9986k8ZP3uGf/zxj/nYngf4zM734vbX0Pqe22luX0dT+zrc1Vens1pRVfWqXGgubNy4UT1w4MBi34YgXJdMTExUxNUTiQQA1dXVeghmqcbV33rrLR5++GH2799PsVhk8+bN7Nu3j/b2dlasWMHx48exWq08/fTTAHzhC18AYNeuXXzpS19iy5YtxGIxtm3bxpEjRy7p2qnxMXqOvkvvkXfpOfIOY0MDANg9VVQvX8nvf+XrdBx5F29t/YKlliqKclBV1Y2zfba0TbggLGGy2WxFXH1kZASYiqtPLjdCXH3Tpk3s3r2bL37xi6TTaT72sY/R3t5OV1cXPp8Pq9UKQF9fH3fccYf+vVAoRF9fHwA+n49sNkskEpkxJ+wDDzxAR0cHoMXVc5kM+Uya7e2rWe3WCqpZ7A5Ca9q57dc+RFP7egKhZv7qK19h90c+gq+u4Wo8hlkRkReE64RisUhfX58egunr66NUKmEymVi2bJmeBVNbW7vk4upz4YknnmDTpk3YbDaeeeYZQIvHT585arbIxXTvOhgM0t/fXyHy+WyG//vxL9Bz5DA9R99l6HQnqlrCZLbQcPMamteuo/mW9dQuX4nBaNS/97Of/YwXXniBn//85wvxc+eMiLwgXKOoqsro6CinT5/mzJkzdHd3k8vlAG0g0tatW2ltbSUUCmE2mxf5bhefaDRKIpEgn8+TyWRwOp3Y7XYymYzeJhQK0dvbq++Hw2EaGqa87Ewmg8Vspu/EMXqOvkPPkXcYOHmCb7+5n5F4ErPVitlmx2KzY7Ja+e/3fITb77t/xr28++67fPKTn+Qf//EfZ7wVXG0kJi8I1xDxeLwirh6PlzMxfD5aW1v1fHWHw7HId3rtsXv3bvbs2UNXVxcDAwM8++yzJJNJ1q5dS3d3NwBHjx5l79697N+/n/7+fnbs2MHJjg4ivWc5e+QdPvDx3+KLH95JKZ8DRSHYsoLm9vU0r11H481rsNgv/tx7enrYvn07L730Elu3bl3gX60hMXlBuEbJZrOcPXtWF/XJkgF2u70iri5zml6Yl156CZPJxN69eykWi2zdupXXX3+d7du309raSmdnJytXrmTt2rX8xm/8Bjevvgm1UOQTv7aD5z71EJlEnN7oGMuDAdZt/080r11PaE07dvelTzn4Z3/2Z0QiET796U8DYDKZWEznVTx5QbhKqKrK2NgYvb29hMNhwuEwg4ODely9ublZF/W6urobMq6+ELz66qv84s03+MSHP0hPeRBSMhYFwB2oobl9Pcva1/O1777Mr99/Pzt2zP+ApIVGPHlBWARyuRz9/f2Ew2Fd2CfLBZjNZhobG9m6dauery5x9flBVVXioyP0nTxO79F3iRx5h7P/vp9/7juFo8pL09p1NLevo7n9VqqCUxUzbzvReV0K/MUQT14Q5gFVVYnFYjO89Mm/r+rqakKhEE1NTYRCIYLBIMZpmRjC5ZNNpRg8fZLBzpMMdHYwcKqD1PgYAFaHk9CaWzRRX7sOf9OyJVkGWTx5QZhnstnsDC89lUoBYLFYaGxs5M477yQUChEKhXA6nYt8x0uDUrHIaO9ZBk51MNDZwWDnSSJ9vVA2pr76RlrW3UbdqtXUr1xNsGVFRVrjjYiIvCBcBFVViUajFV760NCQ7qX7/X5WrVpV4aVLPP3KUVWVeGSUwc4OBjpPMnCqg6GuTgrZLAA2t4f6lTexess26lfeRO3Km7C73It819ceIvKCcA7ZbJa+vr4KLz2d1ioFWiwWQqEQ27Zt0710SWecH3LpFIOnO8seuibskx2kRpOJYEsrt2y/h/qVmpdeVVu3JEMv842IvHBDo6oqkUikQtCHh4d1Lz0QCLB69WrdS6+pqREvfR4olYpEenvKMfSTDHZ2EAn3oqolALx19TSvXUfdytXUr7qJmmUrMEnH9GUhIi/cUGQyGd1Ln1wmvXSr1UooFOLmm2/WvfSlWMxrMYhHRxk8Ve4Y7exg6HQn+aw2EtXmdFG3ajWrbt9K/crV1K286bLy04XZEZEXliylUmlWL32SmpoaXdCbmpoIBALipc8D+UyGwTOnGDjVoWe8JKIRAAxGE8GW5ay9eyf1q1ZTv/ImvHUNEnZZQETkhSWBqqrE43GGhoYqPPXJuiU2m43GxkbWrFlDKBSisbFRvPR5oFQqEg33ah2jnR0MnupgtLdHD7tU1dYRamunfuVN1JWzXUwWyyLf9Y3Fgou8oijdQBwoAoXz5XIKwlwpFAqMjo4yODjI0NCQvp5MYQStmuCkoDc1NeH3+8VLnwcSsagu5gOdJxk6c4rcZLjL6aSu9SZu33SHHnZZqNmOhLlztTz596uqOnqVriUsIZLJ5AwxHxkZoVTSPEWj0UgwGGT16tXU1tZSV1dHXV0dNpttke/8+qZYyBPtCzPa081IT7e+ngq7GKlZtpy2bdupX3kT9atW46trQBFDes0h4RrhmqBYLBKJRCrEfHBwUJ/dCLT5Muvq6li5cqUu5tXV1TJy9ApQVZVENMJITxejPWcZOdvFaE830f4wpWIR0OLo/lATTWvXUbu8VQu7LF+B2WJd5LsX5sLVEHkV+BdFUVTgOVVV/2b6h4qiPAo8CtDc3HwVbkdYbNLp9AwxHxkZoVAoAGAwGKipqaG1tZXa2lrdQ5dRo1dGLpNmtOfsOd55F9lyPR0At7+GmmUtrHjPZgLNLdQ0t+Crb8S4xOeBXcoseO0aRVEaVFXtVxQlCPwr8Huqqr4xW1upXbO0KJVKxGKxGeGW8fFxvY3D4aCurk4X8traWgKBwJKfXHohKZWKjA0O6GI+crab0d5uxocG9TZmm52a5hYCzcuoaV5OoHkZgeYWbE7XIt65cLksau0aVVX7y+thRVFeBTYDs4q8cP2SzWYZGhqqEPOhoSHy+TygTbEWCARoampi06ZNuqi7XC5Jn7sCUhPjeohl0juP9PZQyGszSCmKAV99A7UrVtH+vp0Eli2npnkZnkBQ4uc3CAsq8oqiOAGDqqrx8vY9wJ8t5DWFhWWyJvq54ZZYLKa3sdls1NbWsmHDBl3Ma2pqpJTuFVDI5Yj09WpifraL0V4tfj5ZbRHAUeUl0NzC+ns+qIdaqkNNEju/wVloT74WeLXsqZmAv1VV9Z8W+JrCPKCqKqlUitHRUUZHRytEPVsuEAVaCd36+npuvfVWPdxSVVUl3vlloqoqEyPDjPZqYZZJ7zw20IdazigymS34m5pZfutGapa1EGjSwi5Or8weJcxkQUVeVdUzwPqFvIZwZRQKBaLRKKOjo0QikYr19AmQLRYLtbW13HLLLbqYB4NBrFbxEi8HVVVJxyeI9ofLnaFd5dj5WXLpqXz/qmAtgebl3HT7VgLl2LmvruGGL58rzB3p3boBmBwNeq6IRyIRxsbGmN757na78fv9tLe34/f7CQQC+P1+vF6vDCa6DHLpFLGBfmIDfdp6cHK7ryKrxep0EmhqYc1d7yfQ1FL20JfNaeJoQbgQIvJLiFwuN6uQRyIRcrmc3s5sNuP3+2loaGDdunUVYi6e+aVTyOcZHxqYJuZTgj5ZKncSd6AGX30jN7/3bnx1DfgaGgg0teD2ByTEJSwIIvLXGaVSifHx8VnFfGJioqKt1+vF7/fT3NxcIeQej0cE5RIplYrER0eI9fcRLYv5WNkrnxgZ0Wu1gNYB6q1roGX9Bnz1jVTXN+Krb6Cqrl46QReBgYEBHnnkEV577TUAnn76aV544QWMRiPPPPMMu3btAmDnzp288sor+HxX1rdx4sQJfuu3fotDhw7x1FNP8Yd/+IdX/BuuBBH5a5RMJjNrnDwajeqDhkArjxsIBGhpadFFPBAIUF1dLdksl4iqqiTHYowN9BOd7pEP9DE+NEBx2nO32O346hupX3Uzbdu2U13fgK++EW99g+SaX2N85Stf4ZFHHgHg2LFjvPzyyxw9epT+/n527tzJyZMnMRqNfPzjH+cb3/gGjz/++BVdr7q6mmeeeYYf/OAH83H7V4yI/CJSLBaJxWJ6SGW6mCenxWsVRcHn8xEIBGhtba0Qc6fTKV75JZJJJioEfPo6n0nr7YwmE966BqobGml9z2Z89Y3lEEsjjiqvPPdriLfeeouHH36Y/fv3UywW2bx5M/v27aO9vZ3vf//7PPnkkwD88Ic/ZM+ePVitVpYvX87KlSvZv38/W7ZsYffu3Wzbtu2KRT4YDBIMBvnRj3503jbFRIJ8OEyut5d8uI98OIx11Up8e/Zc0bVnQ0R+AVFVlXQ6zdjY2IwlEokQi8X0Qlugjf70+/3cdNNNFeEVn88nI0AvkXw2w9jggNbR2T9NyAf7SU9MjbhVFAOeYBBfXQONq9fgK3vkvvpG3IEABoNksVwPbNq0id27d/PFL36RdDrNxz72Mdrb2+nq6sLn8+l9TX19fdxxxx3690KhEH19fQD4fD6y2SyRSAS/319x/gceeICOjo4Z133sscd46KGHZhwv5XIUYjGywyPEXn65Qszz4TDFaaO+AQwuF1W7P3zFz2E2RDmukEwmQywWm1XIx8bGKnLKQQuveL1eampqaGtrqxBzmSt07hQLeeKRCBMjQ0yMDDMxOqytR4YZGxokHhmpaO/0VeOrb2Dlpjs0b7ws5FW1dTKt3BLhiSeeYNOmTdhsNp555hlAi8fX1NTobWYr4zL9jSwYDNLf3z9D5Pft21exrxaLFIaHyYfDjL36A128c31h8r1hCsPDxEaGyRoMDL7+OorZjLmxEXMohO2WdiyhEOZQCHOoCUuoEcMCji0Rkb8I2WxWF+zZxHx6Ljlo+eRerxev18uyZcvw+Xz6vtfrlYkq5kg+l2ViZJj4yDAToyOM62I+wsTIEIlYFKb/wSoKLl81nkCQUNtafA2NU+GV+gZJRbwBiEajJBIJ8vk8mUwGp9OJ3W6v+BsNhUL09vbq++FwmIaGBn0/k8lgt9tRVZXi2FjZ++7loSee4FRvL+TzqJOLqvKbvmruraoCRcFUW4slFMK5ZQvmUAj3z9/EU1vLyi98AVNw8cpI3PAin8vlKkT7XCGfnP9zErPZrAt2U1OTvj0p5na7XWK1cyCbSk1539O88Mnt6cP1ARSDAbe/Bk9NDctuuRV3IIinpoaqmlo8gSDuQACjSTzyG5lHH32UL3/5y3R1dfG5z32OZ599lptuuonu7m69ze7du9m7dy+PPfYY4TNnOHnsGG3JJNHvfJdcby99HR2U/uAxTvb3U5rWL/Y0YGxZXva+Q1iaQpgby9uhRkwNDRjOmfHKNjqCxeXCXFd3lZ7A7Cx5kc/n8zO87+lCPn02IQCTyaQLd2NjY4UX7vP5cDgcIuIXQVVVMon4rAI+XvbOM8lExXeMZjOeQA2emlpa37MZTyCIJ1hbPhbE5fPLKE/hvLz00kuYTCb27t1LsVhk69at/PRf/5VtbW20BAIc/NrXaCoU8YbD/KeSyqoqL8ZSic8Hgwz8t98D4Kiqcqvbjb2xEfPtt2MONWJpatKEvbERo2tuWVODg4Ns3LiRiYkJDAYDf/3Xf82xY8fweBZncvIFLzV8KVxuqeFMJkM4HJ5VzKdnqYA2k1BVVdWMMMrkIlURL46qqqTGx5gYGZ4RRpncnp6lAmC22vDUBLUlEKzYrgrW4vBUSVVEYc6UslkKQ0PkBwcpDA1TGBokPzRMYXCQ/PAQhcEhCqOjUCzyk3ico5kMn6mrw1xfj7kppMXEp3ni5lCIx/7sz7j33nvZsWPHYv+8S2ZRSw1fDUZHR/nud78LaBNOVFVV4fV6Wb169awiLsPzz8+kF56IRUlGIyRiURLRiOaRj47oHnmxXEJ4EpvThbsmiLeugWW33DpDzG0utxhP4aKoqkopHj9HvDXR1sV7aIji2NiM7xpcLky1tZhra7FubcVUp8XIf7sxxHfefIOb/+APUC6QpXbLLbdclwJ/MZaEJ5/L5RgYGMDr9eJ2u0XEZ0FVVXLpFIlolEQsQjIW1QQ8FiEZndyOkhyLzhBw0EZxToZTNAGfvh3EKplBwkVQi0UKkcgM8S4MD5Evi3d+aAj1nH4wAKPfj7m2FlNtLaa62vJ2HebaIKa6OkzBWoyuG3fmsCXvyVssFpYtW7bYt7Fo5DMZErHIlFBHp22XhTwRi1I4J50TwGJ34PJV46qupvHmNTi9Plw+P67qapy+alw+P06fT4bjCxeklMtRGJoS6tnEuzAyAtNGDQNgMmEOBjHV1mJtuxnX+96Hqa5SvM3BGpRzOjWFubMkRH6pUsjlSI5Fy953lOR0IY9F9OPTS9NOYrJYcVVrIl27YhWtvmpcvmqc1X5N1H2aiFtsktIpnB81n6cQjVGMjFKIRLXc8KFyKGVwkPywti5OmzRmEoPDoQl1bRDn5s36tnlSvOtqMVZXS1/MAiMivwgUC3mSYzES0UpPWw+hRLVwyrkZKKANtXf6NKEONC1j2frbNM/bN+V5u6qrsdglC0iYnVI2S3F0VAudjEYoREYpRqIUIhFNzEcj5e3IrLFvAGN1tR7/tq9bh7muFlNwWiilrm7O2SjCwiIiPw8UCwXS8QlS42OkJyZITYyRnhgnNbmMj5f3x0iNj8/qeRuMRpxezcP21TfStPaWcqikesrzrvZjc0r2j1CJqqqUkqmytx2hMDpKMRLRPO/IKMWyaE+KeSkx03mAcsel348xEMC6YgXGzZsw+QOYAn6Mfj8mvx9TTQ2mYBCDlKS+bhCRn4VSsUg6PjGLUI9NbcenxHs2jxu0ATwOTxV2TxUOTxW1K1bhqKrC4a7C4fXiqvbrXrjd7ZHXVkFncsTlpFhXeNjRyNR22SNXZ+lvATB6vRgDfkz+APa17dp2tX9KuAMBTNXVGP1+DDbbVf6VwtXghhB5tVQinYhroj0+Rmq6t6172VNLJhGvHDJfRlEM2NxuHGXRrmlZgcPjweHxakJeVTUl6lVebA6nCLcAaJklxYkJirExTbynLbq3HY1WbM/opAQwGjVRLouzdXkLRn9A88D91ZWed3X1BVMGhRuDJfE/IDkWo/Otfyc1PiXUk4KueeQTFZM6TMfm9uBwe3BUefGHmpS7SIUAACAASURBVGjy3DJNsL2aiFdpIm5zuaQq4Q2OqqqoqZQmzmWRLo2P69uVy7i+XTpnQpfpKBaL7m2bg7XY2trOCZMEMPk1YTdWyaAx4dJYEiKfiEb4yfPfALS5Mic9a199I42r1+ComgqZTAq2w1OF3e2RofI3MGo+r3nXswr07GJdHBtDnWUcwSQGl0sLkXi9GKuqsDQ1Te2fu/i0tUHmBBAWkCUh8oHmZfzON1/E7vFIkaobEDWXoxiPU4rHKcbjU6I8Pn5B8T5fByQAZjNGbxUmrxdjlRdLyzKM3vW6eM8q2lVVKFK2WLjGWBIibzSZcVX7L95QuOZQVRU1ndZEemKCYjxBKT5BcSJOKRHX1vFzjk+KeXyCUjyBek6553MxuN3TvGcfluXLy/vniHXV1LbBKSmowtJgSYi8sHioxSKlROIC4hynNBGv8LRLExMUE4myqMehWLzgNRSzGYPHg9Ht1tYuF6b6em3f7cboKa8n96d72h6PdD4KNzTyv/8GRc3nKSWTlJJJiuV1KZnSj527FBNlsU5UivYFQx5lDA5HWaRdGNweTDU1WFasKIvz1HGjx43B5a487vFITrYgXAEi8tcJaqlEKXV+EZ4p2BcWbTWXm9uFTSYMTmfZi3ZjdHswNzdhm0WUDW4Xxsnjk561yyWetCAsIvLXtwCoxSKldAY1naKUTlNKp1HL61I6TSmVppS6sAjrol1up6ZmjpKdFUXB4HTOWMw+HwanQxPsyeMOx6xtpy+KxSKxaUG4jrkhRV4tlVAzmZkCnEpTSqe0z1KTopxCTWembacpTd+fbJfJoKY0UZ+zl1xGsdvLoloWYYcTU00NhpaWCwjwlEAbp4uyTD8oCMI0loTI5/v7iX7ve7MLcCYzU4xnqVd9QRRFE2KbDYPdjsFhR7E7MNhsGGtqUBx2DHaH9pndprW1O7R2Npu+bbDbUWza93Uv2uFAkVx9QRAWiCUh8oVYjNh3v6eJ6KTglgXZ6PNibmjAYLNVivGFBNhe3p9c22ziHQuCcF2yJETevnYtN79zeLFvQxAE4ZpDimAIgrCkOXz4MLW1tfr+rl27MJvNWCwWnnrqKf14dXU1XV1dV3y9UqnE+vXrMZvN2O12vve9713xOa8EEXlBEJY0jzzyCA8//DAAf//3f88bb7xBJBLhJz/5CV/60pfIlRMlfv3Xf51HH330iq/35S9/mf7+frLZLF/72tf4nd/5nSs+55UgIi8IwnXPiy++iN1uZ2xsjOHhYWw2G6+++iqgefKf//znAfj617/OXXfdhcfj4a677qKqqopvf/vbAPzRH/0Rb7755hXfy8svv8yePXswGAx88pOfJJfLcfjw4oWTReQFQbju+cQnPsF73vMedu7cyY4dO3jf+97HfffdxxtvvIHZbMbj8QAwMDDAihUr9O9VV1fT0dEBwPLlyymVSpw6dWrG+Zubm7Hb7TOW2Tz/aDRKW1ubvu9yuXj33Xfn+yfPmQXveFUU5deArwJG4HlVVf98oa8pCMKNxz/90z9RU1OD0Wjk4MGDABw/fhyHw6G3UWeZDMgwrT6/zWbj3XffZdWqVRVtenp65nwfF7vG1WZBr6woihH4OvABYA3woKIoaxbymoIg3JicOXOGQqFAPp9nrDwBucfjoTBthq2GhgbOnDmj70ej0QpBLxQKeL3eGee+FE/e7/dz/PhxfT+RSNDe3j4vv/FyWGhPfjPQqarqGQBFUV4G7gWOzfeFfvmNR3CPHb94Q0EQliTbn97Pf3lPgO7RNHeuW86rv38bjYkciYlxjv6POwH4UE2EP9zXwX88cQcnBpKMxSJsGvw2R//HS5RKJbKZNN43/5ij//EnFef+x081A82zXPWYfu5J3l+f4Hvf+l88WnuQ/31wGBNFzD/+bxz98YXvP+5t445P/79X8ARmZ6FFvhHonbYfBm6f3kBRlEeBR0GzloIgCJfKl37QiUGBL+5uJVco8d4n/4NvvxnmN7eFcFqN/LIzxh0rfWxf42fDMg/vfWo/igKfen8TFpMW0PiHwyP4XRZslisLcPzu9ib+9WiU2/7klxgN8Cf3ts7HT7xslNniR/N2ckX5DWCXqqqfLO9/HNisqurvzdZ+48aN6oEDBxbsfgRBuPH43Oc+x//5P/+Hn//85xdsd+utt/LRj36Uz372s1fpzuYPRVEOqqq6cbbPFro3IAw0TdsPAf0LfE1BEASdv/iLv6C19eLe9Lp1665Lgb8YC+3Jm4CTwA6gD3gL2Kuq6tHZ2osnLwiCcOlcyJNf0Ji8qqoFRVH+G/DPaCmU3zqfwAuCIAjzz4Lnyauq+mPgIv3KgiAIwkIgI14FQRCWMCLygiAIS5glUU/+9NhpPv/m57Gb7PpiM9kq9mdbztfGZrRhNMhsTYIgXP8sCZE3KAbqHHWkC2kS+QQj6RHS+TTpQppMMUO6kKakli7pnFajdYYRsBlt2M12HCbH1L7Jjt089zY2ow2byYZBkZcoQRAWniUh8surlvO1HV877+eqqpIr5UjnNdFPFVKkC+kZhmByP11Mz/55Ic14Zpyh4pC+P7lcKnaTHafZidPsxGFyTG2by9umc/Yv0NZmlOkJBUGYnSUh8hdDURSsRitWo3VBzl9SS2SL2SnRz1cahlQhRaaQqTAKqXyKZD6prQtJkvkkw6lhknlte9IQzQWDYsBpqjQIDrNjTobCaXZWGByn2YnFaFmQ5yQIwtXnhhD5hcagGPRwzXxSLBVJFaYZg3xSNwj6/jSjoG+XP4tmohXt8qX8nK5rMpj0twmH2YHL7MJlceG2uPFYPLjM2vaMxTy1bTVa5e1CEK4BROSvYYwGoy6a80G+mJ+zoTjXsMQyMXomekjkE0zkJiiUChe8lslgwmPx4La4ZxoFs/uiRsNpdkq/hSDMAyLyNxBmoxmv0YuXmfWyLwVVVckUMyRyCeK5OBO5CRJ5bXvGkp/aHkmN6McuFopSUCrE/1LeJCaNi2RICYKIvHAZKIqih6dqHDWXdY58Ka8biemGIJHT3hTiubhuOCZyEyRyCQYSA5zMnSSe19qpnL/ukoKCx+rBa/VSZamiylqlbZfXXquXKtvUttfqxWPxYDfZJcwkLClE5IVFwWww47P58Nl8l/X9kloimU/OMAoTuQkmshOMZccYz44znh1nLDvGaHqU02OnGcuOkSqkzntei8EywwBUGAZrFVWWKry2qeMeiweTQf6UhGuTJfE/czQc53//z0PY3RbsLrO2dptn3Xe4LdhcZowmifdezxgUgx6mqaf+kr6bK+Z08Z9uCCb3p2+fHjuttymo5++HcFvcM4xClbVqhoGY/uYgbw3C1WBJiLzFZqJtaz3peJ50PEc8kmH47ASZeJ5SafZXeqvDhM1l1kXf7jm/QbC7zBiMYhSWChajhRpHzSWFmlRVJZlPzjAEY9kx/c1h8rNYJkbXeBfj2XES+cR5z2k2mPFavVTbqvHb/fhtfvx2/4x9v82P1+bFbDDPx88XbjAWtJ78pTLf9eRVVSWbKpCO50gnNAMwaQjS8TzpRK7iWCaR53yPw+o06aLvcFuwTRoDl7bWjmlrq9OMwSAemqD1PYxnx2cYgulGIpqOEslEiGaiRNIRMsXMrOfyWr0Vwl9tr67cnzQOdv+CjQkRrk0WrZ78YqMoCjanGZvTzFwiv2pJJZPKVxqCWQxEdCBJ+tQYmWSeWfv+FLC7zNhcFhzuaWv3tP0qC84qK84qCyaLZIEsVcwGMwF7gIA9MKf2qqqSKqSIpCOa8JcNwOT+5PpY9BiRdOS8bwous2vqreA8hmDyuMPkkLDREmZJi/ylohgUzTN3WaDeedH2pZJKJjG7IZjazxHpSxBO5MgmZ4/pWh0mnF5N8J1VVhxeq2YAvGVD4LXi8FikH+EGQFEUfeRxs+fiE9tnChn9DSCaqTQIkwaia7yLA0MHGMuOzXoOm9F20ZDRpLGoslaJQbjOEJG/AgwGBYfHgsMztzIAxWJJNwqpiRzJsRzJ8SypsSzJcW07NhgjNZ6btS/B7jbjqJoyBpOGwaFvW3F4pP/gRsJmstHgaqDB1XDRtvlSnrHM2Iw3g0kjEclEGEgOcCRyhFgmRlEtzjiHxaD1Z9Q6arXFWUvQESToCOrHAo6A9B9cQ4jIX0WMRkM5RHPheKlaUkkn8poBKIt/ctIQjGVJjWeJhBOkJnIz+xAUcLgtUwZg8q2gyqIbAqfXit1lRpF+gxsKs8E85w7nklpiPDs+I0w0kh5hKDnEcGqYI5EjvN77OtlituK7Cgp+u59aR+2U+DtrK/aDjiAOs2OhfmoFAwMDPPLII7z22msAPP3007zwwgsYjUaeeeYZdu3aBcDOnTt55ZVX8PkuL633fDz77LP89V//NadPn2ZkZIRAYG6hu/liSXe8LnVKJZV0PFdhAM59M0iOZUnHZ9asMRgUHJNvARUGYPKYtm1zmuX1XDgvqqoynh1nKDXEUEoT/+nrSYMwkZuY8V23xa17/0FHUH8r0N8SHLXzEh767Gc/y5133sm9997LsWPHePDBB9m/fz/9/f3s3LmTkydPYjQaefHFFwmHwzz++ONXdL1zefvtt/H5fNx9990cOHBgQUT+Qh2vS0Lk1WKJUrqAQQRpVoqFkhYeGs+SGpv+ZlBpHGbrMzCYFJweK26/DXe1TVtP3/bZMJolPCRcmFQ+pb8FVBiE5JRBGE2PzhjFbDVaZ4SDphuDoCNIwB7g7YNv8/DDD7N//36KxSKbN29m3759tLe3s2LFCo4fP47VauXpp58G4Atf+AIAu3bt4ktf+hJbtmwhFouxbds2jhw5siDPoKWlZVFEfkmEa/KDKYa/9jaKxYDRZ8Pks2H0WctrGyafFaPPhsFhuiGNgNFk0ES52nbBdoV8UQsPnfNmkBzLEo9m6DsZIzmWnREiclRZ8FQYAbt+Pbffhtkq2UM3Og6zg2XmZSzzLDtvm3wpTyQdqXgDmDQIQ8kh3h15l6HU0IxqqgbFQMAWoLSmxHsfei/mkpmbd97MWcdZzr51FneVG4tF6zfr6+vjjjvu0L8bCoXo6+sDwOfzkc1miUQi+P3+ims88MADdHR0zLjnxx57jIceeuiyn8vVYEmIvNFtoerDKyhGMxRiWYqxDNnucdRMZceRYjXqgl9hAKq1fYN9STyOy8ZkNuIJ2PEEzl8yuVgskYxliUcyxKMZJsrreCTDUPcEp98eoVSstAI2lxl3tQ2P34arbAw8094IrA7ppBO0PoM6Zx11zjo4T7eBqqqMZcf0N4HB5KBuDAYeGuCV//oKJWOJxP+V4PNvfp7UqRTDhWE2fW8TDa4GzvScYSgwRPFokUZXI2PZMdLFqWJ5wWCQ/v7+GSK/b9++hfzpC8qSUDWjx4L7vY0zjpfSBQqxDMVYhkJUE//J/ezpcdTcOUbAZpzh/Zuqp/YNtiXxuK4Io9FwQUNQKqmkxnPEI2lN/MuGIBHJEB1I0n0kQjFfORWjxW6aGQqqtuEJaGubS8JwgoaiKHrNo5urb674bHBwkJ8ZfobVYuVn9/2MOHF+9suf8Vf/8lfsWb2HvkQfXd4u9p/Yz4kDJwDofqebw8sP883SNwm5QnSOdPK3nX/LevN6Gl2NNLobaXQ18tDeh+bsye/atYuhoSE2btzI888/v3APY44siZj85aCqKmq6QCGWpRDNTDMAWW0dzaCeI0YGh0kTfO+U9z89LGSQsMRFUVWVdDyve/9x/U0grRuE/DlvYCaLYcoIVBgDOx6/DYfHIplCArt372bPnj10dXUxMDDAs88+SzKZZO3atXR3dwNw9OhR9u7dy7+++a8c7jzMJz7yCZ76h6foT/cTjof59p5vs+av1pBXKkNCQXtQF/zJJeQOEXKFCDqCcyprLTH5q4yiKCgOMxaHGUuja8bnqqpSSuanRF9fZ8iPpMicjM00Ak7TeUNBRq8Vg4xsRVGmxhbUtnhmfD5ZimKmEdDWw91xbaTxNAwmBZfPpvcLeAJ2vLUOvLV2qmoc0idwA/DSSy9hMpnYu3cvxWKRrVu38vrrr7N9+3ZaW1vp7Oxk5cqVrF27lvvvv58tt23BZDLxree+xQfWfwCAAwcOMHL3CK984hVG06P0JfoIx8OEE2H64n30Jfo4OHSQH3f9mJI69bdvMpiod9bPEP/JN4HvPPcd/vIv/5LBwUHWrVvHBz/4wavq4d+wnvyVoqoqpUR+KhwUy06toxkKYxkoVD5bg8s85f377ZgCdkw1dswBOwaJS8+ZXGYWIzBtPzWRq2jv9Frx1trxBh2a+AcdVAW1kJOMIl76vPrqqxw8eJAnn3zygu0+85nPsHv3bnbs2HHBdvlinsHkIL2JXvoSfboBmDQKsWysor3D5KDB1VAp/mUDEHKF5mW8gHjyC4CiKBjdFoxuCzTP4pGWzjUC5beBaIZcX4L0kQhMG9VqcJox1WjCb65xaNs1dkzVNhQZwVqBxWbC3+DC3zDzDQw0IzA+kmZsKMX4cJqx4RRjQyk6Dw1XpIkqBgWP34a3VhN93QjUOnB5rRICWiLcd999RCKRi7Zrb2+/qMCDNsNak6eJJk/TrJ+n8qkK719/I4iH+Y+B/5gxK1q1rZpGVyPvb3o/j6x7ZG4/6hIQT36RUIslCtEMhZE0hdE0hZE0+ZEUhdE0pcS0cIRBwVRt00XfHJgyADIu4NLJJPKa6JeFf2wozfiItl3ITb2CG80Gqmrsuuc/+SZQFXRgd8tzFy4PVVWJZWO6AQgnNPHvS/SxrmYdv3fb713WeZf8YKilRimVJ18Wfs0ApMiPpClE0hUhIMVmwjzp8dfYMQUc2r7fjiIDlC4JVVVJjuV08R8fTjE2rL0NTIykK2oJWewmvEG77vXrbwFBB5YbPA1XWBxE5JcIakmlOJadEv2yASiMpClOj0MraB2/AfuUESgbAIPHIl7oJVIqlohHM4wNaaI/ZQjSxGOZinLTDo9FE339DaDcB1Bjl5HBwoIhMfklgjIZuqm2YVtd+VkpW6z0+ke1Jdk9jjotDKFYDOUOX8c0I6BtSwro7BiMBqpqHFTVOFjWXjlIppAravH/4XL8v2wEut8drawZpIC72lYR/qkqe/8ev03i/8KCISK/RDBYjVgaXTPSQVVVpTSRKwt/qhz7T5PrjZN+d6TCCzV6LOXQT6UBMEon5HkxWYz4G134Z0nDzaYLWthnqLwMpxkfTnHilwMVYwFMViP+Bqd+Hn+jtm1zSsaVcOVIuOYGRs2XKETSFQZg0giomWnFykwGzLUOzPVOzPVOLPVOzPWuG74MxOUyOSBsbChFbDBJtD9JpC/BaF+iIvvH6bXib3QRCDmpbnARCLnw1jok7VOYgYRrhFlRzAbMdU7MdZWzYE0OBNNFfzhFfjBJ5niU1IEhvZ3Ra9WF31wWflO1hB4uxvQBYQ2rvPpxVdVKQoz2JYhMLuEk4RNRvR6Qwajgq3NM8/o1z9/ptUpfizArIvLCDBRFweiyYHRZsC6v0o+rqkopnic/kCA3kCRfXjIdUSiH/RVL2XBME35znQODVf6rXQxFUbS6/l4ry9ZOxf6LxRJjgyki/ZroR/oS9J8a4+T+KYNrdZhmhHuqG5xYpN7SDc+ChWsURfkS8AgwUj70R6qq/vhC35FwzfWJmi+RH5oS/UkDMD3kY/TbsNRNE/56J0afeJ9XQiaZ10M9U0uSfHYq3u8J2DThD2mDxwIhF54aOwZ521pSLGa45v9RVfV/LvA1hEVGMRuwhNxYQm79mKqqFMez5PvL4j+ordPHInpnr2Iz6l6/pSz85joHilmyfOaCzWmmYZW3MuRTUolHM4yGE0T7E4yGk0T7E3S/O6rPA2AyG6hucFLd6CIwzfO3u+c2V7FwfSHvcsKCoCgKJq8Nk9eGfc1U6KGULVZ4/fmBJKmDwyRzA+UvomX2TPP4LfVOye+fI4pB0UtBr7h1qih7IVckOpAk0jfl+Z/91SgnfjGgt3F4LLrgT3r+vnoHJjG61zULHa75TWACOAD8d1VVY7O0exR4FKC5ufk9Z8+eXZD7Ea5d1JKqVfecFurJDyQoxqYmiDY4TBXCb65zYq51oEimyRWRmsjNCPdE+5MUC1oni2JQ8AbtuujXNLupbfFgc0l657XEgo14VRTlJ0DdLB89DvwSGEV7Of8yUK+q6m9f6HwSkxemU8oUKjz+3ECC/GAKygKEQcFUY9dTOic7e40SdrgiSsUS4yPpcsgnyWhYMwDxSEZvU1Vjp3a5R1taqgiEXDKidxFZ9LIGiqK0AK+pqtp+oXYi8sLFUEsqhdH0NPFPkB9IVpR1MHosWJZ5sDS5tXWDS2r5zAO5dIGRnjhD3RMMnhlnqHuC1Lj23A0mhZomzcufFH9PwC4htqvEooi8oij1qqoOlLf/ALhdVdU9F/qOiLxwuRSTeV30c+EEuZ6JqXCPUcHS4MLS7MbS7MGyzI2xSjJ7rhRVVUnEsgx1TTDUPcFw9wTDZyf0ap42l3lK9Fs8BFs8Mop3gVgskf8OcCtauKYb+J1J0T8fIvLCfFKcyJHrmSDbGyd3doJ8X0KfzcvgtmBtLnv6zW4sjS7J6pkHSsUS0YGkJvxl8Y8OJPWMKm+to8Lb9ze6ZATvPLDo4Zq5IiIvLCRqsaTF9nviZHsmyPXEKUbLcWajgrneibXZo3v8ksc/P+TSBYbPaoI/Kf6Ts3cZTQZqml3UtlTpwu/22+S5XyIi8oJwHorxHLmeOLneCbJn4+TD8WnevhlLkwfrMk30zY0umad3HqgI83Rpsf2Rs3EK5edud08P81QRbHFjlekxL4iIvCDMEbWokh9Mkit7+tmeCYqTWSUGzdu3NLt1j99YLV7nfFAsloj2Jcve/jhDXRPEBlP657666WGeKqobnRhlWkwdEXlBuAKKibK33xPXxD8c12v0G1xmrTO32Y212Y055BZvf57IpgsMd0/F9oe6xvUa/UazgWCzm2C5U7d2uQf3eQzuwMAAjzzyCK+99hoATz/9NC+88AJGo5FnnnmGXbt2AbBz505eeeUVfD7fvP6Oj370oxw4cACz2czmzZt57rnnMJvn981ERF4Q5hG1qJIfSk6Jfk+cwmh5cmYDmOuc5SweD9YmN0aJMc8LqqoSj2QqYvsjvXGKk2Eej6WiU7d2mQeL3cRnP/tZ7rzzTu69916OHTvGgw8+yP79++nv72fnzp2cPHkSo9HIiy++SDgc5vHHH5/X+/7xj3/MBz7wAQD27t3LXXfdxe/+7u/O6zWk1LAgzCPKZEpmgwvuqAe0FM5cOYsn1xsndWiY5C+1ZDKD0zyVvtms1fiRWbguHUWZKtmwamMtoIV5IuEEP/3HN3j8qd/ljx/6GzoO9vKXr/5XHv5Pf8z6W9fxt995mUc/+gfks0V++MMfsmfPHqxWK8uXL2flypXs37+fLVu2sHv3brZt2zbvIv/BD35Q3968eTPhcHhez38xROQFYR4wOs3Yb67GfnM1oA3ayg+ldE8/1zNB5nhUa6xMevtaCqe11YupyrqId3/9YjQaCC7z8OCnPsTR8C85mvkHEsYkD33i4+ze/n4O/ftRjCU7//K/TmAwdvDm4cPcfscd9J8ao3a5h1AoRF9fHwA+n49sNkskEsHvr5zm8YEHHqCjo2PG9R977DEeeuihOd1rPp/nO9/5Dl/96lev/IdfAiLygrAAKAYFS7m4Grdr3n4plddz9nM9cVKHR0j+xyAApho7tlU+rCu9WFdUYZA68JfME088waZNm7DZbPziF7/AaDRS8A/SdmA5H/799fR1xPjBwRJnDg/z6l8dwmQx0HM0QlfDKMObJgg0uQkGg/T3988Q+X379l3x/X3605/mrrvuYtu2bVd8rktB/icJwlXC4DBjX12NffU0b38wSbZzjEznGMm3Bkn8oh8MaKmbK73YVnmxNLlRJJPkokSjURKJBPl8nkwmg9PpxG63k81laV7jp3mNn7tPbKCQK/KBe28h3BFj6MeDDP4qzytPH8DqMDHUG+PsO2M0+pP46h16X8qlePK7du1iaGiIjRs38vzzzwPwp3/6p4yMjPDcc88t/IM4B+l4FYRrBLVQInt2Qhf9fDgOKigWI9YVVbrom4IO6cidhd27d7Nnzx66uroYGBjg2WefJZlMsnbtWrq7uwE4evQoe/fu1Tted+zYweEDRxg8PUH4eJT/8l/fx5ce/FuMBiMOj4XG1T5CN/sIrfbhCdgv676ef/55vvWtb/HTn/4Uu/3yznExpONVEK4DFJMBW6sXW6uXql3l8M6ZcTKdY5rwn4gyjlaSwbbSq4u+0SPx/JdeegmTycTevXspFots3bqV119/ne3bt9Pa2kpnZycrV65k7dq13H///axZswaTycTXv/51PNUOPNUOJpQw77/nLn7zf9xJuCNG+ESMvo4Yp97Spln0BGya6K/20bjah3OO/Sif+tSnWLZsGVu2bAHgP//n/8wTTzyxYM/iXMSTF4TrhEIso3v52c4YpaQ2vaIp6NBEf1U5ni/z6Vbw6quvcvDgQZ588skLtvvMZz7D7t272bFjh35MVVVigyld8PtOxsimtOfuq3cSKot+w03eRS2+Jp68ICwBTD4bpk11ODfVVcbzT8VI7J+M5ytYmtzYVmmevsTz4b777iMSiVy0XXt7e4XAg5a2WV3vpLreybr3hyiVVEZ744Q7NNE//ot+fvVvYVCgpsmtif7NPupXejFfI2my4skLwhJAzU+P58fI9yW0eL51ejzfh6lGarzPJ8VCiaHuCd3THzwzTqmoYjAqWopmWfRrW6oWdFIVGfEqCDcYpVSezOlxsp0xMp1jev0dg6ccz1/lw9bqxeiRWbTmk3yuyEDnGH3lmP5ITxxV1SZPr1/l1eP5Nc1uDIb5M7Yi8oJwg1OIZsh0xsiWO3FL5biyqdahi751eZWMxJ1nsqk8nFo0RgAADTNJREFUfSfLot8RI9qfBMBiN9F4k1fP3qmud17RG5aIvCAIOmpJJT+Q1Lz8U2Nku8ehoGrx/Ga3LvqWkBvFKKGd+SQ1kSt7+VHCHTEmRrU3LLvbzNptjdy+e8VlnVc6XgVB0FEMCpZGF5ZGF+73NaHmi1o8/5SWuTPx0x74SY8ez58ciSvx/CvH4bGwalMtqzZptXcmRtN6J67JsjAxe/HkBUGooJjMkz09pqdrTs6eZayyYFtdjX2tH2urF0Wm7btmEE9eEIQ5Y3SacayrwbGuBoBCJK3l5p+KafV29g+i2IzY2/zY2/1YV/mkhv41jIi8IAgXxOS34/Lbcd1ej5ovkTkVI300QvpYhNTbwyhmA7abNQ/fdnO1FFe7xpB/DUEQ5oxiNmBf48e+xo9aLJE9M64J/tFR0r8aBaOCbZUPe7sfW5sf4yKOAhU0ROQFQbgsFKMB2yoftlU+vLtbyfVMkD4SIX1klMyJKBhOYV3hxb7Wj31tQHLyFwnpeBUEYV5RVZV8X0Lz8I+MUhhJgwKWZg/2dk3wTdW2xb7NJYV0vAqCcNVQFAVLSJvm0HPPMgrDKd3DH/9RF+M/6sLc4MTeHsDeHsAcdCz2LS9pROQFQVgwFEXBXOvEXOvEs6OZQiSte/gT/3KWiX85iylo1wR/bQBzw5WN/BRmIiIvCMJVw+S3474rhPuuEMXxrC748Z/1En+9F2O1TYvhtwe0CprzWN/lRkVEXhCERcFYZcW1tQHX1gaKiRyZ41HSR0ZJ/KKfxJt9GDwWvdPWurxKSixcJiLygiAsOkaXBWe5Vn4pU9AFP3VgiOS/D2BwmLCt0Tx820oZbXspiMgLgnBNYbCZcNwWxHFbkFKuSPZkjNQRLQ8/dWAIxWrUBl+1B/7/9u4+pqr7juP4+8tF4PIMCmLFR2hVYGBbMNa1TauktnbFubjqXKMhNs2SzrqYmrrZuCZr4pbUZHF2NlvXB7PMutV1S1c70ta1LmstlPmInY/opIIPiKKXixT87o97awERDgL3ad9XcpN77vnl3A9Hfl9+nvM75xA3ya627YsVeWNMyIqKcV2bhaPtV2k9egHvvnO0HmjEu+csMiyK2NvSfG0mpxPltpLWne0RY0xYkOgo3JPScU9KRzuUK8cv4t1/Dm9NI601jTS5hNicVOILRhCXl44r0S6+ArsYyhgT5vSq0lZ3yVfw9zf67popEDshhfg7R+L+xoiIP6TT28VQdvbCGBPWJEqIHZtM6pyJZK0sJvOp20m6fwwdzW00/ekQu1e9w+zi+2mru4SqsnbtWnJzc5k0aRIVFRXXtlNaWkpTU9Og51u6dClFRUUUFhYyf/58Ll++POjf0RsbyRtjIpKq0lbbzNM/WsHtcbfywMRvcozTPPnmGiqrKmm4cJbS0lIOHTqEy+Xi9ddfp66ujtWrVw9qjubmZpKTkwFYsWIFmZmZrFq1alC/w0byxpiIVlVVRWFhIa2trXg8HvLz86mpqSF2Ygrv7P8H3/vtMlK/nUtFzUc8PPoeGtftJnnnFSbeMo5Pd34KQFlZGZs3bx70bF8VeFXF6/UG/IpeO/FqjAl7JSUllJWV8eyzz+L1ennssccoKCigtraWtLQ03CkJMD2Bi+OU4kklJIzNomXXWYZfdFOzcQcFrdkkFY/kypUrNDY2Mnz48C7bX7BgAQcPHrzue1esWMHixYv7zFdeXs62bdvIy8tj3bp1g/ZzOzGgIi8i3wWeA6YA01T1s07rfgwsBTqAp1S1oseNGGPMIFizZg0lJSXExcWxfv16AOrr68nIyLjWRlWJTosjbW4uqXMmELs7BVdiDM3vnaD5/ROkkcixj2pIn3s34vr6QMeWLVsGlO3VV1+lo6ODZcuWsWXLFsrLywe0vf4Y6OGa/cB3gB2dPxSRPGAhkA88CPxaRCL79LYxJqjOnz/P5cuXuXTpEq2tvufSut3ua+8BsrOzOXnyJAAyzEWDt5EpS+4ia2UxSfeNobXFS2vFF9SvreTCtlq+PNsC+EbyU6dOve61adOm63LMnj2bqVOn8vjjj3f53OVysWDBArZu3TpUu6BHg3LiVUQ+BJ7+aiTvH8Wjqmv9yxXAc6r6SW/bsROvxpibVVZWxsKFC6mtraW+vp4NGzZcOz5//PhxAGpqali0aBGVlZWcOnWKWbNmcfjwYVwuF6pKdnY2/6nYxZVd52j9/DxcVWLGJ5NQnIW7sP9TMVWVo0ePkpubi6qycuVKAF544YVB/dmDcT/50cDOTst1/s+uIyJPAE8AjB07dojiGGMi2aZNm4iOjmbRokV0dHQwY8YMtm/fzsyZM8nJyeHIkSPk5uaSn5/Po48+Sl5eHtHR0bz44ou4XL7CXV1dzfTp00kqyCSpIJOOS220/PsMnqoGmt48xIW3jxJflEFCSRbDshMdnUBVVZYsWUJzczOqSlFRERs3bhzq3dFFnyN5EXkfyOph1WpV/au/zYd0Hcm/CHyiqr/3L/8O2Kaqvf4/xUbyxpjB9tZbb1FdXc3zzz/fa7vly5dTVlbGrFmzunyuqrSdaMZT2YB33zn0y6sMy4onvjiL+NszQ+I5tgMayatq6U18Zx0wptNyNnDqJrZjjDEDMm/ePBobG/tsV1BQcF2BB9+DT2LHpxA7PoWrZTm07DmLp6qBi387xsV3a3HnDyehJIvYnNSQvP/9UB2Tzwf+AEwDbgE+AG5V1Y7etmMjeWNMuGir99BS1YBn1xnU244rLZaE4izi7xxJdGpsQLMM2TF5EZkH/ArIAN4Rkd2qOltVa0Tkj8ABoB14sq8Cb4wx4SRmVAIxZTmkPDQB74FzeKpOX5uKGXdbGvHFWbinpAf93vd2WwNjjBkk7Y1ePNWnafnsNB3NbUQlDCP+jkzfydohfGB5MGbXGGPM/53o4W5SHhhPcuk4Wg810VLVwOV/+R5nGDMumYTikbgLM4iKDdxlQ1bkjTFmkEmU4J6cjntyum8q5i7/VMyth7nw9jHiizKILxnpe1j5EN/Lxoq8McYMIVdSDEn3ZpN4z2jfVMyq07Ts9hX96JHxJJQM7VRMK/LGGBMAXaZiPjKRlr1n8VSdvjYVM/GuW0j91sRB/14r8sYYE2BRcdEkThtF4rRRfNngwVPVgCttaKZdWpE3xpggGpaVQOojOUO2fXtoiDHGRDAr8sYYE8GsyBtjTASzIm+MMRHMirwxxkQwK/LGGBPBrMgbY0wEsyJvjDERLKRuNSwiZ4ETA9jECODcIMUJJMsdeOGa3XIHXjhkH6eqGT2tCKkiP1Ai8tmN7qkcyix34IVrdssdeOGcHexwjTHGRDQr8sYYE8Eircj/JtgBbpLlDrxwzW65Ay+cs0fWMXljjDFdRdpI3hhjTCdW5I0xJoKFXZEXkQdF5KCIHBGRVT2sFxFZ71+/V0TuCEbOnjjI/n1/5r0i8rGIFAUjZ3d95e7UrkREOkRkfiDz3YiT3CJyn4jsFpEaEfko0BlvxMHvSoqIvC0ie/zZy4ORszsReUVEzojI/husD8n+6SB3SPZNR1Q1bF6ACzgKTARigD1AXrc2c4B3AQGmA58GO3c/ss8A0vzvHwqF7E5yd2q3HdgGzA+H3EAqcAAY61/ODHbufmT/CfAL//sM4DwQEwLZ7wXuAPbfYH2o9s++codc33T6CreR/DTgiKoeU9U24A1gbrc2c4FN6rMTSBWRUYEO2oM+s6vqx6ra5F/cCWQHOGNPnOxzgGXAVuBMIMP1wknuRcCfVfW/AKoaTtkVSBIRARLxFfn2wMa8nqru8Ge5kZDsn33lDtG+6Ui4FfnRwMlOy3X+z/rbJhj6m2spvhFPsPWZW0RGA/OAlwKYqy9O9vdtQJqIfCgi1SKyOGDpeuck+wZgCnAK2AcsV9WrgYk3IKHaP/sjVPqmI+H2IG/p4bPuc0CdtAkGx7lE5H58v0h3D2kiZ5zk/iXwjKp2+AaWIcFJ7mjgTmAW4AY+EZGdqnpoqMP1wUn22cBuYCaQA7wnIv9U1eahDjdAodo/HQmxvulIuBX5OmBMp+VsfCOZ/rYJBke5RKQQeBl4SFUbA5StN05yFwNv+Av8CGCOiLSr6l8CE7FHTn9XzqmqB/CIyA6gCAh2kXeSvRz4ufoOEh8RkVpgMlAZmIg3LVT7Z59CsG86E+yTAv154fujdAyYwNcnpPK7tXmYrid2KoOdux/ZxwJHgBnBztuf3N3av0ZonHh1sr+nAB/428YD+4GCMMm+EXjO/34k8AUwItjZ/XnGc+MTmCHZPx3kDrm+6fQVViN5VW0XkR8CFfhmILyiqjUi8gP/+pfwze6Yg+8fpAXfiCfoHGZfAwwHfu0fFbdrkO9+5zB3yHGSW1U/F5G/A3uBq8DLqtrjFLpAcrjPfwa8JiL78BXMZ1Q16LfDFZHNwH3ACBGpA34KDIPQ7p8Ocodc33TKbmtgjDERLNxm1xhjjOkHK/LGGBPBrMgbY0wEsyJvjDERzIq8McZEMCvyxhgTwazIG2NMBPsfaUmW1pjkcwcAAAAASUVORK5CYII=\n", "text/plain": [ "