RK_ODUSystem_Solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Численное интегрирование систем ОДУ первого порядка методом Рунге-Кутты"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Необходимо найти решение\n",
    "$$ \\frac{d\\vec{x}}{dt} = \\vec{f}(t, \\vec{x}), \\hspace{1cm} \\vec{x} \\in \\mathbb{R}^n $$\n",
    "\n",
    "удовлетворяющее начальному условию\n",
    "$$ \\vec{x}(t_0) = \\vec{x}_0 $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На примере одношагового *метода Рунге-Кутты* порядка $s$ с параметрами: $a_{ij}, b_i, c_i$, рассмотрим его распространение на системы при численном решении полученной выше задачи Коши.\n",
    "$$ x_{m+1}^i = x^i(t_m + h) \\approx x^i(t_m) + h \\sum_{i=1}^s b_i k_{ij}(m) \\hspace{1cm} i = \\overline{1, n} $$\n",
    "\n",
    "Здесь функции $k_{ij}(m)$ вычисляются по правилу:\n",
    "$$ k_{i1}(m) = f_i(t_m, x^1, \\ldots, x^n) $$\n",
    "\n",
    "для $j = 2, \\ldots, s$:\n",
    "$$\n",
    "k_{ij}(m) = f_i \\left(t_m + c_j h, x^1 + \\sum_{l=1}^{j-1}a_{jl}k_{1l}(m),\\ldots, x^n + \\sum_{l=1}^{j-1}a_{jl}k_{nl}(m)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Чтобы записать это решение в векторном виде, предварительно представим параметры метода в векторно-матричной форме:\n",
    "\n",
    "$A$ -- матрица ($s \\times s$) с элементами $a_{ij}$, $\\vec{b}$ и $\\vec{c}$ -- вектора из $\\mathbb{R}^s$ с элементами $b_i$ и $c_i$ соответственно.\n",
    "\n",
    "Тогда:\n",
    "$$ \\vec{x}_{m+1} = \\vec{x}(t_m + h) \\approx \\vec{x}(t_m) + h \\left(K(m) \\cdot \\vec{b} \\right)$$\n",
    "\n",
    "$K(m)$ -- матрица ($n \\times s$) с элементами $k_{ij}(m)$.\n",
    "Столбцы этой матрицы вектора $\\left\\{ \\vec{k}_j \\right\\}_{j=1}^s$:\n",
    "$$ \\vec{k}_j = \\vec{f}(t_m + c_j h, \\vec{x}_n + \\vec{p}(j, m))$$\n",
    "\n",
    "Здесь вектор $\\vec{p}$ отвечает за смещение аргумента функции $\\vec{f}$:\n",
    "$$ \\vec{p}(j, m) = \\left(\\sum_{l=1}^{j-1}a_{jl}k_{1l}(m),\\ldots, \\sum_{l=1}^{j-1}a_{jl}k_{nl}(m) \\right)^T $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Методическая погрешность $\\vec\\sigma (h)$ -- вектор функция. Каждая ее компонента определяется равенством:\n",
    "$$ \\sigma^i (h) = \\frac{\\max\\limits_{i \\in [1, n]}|x_i^{(h)} - x_i^{(h/2)}|}{2^s - 1} $$\n",
    "\n",
    "Чтобы найти решение системы с заданной точностью $\\varepsilon$ необходимо добиться выполнения неравенства:\n",
    "$$ \\sigma^i < \\varepsilon \\hspace{1cm} \\forall i \\in [1, n] $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Реализация представлена ниже"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_odesystem(ode_system: 'cuple of functions-list and 1Darray initial conditions',\n",
    "               RK_method: 'cuple of params-list like [A, b, c] and accuracy order',\n",
    "               integration_segment: '[t_0, t_N]', N0, epsilon: 'accuracy' = 0):\n",
    "    \"\"\"Решает систему ОДУ с заданной точностью\"\"\"\n",
    "\n",
    "    def runge_kutta_odesystem(ode_sys, rk_method, integration_seg, n):\n",
    "        \"\"\"Решает систему ОДУ методом Рунге-Кутты\"\"\"\n",
    "        f, r0 = ode_sys\n",
    "        (A, b, c), accuracy_order = rk_method\n",
    "        sys_order = r0.size\n",
    "        r = [r0]\n",
    "        t_grid, h = np.linspace(*integration_seg, n, retstep=True)\n",
    "        for n, tn in enumerate(t_grid[:-1]):\n",
    "            K = np.array([[f[i](t=tn, r=r[n])] for i in range(sys_order)])\n",
    "            for j in range(1, accuracy_order):\n",
    "                p = np.array([A[j][:j] @ K[i] for i in range(sys_order)])\n",
    "                Kj = np.array([f[i](t=(tn + c[j]*h), r=(r[n] + h*p)) for i in range(sys_order)])\n",
    "                K = np.hstack([K, np.transpose([Kj])])\n",
    "            r.append(r[n] + h*(K @ b))\n",
    "        return t_grid, np.array(r)\n",
    "    \n",
    "    t_grid, r = runge_kutta_odesystem(ode_system, RK_method, integration_segment, N0)\n",
    "    if epsilon:\n",
    "        sigma = 2 * epsilon\n",
    "        k = 1\n",
    "        R = [r]\n",
    "        while np.all(sigma > epsilon):\n",
    "            k *= 2\n",
    "            R.append(runge_kutta_odesystem(ode_system, RK_method, integration_segment, N*k - 1)[1][::k])\n",
    "            diff = np.absolute(R[-2] - R[-1])\n",
    "            sigma = diff.max(axis=0)/(2**rk_method[1] - 1)\n",
    "        return t_grid, np.array(r), sigma\n",
    "    return t_grid, np.array(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Классический 4-стадийный метод РК 4-го порядка"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.array([0, 0.5, 0.5, 1])\n",
    "b = np.array([1/6, 2/6, 2/6, 1/6])\n",
    "A = np.array([[0  , 0  , 0, 0],\n",
    "              [0.5, 0  , 0, 0],\n",
    "              [0  , 0.5, 0, 0],\n",
    "              [0  , 0  , 1, 0]])\n",
    "accuracy_order = 4\n",
    "rk_method = ([A, b, c], accuracy_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Пример:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{cases}\n",
    "    \\dot{x} = -y  \\hspace{1cm} &x_0 = 1\\\\\n",
    "    \\dot{y} = x  &y_0 = 0\\\\\n",
    "\\end{cases}\n",
    "\n",
    "Решением этой системы является:\n",
    "\\begin{cases}\n",
    "    x = \\cos t \\\\\n",
    "    y = \\sin t \\\\\n",
    "\\end{cases}\n",
    "\n",
    "То есть, если построить график y(x), мы должны увидеть окружность"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# r = [x, y]\n",
    "# dx/dt = -y\n",
    "def fx(*, t=0, r=[0, 0, 0]):\n",
    "    r = np.array(r)\n",
    "    return -r[1]\n",
    "# dy/dt = x\n",
    "def fy(*, t=0, r=[0, 0, 0]):\n",
    "    r = np.array(r)\n",
    "    return r[0]\n",
    "\n",
    "# x(0) = 1, y(0) = 0\n",
    "r0 = np.array([1, 0])\n",
    "ode_system = ([fx, fy], r0)\n",
    "\n",
    "# Интегрируем по всей окружности [0, 2pi]\n",
    "integration_segment = [0, np.pi*2]\n",
    "N = 100     # кол-во точек в сетке\n",
    "epsilon = 1e-3 # точность\n",
    "t, r, sigma = solve_odesystem(ode_system, rk_method, integration_segment, N, epsilon)\n",
    "x, y = r.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy x(t): 4.113091413386153e-08\n",
      "accuracy y(t): 5.301527154777386e-08\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAJRCAYAAADMCCpKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABkbUlEQVR4nO3deXxU1fnH8c9Jwg4qCiIoChEIBhcU6oYbiopLcF+wWrRaq9Vat9atrZa6/upW11pX3Iq4VeKGqGgrKgqKIksQoggCigIqsiY5vz+eSTNgAglk5tw79/t+veY1987cyTyHSciTszzHee8RERERkfjKCx2AiIiIiGwYJXQiIiIiMaeETkRERCTmlNCJiIiIxJwSOhEREZGYU0InIiIiEnMFoQMIqV27dr5Lly4ZfY8ff/yRVq1aZfQ9oirJbYdktz/JbYdktz/JbYdkt19tz3zbJ0yY8I33vn1tzyU6oevSpQvjx4/P6Hu88cYb7Lfffhl9j6hKctsh2e1Pctsh2e1Pctsh2e1X2/fL+Ps452bV9ZyGXEVERERiTgmdiIiISMwpoRMRERGJOSV0IiIiIjGnhE5EREQk5pTQiYiIiMScEjoRERGRmFNCJyIiIhJzSuhEREREYk4JnYiIiEjMKaETERERiTkldCIiIiIxp4ROREREJOaU0ImIiIjEnBI6ERERkZhTQiciIiISc0roRERERGJOCZ2IiIhIzCmhExEREYm5SCV0zrkHnHNfO+c+qeN555y7zTk3wzn3sXNul7TnhjjnPk3dhmQvahEREZGwIpXQAQ8BA9fy/CFA99TtTOBuAOfcpsCVwG7ArsCVzrm2GY1UREREJCIildB57/8DLFzLJUcAD3vzLrCJc64jcDAw2nu/0Hu/CBjN2hNDEZENV14OvXpBQYHdl5dn9/UiIikFoQNooC2B2Wnnc1KP1fW4iEjdZsyAww6DmTNhm23guutg441hyRL48Ue7r+1W/dzbb8Py5fa1pkyxpGzHHdn5xx9hs80sUVvbbeRIWLzYXj91Kuy1F9x1F7RrV3Nr2xby82uPv7wcSkqgrAyKiqC0FAoLs/JPJyLR4rz3oWNYjXOuC/C89377Wp57Hrjee/9W6vw14BJgP6C59/7q1ON/ApZ572+s5WuciQ3X0qFDhz7Dhw/PUEvMkiVLaN26dUbfI6qS3HZIdvuz1fbmc+eyw+WX03L2bJZ27syka69leadOALjKSpouXEizBQtoumABzRYsoNk339h96rj5vHm4eryPd47K5s2pbNFitdvGkyat9noPLNx1V6pWrqQgFcPabs3nz1/n+3vnqGjThlUbb8yqjTay+9Rti5dfpsl33+G8p8o5lm29Ne8/9NB6/Vs2liR/30Oy26+2Z77t/fv3n+C971vbc3FL6O4B3vDe/yt1XoYlc/sB+3nvf13bdXXp27evHz9+fKPGv6Y33niD/fbbL6PvEVVJbjsku/1ZafvKlVBcbL1U3oNz0Lq1PTZnDsybB1VVq7+mWTPYaqua2+OP22ur5eXBf/9rXyf91qKFff019eoF06bZ++TlQc+eMHly/du/5uu33RaeeAK++eant2+/Xf18wQL7N1hT9+7QrdtPb126QNOmDfkXXi9J/r6HZLdfbd8v4+/jnKszoYvbkOtI4Fzn3HBsAcR33vt5zrlRwLVpCyEOAi4LFaSINKLKSvjsM/jkk9VvZWVQUVFznffwww+w0UZw4IGrJ27Vt003XT0x+/DDnyZke+5Z/9hKS3865NkQtb2+vkOm3lvyOn26xe+cDfP27m1DyW+9Zf8e1fLyYOutV0/yWreGv/0NPv9cQ7YiMRephM459y+st62dc24OtnK1CYD3/h/Ai8ChwAxgKXBa6rmFzrm/Au+nvtRQ7/3aFleISBSsOQfsn/+0JCQ9cZsyBZYtq3lN166w/fYwaBA8/HBNT1x1QvbKK/V//w1NyAoLYfLkhr2msV7vHLzwQt0JoffWizdzpiV46bcRI2DhGv9FTpkCffvCjTfafXGxzfMTkViI1E+r937wOp73wDl1PPcA8EAm4hKRDJg9G/bZB+bOteRjyhRbFFCtY0dL3M46y+63396SjPR5KmecETYhC21t8TsHm29utz32+OnzixbZoov0YelFi+D00+24RQvr7fvZzyzB69sXevSoe4GGiAQVqYRORHJUVZWt4nzrLZuj9tZbMGvWT6/Ly4MxY2xu2Wabrfvrxj0hC6ltW+vRXHPI+dlnYfz4mtv998Ntt9lrWreGXXapSfA6dIDf/nb9E2oRaTRK6ESk0blVq+Cdd2oSuLFja4b4OnSAvfeGCy+Ev//d5m+lJxT77BM09kSpaw5fjx5w0kl2TWWlJX3pSd5dd9WUa6k2dSocfLANm4tI1imhE5H1lz4Hbsst7fiTT9jrnXdqVmD26AFHHmlJ3F572UrO6oUJhx++YUOmsmHq08OZn289pr16wZDUroqrVtkQ+S671AzZeg8zZrD74MH2uR54IBxwgA3rikjGKaETkYarnvO2//7w9df22BdfwJ13Qt++zB00iM6DB0O/ftYjVxcNmcZTkyaw004/HbLt0IEfunWj+VNP2VCtc7DzzpbcHXigfT80bx46epGcFKmtv0Qkwlatgtdeg/PPt1627bevSeaq5efD++8z85xz4Oij157MSfyVllpSl59v92+9xeShQ61O3jvvwF/+Aq1awU03wYABNm/v4INtJe3zz2vbM5FGpB46EanbokXw0ku2RdXLL8N331lx3gED4JJL7Bf1zJk1PTRFRaEjlmyqrYf1iy8sSdt9d7v96U9WiubNN2H0aLv9/verv2bqVBumnTIle7GL5BgldCJJt2YtuNtvh4kTrfflv/+1SfGbbw7HHGO13wYMsF4XsGE0zYGTdWnTxhK2ww+38zlzbO/c9Pl3U6fCaafB8cfb91iTJuHiFYkhJXQiSVdSYr9Mq+fFHXCAPb799tYLV1ICu+5qPXBr0hw4WR9bbbX6/DvnbIePZ56Bhx6yodmjj7bkbv/9VeBYpB70UyKSVLNmwSOP/HSYKy/PdhPo2jVMXJIMtZVM2XJL2+ljxAi73X+/1SM85hhL7vbdV8mdSB30kyGSJEuWwNNPw7BhVsAXoGVL21rL+5pacErmJNPq6t0tKbHb8uU2b3PECHjsMatvVz30v9decPXVto+t9qAVAbTKVST3VVVZ8nbqqbDFFnY/axYMHWqb3k+aBNttV7NSUfPgJAqaN7f6hY8/bqupn3oK9tvP/hj5+c9tmkB10eOSktDRigSnHjqRXDVjhm1e//DDlsC1aQODB1tx2H79aor7gubBSbS1bGk9c8ccAz/+aPPtqhdUVFXZtIHHHrPnVedOEko9dCK5oLy8pqZXp062oXr37jYsVVRkv+zmz4d777XhqvRkTiROWrWynuTqRTrO2YrYk0+2xRZ/+IP9MSOSMEroRHLBgAHWS1FZCfPmwccfw3XXWU2wUaNsX86WLUNHKdI40gsab7edfe+PHm1DsjffbH/MHHwwPPssVFSEjlYkKzTkKhJXlZVW8PeWW2wu3JrPXXppmLhEMq22BRXdutkfNnPn2urYf/7TSp906gS/+hWccYb14InkKPXQicTNkiVW/LdHD/uF9cUXtsVW9RCUdmyQJOvUyXan+OwzeO4523N26FDo0gWOOsrq3GnLMclBSuhE4mLOHCv027kznHeeJXFPPmnzhd5+e/U9NbVSVZKuoMB2NnnxRfsZ+f3vYexY242ienqCVshKDtGQq0jUjR9vw6ojRtiKvmOOgQsugD32qLlGOzaI1K2w0OaUXnWVzSVNXyE7dartWdy2bdAQRTaUeuhEoqiy0oaL9t3XVqyWlsJvfwszZ1pil57MiUj9NGu2+gpZsILa22wDl11m9e5EYkoJnUhUlJfXFPht0cKKqs6aZav25syx+y5dQkcpEm/pK2SLi+GFF+DQQ+GGG+zn6/zz4csvQ0cp0mBK6ESiYPly2HPPms3KV62yFXkzZtjw6kYbhY5QJDdUT0+oqLD7Qw+F4cNt6PX44+GOO+yas8766epxkQhTQicSUkUF3HefrVj96qvVn5s3TxuRi2RLUZGtgP30U/jlL+HBB62e3ZAh9oeWSMQpoRMJoarKegWKi61GVqdONo9HpUdEwuraFe6+26ZAnHeerSQvLrbeuxdeUMkTiSwldCLZ5L3N4dl5Z9tXtXlzW/zwzjvw+usqPSISFVtuafNWZ82yBRMvvwyHH66SJxJZSuhEsuX112116qBBsHQpPP44TJxo5879dG5PYWHoiEWkfXu45hpL7NL3QK6qgrKycHGJrEEJnUimjRtnWxIdcICtnrv3Xvsrf/Dg1csniEh0tW1rq9DTf2arqizZW7YsXFwiKfptItLYysuhVy/2PeAAaNMGdt8dPv4Ybr3VJlyfcQY0aRI6ShFpqPSSJ926wUEHwR//aPNdH3uspmCxSABK6EQa2yGHwJQpuKoq23d1880tyfvd72zOnIjEU/q0iE8/tXl1b75pP+Mnnwy7787GkyaFjlISSgmdSGOpqLBeuOnTV3/822+hdesgIYlIhu2zD7z3Hjz8MMydy87nnQfHHmu7uohkkRI6kcYwdiz06WNFgFu1UvkRkSTJy4NTToGyMj479VR46SUrdfL738PixaGjk4RQQieyIb7+Gk49Ffbayzb4fuYZ+Ogj6NkTn5en8iMiSdKqFbOGDLHh2JNPhptusrl2V11lCZ7q10kGKaETWR+VlXDXXdb79vjjcOmltnXQUUfBttvC5Mm8+dprKj8ikkSdOsH998MHH8COO8Jf/mL/P6h+nWSQ9hUSaahx4+A3v7H/rPffH+6803riRETS9e4Nr71mPXPVK2BVv04yRD10IvX17bdw5plWHHj+fNu669VXlcyJSN2cs/8j0uvXOaepGNLolNCJ1CVVT46CAhtC6dYNHngALrzQhk1OOGH1yvEiIrVJr1/XtatNyxg0yIqLf/116OgkRyihE6lLSYklbpWVMG8erFxpW3XdeKMVDBYRqY/0+nXl5VZofOhQePppWyzx2GO2z7PIBlBCJ1KbykqbxJxe+X3FCth++3AxiUhuaNoU/vQn+PBD6N7dVsQefjjMnh06MokxJXQia/rsM+jff/W/mFVPTkQaW69e8NZbVpD8jTest+6uu7SFmKwXJXQi1byHBx+0MgMffWRDq8XFNu9F9eREJBPy821bwE8+sX2fzzkH9ttPK2GlwZTQiYBNTD7qKPjlL6FvX5vjctFFNfNeVE9ORDKpa1d45RVbeDVpkv1h2aGDihFLvSmhExk5EnbYwbbruekmqxu1zTahoxKRpHEOTjsNpkyBZs3sD00VI5Z6UkInyfXDD/CrX8ERR0DHjjBhgpUkydOPhYgE1LEjLF1ac15VZUmdyFroN5ck09ixsNNOtj3PpZfa7g9awSoiUVFUtPofl1VVNr9u+fJwMUmkKaGTZEgvEtyuHey9tz3+n//AddfZ8IaISFSkFyPebjs4/XRbAbv77jB9eujoJIK0l6skQ3WR4Koq28Jrk01sJasKBItIFFUXI0531FEwZAjssgv84x9Wv04kRT10kgzVyVy1H35QMici8XLYYbZbzS67wCmn2Kr8H38MHZVEhBI6yW0VFfCHP6yezKlIsIjE1VZbweuvwx//CA89BD/7mdWwk8RTQie565tvYOBA+Nvf4Oc/t3koKhIsInFXUAB//avVrVu40JK6e+/VfrAJp4ROctOECdCnj22r88AD8OijVttJRYJFJFcMGGBzgffaC848EwYNsj9cVYw4kZTQSe4ZNgz69bO/Vt96ywp1iojkog4dYNQouOYaeP55my+sYsSJpIROcsfKlVan6dRTLaGbMMG28RIRyWV5eXD55T+tW6f9YBNFCZ3khrlzoX9/q9N08cX2F2v79qGjEhHJnp49V0/q2rXTvLoEUUIn8Td2rM2X++gjeOIJWwRRoBKLIpIw6cWI27SBr76yuXUrV4aOTLIgUgmdc26gc67MOTfDOXdpLc/f4pybmLpNd84tTnuuMu25kVkNXLKveueH/HybENysGbz7Lhx/fOjIRETCqC5GXFEBixfDFVfAfffZav+FC0NHJxkWmW4M51w+cCdwIDAHeN85N9J7P6X6Gu/9BWnX/xbYOe1LLPPe985SuBLa4YfD1Kk1582bay9WEZFqeXlw9dVWc/OMM2CPPWzRRPfuoSOTDIlSD92uwAzvfbn3fiUwHDhiLdcPBv6VlcgkWr7/fvVkDmDGjDCxiIhE2SmnwGuvWQ/dbrvBG2+EjkgyJEoJ3ZbA7LTzOanHfsI5tw3QFXg97eHmzrnxzrl3nXNHZixKCWvePNhnHzt2zu6184OISN322gvGjYOOHeHAA+H++0NHJBngfERWwDjnjgUGeu/PSJ2fAuzmvT+3lmsvAbby3v827bEtvfdfOucKsUTvAO/9zFpeeyZwJkCHDh36DB8+PDMNSlmyZAmtW7fO6HtEVWO3vcXs2ez4hz/QdPFipp93Hls/8QQtZ89maefOTLr2WpZ36tRo79UY9Nkns+2Q7PYnue0Q7fbnL1lCr7/8hU3Hj+eLE06g/Fe/snnIjSTKbc+0bLW9f//+E7z3tdfj8t5H4gbsAYxKO78MuKyOaz8E9lzL13oIOHZd79mnTx+faWPGjMn4e0RVo7Z93Djv27Xzvn17799/v/G+bgbps0+uJLc/yW33PgbtX7XK+3PO8R68HzDA+549vc/P97642PuZMzfoS0e+7RmUrbYD430dOU2UhlzfB7o757o655oCJwI/Wa3qnOsJtAXeSXusrXOuWeq4HdAPmLLmayWmXnrJasxttBG8/baKBYuIrK+CArjjDrj9dnj1Ve0skUMik9B57yuAc4FRwFRghPd+snNuqHNuUNqlJwLDU5lqte2A8c65j4AxwPU+bXWsxNiwYfafTFGR1Zvr1i10RCIi8XfuudpZIsdEpmwJgPf+ReDFNR778xrnV9XyureBHTIanGSX93DDDXDZZbYB9TPPWKFMERFpHD17Ws9cVZWdd+kSNBzZMJHpoRP5n6oqOP98S+YGD4YXXlAyJyLS2NJ3lsjPt9ImH34YOipZT0roJFpWrLAk7rbb4IIL4NFHoWnT0FGJiOSe9J0lpk6F1q1h//2txInEjhI6iYbycthuO9vxYcQI6527+ebV53iIiEhmdO8O//0vbLaZTXN5883QEUkD6belRMOhh9pcDrCCwc89FzYeEZGk2WYb+M9/oHNnOOQQGDUqdETSAEroJLxFi1ZfXeW9VluJiITQqZP1zvXoAYMG6Y/rGFFCJ2F99x0cfLAdaysvEZHw2reHMWOgd2845hh44onQEUk9KKGTcH74wbr1J06Ef/7T5tDl59uqq9LS0NGJiCRX27YwejTsuSecdBI89FDoiGQdlNBJGD/+aPPm3nvP/vr71a9qVltNnmyrr0REJJyNNoKXX4YDDoDTToOOHW2niV69bCGbRIoSOsm+pUvh8MNtG6/HH4ejjgodkYiI1KZlSxg50kqazJ+vbcIiLFI7RUgCLF8ORxxhk24ffRSOPz50RCIisjbNm8OyZTXn2iYsktRDJ9mzYoX1xr32Gjz4oM3LEBGR6CsqWr0u6BZbhItFaqWETrJj5Uo47jibj3HvvTBkSOiIRESkvtK3CWvZ0oZfX345dFSSRgmdZN6qVXDiifYfwt13w+mnh45IREQaIn2bsLlzYYcd4Oij4Z13QkcmKUroJGOaz50LxcXQrBk8+yz8+c9w1lmhwxIRkQ2x8cbWO7fllnDYYZboSXBK6CRjdrj8ctvw2XsrGvzUU6FDEhGRxtChA7zyii2YOOgg+Pzz0BElnhI6yZiWX3xRc6LtvEREckvXrrbf69KlcNBBNFm0KHREiaaETjLjnntw1T1zoO28RERy0Q47wPPPw5w57HjJJfD996EjSiwldNL4Xn4ZzjmHRb1716yK0nZeIiK5qV8/eOopWpWXw5FHWr1RyToldNK4Pv7YigXvsAOfXHstTJmi7bxERHLdoYcy7ZJLYMwY+PnPbUcJySoldNJ45s61FU8bbQTPP09lixahIxIRkSz5+sAD4dZb4ZlnoF077fuaZUropHEsWWL7sy5eDC+8YMvZRUQkWX73O0vmFi/Wvq9Zpr1cZcNVVsLgwfDRRzY5dqedQkckIiKhpK921b6vWaMeOtlwF1xgidwdd8Ahh4SORkREQlpz39dttgkXS4IooZMN8/e/w+23w0UXwdlnh45GRERCS9/3taAAli2Dr78OHVXOU0In6++556x37qij4P/+L3Q0IiISBen7vr77rg3BHnec7estGaOEThqmvNxWLeXnWyK3/fbw6KOrd6+LiIgA9OkD998P//mPdQBIxmhRhDRMSYmtWqqqsvMVK6Bly7AxiYhIdJ10Enz4Idx4I+y8M5x+euiIcpK6VaRhyspqkjmAmTPDxSIiIvFw3XVw4IE21/qdd0JHk5OU0EnDbLxxzbH2ZxURkfooKIDhw6FzZzj6aCtEL41KCZ3U37PPwsKFsOmm2p9VREQaZtNNbTHdDz9YUqc9XxuVEjqpn/JyOO00+NnP7C8r7c8qIiINtf328PDDMG4c/OY34H3oiHKGEjpZt+XLbcm5czBiBDRrFjoiERGJq6OPhj/9CR58EDp10p6vjUQJnazbRRfBBx/AsGHQpUvoaEREJO6uugpat4b587XnayNRQidr98QTcNddltQNGhQ6GhERyQV5ebaDRDXt+brBlNBJ3aZPhzPOgD32sCXnIiIijWXNPV9VNWGDKKGT2i1bZvPmmjWzXromTUJHJCIiuaR6z9fqpG7w4LDxxJwSOqndeefBxx/DI49Y3SAREZHGVL3n66pVVnT42ms17LoBlNDJTz3yCNx3H1x2GRxySOhoREQkl+XlwUMP2TaSP/85rFwZOqJYUkInNcrLoVs3+MUv7Afr1FNDRyQiIknQqRPcey9MmGArYKXBlNBJjcMPr9mbdflyOOqosPGIiEhyHHUUnH46XH89/Oc/oaOJHSV0UmPatJpjLSEXEZFsu/VW2HZbOOUUWLw4dDSxooROzPjxq2/BkpenJeQiIpJdrVvDo4/Cl1/COeeEjiZWlNAJrFhh+7Ruvrklcfn5tpS8tDR0ZCIikjS77QZXXgmPP243qRcldGJLxT/5BO6/34ZdKypsKXlhYejIREQkiS67DPbcE84+G2bNCh1NLCihS7qJEy2hO/lkWxQhIiISWkGBDb1WVkJxsZ336mXVGKRWSuiSbNUqG2rdbDP4+99DRyMiIlKja1fYeGNYutQSu2nToKQkdFSRVRA6AAnohhush+6ZZ2DTTUNHIyIisrqvvqo5VvWFtVIPXVJ98gkMHQonnKB6cyIiEk1FRTV7vVafS62U0CVRRQX88pewySZw++2hoxEREaldaalVXahO6k4/PWw8EaaELoluvhnefx/uuAPatw8djYiISO0KC63qwqpVsPvuNlVo0aLQUUWSErqkmTYN/vxnG2Y97rjQ0YiIiKxbXh7cfTd88w1cfnnoaCJJCV1SlJfb0u/ttrO/dC65BJwLHZWIiEj99O4N550H99wD48aFjiZyIpXQOecGOufKnHMznHOX1vL8qc65Bc65ianbGWnPDXHOfZq6Dclu5DFQUlKzV6v3NodOREQkToYOhY4d4ayzbD64/E9kEjrnXD5wJ3AIUAwMds4V13LpE9773qnbfanXbgpcCewG7Apc6Zxrm6XQ46GsrGavVu+19FtEROKnTRurmzpxItx5Z+hoIiUyCR2WiM3w3pd771cCw4Ej6vnag4HR3vuF3vtFwGhgYIbijKc2bWqO8/K09FtEROLpmGNg4ED405/gyy9DRxMZUUrotgRmp53PST22pmOccx87555yznVu4GuT6aOP4LvvrHhwfr4tAS8tDR2ViIhIwzlnVRpWrYILLggdTWQ4Xz0MF5hz7lhgoPf+jNT5KcBu3vtz067ZDFjivV/hnPs1cIL3fn/n3MVAc+/91anr/gQs897fWMv7nAmcCdChQ4c+w4cPz2i7lixZQuvWrTP6HmvlPTtddBGtZ85k3KOPUpHeU5dhwdseWJLbn+S2Q7Lbn+S2Q7Lbn+22b/PII3R94AE+vuEGFu66a9betzbZanv//v0neO/71vqk9z4SN2APYFTa+WXAZWu5Ph/4LnU8GLgn7bl7gMHres8+ffr4TBszZkzG32Ot/v1v78H7O+7I+lsHb3tgSW5/ktvufbLbn+S2e5/s9me97cuXe19U5H1hofdLl2b3vdeQrbYD430dOU2UhlzfB7o757o655oCJwIj0y9wznVMOx0ETE0djwIOcs61TS2GOCj1WLKtWAEXXWTlSn7969DRiIiINJ5mzeCPf7SyXK1aQa9edpxQBaEDqOa9r3DOnYslYvnAA977yc65oVhGOhI4zzk3CKgAFgKnpl670Dn3VywpBBjqvV+Y9UZEze23w8yZMGoUFETmoxYREWkc111n995baa6SEttZIoEi9Vvee/8i8OIaj/057fgybCi2ttc+ADyQ0QDj5Ouv4a9/hcMOg4MOCh2NiIhI40svwVVVleiSXFEacpXG9Oc/w9KlcNNNoSMRERHJjKIiK8VVbZttwsUSmBK6XPTxx3DvvXDuuao3JyIiuau01Epx5edbYrf11qEjCkYJXa7xHs4/HzbZxHrpREREclVhoc2Zq6iAa6+FN96A//43dFRBKKHLNSNHwpgxtt9dW+1+JiIiCfHb39o+r5ddVrPVZYIoocsV5eVWnuTII6FpUzjwwNARiYiIZE/LljYyNXYsvPBC6GiyTgldrigpsSXbYNuhHHVU2HhERESy7fTTYdtt4fLLbdVrgiihyxVlZTVdzN4neum2iIgkVJMmVrJr0iT4179CR5NVSuhyRfv2Ncd5eVrdKiIiyXTCCbDTTjb8unJl6GiyRgldLlixwnrlWra0pds9e9pSbhERkaTJy7MdJMrL4b77QkeTNUrocsFDD8FXX8G//21LtydPtqXcIiIiSTRwIOy9tw2//vhj6GiyQgld3K1aBddfD7vtBgMGhI5GREQkPOesl27+fLjtttDRZIUSurh79FH4/HP405/sG1hERESgXz/o3x/++EcoKIBevWwYNkcpoYuz6srYu+wChx4aOhoREZFomTXLypdUVlppr5KS0BFlTEHoAGQDPPEEzJgBzzyj3jkREZE1zZpVc1xVldMlvdRDF1dVVXDNNbD99nDEEaGjERERiZ6iopoOD+dyuqSXErq4evppmDrV5gbk6WMUERH5idJS2G47Oy4osGoQOUqZQBxVVcHVV9tfGsceGzoaERGRaCostFJeTz1lVSE+/jh0RBmjhC6ORo60b8orrrBCwiIiIlK3I4+05O6mm0JHkjFK6OKkvNyWXR91lO1Xt9tuoSMSERGJvvx8OP98eOcdu+UgJXRxUlJi8+bASpYcdVTYeEREROLitNNgk01ytpdOCV2clJXZnq1g9zm8/FpERKRRtW4NZ50Fzz6bkwWGldDFSfr+rHl5Ob38WkREpNH99rc2/HrrraEjaXRK6OJk//3tPj8feva05dgiIiJSP506weDB8MADsGhR6GgalRK6uFixwmrPHXmkzZ+bPHn1HjsRERFZtwsvhB9/hHvuCR1Jo1JCFxdPPw3ffANnnx06EhERkfjaaScYMABuvx1WrgwdTaNRQhcXd90F3brZN6GIiIisv5NOgrlzoUULKweWA4sklNDFwaRJMHasrc7RNl8iIiIb5sYb7b6qCqZNs7JgMafsIA7uvhuaNYNTTw0diYiISPyll/2qqsqJMmBK6KLuhx/gkUfghBNgs81CRyMiIhJ/RUU1I17O5UQZMCV0Uffoo7BkCfzmN6EjERERyQ2lpVb+yzm7PfVU6Ig2mBK6KPPehlt33hl23TV0NCIiIrmhsNDKf73yig25TpoUOqINpoQuyt5+277Jzj7b/oIQERGRxtO/P3TuDA89FDqSDaaELorKy20Z9d572xj/HnuEjkhERCT35OfDKafAqFFWxiTGlNBFUUmJLaP23m4nnBA6IhERkdw0ZIgNuz72WOhINogSuigqK7NvLrCELgeWU4uIiERSjx6w55427Op96GjWmxK6KEpfPp2XlxPLqUVERCJryBCYMgUmTAgdyXpTQhdFd99t987ZsurS0rDxiIiI5LLjj4fmzWO9OEIJXRS9957dz5hhy6oLC8PGIyIikss22QSOOgoefxxWrAgdzXpRQhdFTzwBffsqkRMREcmWIUNg0SJ4/vnQkawXJXRRU14O48db96+IiIhkx4AB0KlTbIddldBFzZNP2v1xx4WNQ0REJEny82HQIOuhKyiwerDl5aGjqjcldFEzYgTstht06RI6EhERkWQZPdruKyutHmxJSdh4GkAJXZTMmAEffKDhVhERkRA+/7zmuKoqVnVgldBFyYgRdq/hVhERkewrKqrZOz1mdWCV0EXJiBFWrbpz59CRiIiIJE9pac2Upy22iFUdWCV0UVFWBh99pOFWERGRUAoLbfrTFltAv36xKh+mhC4qRoywbt5jjw0diYiISHLl5dliiJdfhpUrQ0dTb0roQisvt6XRf/4ztGgR2wrVIiIiOaOkBH74Ad58M3Qk9aaELrSSElsaDbBsWayWSIuIiOSkAw6wvV01h07qrazMlkYDeB+rJdIiIiI5qWVL2zmitNR+N8eAErrQ0pdEx2yJtIiISM4qKbG6dJMnh46kXpTQhVZde8456NkzVt27IiIiOevww+0+Jr+XldCF9uWXdv/ii/ZXQIyWSIuIiOSsTp2gb18ldOvDOTfQOVfmnJvhnLu0lucvdM5Ncc597Jx7zTm3Tdpzlc65ianbyOxGvgFefRWaNoW99w4diYiIiKQrKYF334Wvvw4dyTpFJqFzzuUDdwKHAMXAYOdc8RqXfQj09d7vCDwF/F/ac8u8971Tt0FZCboxjB5tu0O0ahU6EhEREUlXUmKLIl54IXQk6xSZhA7YFZjhvS/33q8EhgNHpF/gvR/jvV+aOn0X2CrLMTauBQtg4kRbSSMiIiLRstFGUFAAp59uNWPLy0NHVKcoJXRbArPTzuekHqvL6cBLaefNnXPjnXPvOueOzEB8je/11+1eCZ2IiEj0DBoElZXWSzdtWqRrxTofkfoqzrljgYHe+zNS56cAu3nvz63l2pOBc4F9vfcrUo9t6b3/0jlXCLwOHOC9n1nLa88EzgTo0KFDn+HDh2esTQBLliyhdevWtT7X48Yb2fyNNxj73HP4/PyMxhHC2tqeBEluf5LbDsluf5LbDslufy62fd8DDsBV14oFfF4eb7722k+uy1bb+/fvP8F737fWJ733kbgBewCj0s4vAy6r5boBwFRg87V8rYeAY9f1nn369PGZNmbMmNqfqKryfpttvD/yyIzHEEqdbU+IJLc/yW33PtntT3LbvU92+3Oy7cXF3ufleQ/eO2fntchW24Hxvo6cJkpDru8D3Z1zXZ1zTYETgdVWqzrndgbuAQZ5779Oe7ytc65Z6rgd0A+YkrXI10d5OcyapeFWERGRqCottRqxYIsXI1zCJDIJnfe+AhtGHYX1wI3w3k92zg11zlWvWv0b0Bp4co3yJNsB451zHwFjgOu999FO6EaPtnsldCIiItFUWGg1Yk87DZo1g65dQ0dUp4LQAaTz3r8IvLjGY39OO641+/Hevw3skNnoGtmrr0LnztCjR+hIREREZG322gsefND2W6/usYuYyPTQJUplpa1wHTDAtvwSERGR6OrXz+7feitsHGuhhC7bysuhe3dYtAheeSXSNW1EREQEG01r1w7Gjg0dSZ2U0GVbSQl8/rkdz5sX6Zo2IiIigo2m9eunhE7SlJVZgUKAqio7FxERkWjr1w8+/RS++ip0JLVSQpdtRUU1x3l5q5+LiIhINFXPo3v77bBx1EEJXbY99pjdO2crZSJc00ZERERS+vSx0iURXRgRqbIlibBokd2/+CIMHBg2FhEREamfZs3gZz+L7Dw69dBl24QJdt+nT9g4REREpGH69YMPPoClS0NH8hNK6LLtgw+soHD79qEjERERkYbo1g1WrYI2baBXr0iVHlNCl20TJqh3TkREJI5uvNHuq6pg2rRIlR5TQpdN338P06fDLruEjkREREQaasaMmuOIlR5TQpdNEyfavXroRERE4ifCpceU0GWTFkSIiIjEV2mpbQEGlsxFqPSYErps+uAD6NQJOnQIHYmIiIg0VGEh3HabHY8YYecRoYQum7QgQkREJN6Ki+1+ypSwcaxBCV22LFliK2KU0ImIiMRXjx42f04JXQKVl8MOO4D3cP/9kapbIyIiIg3QooUNtSqhS6CSEpg1y46//DJSdWtERESkgYqLldAlUlmZ9c5B5OrWiIiISAMVF1td2VWrQkfyP0rosiHCdWtERESkgYqLLZmbOTN0JP+jhC4bSkuhaVNwDnr2jFTdGhEREWmgCK50VUKXDV262P3FF8PkyZGqWyMiIiIN1LOn3SuhS5h582DlSiVyIiIiuaBVK+usUUKXMNVlSrp2DRuHiIiINI7iYht1iwgldNlQndCph05ERCT+ysth3Dj4+GMoLqb53LmhI1JClxWffWYLIrbeOnQkIiIisqFKSmDhQjsuK2OHyy8PGw9K6LKjvBy22gqaNQsdiYiIiGyoNerLtpw9O2w8KKHLjs8+0/w5ERGRXFFUZHVlAZxjaefOYeNBCV12lJdr/pyIiEiuKC2F7t3tePPNmXTttWHjQQldxuWtXAlz56qHTkREJFcUFsLUqdC8OZxyCss7dQodkRK6TGs+f74dqIdOREQkdzgHHTtardkIUEKXYf9byqweOhERkdzSsaONwkWAErpMKi+n6G9/s+Nf/rKmHp2IiIjEn3roEqKkhKbVdWpmzLC6NSIiIpIblNAlRFkZrvq4qsrq1oiIiEhu6NgRvvuOvBUrQkeihC6jiorw1cd5eVa3RkRERHJDx44ANP3228CBKKHLrNJSqpo3t+OePa1ujYiIiOSGVLmSKCR0BaEDyGmFhSzdaivaKJkTERHJPakeumbV8+UDUg9dhhUsWQKbbBI6DBEREWlsGnJNjiY//ABt24YOQ0RERBrbd98B0O2OO6BXr6DlyZTQZVJlJQU//qiETkREJBcdcQQAznuYNi1oeTIldJn0/fd2r4ROREQk96SXIwtcnkwJXSYtWmT3mkMnIiKSe9LLkQUuT6aELpOqEzr10ImIiOSe0lJo1cpqzgauaKGyJZm0eLHdK6ETERHJPYWF0L8/S8rKaDN5ctBQ1EOXSeqhExERyW0tWpC3cmXoKJTQZZTm0ImIiOS25s2V0OW08nL44x/t+KCDgtamERERkQxRQpfjSkpgwQI7nj49aG0aERERyZDmzclbsSJ0FEroMqasDLy348C1aURERCRDNIcuxxUVgXN2HLg2jYiIiGRI8+bkr1xZ04kTiBK6TCkthbZtI1GbRkRERDKkeXO7DzzsqoQuUwoL4ZRTqGzVCiZPtnMRERHJLS1a2P3y5UHDiFRC55wb6Jwrc87NcM5dWsvzzZxzT6SeH+ec65L23GWpx8uccwdnNXARERFJpuoeurgkdM65fzvnDnfOZSQJdM7lA3cChwDFwGDnXPEal50OLPLedwNuAW5IvbYYOBHoBQwE7kp9PREREZHMqU7oli0LGkZDkrMfgSeAOc65a51z3Rs5ll2BGd77cu/9SmA4cMQa1xwBDEsdPwUc4JxzqceHe+9XeO8/A2akvp6IiIhIZrzxBvzqV3ZcVGTngdQ7ofPe/xzoCPwVGACUOef+45z7hXOuRSPEsiUwO+18TuqxWq/x3lcA3wGb1fO1IiIiIo3n4IOhosKOV62y80AKGnKx9/574G7gbudcL+AM4B7gNufcE8Ct3vupjR9m43HOnQmcCdChQwfeyGA23W3OHDpARt8jypYsWZLYtkOy25/ktkOy25/ktkOy25/Etu+7ciUu7dyvXMmbgf4NGpTQVXPOdcKGOQ8HKoCngc7Ax865y7z3N67Hl/0y9TWqbZV6rLZr5jjnCoCNgW/r+VoAvPf/BP4J0LdvX7/ffvutR6j19O9/UwFk9D0i7I033khs2yHZ7U9y2yHZ7U9y2yHZ7U9k25s2hbSiwq5p02D/Bg1ZFNHEOXesc+5FYBZwJPB/QEfv/ene+0OBY4A/rmcs7wPdnXNdnXNNsUUOI9e4ZiQwJHV8LPC6996nHj8xtQq2K9AdeG894xARERFZt1GjoCDVN9akiZ0H0pAeunmAAx4HLvXef1zLNf8BFq1PIN77CufcucAoIB94wHs/2Tk3FBjvvR8J3A884pybASzEkj5S140ApmA9hud47yvXJw4RERGRetlvP7j9djj7bPjiC9hii2ChNCShuwB40ntfZ6EV7/1ioOv6BuO9fxF4cY3H/px2vBw4ro7XXgNcs77vLSIiItJg1fXnqsuXBFLvhM57/0gmAxERERGJner6c4ETukjtFJFTysvhkUfI//FH6NXLzkVERCS3VPfQNWsWNAwldJlSUgKLFtly5mnT7FxERERyy/LlVDZtCs6t+9oMUkKXKWVl4L0dV1XZuYiIiOSWZcuoato0dBRK6DKmqKgmW8/Ls3MRERHJLcuXUxV4uBWU0GVOaWnN8uVtt7VzERERyS3Ll6uHLqcVFsJtt9nx00/buYiIiOQWJXQJ0Lat3S9ar1rLIiIiEnWaQ5cAm2xi90roREREcpN66BJAPXQiIiK5q7wcxo5lo8mTg9ecVUKXSdUJ3eLFQcMQERGRDCgpgWXLIlFzVgldJm20kd2rh05ERCT3pNeYDVxzVgldJuXnU9GqlRI6ERGRXJReYzZwzVkldBm2qk0bJXQiIiK56KGHAPDOQc+eQWvOKqHLsIrWrTWHTkREJBeldoiYcuWVMHly0JqzSugyrEI9dCIiIrlp3jwAVmy6aeBAlNBlVnk5bcrKYOzY4MuZRUREpJHNnQvAys02CxyIErrMKikhf+lSOw68nFlEREQaWaqHTgldrisrs9o0EHw5s4iIiDSyefNg442pSs2lC0kJXSYVFdnKFwi+nFlEREQa2bx50LFj6CgAJXSZVVrKivbt7XjrrYMuZxYREZFGpoQuIQoL+fhvf7Pjv/416HJmERERaWTz5kGnTqGjAJTQZdzyLbawA61wFRERyR3eq4cuSaqaNrXs/bPPQociIiIijeW772D5ciV0iVJYqB46ERGRXFFeDn372vHNN9M8VY8uJCV02dC1q3roREREckVJSU1Hzdy57HD55WHjQQlddhQWwpw5sGJF6EhERERkQ5WV2Rw6AO9pOXt22HhQQpcdXbvaB//FF6EjERERkQ1VVARpdWaXdu4cNh6U0GVHdbkSzaMTERGJv9JSaNXKjnv2ZNK114aNByV02VGd0GkenYiISPwVFkK7djB4MEyezPII1KJTQpcNHTtCs2bqoRMREckFP/4In38OxcWhI/kfJXTZkJcHXbqoh05ERCQXTJtm90roEqa83Fa5PvUU9OqlnjoREZE4mzLF7pXQJUxJCSxdasfTptm5iIiIxNOUKdCkCWy7behI/kcJXTak16upqrJzERERiacpU6BHD0vqIkIJXTak16txzs5FREQknqZMidRwKyihy47S0pokrl07OxcREZH4WbbM5sJHLKErCB1AIhQWwtSp9uF361ZTl05ERETiZfp0mz4VsYROPXTZ1KcPfPBB6ChERERkfUVwhSsoocuuXXaBL7+Er74KHYmIiIisjylTID8funcPHclqlNBlU58+dj9hQtg4REREpOHKy+Hvf4fKSuukiVBdWSV02dS7t91r2FVERCR+Skrghx/sOGJ1ZZXQZdNGG1ndGvXQiYiIxE96HdmI1ZVVQpdtffoooRMREYmjTp1qjvPyIlVXVgldtu2yC8yeDQsWhI5EREREGuLQQ+0+Px969oxUXVkldNlWvTBC8+hERETiZfJk2G03qKiw4wjVlVVCl20772z3SuhERETiY8UKeP996NcvdCS1UkKXbQsX2ma+V1wBvXpFasmziIiI1GHCBEvq9tordCS1UkKXbSUlsGoVeB+5Jc8iIiJSh7Fj7X7PPcPGUQcldNkW4SXPIiIiUoexY213iA4dQkdSKyV02VZUBM7ZsXORWvIsIiIitfDeErqIzp8DJXTZV1pqS50BNt44UkueRUREpBbTp8M33yihkzSFhbax7zHHQJs20LVr6IhERERkbarnz0V0QQREJKFzzm3qnBvtnPs0dd+2lmt6O+fecc5Nds597Jw7Ie25h5xznznnJqZuvbPagPUxYIAVGP7009CRiIiIyNq89RZstlmkp0lFIqEDLgVe8953B15Lna9pKfAL730vYCBwq3Nuk7Tnf++97526Tcx0wBvswAPt/tVXw8YhIiIiazd2rK1urZ4DH0FRSeiOAIaljocBR655gfd+uvf+09TxXOBroH22Amx0hYXQpYsSOhERkagqL7deuenTYdy4SNeOjUpC18F7Py91PB9Y65pg59yuQFNgZtrD16SGYm9xzjXLUJyNxzkbdn39daisDB2NiIiIrKmkpGZq1IIFka4d67z32Xkj514FtqjlqSuAYd77TdKuXeS9/8k8utRzHYE3gCHe+3fTHpuPJXn/BGZ674fW8fozgTMBOnTo0Gf48OHr26R6WbJkCa1bt671ufavv06vv/6VCXfeyQ/FxRmNI4S1tT0Jktz+JLcdkt3+JLcdkt3+XGz7vgccgKuq+t+5z8vjzdde+8l12Wp7//79J3jv+9b6pPc++A0oAzqmjjsCZXVctxHwAXDsWr7WfsDz9XnfPn36+EwbM2ZM3U9+/bX34P3VV2c8jhDW2vYESHL7k9x275Pd/iS33ftktz8n215cbL+nwfu8PDuvRbbaDoz3deQ0URlyHQkMSR0PAZ5b8wLnXFPgWeBh7/1TazzXMXXvsPl3n2Qy2EbTvj307q15dCIiIlF0/fV2n5dnNWQjXDs2Kgnd9cCBzrlPgQGpc5xzfZ1z96WuOR7YBzi1lvIkjznnJgGTgHbA1VmNfkMceCC8/Tb8+GPoSERERCTdhx/anPd582DyZFvQGFEFoQMA8N5/CxxQy+PjgTNSx48Cj9bx+v0zGmAmDRgAf/ub1bg5+ODQ0YiIiEi10lLYfXfYfPPQkaxTVHrokmuvvaBpUxg9OnQkIiIiUm3uXBg/PtIrW9MpoQtt/nwoKICbboJevSJd40ZERCQxnn/e7pXQSb2UlMCyZXY8bVpsvnFERERyWmmpbQDQq1foSOpFCV1oZWW2IBqgqsrORUREJJylS60CRUlJpLf7SqeELrSiIlsOnX4uIiIi4bz2GixfHqtRMyV0oZWWWm2b6r8A7r03bDwiIiJJV1oKbdrAvvuGjqTelNCFVlhotW2mTrXz998PG4+IiEiSVVVZQjdwoFWhiAkldFFRVAQ77QQjRoSOREREJLkmTLAKFDEabgUldNFy/PG2a8Ts2aEjERERSZ7ycjj0UDu+5ppYlRJTQhclxx9v908+GTYOERGRJCopgW++seNPP41VL50Suijp1g122UXDriIiIiGklw6LWSkxJXRRc/zxMG4cfP556EhERESSZZNNao7z8mJVSkwJXdQcd5zda9hVREQkeyorIT8fWre2+549bbVrTCihi5rCQujbV8OuIiIi2fTqq/D11/DQQ1BRYSXFCgtDR1VvSuii6IQTYPx4mDkzdCQiIiLJMGwYtG0Lhx8eOpL1ooQuijTsKiIikj2LF8Ozz8JJJ0GzZqGjWS9K6KKoshKaN4fLLoNevWJVB0dERCR2RoywvVtPPTV0JOtNCV0UlZTAihV2PG1arOrgiIiIxM6wYVBcDH36hI5kvSmhi6KyMvDejmNWB0dERCRWpk+3XZpOPRWcCx3NelNCF0VFRVb/plr37uFiERERyWXDhtnv3JNPDh3JBlFCF0WlpVb/pjqpO+WUsPGIiIjkospKePhhOPhg6NgxdDQbRAldFBUWWv2bigrYYQd46qmaIVgRERFpHGPGwJw5sV4MUU0JXZQ5B2efDR9+CO+9FzoaERGR3DJsmG33NWhQ6Eg2mBK6qDv5ZNuG5K67QkciIiKSG8rLYbvt4NFHrfNk7tzQEW0wJXRR16aNzaF74gn49tvQ0YiIiMRfSUlNBYnFi3OiPJgSujg4+2yrS/fQQ6EjERERib/08mDe50R5MCV0cbDDDtCvH/zjH1aXTkRERNbfllvWHOflWbmwmFNCFxe/+Q3MmAGvvho6EhERkXjbaisoKID8fCsTVloaOqINpoQuLo45Btq3h7vvDh2JiIhIfH30ke0M8de/WnmwyZOtXFjMKaGLi2bN4Je/hJEjrWaOiIiINNzNN0OrVvDrX4eOpFEpoYuTQw6xOXRbbw29etmyaxEREamfuXPhX/+yDpK2bUNH06iU0MXJb35j997DtGk5scxaREQka26/3bb7Ov/80JE0OiV0cZK+rLqqKieWWYuIiGTFkiVWLeKoo3JiztyalNDFSVGRLa+u1qNHuFhERETi5MEHrYjwRReFjiQjlNDFSWmpLa+uTupOPz1sPCIiInFQWQm33gp77GG3HKSELk4KC2159apVsOOOcO+99k0qIiIidfv3v20hYY72zoESunjKy4M//tHm0D31VOhoREREou2mm6xT5MgjQ0eSMUro4uqYY2C77eDqq7UdmIiISG3Ky6FrV3jnHVi6FGbNCh1Rxiihi6u8PLjiCvjkE3juudDRiIiIRE9JCXz+uR1//XVOl/tSQhdnJ5wA3brZ9iXeh45GREQkWhJU7ksJXZwVFMDll8OHH8KLL4aORkREJFpatKg5zsuz8l85Sgld3J18MnTpol46ERGRdGPHWjHhzTeH/Hwr+1VaGjqqjFFCF3dNmsCll8K4cfDqq6GjERERCc97uOwy2GILWxhRUWFlv3Jwh4hqSuhywamnwlZbWS+diIhI0r38Mvz3v/CnP0GrVqGjyQoldLmgWTPbNeK//7Vu5V697C8SERGRpKmqst65wkI444zQ0WSNErpcMWKE3VdVwbRpOb00W0REpE5PPAEffQRDh0LTpqGjyRoldLli+vSa4xxfmi0iIlKrVatsmHXHHWHw4NDRZFVB6ACkkRQVWc9c9a4RPXqEjUdERCTb7r8fZs601ax5yeqzSlZrc1lpqS3Jrv4GPuGEsPGIiIhk09KlNszarx8cdljoaLJOCV2uKCy0JdkVFbD//nDbbbBoUeioREREsuP222HePLjuOnAudDRZp4Qu1zgHt9wCixfbXyoiIiK5btEiuP56OPRQ2Hvv0NEEoYQuF+24I/zqV3DHHTavTkREJBeVl1uprs02s46Ms88OHVEwSuhy1dCh0LIlXHxx6EhEREQyo6TEOi6qt7685JKw8QQUiYTOObepc260c+7T1H3bOq6rdM5NTN1Gpj3e1Tk3zjk3wzn3hHMuOYVn6rL55rZ0+4UXYNSo0NGIiIg0vrKymuoO1ecJFYmEDrgUeM173x14LXVem2Xe+96p26C0x28AbvHedwMWAadnNtyY+O1vYdtt4cILbbGEiIhILtlqq5rjvDwr4ZVQUUnojgCGpY6HAUfW94XOOQfsDzy1Pq/Pac2awU03wZQpcM89oaMRERFpPCtW2ELAJk1s28uePa2EV0JFJaHr4L2flzqeD3So47rmzrnxzrl3nXNHph7bDFjsva/ugpoDbJm5UGNm0CArY/LnP6uMiYiI5I6//Q0+/xxGjrRRqMmTrYRXQjlfPZEw02/k3KvAFrU8dQUwzHu/Sdq1i7z3P5lH55zb0nv/pXOuEHgdOAD4Dng3NdyKc64z8JL3fvs64jgTOBOgQ4cOfYYPH75hDVuHJUuW0Lp164y+x7q0mjmTvmeeyZyjj2bmOedk7X2j0PaQktz+JLcdkt3+JLcdkt3+bLa9+Zdfsutpp/HNnnsy5aqrsvKea5Ottvfv33+C975vrU9674PfgDKgY+q4I1BWj9c8BBwLOOAboCD1+B7AqPq8b58+fXymjRkzJuPvUS8nnug9eJ+f731xsfczZ2b8LSPT9kCS3P4kt937ZLc/yW33Ptntz1rbq6q8HzjQ+9atvZ8zJzvvuQ7Zajsw3teR00RlyHUkMCR1PAR4bs0LnHNtnXPNUsftgH7AlFQDx2DJXZ2vT7wPPrD7ykpb4l1SEjYeERGR9fH00/Dyy3D11bClZlhVi0pCdz1woHPuU2BA6hznXF/n3H2pa7YDxjvnPsISuOu991NSz10CXOicm4HNqbs/q9HHwcyZNcdVVYle2i0iIjH1/ffwu99B796QxSlEcVAQOgAA7/232Hy4NR8fD5yROn4b2KGO15cDu2YyxtgrKrKeuep6PdtuGzYeERGRhrryStuv9ZlnoCASKUxkRKWHTjKttNSWdOfn2zLvLl1qKmuLiIhE3cSJcNtt8Otfw267hY4mcpTQJUVhoS3prqiwDYxfeQWefDJ0VCIiIutWVQVnnQXt2sG114aOJpKU0CXRhRfCz35m8w8WLAgdjYiIyNrdey+MG2fF8tvWujto4imhS6KCAnjgAfjuO9seTEREJIrKy20O+FlnQcuWsMceoSOKLCV0SbX99rZ7xBNPwLPPho5GRETkp0pKYPp0O16+3HY/klopoUuySy6BnXeGs8+Gb78NHY2IiMjqpk2rOVbJrbVSQpdkTZrAgw9aMnf++aGjERERqfHZZ6uf5+XZ8KvUSgld0u20E1x+OTz6KDz/fOhoRERErCLDySdDq1bQvbuV3OrZ00pwSa2U0AlccYXNqfv1r2Hx4tDRiIhI0l13Hbz9NvzjHzaHrqLCSm8VFoaOLLKU0Ak0bWpDr/PnQ9eutgq2Vy9bXSQiIpJN48bBX/4CJ51kN6kXJXRi+va12j6LF0NlpU1ELSkJHZWIiCTJkiU21LrllnDnnaGjiRVthCY10odbtZpIRESy7fzzYeZMeOMN2GSTwMHEi3ropEZRka0iqtajR7hYREQkWZ59Fu6/Hy69FPbZJ3Q0saOETmqUltoqouqkbt99w8YjIiLJMHcunHEG9OkDV10VOppYUkInNQoLbRVRZaX9YP3jH/DSS6GjEhGRXFZVBaeeajtBPPaYLdSTBlNCJ7W77TbYcUebnDp7duhoREQkV912G4weDTffrMLBG0AJndSuRQt48klYuRJOOAFWrQodkYiI5JLycth2W7jgAmjdGgYMCB1RrCmhk7r16AH33QfvvAOXXRY6GhERySWHHVZT73TpUhg0KGw8MaeETtbuhBPgN7+Bm26C554LHY2IiOSCNUtjqVTWBlNCJ+t288228ujUU3+6WbKIiEhDXXUVeA/O2XlenubPbSAldLJuzZrBiBH2w3fCCbBiReiIREQkrp55Bv76Vzj2WNhuO8jPt5JZpaWhI4s1JXRSP4WFtt/r++/D738fOhoREYmjTz6BX/wCdtsNHn3USmVVVNh9YWHo6GJNCZ3U31FHwWmnwe23219UvXrVTGgVERFZm4UL4YgjoE0b66Vr1ix0RDlFCZ00zLvv2n1VFUybBiUlYeMREZHoq6iAE0+EOXMsmevUKXREOacgdAASM9On1xxrVZKIiNTHZZdZ8eD77oM99ggdTU5SD500TFFRzV6vAE2aWP0gERGR2jz2GNx4I5xzDpx+euhocpYSOmmY0lJbjZSfD5072957J59s+7+KiIikmzDB9gbfZx+45ZbQ0eQ0JXTSMIWFNauSvvgCbr0Vnn0WLrkkdGQiIhIlX39ti+nat7etJJs0CR1RTtMcOtkwv/sdzJxpO0lsuy2cfXboiEREJKTycjj8cJg61QoH//vfsPnmoaPKeeqhkw13yy32w3vuufDii6GjERGRkEpKLJmrpr3As0IJnWy4/Hz4179gp51sJ4mPPgodkYiIhJKezHmvaghZooROGkfr1vD887DJJnDYYfDll6EjEhGRbPv73y2Jq6Y9WrNGCZ00nk6d4IUX4Lvv4PDDyVc5ExGRxOgwejScfz4cfLD2aA1AiyKkce24o61mOvxwiocOhYMOggJ9m4mI5LQXXqDohhugf39bBNG8eeiIEkc9dNL4Bg6EO+9ks3HjbGVTQYH2fRURyVVjx8Jxx/FjYaGSuYCU0Elm/PrXrNxoI1i0yIoOa99XEZHcM2mSVTno3JmPb7gBNtoodESJpYROMqbJkiU1J9r3VUQkt3z2mc2Xa9kSXnmFVW3bho4o0ZTQScYs7dx59X1f27cPF4yIiDSer76CAw+07R9feQW22SZ0RImnhE4yZtK119bs+9qmDcyfD3ffHTosERHZEN99Z3Ol582zyga9eoWOSNAqV8mg5Z062b6vACtXwrHHwm9+Y/v5nXFG2OBERKThli2DQYPgk0+sHMkee4SOSFLUQyfZ0bSplTMZOBDOPBOGDQsdkYiI1Fd5ORQXQ6tW8J//wN/+Zv+fS2QooZPsadYMnnkGDjgATjsNHnssdEQiIlIfhx9uW3p5D87BvfeGjkjWoIROsqtFC3juOdh3X/jFL2DEiNARiYjI2ixfbqWnqml/1khSQifZ17Klzb3Yc0846SR49tnQEYmISG2WLrU5c9U9c6D9WSNKCZ2E0bo1vPgi7LornHCC9voTEYma77+3eXKvvQY33KD9WSNOCZ2E06YNvPQS9O4NxxxjdYy0TZiISHgLF1qduXfegccfhz/8waoWVFTYfWFh6AhlDUroJKyNN4ZRo+yvvi++0DZhIiKhff017L8/TJwITz9toygSeapDJ+G1bQurVtWca5swEZEw5s61SgSzZtmw6kEHhY5I6kk9dBINRUWrbxPWrp1NwhURkeyYNQv22QfmzIGXX1YyFzNK6CQaSktrtgnbaCPbJ/Cii6y3TkREMuvTT2HvveHbb+HVVy2xk1hRQifRUFhYM+F20SI47zy45RY4+WRYsSJ0dCIiuWvyZEvgli2DMWNgt91CRyTrQXPoJHry8uDWW6FTJ7j0Upug+8wz1nMnIiKNo7zcVrKWl1uFgeeft6oDEkuR6KFzzm3qnBvtnPs0dd+2lmv6O+cmpt2WO+eOTD33kHPus7Tneme7DdLInINLLoGHHoI33oD99oP58wMHJSKSQ/r3rykRVVUFF14YNh7ZIJFI6IBLgde8992B11Lnq/Hej/He9/be9wb2B5YCr6Rd8vvq5733E7MQs2TDkCE2v66szHaW+PTT0BGJiMTf7bdbqahqqi4Qe1FJ6I4AhqWOhwFHruP6Y4GXvPdLMxmURMQhh9i8jh9+gH794P33Q0ckIhJPFRVwzjk2T7lNm5rqAtrOK/aiktB18N7PSx3PBzqs4/oTgX+t8dg1zrmPnXO3OOeaNXqEEtauu8LYsdCqlQ0TvPxy6IhEROJl8WI47DC46y7b+eGDD2qqC2g7r9hzPku1vpxzrwJb1PLUFcAw7/0madcu8t7/ZB5d6rmOwMdAJ+/9qrTH5gNNgX8CM733Q+t4/ZnAmQAdOnToM3z48PVuU30sWbKE1q1bZ/Q9oioTbW+6cCE7XHIJrcrLWbnppjRbuJClnTsz6dprWd6pU6O+14bSZ5/MtkOy25/ktkN029987lx2uPxyWsyZw/QLLmD+YYc1+ntEte3ZkK229+/ff4L3vm+tT3rvg9+AMqBj6rgjULaWa38H/HMtz+8HPF+f9+3Tp4/PtDFjxmT8PaIqY23/7jvvW7b03koPe5+X531xcWbeawPos0+uJLc/yW33PqLt/+9/vW/XzvtNN/U+g/FFsu1Zkq22A+N9HTlNVIZcRwJDUsdDgOfWcu1g1hhuTfXQ4Zxz2Py7Txo/RImMjTZavTadJvOKiNTu4YdtK69NN4V337WKAZKTopLQXQ8c6Jz7FBiQOsc519c5d1/1Rc65LkBn4M01Xv+Yc24SMAloB1ydjaAloDW3CsvLg0+Ux4uIAPaH7hVXWKWAfv3gnXege/fQUUkGRaKwsPf+W+CAWh4fD5yRdv45sGUt1+2fyfgkgkpLoaTEeuY6d4YlS6y6+QMPwAknhI5ORCScpUvhF7+Ap5+GM86wRRBNmoSOSjIsKj10Ig2TvlXYZ5/Bxx/DzjvDiSfCxRfb4yIiSVJebqMXrVpZMnf55fDPfyqZSwgldJIbOnaE11+3+ko33QQHHQQLFoSOSkQkew44AKZPt2Pn4N//tntJBCV0kjuaNoU77rDtwt55B/r0gfHjQ0clIpJZlZVwzTXw+ec1j3mvxWIJo4ROcs+QIVaE2DnYay+bVycikovmz4eBA+GPf7QKANr5IbGU0Elu2mUXmDDBErrTT4ezz4aVK0NHJSLSeF59FXr3tj9g771XOz8knBI6yV3t2tkWYX/4A/zjH7DJJlBQAL162eRhEZE4qqiwHrmDDoLNNoP33rPVrNtuW7NYbPJkWzwmiaGETnJbQQHccANstRUsW2ZzTaZOtZInIiJxM2eO7Wd9zTVw2mmWzG2/feioJAKU0EkyzJtXc+y9JXU//BAuHhGRhnr+eRtinTgRHn0U7r/fSpSIoIROkiJ9ZwnnLKnbaSd4662wcYmIrMvKlXDRRTay0LmzzQ/++c9DRyURo4ROkqG0tGay8HbbwRNPWGK3zz5w6aWr7w0rIhJaebnN9y0ogI03hptvtjqb77wDPXqEjk4iSAmdJEP6zhKTJ8Pxx9uwxRln2By73XbTXrAiEh0lJTY1pLISli+3ecB33AHNm4eOTCJKCZ0kV5s2ti3OyJE2x65PH9tloqoqdGQikmTz5lky5/3qj4mshRI6kZIS65079FDbB3b//WHWrNBRiUjSeG+F0IuL7bx62y4VCZZ6UEInAtC+PTzzDDz4oBXn3HFHePjh1f9CFhHJlPJyqyt3+umwww4werTN91WRYKknJXQi1ZyDU0+Fjz+2FbBDhthkZBUjFpFMqayEW2+1JG7cOLjrLnjjDTjgABUJlgZRQieypi5dYMwY2Hxzq1WnYsQikgmTJ0O/fnDBBbDffnZ+9tk1JZZEGkDfNSK1yc+Hb7+tOfcepkyBSZPCxSQiuWHlShg6FHbeGWbMgMces6LBnTuHjkxiTAmdSF3WLEacl2f/AV90EXz/fdjYRCSe3nvPVtRfeSUce6z1/p90Us0CCJH1pIROpC5rFiN+7z2bsHzLLfb48OFaNCEi61Zebv+H5OVZzcsFC6xc0uOP24IskUaghE6kLmsWI+7TB+65B959Fzp1gsGDYcAA+wtbRKQ23tv8uGnT7Ng5aNtWc3Kl0SmhE2moXXetWY32wQe2IvbSS+HHH0NHJiJR8uGHtlp19uyax7yHTz8NF5PkLCV0IusjP99Wo5WV2SbZN9xgQyrPPKNhWJGkmzsXfvlL69X/+GPYYoua+bgqEiwZooROZENsvrkVI37rLRtGOeYY2Hdf6N6dfQ84QPXrRJLkxx/ZZtgw6N7dVq5efLGtYh07tmY+rooES4YUhA5AJCf06wcTJsCdd8KFF0JVFQ5s3kxJic3BE5HcVFUFjz4Kl19O1y+/tNWr118P225rz2+yif4PkIxTD51IYykogN/9bvXyA1VVtmhiyZJwcYlI5rz5JvzsZ7azTKdOfHjbbfDkkzXJnEiWKKETaWzp9evA5tQVFtr2PsuXBwtLRBrRp5/C0UfbCtYFC2yI9d13+W6HHUJHJgmlhE6ksaXq1/m8PCguhqeftpWwF1xgc2vuvRdWrQodpYg0VHm5zYHLy4MePWDUKLjmGlscddJJ2rJLgtJ3n0hjS9Wve/O112zezNFHw+jR8NprsNVWcOaZtiL28cdtSFZEom/RIth9d0veqleyb7UVXH45tGgRNjYRlNCJZM/++8Pbb1sPXqtWVu6kd2947jmVOhGJqgULLGnbZhs7TjdzZpiYRGqhhE4km5yDww+3gqPDh9ucuiOPtL/8H3nEypwUFKjciUhoX35p0yS22cZWrB5yiPW+q56cRJQSOpEQ8vLghBNgyhS47z6YNw9+8Qs7r6ysKXciItn12Wdw1lmWvN1+Oxx/vP1cPvGETZ1QPTmJKNWhEwmpoABOPx1OPhlatqyZU1dVZUldRYVdIyKZNW0aXHedrVbNz7edHv7wB+jateaa6v2dRSJIPXQiUdCsWc3quWpVVVbL6uab4bvvwsUmkss++sh6y4uLrX7cb39r0x3uvnv1ZE4k4pTQiURFqtwJ+fm2Cvaee+wXykUXQefONp/ns89CRykSb+XlNkc1Px/atLGFSS+9BJdeCp9/DrfcAltuGTpKkQbTWI5IVNQ2nHPmmbal2C23wB13wG23WRmUCy+EPfYIE6dIXK1YAfvsYwsewHZwad/eSpG0bRs2NpENpB46kajr08f2ifzsM/j97+HVV2HPPS2he/JJm2cnInUrL7ceuM6da5K5agsXKpmTnKCETiQuttrKyifMnm2r7xYssBV43brBH/9ow7QqeSJiKiqsxuMhh9jPyI03wl57wdZbq/SI5CQldCJx07o1nHuuDRM9+6zVybrmGlulV1kJU6eq5Ikk15dfwtChNv/0yCNh0iS48kqYNQueeQbGjFHpEclJmkMnElf5+fYL68gj7bi65In3VjfruuvglFOsZ08kl1VV2dZ6//iH9cpVVsLBB1tP9uGHr176R6VHJEeph04kF6SXPHHOatpdfrkNLx10kNXWWro0bIwijaV6pWpBAWyxhSVpBx0E//mPrQqfMQNeftn+2FEdR0kIJXQiuWDNkieTJtkvtT/9CaZPt8LFW2wBZ5wB//2v9o6V+PrxR9h335pdVb76Cr7+2v5omTMHbrjB6jeKJIwSOpFcUD2MVFFh94WF9kvtL3+x3owxY+CYY2z/2H32sUniQ4da3S2RqFu6FJ56yhYBtW9viVu6lSvhpJOsQLdIQimhE8l1eXmw337w4IMwfz4MGwZdusBVV9nE8d12s0KqWiErUbJsmS36GTwYNt8cjjsO3nwTTjvNFgJpparIapTQiSRJ69bwi1/YBPLPP4err4YPP4S5c234asoUq3GnpE5CWLECRo60KQKbb25FtF991c5fe82+T++8E15/XStVRdag2aIiSbX11nDFFVbSId1XX9lwba9eMGiQlUDZdVf75SnSWMrL7XurrAw6dbIC2q+/Dt9/D5tuCieeaEOs/fv/dGGDVqqK/IR66ESSrqho9eGrbt3g5puth+T//s967Dp1gl/+Ev79b5uULrIh5syBvfeuWdgwe7b1zB1zjO2rOn8+3HsvHHigVqmK1JMSOpGkS18h27MnjBoFF1xgvSULFsDjj8P++1tR1qOOgs02g8MOs5pfb79dUz5C8++kLj/8YN9n551nq7A7d7bh03TOwQMPwMCB0KRJmDhFYkx/+ogk3dqGr9q2tUnpgwfDqlVW8qS01HpTXnxx9WunTrVEb+rUzMcs0VZRAe+/D6NH2+3dd+2xFi1slfUZZ8Bdd9k8zqoqLWwQaQRK6ESkfpo0sZ66/fe3IdmpU2H77Wtq2nlv24/17Uu3Ll3g229t78wOHYKGLRmUPg+ua1c45RR6jR5tdRC/+8563XbZBS6+2IZP99wTmje31x51VM1ri4q0sEFkAymhE5GGcw6Ki234bNo062VxzoZj27ShY2kpPP20Xdu9uyV2e+9t99262bUSX6tWWa/uoYfCvHn22IwZcOWVtOnQwUqMHHigJf/t2tX+NbSwQaRRKaETkfVXWvrTXpbCQt4aPZp927SBt96yYdrnnrM6eGA9dtUJXpcucNlltptF2uslQiorrTd2/Pia28SJVmJkTfn5vPuvf7Ff//5ZD1Mk6ZTQicj6q6OXxTdpArvvbreLL7YevGnTahK8t96q6cGrNnUq9Otnu1lsv7319klmpQ+ZFhXVbGyfnrx98EHNPsCtW1t5kXPPhb59bWu58vLV58Gp91UkCCV0IpJ5eXk2RFtcDGeeaY/NmWMV/6uq7Nx7K1ex3352vsUWltil34qLoU2bIE3IOYsWwYABtjDBeyshUlRU83m0aAE77wy/+pUlb337Qo8eNSVuwOoTrtlD+8UXQZojknSRSOicc8cBVwHbAbt678fXcd1A4O9APnCf9/761ONdgeHAZsAE4BTv/coshC4i62urraxMSvUcvLw8K2h8++3wySc1t3vusW2gqnXpUpPgbb453HEHzJqlIds1eW+b1s+cafPb1rwtWlT7ax54wJK37bZbdw242npoldCJBBGJhA74BDgauKeuC5xz+cCdwIHAHOB959xI7/0U4AbgFu/9cOfcP4DTgbszH7aIbJA65uBx8ME111RWWi9SepL3ySfw8stWCqPalCmWhOyzj+1Nu9VWP71tttnqQ4JrDjnGKSH03laTHnusJW0dOsARR8A339QkbUuW1Fyfl2c9ot262S4M3brBrbfCl1/WJNQ9e9peqSISO5FI6Lz3UwHc2ude7ArM8N6Xp64dDhzhnJsK7A+clLpuGNbbp4ROJOrqs9IxP9967rbd1hKWaitXQsuWlvBVW7XKkpjqfT+rhw+rNWu2erI3ahQsXGjJ0dSpts3U44/bXLH0W4sWqw81VtvQhDD99d272+4ILVtaUrbm7dtvf/rYyrSBiHnzrNhzjx6WrO2zj91X37bZBpo2Xf39jzxSpUNEckQkErp62hKYnXY+B9gNG2Zd7L2vSHt8yyzHJiLZ1rSpJSHpQ7Y9e8I779jzFRW2L+2cOdYLNWfO6rd337UkqZr3Nly4114/fS/nLNFaM9GbOLFmK7QpU2CnnaBfP3ZasMDm+lVUrP325Zc1vYzTptnK39ree7PN7NaundV7+9nP7Pj//q+mDiBY8ltWVv9/Q5UOEckZzqf/Z5DJN3LuVWCLWp66wnv/XOqaN4CLa5tD55w7FhjovT8jdX4KltBdBbzrve+Werwz8JL3fvs64jgTOBOgQ4cOfYYPH76BLVu7JUuW0Lp164y+R1Qlue2Q7PZnq+3N585lh8svp+Xs2Szt3JlJ117L8k6d6v36n516Ki2++II876lyjhUdOjD9wgvJX7aM/OXL7b62W+q5TT74gPRxBQ/8sN12VAB5TZvi8/PXeuvwyiu4tP+DvXN8MnQoqzbe+H+3itatLVGrR/zLtt6a9x96aL3+LRtLkr/vIdntV9sz3/b+/ftP8N73rfVJ731kbsAbQN86ntsDGJV2flnq5oBvgILarlvbrU+fPj7TxowZk/H3iKokt937ZLc/Nm2fOdP74mLv8/PtfubMhr2+uNj7vDzvwe6Li733DWh/Ha/PWvwZEJvPPkOS3H61PfOA8b6OnKaWSSGR9T7Q3TnX1TnXFDgRGJlq4Bjg2NR1Q4DnAsUoInFSPeRYUWH3DV0QUVpqw7z5+Xbf0DloG/r6DY1fRHJGJObQOeeOAm4H2gMvOOcmeu8Pds51wsqTHOq9r3DOnQuMwsqWPOC9r578cQkw3Dl3NfAhcH+AZohI0mzoHDTNYRORRhKJhM57/yzwbC2PzwUOTTt/EXixluvKsVWwIiIiIokTpyFXEREREamFEjoRERGRmFNCJyIiIhJzSuhEREREYk4JnYiIiEjMKaETERERiTkldCIiIiIxp4ROREREJOaU0ImIiIjEnBI6ERERkZhTQiciIiISc0roRERERGJOCZ2IiIhIzCmhExEREYk5JXQiIiIiMaeETkRERCTmlNCJiIiIxJwSOhEREZGYU0InIiIiEnNK6ERERERiznnvQ8cQjHNuATArw2/TDvgmw+8RVUluOyS7/UluOyS7/UluOyS7/Wp75m3jvW9f2xOJTuiywTk33nvfN3QcISS57ZDs9ie57ZDs9ie57ZDs9qvtYduuIVcRERGRmFNCJyIiIhJzSugy75+hAwgoyW2HZLc/yW2HZLc/yW2HZLdfbQ9Ic+hEREREYk49dCIiIiIxp4SuETjnjnPOTXbOVTnn6lzl4pwb6Jwrc87NcM5dmvZ4V+fcuNTjTzjnmmYn8g3nnNvUOTfaOfdp6r5tLdf0d85NTLstd84dmXruIefcZ2nP9c52GzZEfdqfuq4yrY0j0x7P9c++t3PundTPx8fOuRPSnovdZ1/Xz3Da881Sn+OM1OfaJe25y1KPlznnDs5q4I2kHu2/0Dk3JfVZv+ac2ybtuVp/BuKiHm0/1Tm3IK2NZ6Q9NyT1c/Kpc25IdiNvHPVo/y1pbZ/unFuc9lzcP/sHnHNfO+c+qeN555y7LfVv87Fzbpe057L32XvvddvAG7AdUAS8AfSt45p8YCZQCDQFPgKKU8+NAE5MHf8DODt0mxrQ9v8DLk0dXwrcsI7rNwUWAi1T5w8Bx4ZuR6bbDyyp4/Gc/uyBHkD31HEnYB6wSRw/+7X9DKdd8xvgH6njE4EnUsfFqeubAV1TXyc/dJsy0P7+aT/bZ1e3P3Ve689AHG71bPupwB21vHZToDx13zZ13DZ0mxq7/Wtc/1vggVz47FPx7wPsAnxSx/OHAi8BDtgdGBfis1cPXSPw3k/13pet47JdgRne+3Lv/UpgOHCEc84B+wNPpa4bBhyZsWAb3xFYzFC/2I8FXvLeL81kUFnU0Pb/TxI+e+/9dO/9p6njucDXQK1FMWOg1p/hNa5J/zd5Cjgg9TkfAQz33q/w3n8GzEh9vThZZ/u992PSfrbfBbbKcoyZUp/Pvi4HA6O99wu994uA0cDADMWZKQ1t/2DgX1mJLAu89//BOiLqcgTwsDfvAps45zqS5c9eCV32bAnMTjufk3psM2Cx975ijcfjooP3fl7qeD7QYR3Xn8hPf9CvSXVT3+Kca9boEWZWfdvf3Dk33jn3bvVwMwn77J1zu2J/3c9MezhOn31dP8O1XpP6XL/DPuf6vDbqGtqG07Fei2q1/QzERX3bfkzq+/kp51znBr42yurdhtQwe1fg9bSH4/zZ10dd/z5Z/ewLMvWFc41z7lVgi1qeusJ7/1y248mmtbU9/cR7751zdS6bTv3FsgMwKu3hy7BkoCm27PsSYOiGxtyYGqn923jvv3TOFQKvO+cmYb/sI62RP/tHgCHe+6rUw5H/7GX9OOdOBvoC+6Y9/JOfAe/9zNq/QiyVAv/y3q9wzv0a66ndP3BMIZwIPOW9r0x7LNc/+0hQQldP3vsBG/glvgQ6p51vlXrsW6x7tiD1F33145GxtrY7575yznX03s9L/dL+ei1f6njgWe/9qrSvXd3Ds8I59yBwcaME3Ygao/3e+y9T9+XOuTeAnYGnScBn75zbCHgB++Pn3bSvHfnPfg11/QzXds0c51wBsDH2M16f10ZdvdrgnBuAJfz7eu9XVD9ex89AXH6pr7Pt3vtv007vw+aYVr92vzVe+0ajR5hZDfn+PRE4J/2BmH/29VHXv09WP3sNuWbP+0B3Z6sam2Lf9CO9zZwcg80tAxgCxKnHbyQWM6w79p/Mq0glAtXzyY4Eal1FFGHrbL9zrm31cKJzrh3QD5iShM8+9b3+LDa/5Kk1novbZ1/rz/Aa16T/mxwLvJ76nEcCJzpbBdsV6A68l6W4G8s62++c2xm4Bxjkvf867fFafwayFvmGq0/bO6adDgKmpo5HAQel/g3aAgex+ihFHNTnex/nXE9s8v87aY/F/bOvj5HAL1KrXXcHvkv9wZrdzz5Tqy2SdAOOwsbGVwBfAaNSj3cCXky77lBgOvaXyRVpjxdi/7nPAJ4EmoVuUwPavhnwGvAp8CqwaerxvsB9add1wf5ayVvj9a8Dk7Bf5o8CrUO3qbHbD+yZauNHqfvTk/LZAycDq4CJabfecf3sa/sZxoaJB6WOm6c+xxmpz7Uw7bVXpF5XBhwSui0Zav+rqf8Dqz/rkanH6/wZiMutHm2/DpicauMYoGfaa3+Z+p6YAZwWui2ZaH/q/Crg+jVelwuf/b+wFfqrsN/1pwNnAWelnnfAnal/m0mkVbvI5mevnSJEREREYk5DriIiIiIxp4ROREREJOaU0ImIiIjEnBI6ERERkZhTQiciIiISc0roRERERGJOCZ2IiIhIzCmhExEREYk5JXQiIhvIOdfeOTfPOXdl2mM7OueWO+eOCxmbiCSDdooQEWkEzrmDgVJgX2zbq/HAe97700LGJSLJoIRORKSROOduxTZmfxPYG9u3dknQoEQkEZTQiYg0EudcM2wT8u7Ant77cYFDEpGE0Bw6EZHG0wXoDHigMGwoIpIk6qETEWkEzrkmwLvAdGAccCWwk/f+i6CBiUgiKKETEWkEzrnrgZOAHYHvgJeA5sD+3vuqkLGJSO7TkKuIyAZyzu0LXAT8wnu/2NtfyqcCxcAlIWMTkWRQD52IiIhIzKmHTkRERCTmlNCJiIiIxJwSOhEREZGYU0InIiIiEnNK6ERERERiTgmdiIiISMwpoRMRERGJOSV0IiIiIjGnhE5EREQk5v4f0T8+7R7oSHIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 10))\n",
    "plt.plot(x, y, 'o-r', markersize=4)\n",
    "plt.xlabel('x', fontsize=14)\n",
    "plt.ylabel('y', fontsize=14)\n",
    "plt.grid()\n",
    "\n",
    "print(f'accuracy x(t): {sigma[0]}')\n",
    "print(f'accuracy y(t): {sigma[1]}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}