{ "cells": [ { "cell_type": "markdown", "id": "05cedd92", "metadata": {}, "source": [ "## Ordinary Differential Equations (ODEs)\n", "\n", "Differential equations describe lots of things in pretty much every modern field including physics, chemistry, engineering, finance etc....\n", "\n", "They describe the rate of change generally of some quantity be it position, momentum, energy and so on. This is expressed as a derivative of some function of that quantity.\n", "\n", "The simplest type of this is when your equation is a function of one variable. This is called an ordinary differential equation.\n", "\n", "$$\\frac{d^n}{dt^n}f(t) = g(t)$$" ] }, { "cell_type": "markdown", "id": "3a68f32f", "metadata": {}, "source": [ "The simplest type of ODE is when the derivative is a first derivative, for example\n", "\n", "$$f'(t) = t$$ \n", "\n", "Does this equation have a unique solution? No, we need an initial condition\n", "\n", "$$f(t_0) = f_0$$\n", "\n", "In general though, a good question to ask is whether a solution actually exists for any right hand side $g(t)$? Furthermore, if a solution exists is it unique?\n", "\n", "Here are some cases:\n", "\n", "$$f'(t) = \\sqrt{f(t)}, f(0) = 0$$\n", "\n", "$$f'(t) = f(t)^2, f(0) = 1$$\n", "\n", "Do these have unique solutions? Turns out the first case does not have a unique solution, the second case does not have a global solution." ] }, { "cell_type": "markdown", "id": "27d28382", "metadata": {}, "source": [ "For the first case, can solve by just good old integration\n", "\n", "$$\\int \\frac{1}{\\sqrt{f}}df = \\int t dt$$\n", "$$2\\sqrt{f(t)} = t$$\n", "$$f(t) = \\frac{t^2}{4}$$\n", "\n", "Which is a solution. However so is the solution $f(t) = 0$. So which one is correct?\n", "\n", "For the second case\n", "$$\\int \\frac{1}{f^2} df = \\int t dt$$\n", "$$-\\frac{1}{f(t)} = t - 1$$\n", "$$f(t) = \\frac{1}{1-t}$$\n", "so the solution does not exist at $t=1$. So then what do we do?" ] }, { "cell_type": "markdown", "id": "0056c453", "metadata": {}, "source": [ "We are thus interested in knowing when we can expect the ODE to be uniquely solvable. This is called well-posedness. There is a very well known theorem for well-posedness of ODEs called the Picard-Lindelof theorem" ] }, { "cell_type": "markdown", "id": "cdf9af80", "metadata": {}, "source": [ "### Picard-Lindelof Theorem\n", "\n", "If the right hand side $g(t)$ is Lipschitz continuous, then a unique solution to the ODE exists.\n", "\n", "Lipschitz continuity: A function f is Lipschitz continuous if there exists some constant $L$ such that for all $x,y$, $|f(x)-f(y)| < K|x-y|$." ] }, { "cell_type": "markdown", "id": "74147c16", "metadata": {}, "source": [ "The fact that such a theorem exists is actually pretty amazing - for context there doesn't exist anything like this for PDEs (when f is a function of more than one variable). A lot of the problems in modern mathematics are on the topic of showing well-posedness of various PDEs (such as the Navier-Stokes equations)." ] }, { "cell_type": "markdown", "id": "ced1a8dd", "metadata": {}, "source": [ "For problems in this class, we will assume that they are well-posed so that we can actually even attempt to solve them computationally.\n", "\n", "The idea for solving these computationally lies in basic calculus, namely the FTC\n", "\n", "$$f'(t) = g(t)$$\n", "$$f(t) = f(t_0) + \\int_{t_0}^t f'(t) dt$$\n", "$$f(t) = f(t_0) + \\int_{t_0}^t g(t) dt$$" ] }, { "cell_type": "markdown", "id": "f18e662f", "metadata": {}, "source": [ "If we knew the value of that integral, then we could calculate $f(t)$ exactly. Obviously we can't though in general which is where numerics come in. The goal then is to approximate that integral.\n", "\n", "The first thing you might say is that the integral is roughly approximated by $g(t_0) * (t-t_0)$ given our knowledge of left Riemann sums. This gives\n", "\n", "$$f(t) \\approx f(t_0) + g(t_0)(t-t_0)$$\n", "\n", "This is known as the forward Euler method.\n", "\n", "The idea then is to use small timesteps $h$ such that if we know the solution at $t_0$, we can then use it to get the solution at $t_0 + h$, then use that to get the solution at $t_0+2h$, and so on." ] }, { "cell_type": "code", "execution_count": 3, "id": "955fce95", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "euler (generic function with 2 methods)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using PyPlot\n", "\n", "function euler(f, y0, h, N, t0=0.0)\n", " t = t0 .+ h*(0:N)\n", " y = zeros(N+1, length(y0))\n", " \n", " y[1,:] .= y0\n", " for n = 1:N\n", " y[n+1,:] = y[n,:] + h * f(t[n], y[n,:])\n", " end\n", " \n", " return t,y\n", "end" ] }, { "cell_type": "code", "execution_count": 4, "id": "e50b02ed", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABIKUlEQVR4nO3de1xUdf4/8NeZGYaR2yh3UERUvCB5AUNRydREzdhsu5ilVnbTr2Xq1qb5/a2r23fZ+lbbZoqVWlua2Te1tIxk07xiJIKKeBcFZUYEdLjJbeb8/kBIZIAZLnPmDK/n4zF/ePiMvI+nmBefqyCKoggiIiIiiSikLoCIiIg6NoYRIiIikhTDCBEREUmKYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUiqpC7CEyWRCbm4u3N3dIQiC1OUQERGRBURRRHFxMQIDA6FQNN7/IYswkpubi6CgIKnLICIiohbIyclBt27dGv26LMKIu7s7gJqb8fDwkLgaIiIiskRRURGCgoLqPscbI4swUjs04+HhwTBCREQkM81NseAEViIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiIiISFIMI0RERCQpWWx61h6MJhEpWYXIKy6Hr7sGUSGeUCp47g0REZGtWd0zsnfvXsTFxSEwMBCCIODbb79t9j179uxBZGQkNBoNevbsidWrV7ek1jaTmKHDqLd2Ydonh/DKV+mY9skhjHprFxIzdJLWRURE1BFZHUZKS0sxaNAgfPjhhxa1z8rKwv3334+YmBikpaXhjTfewLx587B582ari20LiRk6zFl/BDpDeb3rekM55qw/wkBCRERkY1YP00yaNAmTJk2yuP3q1avRvXt3vP/++wCA/v374/Dhw3jnnXfw8MMPW/vtW8VoErFseyZEM18TAQgAlm3PxPgwfw7ZEBER2Ui7T2BNTk5GbGxsvWsTJkzA4cOHUVVVZfY9FRUVKCoqqvdqCylZhQ16RG4nAtAZypGSVdgm34+IiIia1+5hRK/Xw8/Pr941Pz8/VFdXIz8/3+x74uPjodVq615BQUFtUkteceNBpCXtiIiIqPVssrT3zqODRVE0e73W4sWLYTAY6l45OTltUoevu6ZN2xEREVHrtfvSXn9/f+j1+nrX8vLyoFKp4OXlZfY9zs7OcHZ2bvNaokI8EaDVQG8oNztvRADgr61Z5ktERES20e49I9HR0UhKSqp3befOnRg6dCicnJza+9vXo1QIWBoXBqAmeJizNC6Mk1eJiIhsyOowUlJSgvT0dKSnpwOoWbqbnp6O7OxsADVDLDNnzqxrP3v2bFy6dAkLFy7EyZMnsW7dOqxduxavvvpq29yBlSaGByBhegT8tfWHYpQKAauejMDE8ABJ6iIiIuqorB6mOXz4MMaMGVP354ULFwIAnnrqKXz22WfQ6XR1wQQAQkJCsGPHDixYsAArV65EYGAgPvjgA5sv673dxPAAjA/zR0pWIbILS7FkawaqTSJCfFwlq4mIiKijEsTa2aR2rKioCFqtFgaDAR4eHm3+9z//+WEkZV7Fy2N740+xfdv87yciIuqILP385kF5AB4YWDM088MxHWSQzYiIiBwKwwiAcf39oFYpcCG/FCd1xVKXQ0RE1KEwjABwc1ZhTF8fAMAPx3MlroaIiKhjYRi5ZfLAQADAjuN6DtUQERHZEMPILeP6+cJZpUBWfikydW1zFg4RERE1j2HkFldnFcb09QVQM5GViIiIbINh5DaTa1fVHOeqGiIiIlthGLnN2H6+0DgpcKmgDCdyOVRDRERkCwwjt3F1VmFsv5qhmu85VENERGQTDCN3mHxXzaqaH47ncqiGiIjIBhhG7jCmnw86OSmRU3gTx68YpC6HiIjI4TGM3MFFrcLY/rdW1RznUA0REVF7YxgxY/JdPKuGiIjIVhhGzBjT1xednJS4fP0mjl3mUA0REVF7Yhgxo5NaiXEcqiEiIrIJhpFGPDCQQzVERES2wDDSiHv7+sJFrcSVGzeRnnND6nKIiIgcFsNIIzROStzX3w8Az6ohIiJqTwwjTag9q2bHcR1MJg7VEBERtQeGkSaM7uMDV7USuYZypF++IXU5REREDolhpAkaJyXuC+NQDRERUXtiGGlG7QZoHKohIiJqHwwjzbinjw/cnFXQGcqRlnNd6nKIiIgcDsNIMzROSoy/NVTzPYdqiIiI2hzDiAU4VENERNR+GEYsENPHG+7OKlwtqkBqNodqiIiI2hLDiAWcVUqMH8BVNURERO2BYcRCD3ADNCIionbBMGKhUb194K5RIa+4AocvcaiGiIiorTCMWEitUiA2zB8A8MOxXImrISIichwMI1aoG6rJ0MPIoRoiIqI2wTBihZG9veGhUeFacQV+u1godTlEREQOgWHECmqVAhMG1A7VcFUNERFRW2AYsdLkW0M1P2boOFRDRETUBhhGrDSytze0nZyQX1KJX7MKpC6HiMwwmkQkny/Ad+lXkHy+gL84ENk5ldQFyI2TUoGJA/yx6XAOdhzXYUQvb6lLIqLbJGbosGx7JnSG8rprAVoNlsaFYWJ4gISVEVFj2DPSAvffGqpJzNCj2miSuBoiqpWYocOc9UfqBREA0BvKMWf9ESRmcK4XkT1iGGmBEb280NmlZqgmJYuraojsgdEkYtn2TJgbkKm9tmx7JodsiOwQw0gL1A7VAMD3x/mbFlFLteXcjpSsggY9IrcTAegM5a3+BYLzUYjaHueMtNDkgQH46rccJGbosfwPA6BSMtcRWaM1czuul1bizNXiW68SnL5ajBNXDBZ93+c+/w0DArXo4+eGUF93hPq6obefG3zcnCEIQrvVTESNE0RRtPtYX1RUBK1WC4PBAA8PD6nLAQBUG024+3/+g+tlVVj/7DCMCuVEViJL1c7tuPOHT20USJgegYnhASgur8KZqyU4e7UYp68W4+yt4HGtuKLNa9J2ckKorxtC/dzQ+1ZICfVzg7+HBoIgWFwzEf3O0s9v9oy0kEqpwMTwAGxMycYPx3MZRogsZMncjle+SoenywnoihoPHd26dEIfP3f08XNHX3839PR2w4tfHMbVogqzf7cAwM9Dg49mROJCfgnOXi3B2byaoJNdWAbDzSocvnS9wUGY7s4q9PJ1xWl9SaM1C6iZjzI+zB9KRdO9K0TUEMNIKzwwsCaMJGbosfzBcDhxqIaoWSlZhU3O7QCAimpTXRDx83D+PXT4uaOPvzt6+7rBzbnhj6+//mEA5qw/AgGoFxyEuq+HYVBQZwwK6lzvfeVVRly4VoqzecU4l1cbVIpxsaAMxRXVSM9pegjo9vko0b28mv4HIKIGGEZaYViIJ7xc1SgorcShCwWICfWRuiQiu5eVX2JRu1fGhWLWyBBoXZws/rsnhgcgYXpEg3kd/s3M69A4KREW6IGwwPrdyJXVJlwsKMXGX7Px6cGLzX7/vOKmQxYRmccw0goqpQITwv3x5a/Z+OGYjmGEqAkX80uxdn8Wvvot26L2w3t6WRVEak0MD8D4MH+kZBUir7gcvu4aRIV4tmj4RK1SoI+fO2IH+FsURgpLK63+HkTEMNJqD9wVgC9/zUbiCT3+NoVDNUR3Sr10HZ/svYCfMvWonS6vUgiobmRJrICanoyoEM8Wf0+lQmjT4ZKoEE8EaDXQG8rNzhuptWx7JvaeuYZXJ/TFgEBtm31/IkfHT85WigrxhLebGjfKqnDwPM+qIQJqJqkmZujxcMJBPJxwEIknaoLImL4++PL5YVgxbQgE/D6Xo1btn5fGhdnVRFClQsDSuDAAjdccE+oNpULA7tPXMPmD/Xh5Yxou5pfatE4iueLS3jbw398ex/pD2XhsaDe8/cggqcshkszNSiO+OXIZa/ddwMWCMgCAWqnAlCGBeC6mJ/r4ude1leOeHc3VnJVfiveSzmD70VwANSFm6t1BmDc2FP5ajVRlE0nG0s9vhpE2kHy+ANM+OQRtJyf8tuQ+qFXscCLHYTSJzc6/yC+pwOfJl/BF8kVcL6sCULNvx/Th3fFUdA/4epj/ILbk77Y3ltR8IteAd346jd2nrwEAnFUKPD2iB2aP7oUurmopyiaSBMOIDRlNIob9/Wfkl1Tg02fuxpi+vlKXRNQmmusJOH+tBGv2ZWHzkcuorK45NDLIsxOeHRmCR4cGwdXM8tuOJCWrEP/70yn8drFm7xJ3ZxVeuKcnZo0K6fD/NtQxMIzY2F++y8DnyZfwSGQ3vPMoh2pI/pracVQEMLCrFsdu24J9UDctXrinFyYM8OPxCLcRRRG/nL6GtxJP4ZS+GADg7abGS2N6Y9qw7nBWKQHIs5eIqDkMIzb264UCTP34EDw0Khz+7/EcqiFZM5pEjHprV7ObkwHAff398MI9PXF3jy7Nnu3SkZlMIrYfy8V7SWdw6dZ8mq6dO2HB+D5wcVLibz/Ia/4MkSUYRmzMaBIxPP5nXCuuwKdP340x/ThUQ/JVOw+qOe88OgiPRHazQUWOo8powteHc/Cv/5xFXhNn7PDMG3IEln5+89f3NqJUCLg/3B8A8P0xncTVELWOpTuJOinZE2ItJ6UCTw4Lxp7XxuDPE/s2WCpcq/a3xGXbM2FsZE8WIkfBMNKGJg8MBADszNSjotoocTVELefr7mxhOy5XbalOaiWGBHVpchO128+8IXJkDCNtaGhwF/i6O6O4vBr7z+ZLXQ5Ri+gMN7F6z/km2wiomdPQml1SyfIeKJ55Q46Oa8vakEIh4P67AvDZwYv47MBFlFRUc1Y8yYYoivj6cA7e/P4kiiuq67Zsb+wEXHvbJVWOLO1ZYg8UOboW9YysWrUKISEh0Gg0iIyMxL59+5psv2HDBgwaNAguLi4ICAjAM888g4ICx9w63dutZkOjfefy8cpX6Zj2ySGMemsXEjM4j4Ts15UbNzFzXQpe33wcxRXVGNK9MxLnx2D19IgGO4f6azWcVNlGas+8aSrSqRQCgr1cbFYTkRSsXk2zadMmzJgxA6tWrcLIkSPx0UcfYc2aNcjMzET37t0btN+/fz9Gjx6Nf/7zn4iLi8OVK1cwe/ZshIaGYuvWrRZ9TzmspgGa3pcB4Kx4sj+iKOLLlGzE7ziFkopqOKsUeDW2L2aNCqnr9eD+F+2r9ucGgEbnj3Tt3An/nhWF3r5utiuMqA2029LeYcOGISIiAgkJCXXX+vfvjylTpiA+Pr5B+3feeQcJCQk4f/73MegVK1bg7bffRk5OjkXfUw5hpLl9GWpPIt3/+lj+ICe7kFNYhkVbjuHAuZpeyqHBXfD2IwPR04cfeLbW2E63c8f0xtr9WcjKL0VnFyesfepuRAZ3kbBSIutY+vlt1ZyRyspKpKamYtGiRfWux8bG4uDBg2bfM2LECCxZsgQ7duzApEmTkJeXh2+++QaTJ09u9PtUVFSgouL39fdFRUXWlCmJlKzCJjeIun1WfFsebU5kLZNJxIZfLyH+x1MoqzRC46TAaxP64ekRPRiUJTIxPADjw/zN9kBNCvfHrH8fxtGcG3hyzSGsfCIC4/r7SV0yUZuyas5Ifn4+jEYj/Pzq/4/g5+cHvV5v9j0jRozAhg0bMHXqVKjVavj7+6Nz585YsWJFo98nPj4eWq227hUUFGRNmZLgrHiSg0sFpXhizSH8v+9OoKzSiKgenkh85R48e9uwDElDqRAQ3csLDw7uiuheXnXPw8vNGRufH4Z7+/qgvMqEF75IxabfsiWulqhttWgC651bPoui2Og20JmZmZg3bx7+8pe/IDU1FYmJicjKysLs2bMb/fsXL14Mg8FQ97J0OEdKnBVP9sxkEvHZgSxMfH8fDl0oRCcnJZb9YQC+emE4eni7Sl0eNcNFrcInM4fikchuMJpEvL75OFb8fBYy2ECbyCJWDdN4e3tDqVQ26AXJy8tr0FtSKz4+HiNHjsRrr70GABg4cCBcXV0RExODN998EwEBDSd0Ojs7w9nZsk2X7EXtrHi9odzsJLTaOSPcl4HaS2MTTbPyS/H6N8eQcrFm46zhPT3x9sOD0J0rNGTFSanA/z4yEH4ezli5+zzeTTqDq8XlWPaHcPZqkexZFUbUajUiIyORlJSEhx56qO56UlISHnzwQbPvKSsrg0pV/9solTWnVDpSqlcqBCyNC8Oc9Uca7MuAW3/mvgzUXsxNgPT30GBUby98f1yH8ioTXNRKLL6/P56M6g4F/zuUJUEQ8NqEfvBxc8ay7zOx/lA28osr8f7jg6FxUkpdHlGLWT1Ms3DhQqxZswbr1q3DyZMnsWDBAmRnZ9cNuyxevBgzZ86sax8XF4ctW7YgISEBFy5cwIEDBzBv3jxERUUhMDCw7e7EDkwMD0CCmX0ZgJqleePD/CWoihxd7dLQOydQ64vK8c2RKyivMmFkby/8NP8ezBgezCDiAJ4eGYIPp0VArVQg8YQeM9elwHCzSuqyiFrM6h1Yp06dioKCAixfvhw6nQ7h4eHYsWMHgoODAQA6nQ7Z2b9Prnr66adRXFyMDz/8EH/605/QuXNnjB07Fm+99Vbb3YUduXNWfCcnJRZ+nY4rN25iy5HLeHSo/U/GJfkwmkQs257Z5Pkm2k5O+PczUVApefqDI5k8MACermq88PlhpGQV4rHVyfhs1t0I0HaSujQiq1m9z4gU5LDPSFM+3nsef99xCn4eztj96r1wUXMXfmobyecLMO2TQ8222/j8cC4pd1AndUV4al0K8oorEKjV4PNno9Db113qsogAWP75zV+VbGBmdA9069IJV4sq8MneLKnLIQfCJeXUP8ADW/5rBHr6uCLXUI6HE5KReomn/JK8MIzYgMZJidcn9gMAfLT3PPKK+MFAbcPbzbJVZ1xS7ti6dXHBN7NHYHBQZxhuVuGJT35FUuZVqcsishjDiI08MDAAQ7p3RlmlEe/uPCN1OeQACksrkfDLuSbbCKjZVpxLyh2fp6saXz4/DGP7+aKi2oQXvziMjSk18/eMJhHJ5wvwXfoVJJ8vgNFk96Pz1MFwzogNpV4qxMMJyRAEYMe8GPQPkO+9kLSO5tzAf204gis3bsJJKaDKKDZYUs4DGjumaqMJb2w9jq8PXwZQ84vQ4UvXob/j3JulcWH874LaHeeM2KHIYE9MvisAogj8fcdJqcshGRLFmnNlHl2djCs3bqKHlwu2vzwKq80sKffXahhEOiCVUoG3Hh6Il8b0BgB8f0xXL4gAgN5QjjnrjyAxQydFiUQNsGfExrILyjDuvV9QZRTx2TN3496+vlKXRDJxs9KIJd8ex5YjVwAAsWF+eOexQfDQOAFofAdW6piMJhFDlu9EUXm12a/zJHGyhXY5tZdar7uXC56K7oE1+7Pw9x0nMaq3N/d/oGZdzC/F7PWpOKUvhkIA/jyxH168p2e9M6FqD1ojAmpOEm8siAA8SZzsCz8FJfDy2FB0dnHCmasldeO6RI3ZeUKPuBX7cUpfDG83NTY8NxyzR/dq9HBKIoDLvkleGEYkoHVxwryxoQCA95JOo6Si8d9eqOOqNprwVuIpvPBFKoorqhEZ3AXfvxzD32LJIjxJnOSEYUQi04cHo4eXC/JLKrH6l/NSl0N2Jr+kAjPXpSDh1n8bs0aG4KsXhps994jInNqTxJvqP+Oyb7IXDCMSUasUWDSpPwDgk30XoDPclLgishepl67jgQ/24+D5AriolVgxbQj+EhcGJ84tIivUniQOoNFAMn9cKCevkl3gTzcJTRjgh6genqioNuF/fzotdTkkMVEU8dmBLEz9KBn6onL08nHFd3NHIm6QY51uTbbT2EniqlsB5OvUy6isNklRGlE9XNorsaM5N/DgygMAgO9fHoXwrlqJKyIplFZUY/GW49h2NBcAMPmuALz1yEC4OXPBG7Xencu+fd2dMWXVARSXV2P68O54c8pdUpdIDoqbnsnEoKDOeHBwzW++b/6QCRlkQ2oFc9tyn8srwZSVB7DtaC5UCgH/74EwfPjEEAYRajO1y74fHNwV0b280MvXDf96fDAEAVh/KBtf/5YjdYnUwfGnnR14bUJf/Jihx6ELhfjPyTyMD/OTuiRqB4kZOizbngndbbthdnZxws1KIyqqTfB1d8bKJyNwdw9OKKT2N7afHxbc1wfvJZ3Bf3+bgT7+7hgc1FnqsqiDYs+IHejWxQXPjgoBAMT/eBJVRo7hOprEDB3mrD9SL4gAwI2yKlRUmxDq64bv541iECGbemlMb4wP80Ol0YQ561ORX1IhdUnUQTGM2In/urcXvFzVuHCttO6kTXIMRpOIZdsz0dQAXElFNbxcnW1WExEAKBQC3ntsEHr6uEJnKMfcDUf4yxBJgmHETrhrnDD/vpqN0N7/z1kUlVdJXBG1lZSswgY9Ineq3ZabyNbcNU74eMZQuDmr8GtWIQ/xJEkwjNiRaVHd0cvHFYWllVi5+5zU5VAb4bbcZO96+7rh3ccGAQA+PXARW9N4TAXZFsOIHVEpFXjj/pqN0D7dfxE5hWUSV0RtwdJVMdyWm6Q0YYA/Xh7bGwCwaPNxZFwxSFwRdSQMI3ZmbD9fjOjlhUqjCW9zIzTZy8wtwrLtJ5psI4DbcpN9mH9fH4zp64OKahNe/CIVhaWVUpdEHQTDiJ0RBAFLJveHIADbj+YiLfu61CVRC/3f4Rw8tOoAsgtvwtNFDaDhtty1f14aF8ZtuUlySoWA96cOQbCXC67cuImXNx5BNSe0kg0wjNihAYFaPBzRDQDwPz+c5EZoMlNeZcSizcfw2jfHUFFtwr19ffDzn0ZjtZltuf21GiRMj8DE8ACJqiWqT+tSM6HVRa3EgXMFPKqCbILbwdspvaEc976zG+VVJiQ8GYFJd/HDSg6yC8owZ0MqTuQWQRCAhff1wdwxvaG41etx57bcUSGe7BEhu/T9sVy89GUaAGDFtCE8I4lahNvBy5y/VoMXYnoCAP6ReIqHWcnAfzKv4oEV+3Aitwiermp8PisKL48LrQsiQMNtuRlEyF49MDAQL46u+Rn052+O4ZS+SOKKyJExjNixF0f3go+7My4VlOHz5ItSl0ONqDaa8HbiKTz3+WEUlVdjSPfO+P7lUYgJ9ZG6NKJWeS22L0b19sbNKiNe+DwVhjLuf0Ttg2HEjrk6q/Cn8X0AACt2ncONMs5stzfXiiswY20KVv1yHgDw9Ige2PRCNAI7d5K4MqLWUykVWDFtCLp16YTswjK8sikNRpPdj+yTDDGM2LlHhwahr587DDersGIXN0KzJ79dLMTkD/Yh+UIBXNRKrJg2BH/9wwCoVfzfihxHF1c1Vk+PhLNKgV9OX8M/k85IXRI5IP7UtHNKhYA3JtdshPbvg1nYeuRyvePnyfZEUcSafRfw+MeHkFdcgd6+btj20khO8COHFd5Vi388fBcA4MPd55CYoZe4InI0lm0NSZIa3ccHYQEeyNQVYcHXR+uuB2g1WBoXxmWhNlRcXoXXNx/DjuM1P4z/MCgQ8X+8C64W7rJKJFcPDemG45eLsO5AFv70dTp6+45Eb193qcsiB8GeERlIzNAhU9dwJrveUI45648gMUMnQVWOzWgSkXy+oF4v1Cl9ER788AB2HNfDSSlg+YMD8K/HBzOIUIex+P5+GBbiidLKmgmtPNCT2gr3GbFzRpOIUW/tavTUVwE1y4D3vz6Wy0TbSGKGDsu2Z9b7N9d2ckJZZTWqjCICtRqsfDICQ7p3kbBKImnkl1QgbsV+6AzluK+/HxKejMDhS9e5dw6ZZennN3+ls3PNHT8v4vfj56N7edmuMAeVmKHDnPVHcGdCN9ys+Q2wf4A7Njw3HJ6uatsXR2QHvN2csXp6JB79KBn/OXkVQ/6WhJKK6rqvc/iYWoLDNHaOx8/bjtEkYtn2zAZB5HY3yqqg7eRks5qI7NGgoM6YencQANQLIgCHj6llGEbsnKXHyvP4+dZrrhcK+L0XiqgjM5pE/Cfzqtmv1Yb5ZdszueKPLMYwYueiQjwRoNU0OO31dk5KAb18XG1Wk6NiLxSRZawZPiayBMOInVMqBCyNCwPQ8Pj5WlVGEY+sTsb5ayW2K8zBVFQbsf9svkVt2QtFHR2DO7U1hhEZmBgegAQzx88HaDX46x/C0N3TBdmFZfjjqoP49UKBRFXKV0pWISZ/sB//l3q5yXYCav7No0I8bVMYkZ3i8DG1Na6mkYmJ4QEYH+Zv9vj5BwYG4rl/H0Z6zg3MWJuC/310IB4c3FXqku3ejbJKxO84hU2HcwAA3m5qPDi4K9btzwKAehNZa3ullsaFcdkidXi1w8d6Q7nZCd+1Ww4wuJOluM+IgyivMmLBpnT8eGub5j+N74OXxvaGIPCD806iKOLb9Ct48/uTKCitOXxwWlR3LJrYD1oXJ7P7jHC5IlF9tcvgAZgNJKunR/D/F7L485thxIGYTCL+kXgKH++9AAB4NLIb/v7Hu+Ck5Ghcraz8Uvz3t8dx4FzNcFYfPzf8/aG7MLRH/d/gjCbRbC8UEf3OXHAHgM6dnLDv9TFw13AZfEfHMNKBfXHoEpZ+lwGTCIzs7YVVT0Z2+L0xKqtN+GjPeazYfQ6V1SY4qxSYNy4Uz8f05Cm7RK1we3Dv7OKE//dtBrILb+LpET3w1z8MkLo8khjDSAe3+1Qe5n55BGWVRvTxc8O6p+9Gty4uUpcliZSsQryx9TjO5dWsNooJ9cabU8IR7MXl0ERtbf/ZfExf+ysEAdgyZwSPTejgLP385q+EDmpMP198/WI0/DycceZqCR5adRDHLt+QuiybulFWide/OYbHPkrGubwSeLup8a/HB+PzWVEMIkTtZFSoN/44pCtEEVi85TiqjCapSyIZYBhxYOFdtfh27kj083fHteIKTP3oEJIa2TVRjsydrAvcmqCadgXj3t1Tt1JmWlR3/LzwXjw4uCsn9RK1syWT+6OLixNO6YuxZl+W1OWQDHCYpgMoLq/C3C/TsPfMNQgCsPSBMDw9MkTqslqlsRUvc0b3ws7Mq9h/rmYDs1BfN8T/seEEVSJqX9+kXsar/3cUzioFdi64h72RHRTnjFA9VUYT/vJdBjam1PQUzBoZgiWT+8tyhUhjJ+vejhNUiaQliiKeXPMrDp4vQEyoNz6fFcVeyQ6Ic0aoHielAn9/6C68PrEfAGDdgSzMXp+KssrqZt5pXyw5WddZpcCPr8Rg7pjeDCJEEhEEAf/z0F1QqxTYdzYf36ZfkboksmP8Sd2BCIKAOff2woppQ6BWKZCUeRWPf3wIecXljc6/sDeWnKxbUW3C1aIKG1VERI0J8XbFK+NCAQB/+/4kCm9tMkh0J24H3wHFDQpEgFaD5z8/jGOXDZj4z31QKID8kt9/UNjjjqOFpZXYdtSy3654QBeRfXg+pie2pefi9NVi/H3HSbzz6CCpSyI7xJ6RDmpoD09s+a+R8HFTo7Cssl4QAQC9oRxz1h9BYoZOogpr5N64ic8OZOHxj5Mx9M2kujkvzeEBXUT2Qa1S4O9/vAuCUDOp9eA5y07Hpo6FPSMdWHdPl0YnlImoOexq2fZMjA/zb9VEV2u3Vr9wrQSJJ/T46cRVHM25Ue9rYQHuyC68iZIK83NdeEAXkf2JDO6C6cOC8cWhS3hj63Ekzr8HGiel1GWRHWEY6cBqAkLjcytEADpDORZsSsfYfr4I9XNDLx83q36IWHLonCiKOJFbhJ9O6PHTCT3OXC2paysIwNDgLpgwwB8TBvgjyNOl0QO6eLIukf16bWJf7MzU42JBGT7cdQ6vTugrdUlkR7i0twP7Lv0KXvkq3ar3KASgh5cr+vi5o4+/O/r4uaGvnzt6eLs2OJCvsSW4tTHhT7F9cKOsCokn9Lh8/Wbd11UKASN6e2PCAD+MD/MzO+TCk3WJ5CcxQ4fZ649ApRDww7wY9PV3l7okamfcZ4SalXy+ANM+OdRsu/v6+6LoZjVOXy2G4WaV2TZOSgE9vd0Qeiuc9PJxw9JtJ3CtxLJVLRonBUb38cHEcH+M7esHrUvzB/vxZF0ieRFFES98kYqkzKuI6N4Z38weAQX/n3Voln5+c5imA4sK8USAVgO9odzsvh218y8+mjEUSoUAURRxrbgCp68W48zVEpzRF+NMXjHOXi1BSUVNWDl9tRjfw/JJr6N6e2P68GCM7uODTmrrxpCVCgHRvbyseg8RSUcQBCx/cAAOnsvHkewb2JCSjRnDg6Uui+xAi1bTrFq1CiEhIdBoNIiMjMS+ffuabF9RUYElS5YgODgYzs7O6NWrF9atW9eigqntKBUClsaFAfh96KSWufkXgiDA10ODmFAfPDsqBG89MhBb/2skjv81FgcWjcWnT9+NxZP64Y8RXRHUpZNFNTw6tBsmhvtbHUSISJ4CtJ3w2q35Im//eApXi7gMn1oQRjZt2oT58+djyZIlSEtLQ0xMDCZNmoTs7OxG3/PYY4/h559/xtq1a3H69Gls3LgR/fr1a1Xh1DYmhgcgYXoE/LX152X4azVImB5h0fwLQRDQtXMnjOnnixdH98J7jw3G249YtpcAl+ASdTwzontgUFBnFFdU46/bTkhdDtkBq+eMDBs2DBEREUhISKi71r9/f0yZMgXx8fEN2icmJuLxxx/HhQsX4OnZsuWWnDPS/tp6/oXRJGLUW7uaHQLa//pYzvMg6oBO6orwwIr9MJpEfDwjErED/KUuidpBu5xNU1lZidTUVMTGxta7Hhsbi4MHD5p9z7Zt2zB06FC8/fbb6Nq1K/r06YNXX30VN2/eNNseqBnWKSoqqvei9lU7/+LBwV0R3cur1QHB2iEgIupY+gd44PmYngCApdtONLp3EHUMVoWR/Px8GI1G+Pn51bvu5+cHvV5v9j0XLlzA/v37kZGRga1bt+L999/HN998g7lz5zb6feLj46HVauteQUFB1pRJdqIthoCIyHG9Mi4U3T1doDOU452fTktdDkmoRatp7ty1UxTFRnfyNJlMEAQBGzZsgFarBQC89957eOSRR7By5Up06tRwouPixYuxcOHCuj8XFRUxkMjUxPAAjA/z5xJcImqgk1qJ/3koHDPWpuDfyRcxZUhXDA7qLHVZJAGreka8vb2hVCob9ILk5eU16C2pFRAQgK5du9YFEaBmjokoirh8+bLZ9zg7O8PDw6Pei+SrrYeAiMhxxIT64KEhXSGKwOItx1FlNEldEknAqjCiVqsRGRmJpKSketeTkpIwYsQIs+8ZOXIkcnNzUVLy+xbfZ86cgUKhQLdu3VpQMhEROZL/ntwfnV2ccFJXhLX7s6QuhyRg9dLehQsXYs2aNVi3bh1OnjyJBQsWIDs7G7NnzwZQM8Qyc+bMuvZPPPEEvLy88MwzzyAzMxN79+7Fa6+9hlmzZpkdoiEioo7Fy80ZS+7vDwB4/z9nkF1QJnFFZGtWh5GpU6fi/fffx/LlyzF48GDs3bsXO3bsQHBwzS56Op2u3p4jbm5uSEpKwo0bNzB06FA8+eSTiIuLwwcffNB2d0FERLL2SGQ3RPf0QnmVCUu+PQ4ZnFRCbYhn0xARkV3Iyi/FhPf3orLahPceHYSAzp048V3meDYNERHJSoi3K+aN7Y13dp7Bq98chem2X5V5Krdja9HZNERERO0h2MsVAOoFEQDQG8oxZ/0RJGZYfhAnyQfDCBER2QWjScTfd5w0+7XabLJseyaMdyYVkj2GESIisgspWYXQGRo/xVcEoDOUIyWr0HZFkU0wjBARkV3IK248iLSkHckHwwgREdkFX3dN842saEfywTBCRER2ISrEEwFaTYOTvmsJqFlVExXiacuyyAYYRoiIyC4oFQKWxoUBgNlAIgJYGhfG/UYcEMMIERHZjYnhAUiYHgF/bcOhmHv6eHOfEQfFTc+IiMiuTAwPwPgwf6RkFSKvuBz5JRX42/cnceBcAbLySxHi7Sp1idTG2DNCRER2R6kQEN3LCw8O7opnR/XEmL4+MJpEvLPztNSlUTtgGCEiIrv354n9IAjAD8d0OHb5htTlUBtjGCEiIrvXP8ADDw3uCgD4x4+neKqvg2EYISIiWVgwvg/USgUOni/AvrP5UpdDbYhhhIiIZCHI0wUzooMBAG8lnoKJZ9Q4DIYRIiKSjbljesPdWYUTuUXYfixX6nKojTCMEBGRbHi6qvHi6J4AgHd3nkFltUniiqgtMIwQEZGszBoVAh93Z2QXlmFjSrbU5VAbYBghIiJZcVGr8Mq4UADAil1nUVJRLXFF1FoMI0REJDtT7w5CiLcr8ksqsWbfBanLoVZiGCEiItlxUirwamxfAMAney8gv6RC4oqoNRhGiIhIlu6/yx8Du2lRWmnEh7vOSV0OtQLDCBERyZIgCFg0sR8AYMOvl5BdUCZxRdRSDCNERCRbI3p7IybUG1VGEe8m8RA9uWIYISIiWXv9Vu/Id+m5yLhikLgaagmGESIikrXwrlr8YVAgAODtn9g7IkcMI0REJHuvxvaFk1LA3jPXcPAcD9GTG4YRIiKSve5eLngiqjsA4B+JpyCKPERPThhGiIjIIbw8LhSuaiWOXTZgx3G91OWQFRhGiIjIIXi7OeO5mJpD9N7ZeRpVRh6iJxcMI0RE5DCev6cnvFzVyMovxabfcqQuhyzEMEJERA7DzVmFebcO0fvXz2dRVslD9OSAYYSIiBzKtKju6O7pgmvFFVi3P0vqcsgCDCNERORQ1CoF/hTbBwDw0Z4LKCytlLgiag7DCBEROZy4gYEYEOiB4opqrNzNQ/TsHcMIERE5HIVCqNsm/ovkS7h8nYfo2TOGESIickgxod4Y0csLlUYT3ks6I3U51ASGESIickiC8HvvyNa0KzilL5K4ImoMwwgRETmsQUGdMfmuAIgi8HYiD9GzVwwjRETk0F6d0BdKhYBdp/Lw64UCqcshMxhGiIjIoYV4u+Lxu4MA8BA9e6WSugAiIqL29sq4UGw5cgVp2TeQmKFHZxc18orL4euuQVSIJ5QKQeoSOzSGESIicni+Hho8OyoEH+4+h5c2psFo+r13JECrwdK4MEwMD5Cwwo6NwzRERNQh9PRxBYB6QQQA9IZyzFl/BIkZOinKIjCMEBFRB2A0ifjfn8yvpqmNJsu2ZzYIKmQbDCNEROTwUrIKoTOUN/p1EYDOUI6UrELbFUV1GEaIiMjh5RU3HkRa0o7aFsMIERE5PF93TZu2o7bFMEJERA4vKsQTAVoNGlvAK6BmVU1UiKcty6JbGEaIiMjhKRUClsaFAUCjgWRpXBj3G5EIwwgREXUIE8MDkDA9Av7a+kMxSgFY8cQQ7jMiIW56RkREHcbE8ACMD/NHSlYhrlwvw/LvM1FUXo2KKpPUpXVo7BkhIqIORakQEN3LC48MDcLse3sBAFb9cg4m7jEiGYYRIiLqsGYMD4aHRoXz10qReEIvdTkdFsMIERF1WO4aJzw9MgQA8OGuczzRVyIMI0RE1KE9M6IHXNRKZOqKsPt0ntTldEgtCiOrVq1CSEgINBoNIiMjsW/fPoved+DAAahUKgwePLgl35aIiKjNdXFVY/rwYADsHZGK1WFk06ZNmD9/PpYsWYK0tDTExMRg0qRJyM7ObvJ9BoMBM2fOxLhx41pcLBERUXt4LiYEapUCR7JvIPlCgdTldDhWh5H33nsPzz77LJ577jn0798f77//PoKCgpCQkNDk+1588UU88cQTiI6ObnGxRERE7cHXXYPH7w4CAKzcfU7iajoeq8JIZWUlUlNTERsbW+96bGwsDh482Oj7Pv30U5w/fx5Lly616PtUVFSgqKio3ouIiKg9vTi6F1QKAQfOFeBI9nWpy+lQrAoj+fn5MBqN8PPzq3fdz88Per35JVFnz57FokWLsGHDBqhUlu2xFh8fD61WW/cKCgqypkwiIiKrde3cCQ8N6QoAWLmLvSO21KIJrIJQf+9+URQbXAMAo9GIJ554AsuWLUOfPn0s/vsXL14Mg8FQ98rJyWlJmURERFaZc28vKATg51N5yMxlr7ytWBVGvL29oVQqG/SC5OXlNegtAYDi4mIcPnwYL730ElQqFVQqFZYvX46jR49CpVJh165dZr+Ps7MzPDw86r2IiIjaW08fN0weGAgAWPkLe0dsxaowolarERkZiaSkpHrXk5KSMGLEiAbtPTw8cPz4caSnp9e9Zs+ejb59+yI9PR3Dhg1rXfVERERtbO6Ymi3idxzX4VxeicTVdAxWH5S3cOFCzJgxA0OHDkV0dDQ+/vhjZGdnY/bs2QBqhliuXLmCzz//HAqFAuHh4fXe7+vrC41G0+A6ERGRPejn74H7+vvhPyevIuGX83j3sUFSl+TwrA4jU6dORUFBAZYvXw6dTofw8HDs2LEDwcE1G8bodLpm9xwhIiKyZy+N7Y3/nLyKb9OvYP59oQjydJG6JIcmiDLYaq6oqAharRYGg4HzR4iIyCZmrP0V+87mY/rw7nhzyl1SlyNLln5+82waIiIiM+aO6Q0A+Pq3y7haVC5xNY6NYYSIiMiMYSGeGBrcBZVGEz7Ze0HqchwawwgREZEZgiBg7tia3pENv2ajsLRS4oocF8MIERFRI+7t44Pwrh64WWXEpweypC7HYTGMEBERNUIQBLx0a+7IZwcvoqi8SuKKHBPDCBERURNiw/wR6uuG4vJqfJF8SepyHBLDCBERURMUCgH/dWtX1rX7s1BWWS1xRY6HYYSIiKgZcQMD0d3TBYWlldiYwsNb2xrDCBERUTNUSgXm3FvTO/Lx3vOoqDZKXJFjYRghIiKywB8jusLfQ4OrRRX4JvWy1OU4FIYRIiIiCzirlHjhnp4AgNV7zqPaaJK4IsfBMEJERGShaVHd4eWqRk7hTWw7mit1OQ6DYYSIiMhCndRKPBsTAgBYufscTCa7P2tWFhhGiIiIrDBjeDA8NCqcv1aKxBN6qctxCAwjREREVnDXOOHpET0A1PSOiCJ7R1qLYYSIiMhKz4wMgYtaiRO5Rfjl9DWpy5E9hhEiIiIrdXFVY/rwYADAil1n2TvSSgwjRERELfDcqBCoVQocyb6B5AsFUpcjawwjRERELeDrocHUoUEAauaOUMsxjBAREbXQi6N7QqUQcOBcAY5kX5e6HNliGCEiImqhbl1c8NCQrgCAlbvYO9JSDCNEREStMOfeXlAIwM+n8vBVSja+S7+C5PMFMHJDNIuppC6AiIhIznr6uGFI985IvXQDi7Ycr7seoNVgaVwYJoYHSFidPLBnhIiIqBUSM3RIvXSjwXW9oRxz1h9BYobO9kXJDMMIERFRCxlNIpZtzzT7tdpBmmXbMzlk0wyGESIiohZKySqEzlDe6NdFADpDOVKyCm1XlAwxjBAREbVQXnHjQaQl7ToqhhEiIqIW8nXXtGm7jophhIiIqIWiQjwRoNVAaOTrAmpW1USFeNqyLNlhGCEiImohpULA0rgwAGg0kCyNC4NS0dhXCWAYISIiapWJ4QFImB4Bf23DoZinRwZznxELcNMzIiKiVpoYHoDxYf5IySpEXnE5Dp4vwKbfcnDownWIoghBYM9IUxhGiIiI2oBSISC6lxcAYHQfH2xLz8VJXREOni/AyN7eEldn3zhMQ0RE1MY6u6jx2NBuAIBP9l2QuBr7xzBCRETUDmaNCoEgAL+cvoYzV4ulLseuMYwQERG1g2AvV0wI8wcArGHvSJMYRoiIiNrJ8/f0BAB8m5bLXVibwDBCRETUTiKDuyCie2dUGk34/OAlqcuxWwwjRERE7ej5mJrekfW/XkJZZbXE1dgnhhEiIqJ2FDvAH909XXCjrArfpF6Wuhy7xDBCRETUjpQKAc+OCgEArN2fBaNJlLgi+8MwQkRE1M4eHdoN2k5OuFRQhqTMq1KXY3cYRoiIiNqZi1qF6cO7A+AmaOYwjBAREdnAU9E9oFYqkHrpOo5kX5e6HLvCMEJERGQDvh4aPDg4EAA3QbsTwwgREZGNPHdrmW9ihh7ZBWUSV2M/GEaIiIhspK+/O+7p4wOTCKw7kCV1OXaDYYSIiMiGXrjVO/L14RzcKKuUuBr7wDBCRERkQyN7e6GfvzvKKo3Y8Gu21OXYBYYRIiIiGxIEAS/cOkDv3wcvorLaJHFF0mMYISIisrEHBgbCz8MZecUV2HY0V+pyJMcwQkREZGNqlQJPj6jZIn7NvgsQxY69RTzDCBERkQSeGNYdrmolTumLse9svtTlSIphhIiISALaTk547O4gANwinmGEiIhIIrNGhkAhAPvO5uOkrkjqciTDMEJERCSRIE8XTLorAACwZl/H3QStRWFk1apVCAkJgUajQWRkJPbt29do2y1btmD8+PHw8fGBh4cHoqOj8dNPP7W4YCIiIkfy/K1N0LYdvYKrReUSVyMNq8PIpk2bMH/+fCxZsgRpaWmIiYnBpEmTkJ1tfuOWvXv3Yvz48dixYwdSU1MxZswYxMXFIS0trdXFExERyd3goM6I6uGJKqOIzw5elLocSQiileuJhg0bhoiICCQkJNRd69+/P6ZMmYL4+HiL/o4BAwZg6tSp+Mtf/mJR+6KiImi1WhgMBnh4eFhTLhERkd3beUKPF75IhYdGheTF4+DqrJK6pDZh6ee3VT0jlZWVSE1NRWxsbL3rsbGxOHjwoEV/h8lkQnFxMTw9PRttU1FRgaKionovIiIiR3Vffz+EeLuiqLwaXx/Okbocm7MqjOTn58NoNMLPz6/edT8/P+j1eov+jnfffRelpaV47LHHGm0THx8PrVZb9woKCrKmTCIiIllRKAQ8O6pmE7R1B7JgNHWsTdBaNIFVEIR6fxZFscE1czZu3Ii//vWv2LRpE3x9fRttt3jxYhgMhrpXTk7HS4lERNSxPBzRDV1cnJBTeBM/nbDsF3xHYVUY8fb2hlKpbNALkpeX16C35E6bNm3Cs88+i6+//hr33Xdfk22dnZ3h4eFR70VEROTIOqmVmDE8GADw8d6OtUW8VWFErVYjMjISSUlJ9a4nJSVhxIgRjb5v48aNePrpp/Hll19i8uTJLauUiIjIwc2I7gG1SoH0nBtIvXRd6nJsxuphmoULF2LNmjVYt24dTp48iQULFiA7OxuzZ88GUDPEMnPmzLr2GzduxMyZM/Huu+9i+PDh0Ov10Ov1MBgMbXcXREREDsDH3Rl/HNIVQMfaIt7qMDJ16lS8//77WL58OQYPHoy9e/dix44dCA6u6VrS6XT19hz56KOPUF1djblz5yIgIKDu9corr7TdXRARETmI52JqJrLuzLyKrPxSiauxDav3GZEC9xkhIqKOZNZnv2HXqTzMGB6Mv00Jl7qcFmuXfUaIiIio/dX2jvxfag6ul1ZKXE37YxghIiKyM9E9vRDe1QPlVSasP3RJ6nLaHcMIERGRnREEoe4AvX8nX0J5lVHiitoXwwgREZEduv+uAARoNcgvqcC29Fypy2lXDCNERER2yEmpwKyRNXNHPt57Hsnn8/Fd+hUkny9wuO3iuZqGiIjIThWVVyHqzf+gvNpU73qAVoOlcWGYGB4gUWWW4WoaIiIimTt4Lr9BEAEAvaEcc9YfQWKGToKq2h7DCBERkR0ymkQs255p9mu1QxrLtmc6xJANwwgREZEdSskqhM5Q3ujXRQA6QzlSsgptV1Q7YRghIiKyQ3nFjQeRlrSzZwwjREREdsjXXdOm7ewZwwgREZEdigrxRIBWA6GRrwuoWVUTFeJpy7LaBcMIERGRHVIqBCyNCwOARgPJ0rgwKBWNfVU+GEaIiIjs1MTwACRMj4C/tv5QjEohYNWTEXa/z4ilVFIXQERERI2bGB6A8WH+SMkqxKXCUvz1uxMorzbBTeM4H+HsGSEiIrJzSoWA6F5eePzu7ng8qjsAYO3+LImrajsMI0RERDLyzMgeEATgl9PXcC6vROpy2gTDCBERkYwEe7nivv5+AIDPDjpG7wjDCBERkczUnua7OfUKbpRVSlxN6zGMEBERyczwnp4IC/DAzSojNqbkSF1OqzGMEBERyYwgCJg1qqZ35PPki6gyNjzZV04YRoiIiGQoblAAvN3U0BnKkZihl7qcVmEYISIikiFnlRLThwcDkP8yX4YRIiIimZo+PBhqpQLpOTdwJPu61OW0GMMIERGRTHm7OePBwYEAgHUy7h1hGCEiIpKxZ24t8/0xQ48rN25KXE3LMIwQERHJWFigB6J7esFoEvF58kWpy2kRhhEiIiKZe/bWMt+Nv2ajrLJa4mqsxzBCREQkc2P7+SLYywVF5dXYnHpZ6nKsxjBCREQkcwqFgGdG9AAAfHrgIkwmUdqCrMQwQkRE5AAeGRoEd2cVLuSXYs+Za1KXYxWGESIiIgfg5qzC41FBAIB1B+S1zJdhhIiIyEHMjO4BhQDsO5uP0/piqcuxGMMIERGRgwjydMGEAf4AgE9l1DvCMEJERORAapf5bkm7goKSComrsQzDCBERkQOJDO6Cgd20qKw24ctfs6UuxyIMI0RERA5EEATMurVF/OeHLqGy2iRxRc1jGCEiInIw998VAF93Z1wrrsAPx3OlLqdZDCNEREQORq1S4Klbm6Ct3Z8FUbTvTdAYRoiIiBzQtKjucFYpkHGlCL9dvC51OU1iGCEiInJAnq5q/DGiKwBg3X77XubLMEJEROSgnrk1kXVnph45hWUSV9M4hhEiIiIH1cfPHTGh3jCJwGcHL0pdTqMYRoiIiBzYrFuboG36LQfF5VUSV2MewwgREZEDGx3qg54+riipqMY3qZelLscshhEiIiIHplD8vgnapwcuwmiyv2W+DCNEREQO7o8RXaHt5ITswjL8fPKq1OU0wDBCRETk4FzUKkyL6g4AWGeHp/kyjBAREXUAM6ODoVQIOHShECdyDVKXUw/DCBERUQcQ2LkT7r8rAACwbv9FaYu5A8MIERFRBzFrZA8AwPajucgrLpe2mNswjBAREXUQQ7p3wZDunVFpNGHDoWypy6nDMEJERNSB1C7zXX/oEsqrjBJXU4NhhIiIqAOZFO6PQK0GBaWV2HY0V+pyADCMEBERdSgqpQIzR/QAUHOaryhKvwkawwgREVEHM+3u7ujkpMQpfTHW7svCd+lXkHy+QLLdWVsURlatWoWQkBBoNBpERkZi3759Tbbfs2cPIiMjodFo0LNnT6xevbpFxRIREVHraV2cEBXiCQB4c8dJvPJVOqZ9cgij3tqFxAydzeuxOoxs2rQJ8+fPx5IlS5CWloaYmBhMmjQJ2dnmZ+VmZWXh/vvvR0xMDNLS0vDGG29g3rx52Lx5c6uLJyIiIuslZuiw58y1Btf1hnLMWX/E5oFEEK0cLBo2bBgiIiKQkJBQd61///6YMmUK4uPjG7R//fXXsW3bNpw8ebLu2uzZs3H06FEkJydb9D2Lioqg1WphMBjg4eFhTblERER0G6NJxKi3dkFnML/PiADAX6vB/tfHQqkQWvW9LP38tqpnpLKyEqmpqYiNja13PTY2FgcPHjT7nuTk5AbtJ0yYgMOHD6OqqsrseyoqKlBUVFTvRURERK2XklXYaBABABGAzlCOlKxCm9VkVRjJz8+H0WiEn59fvet+fn7Q6/Vm36PX6822r66uRn5+vtn3xMfHQ6vV1r2CgoKsKZOIiIgaYenOq7bcobVFE1gFoX63jSiKDa41197c9VqLFy+GwWCoe+Xk5LSkTCIiIrqDr7umTdu1BZU1jb29vaFUKhv0guTl5TXo/ajl7+9vtr1KpYKXl5fZ9zg7O8PZ2dma0oiIiMgCUSGeCNBqoDeUw9yk0do5I7WrbWzBqp4RtVqNyMhIJCUl1buelJSEESNGmH1PdHR0g/Y7d+7E0KFD4eTkZGW5RERE1BpKhYClcWEAaoLH7Wr/vDQurNWTV61h9TDNwoULsWbNGqxbtw4nT57EggULkJ2djdmzZwOoGWKZOXNmXfvZs2fj0qVLWLhwIU6ePIl169Zh7dq1ePXVV9vuLoiIiMhiE8MDkDA9Av7a+kMx/loNEqZHYGJ4gE3rsWqYBgCmTp2KgoICLF++HDqdDuHh4dixYweCg4MBADqdrt6eIyEhIdixYwcWLFiAlStXIjAwEB988AEefvjhtrsLIiIissrE8ACMD/NHSlYh8orL4eteMzRjyx6RWlbvMyIF7jNCREQkP+2yzwgRERFRW2MYISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJKyejt4KdRuEltUVCRxJURERGSp2s/t5jZ7l0UYKS4uBgAEBQVJXAkRERFZq7i4GFqtttGvy+JsGpPJhNzcXLi7u0MQ2u4An6KiIgQFBSEnJ8dhz7xx9Hvk/cmfo9+jo98f4Pj3yPtrOVEUUVxcjMDAQCgUjc8MkUXPiEKhQLdu3drt7/fw8HDI/8Bu5+j3yPuTP0e/R0e/P8Dx75H31zJN9YjU4gRWIiIikhTDCBEREUmqQ4cRZ2dnLF26FM7OzlKX0m4c/R55f/Ln6Pfo6PcHOP498v7anywmsBIREZHj6tA9I0RERCQ9hhEiIiKSFMMIERERSYphhIiIiCTl8GFk1apVCAkJgUajQWRkJPbt29dk+z179iAyMhIajQY9e/bE6tWrbVRpy1lzj7/88gsEQWjwOnXqlA0rttzevXsRFxeHwMBACIKAb7/9ttn3yOkZWnt/cnt+8fHxuPvuu+Hu7g5fX19MmTIFp0+fbvZ9cnmGLbk/uT3DhIQEDBw4sG5DrOjoaPz4449Nvkcuzw+w/v7k9vzuFB8fD0EQMH/+/Cbb2foZOnQY2bRpE+bPn48lS5YgLS0NMTExmDRpErKzs822z8rKwv3334+YmBikpaXhjTfewLx587B582YbV245a++x1unTp6HT6epeoaGhNqrYOqWlpRg0aBA+/PBDi9rL7Rlae3+15PL89uzZg7lz5+LQoUNISkpCdXU1YmNjUVpa2uh75PQMW3J/teTyDLt164Z//OMfOHz4MA4fPoyxY8fiwQcfxIkTJ8y2l9PzA6y/v1pyeX63++233/Dxxx9j4MCBTbaT5BmKDiwqKkqcPXt2vWv9+vUTFy1aZLb9n//8Z7Ffv371rr344ovi8OHD263G1rL2Hnfv3i0CEK9fv26D6toWAHHr1q1NtpHjM6xlyf3J+fmJoijm5eWJAMQ9e/Y02kbOz9CS+5P7MxRFUezSpYu4Zs0as1+T8/Or1dT9yfX5FRcXi6GhoWJSUpI4evRo8ZVXXmm0rRTP0GF7RiorK5GamorY2Nh612NjY3Hw4EGz70lOTm7QfsKECTh8+DCqqqrardaWask91hoyZAgCAgIwbtw47N69uz3LtCm5PcOWkuvzMxgMAABPT89G28j5GVpyf7Xk+AyNRiO++uorlJaWIjo62mwbOT8/S+6vltye39y5czF58mTcd999zbaV4hk6bBjJz8+H0WiEn59fvet+fn7Q6/Vm36PX6822r66uRn5+frvV2lItuceAgAB8/PHH2Lx5M7Zs2YK+ffti3Lhx2Lt3ry1Kbndye4bWkvPzE0URCxcuxKhRoxAeHt5oO7k+Q0vvT47P8Pjx43Bzc4OzszNmz56NrVu3IiwszGxbOT4/a+5Pjs/vq6++wpEjRxAfH29ReymeoSxO7W0NQRDq/VkUxQbXmmtv7ro9seYe+/bti759+9b9OTo6Gjk5OXjnnXdwzz33tGudtiLHZ2gpOT+/l156CceOHcP+/fubbSvHZ2jp/cnxGfbt2xfp6em4ceMGNm/ejKeeegp79uxp9ANbbs/PmvuT2/PLycnBK6+8gp07d0Kj0Vj8Pls/Q4ftGfH29oZSqWzQQ5CXl9cg8dXy9/c3216lUsHLy6vdam2pltyjOcOHD8fZs2fbujxJyO0ZtgU5PL+XX34Z27Ztw+7du9GtW7cm28rxGVpzf+bY+zNUq9Xo3bs3hg4divj4eAwaNAj/+te/zLaV4/Oz5v7Msefnl5qairy8PERGRkKlUkGlUmHPnj344IMPoFKpYDQaG7xHimfosGFErVYjMjISSUlJ9a4nJSVhxIgRZt8THR3doP3OnTsxdOhQODk5tVutLdWSezQnLS0NAQEBbV2eJOT2DNuCPT8/URTx0ksvYcuWLdi1axdCQkKafY+cnmFL7s8ce36G5oiiiIqKCrNfk9Pza0xT92eOPT+/cePG4fjx40hPT697DR06FE8++STS09OhVCobvEeSZ9huU2PtwFdffSU6OTmJa9euFTMzM8X58+eLrq6u4sWLF0VRFMVFixaJM2bMqGt/4cIF0cXFRVywYIGYmZkprl27VnRychK/+eYbqW6hWdbe4z//+U9x69at4pkzZ8SMjAxx0aJFIgBx8+bNUt1Ck4qLi8W0tDQxLS1NBCC+9957Ylpamnjp0iVRFOX/DK29P7k9vzlz5oharVb85ZdfRJ1OV/cqKyurayPnZ9iS+5PbM1y8eLG4d+9eMSsrSzx27Jj4xhtviAqFQty5c6coivJ+fqJo/f3J7fmZc+dqGnt4hg4dRkRRFFeuXCkGBweLarVajIiIqLfk7qmnnhJHjx5dr/0vv/wiDhkyRFSr1WKPHj3EhIQEG1dsPWvu8a233hJ79eolajQasUuXLuKoUaPEH374QYKqLVO7jO7O11NPPSWKovyfobX3J7fnZ+7eAIiffvppXRs5P8OW3J/cnuGsWbPqfr74+PiI48aNq/ugFkV5Pz9RtP7+5Pb8zLkzjNjDMxRE8dasFCIiIiIJOOycESIiIpIHhhEiIiKSFMMIERERSYphhIiIiCTFMEJERESSYhghIiIiSTGMEBERkaQYRoiIiEhSDCNEREQkKYYRIiIikhTDCBEREUmKYYSIiIgk9f8BHL43CmNPpAQAAAAASUVORK5CYII=", "text/plain": [ "Figure(PyObject