{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Step 2 : build a Runge-Kutta type solver\n", "\n", "📜 _Once obtained the following_ $Q$_-coefficients :_\n", "\n", "$$\n", "\\begin{array}\n", " {c|c}\n", " \\tau & Q \\\\\n", " \\hline\n", " & w^\\top\n", "\\end{array}\n", "$$\n", "\n", "_we can use those to build the associated time-stepping scheme (solver) and for time-dependent problems._\n", "\n", "> 📣 Remember, this is exactly the same as Butcher tables for Runge-Kutta methods ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the following simple ODE, usually named _Dahlquist problem_ :\n", "\n", "$$\n", "\\frac{du}{dt} = \\lambda u, \\quad t \\in [0, T], \\quad u(0)=u_0.\n", "$$\n", "\n", "We choose here $\\lambda=i$ (imaginary unit), $T=4\\pi$, and $u_0=e^{\\frac{i\\pi}{6}}$ :" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "lam = 1j\n", "tEnd = 4*np.pi\n", "u0 = np.exp(1j*np.pi/6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let say we want to solve it numerically, using $12$ time-steps of the famous RK4 method. This is how we do it with `qmat` :" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from qmat import Q_GENERATORS\n", "\n", "rk = Q_GENERATORS[\"RK4\"]()\n", "nodes, weights, Q = rk.genCoeffs()\n", "\n", "nSteps = 12\n", "uNum = np.zeros(nSteps+1, dtype=complex)\n", "uNum[0] = u0\n", "\n", "dt = tEnd/nSteps\n", "A = np.eye(nodes.size) - lam*dt*Q # all-at-once system\n", "\n", "for i in range(nSteps):\n", " b = np.ones(nodes.size)*uNum[i] # ... with its RHS\n", " uNodes = np.linalg.solve(A, b) # ... and its solution\n", " uNum[i+1] = uNum[i] + lam*dt*weights.dot(uNodes) # step update\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To explain a bit, we simply computed $N=12$ time-step solutions $u_1:=u(t_1), \\dots, u_n:=u(t_n) \\dots, u_{N}:=u(t_N)=u(T)$,\n", "with $t_{n+1}-t_{n} = \\Delta{t} = T/N$. One time-step consists on :\n", "\n", "1. solving the **all-at-once system** :\n", "\n", "$$\n", "(I - \\lambda\\Delta{t} Q) u_\\tau = [u_n, \\dots, u_n]^T,\n", "$$\n", "\n", "where $u_\\tau$ stores the numerical approximation at $t_n + \\Delta{t}\\tau_m$, also called the **nodes solution** (for collocation methods) or **stage solutions** (for RK methods).\n", "\n", "2. updating the step solution with the **step update** :\n", "\n", "$$\n", "u_{n+1} = u_{n} + \\lambda\\Delta{t} w^T u_\\tau\n", "$$\n", "\n", "... and that's it\n", "\n", "> 💡 The code is independent from the fact that we used the RK4 scheme $\\Rightarrow$ we can use it for any other scheme ... \n", "\n", "We show the time solution below, starting from the initial solution (orange square), and with the exact analytic solution in dashed line :" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAoARJREFUeJztnQVYVHkXxl+6uxUVkFTs7i7sdo21Nuxu1zU2XNu143PV1V1j7e7uFpMQUUQQQQlBmu85f5iRGEoZnIHze55ZnMud2Zn/XO6898R7VFJSUlLAMAzDMAxThFD92i+AYRiGYRimoGGBwzAMwzBMkYMFDsMwDMMwRQ4WOAzDMAzDFDlY4DAMwzAMU+RggcMwDMMwTJGDBQ7DMAzDMEUOFjgMwzAMwxQ51FEMSU5OxuvXr2FgYAAVFZWv/XIYhmEYhskD5E0cFRWFEiVKQFU15xhNsRQ4JG5KlSr1tV8GwzAMwzCfQUBAAGxtbXPcp1gKHIrcSBbI0NDwa78cpSMhIQEnTpxAy5YtoaGh8bVfjtLD68nrqcjw8cnrqUhERkaKAIXkezwniqXAkaSlSNywwPm8E56urq5YOxY4Xw6vZ8HC68nrqcjw8Vkw5KW8hIuMGYZhGIYpcrDAYRiGYRimyMECh2EYhmGYIgcLHIZhGIZhihwscBiGYRiGKXKwwGEYhmEYpsjBAodhGIZhmCIHCxyGYRiGYYocLHAYhmEYhilyyFXgXLhwAe3btxdDsch1cN++fbk+5vz586hWrRq0tbXh4OCANWvWZNln9+7dKFeuHLS0tMTPvXv3yukdMAzDMAyjjMhV4ERHR6NSpUpYsWJFnvZ//vw5PDw80KBBA9y9exfTpk3DqFGjhKCRcPXqVfTs2RP9+vXD/fv3xc8ePXrg+vXrcnwnDMMoC0nJKfCJUMFBzyBcfRYm7jMMU/yQ6yyqNm3aiFteoWhN6dKlsXTpUnHfzc0Nt27dwsKFC9G1a1exjX7XokULTJ06VdynnxT1oe3btm2T0zthGEYZOPYwCLMOPEJwpBrw+IHYZmOkjZnty6G1u83XfnkMwxQiCjVsk6IzNKE6Pa1atcKGDRvEgDIa7Ej7jB07Nss+ElEki7i4OHFLP42UoOekG5M/JGvGa1cw8HoWDMcfvcHI7ffxMeAhEt76Q9elHtT0TBAcEYuhW+9gea9KaFXeqoD+b8UHPj55PRWJ/HzvKJTACQ4OhpVVxhMQ3U9MTERoaChsbGyy3Ye2Z8fcuXMxe/bsLNtPnDghpmIzn8fJkyd56QoQXs/Ph7JQs++oIfTEany4ewTaZSrh/bmNMKrXB0a1ugBIwU977iHBPwmquQ8hZvj4lDv89/55xMTEKKfAkTUCPSUlJct2WfvkNDqd0ljjxo3LEMEpVaqUiBYZGhoW4KsvPgqa/jgpVUhRNYbX82tCtX5jps/B6xBDaFo6AKpqSHgfhJSEOKQkxCIu2BeJ4cFIcakHi3K1UcvelA/ZfMB/7wULr+eXIcnAKJ3Asba2zhKJCQkJgbq6OszMzHLcJ3NUJz3UbUW3zNCXM39Bfz68fgULr2f+oAsbqrvbsWMHDhw4AHVjG9gMWgGtkm7QMC+NGK/L0HGojjc7piP+tRf0K7XGzapjUd85Yxqc4ePza8B/759Hfr6zFcoHp06dOlnCdpRGql69uvRNZbdP3bp1C/W1Mgzz9UhKSsLvv/+OPn36wNvHB1aOFWDceABU1DWhaVFGRHT1XOtDRU0dOvZVoaprhPgQP4z7pjWW/LmMPzqGKQbIVeB8+PAB9+7dEzdJGzj9++XLl9LU0bfffivdf8iQIXjx4oVIJz158gR//fWXKDCeMGGCdJ/Ro0cLQTNv3jw8ffpU/Dx16hTGjBkjz7fCMIwCQNHbgQMHCmEzYOAgGJtZILJEbWh1+gV6LvWypKpJ4BjX74OSP6yDupEVoKKKTc/10eOHMVi8eDHi4+O/2nthGEa+yDVFRS3eTZo0kd6X1MH0798fmzZtQlBQkFTsEPb29jhy5Ijoklq5cqUwCFy2bJm0RZygSM327dvx008/YcaMGShbtqwIUdeqVUueb4VhmK8IdUFSp+SlS5fEOSI5ORmBDu1hOGAdVNQ14GChhzbu1lh19pnYP73zDUkeNS09TF+4Brsv3EXw+2h4/m85/ktJxunLNzB88LfCf4thmKKFXAVO48aNpUXCsiCRk5lGjRrhzp07OT5vt27dxI1hmKINnT/8/f3FhdCUKVOgpqYG18ad8a5EHQQkG8NIXx1jmjujX50y0FBTRYWSRmk+OJ9sIazTfHCIA0ZWUE+JgWmrEfjofQVHD+7HkT078O/uA+jStqXMWj2GYZQThSoyZhiGkRAQECDSURQJvuP5CLU8esJfzRYf3BpBR1UVfWqVwdgWzjDV05Q+hsz8GjuZYcWOY3AoXxk2xnqoaW+Kk4+DhRcOXW6pqKrBoFJLUaMTcWU74l574aerCRg5wgnfdOuM33/7FQYGBvxBMIySwwKHYRiFawOlMS3UGUU+VTEfY9Fy6gYkVugHbQD1Hc0xo105uFjLFiFqqipwMkqBR0Ub0ZxAoxpmH3ycIW1FqGrpwqTJIKQkJyHiyQWEBQVg3eZ/cfrSNYz6YSC+//57ETFiGEY5YYHDMIzC+IPs3LkTnTt3xqFDh0QEp2zzvrBoPBGJxtawN9fDdA83NHOzzNH3KjM3nr9DUERstr+niI5++SYwMTFD0J2TeHL3LEaNewyXOi1gphKNihUrFtA7ZBimMGGBwzDMV4fSUDRYl0axLFq+BhV7TULsiwgklqkIU211jG7mhG/r2EFTPf+Nn/cDwvO039zRffH6XSf8vmgZElQ00eu3fxDy30x07tYD/1u7GqambBDIMMqEQvngMAxT/IqIu3fvjho1akBTSxu6hsaYf9IXnigDXbuK6F2rNM5NaIzvGjjkW9y8CIvB6O138cexp3na30JfC6NauOH+ziUYOGiQ8M2htvJjN57ApqQtJk+dlmPTBMMwigULHIZhvkqdzaRJk/Do0SOULl0aamrq8FW3g+nANdB2a4w6DmY4PKoBfu9cAWb6+etsonTU9meqaLXsMvbfey22aWvkfqqbvu8Bznu/hZWhNhZ2r4TzW5ei+bRNUDO2QXzsR2w4cg2bTt7Bli1bRJs6wzCKDaeoGIYpNEgYkD0E1docP34c565ch1nHqbAa4AR181IobaqL6W3d0LKcVb7qbIjQD3HCB2fr9ReITyRBk4LGLhaY0NIFr97HiC4qWR45dN9YVwOv3sei/1830KFSCVHEXNHWGCd+6YuD3Zpi0qIN+GBQGiPHTUT0o7PYsfcQfp4yHjVr1izgFWIYpqDgCA7DMIUCpXeGDRuGwYMHI+x9BMxKO8HfpgkehybCpKQ9prRxxclxDdGqvHW+xE1ETAIWHH+KhvPP4q/LzxGfmIyyBinY9l0NbBpYE+4ljUT7+Oq+VYUnTnro/pq+VXFpclMMqmcvJo0fuP8azRadw7YbL0EZqQ6VS+LeX9MxqUtd6FraQUVTB6duPxXmouMmT5PDSjEMUxBwBIdhGLlC3VCUjnJ0dES/AQPx9z/b8MLQHXqNO0JVVQU9q5fC+JYusDDIXyoqOi4Rm674Y+35Z4iMTRTbKtoaYUyzsoj0uoHqZUwy7E8ip0U5a9FVFRIVC0sDbeGRQ23lxM/ty6FzlZKYutcTDwMjMXXPA+y580qkyZysDDCquTN61FiFX/YOwt9Lf0XcS0/seqmNNxPnwtUoGRMnjIe2dkYBxTDM14MFDsMwciEmJgYLFizA9evXcfToUWhp6+AoqsH8+7+gqpkqLn5uV05EWPJDbEIS/rn+EqvP+SL0Q+osKWcrfSGSKLWVmJiII96yH0tipk5Zs2yfu4KtEfYNqyeE0+KT3rjp/x4eyy5iSKOyGN7EUUR8Vg6oj++bbceE9YfhHamO7St/QPLHSJy4eB1jfxwg2tzzm15jGKbgYYHDMEyBp6JoKCaJmlmzZsHQyBhl6rZHvHNzhCZooLSVofCzae2ev1RUQlIydt1+hWWnfaS+NmXMdDG2uTPaVyohjcR8KepqqqJrq00FG8zc/xCnnoRg+RlfHLz/Gr91roB6juaoXMoYJ2f3xr67rzDx9TC8un4YV69ew6UTBzFz/jL8NHYo1NX59MowXxP+C2QYpsDw8fHBoEGDEBISgqPnr8K+6iZElaqHFJd6MNZSF1GQwfXtoa2Rd4fg5OQUHPR8jSUnveEfFiO2WRtqY3RzJ3SrZitmUMmDksY6WP9tdRx7GIyZBx6J/3ef/11Hl6olhUCj7q7OVUuh9dZfseJkb8z743dEPDyPjcE22FCxNprXrYZF8+eyfw7DfCVY4DAM88WEhoZi6tSpuHHjhojevI+IRJOpm6HSYjL0VYBuVW0xsZULLA218xUJOvn4DRad8IbXmyixzUxPE8OaOKJPrdL5EkmfC0WYKJJTz8kcC497Ycu1F9hzJxBnn4ZgmoebEFg6mmqY2LYC+tb/C78ffID/jp7Fmye3sdn3IW48eIqBPTth9MgRYmwEwzCFBwschmE+m/j4eNHuXa1aNWzbtg3R0dEo1XwALJ0aQcXQQhT6zmxfXtS25EfYXPINxcIT3lIXYgNtdfzY0AED69lDT6vwT1uG2hqY09EdnaqUxLQ9D/A0OAoTd3li951XIm1V1kIfNkY6WN63JgY2dMIoCz3cP3sQj2+cwMRbl6DvUBVNy9nA2dm50F87wxRXWOAwDJNvSITcvXsXvXv3hpeXFzbuPgrXrmPxKskIqrZuKGWsI9q+21W0yVedzS3/d1hw3AvXn78T93U01DCwnh1+bFgWRrpfPwJStbQJDo6sjw2XnmPpKW9c83uHNksvitTbkMYO0FJXE/tcWDwce+90wMRfl+Dd22D8cvgphnVthoZNm2PXtq0wNzf/2m+FYYo87IPDMEy++PjxI1q3bi18YEzMzKFrbI5J/15FqE1tmNi7Y3wLZ5we30gU/uZV3DwMjMDAjTfQbc1VIW401VSFsLkwqQkmtXZVCHEjgWp+qKvq5NhGaORsgfikZCw55Q2PPy/iul+Y2Ifa37tWL42HOxdg5qzZSA5+KgwFrz70Q2kHJwwfPRaxsdkPAGUY5sthgcMwTJ549+4dJk6cKNJQGhqaVKCC5/ruMBuwCrqONUXx7dkJjTGymVOe62N8Q6Iw7J/baLf8Es56vRWdUL1qlMLZiY1Faiu/3jiFSSlTXWwaWAPLvqkCc30tPHsbjZ7rrmHyLk+Ex6S2r+tqqmNsC2fc3vI7+s7fCa2SbvgYFY6/dh3GmrM++GvjZiQlJX3tt8IwRRJOUTEMkyPkK7N582Zs2LBBTPt+HPAWweW/gWXJjtAwsUHV0sb4uX150TqdVwLexWDpKR/svfsKySlCK4kRCWOaO8PeXE9pPhGKUNHrbuRkIYZ6kvvxjlsBOPXkjRj30LFyahSrhLEOtozvgtvdmmD4H//Dixh1zJo7HxEXt2Lt5q1YNu9XERFjGKbgYIHDMEyOtTadOnXC4cOHUbNeQxiWcMAtOEFHxQil7axEnU2HfKSi3kTGYsUZX2y/+RIJSalToVqUs8L4ls5wtTZU2k+CUmhzu1QQUSwqQvYJ+YAxO+6JIuRfO7mjjFmqaKtWxgRXVk7AvnuBGDfbH5FaevD0f4vatWuje58B2Ll149d+KwxTZOAUFcMwWfDz80OXLl2wdetWdOjWE9r6RnhmWAnGfZfCxKkqxjR3wpnxjdGxcsk8iZv30fGYe+SJmBdFrdYkbuo7mmPvsLrCa0aZxU16atiZiinoE1o6Q1NdFRd9QtFyyQWsPOsrZmRJ6nO6VLXFo52LMHvrKeiWdBFjP08Fa6DXzPWYMGW6SAMyDPNlcASHYRgpUVFRmDt3rvCzOX36NM5dvgbr79bBYvBaqGrri5TL5NauIuWSF6JiE/C/i89F19GHuERpFIMmfOc0MkGZIWEzoqkT2lUsgen7HuCyb5joDDtw7zV+7+KOamVMxX7U7j6jW20MbLYRE9d2weW3GtizaiwSQl/i5KXrGPfjAHzbtw+PfWCYz4QFDsMwSE5ORkREBBYtWiQEjrVtGVhUbQH1yp0Qk6SCKk62Ym4UiZO88DE+CX9f9cfq888QHpMgtpWzMcSEVs5o4mJZLL607cz1sHVwLey9G4hfDz8RZoVdV18VJoWiM0wntTPM1kQX26b0xM3nYRga/j0eHNmMJz5+GPBtP5y4cgd/r1gANTX5mxoyTFGDU1QMU8y5f/++qAHp1asXPHp/DzP78kiq+S10mo9CqbLOWNyjEvYOrZsncUNpGBI2DRecxdyjT4W4cbDQw8reVXFoZH00dbUqFuJGAr1XSkedHtcI3avZim00KLT54vM45Pla1DhJqGFvhhvrp+PvfSdhWbERVHUMcR4V4NSoE7r06iccohmGyTscwWGYYkpgYCAmT56MFy9eCJEDNQ08Wn4K+j3mQUtdVTgH/9iobJ6cgxOTkkWk4s/TPnj1/qN0lhPV6nSuUlIMsCzOmOhpYkH3SkLsTN/7AH6h0Rjx713sdnklHJKp5VxSn9Ozlh3aHVqPZcfHYN2JO3h++TCeIwUPngWgd4eWmDZpPLS0FLd9nmEUBRY4DFMMjfqoxkZXVxf//PMPVFRVYdlsMDScG0JN30S4D1N3FKVO8jII8+jDYCw66QW/t6mFseRdM7KpI3rWKCWcfZlPUN3R0TENsOrsM6w+90x4/1AR8tgWThhUz14qBElUTu1QCX3rO2GMhQpO7N0B37unMOfWebzXssLo7s3hYFe6WEXDGCa/sMBhmGICpUOuXLmCPn364M2bN1h38CLsPX7AR4ty0LJ2RIWSRvi5fTnRCZSX5zrrFYKFx73xOChSbDPW1cDQRmXxbR07MYCSkQ2JPjL/I6dniuaQc/PvR55i793XotU8vZ8QRXZ2zxyAa33bYcisJfC9exX7XutjTYVKcK9YCfu2b0Hp0qV5qRlGBixwGKYYQIKmR48eePToEcwsraGq8xETN52BdoUOsDXQwqRWLuha1VakSHLj6rMwLDzhhdsv3ov7+lrqGFzfHoMb2IuhlEzecLTUx/YfauO/W6/w25EneBIUic6rLqN/HTvhC2SQbi1rlzXHnc2/YtedV5i+7G8kxEbjwRNvlK9eDx3atcGqJQtgZJT3gaYMUxxggcMwRZi3b99i2bJlmDZtGkLehiIiKhoqNZvCrE0jaOvo4rv69hjWxFGIlNy4FxCOhce9xKRvgup0BtS1E3U6pnqahfBuih6UYupRoxSaulnit8NPRB3Tpiv+OPYwGLM6lEdrd2vpviQ+e1QvBY/1k/FLk5pY/78NeH9lB3b8txs2zQfBLfEZ+vfuCXV1Pq0zDMF/CQxTBImPj8eOHTswZ84c+Pr6wicCSGgwDFaNdaFuaAmPCtaY2sZNWtyaE0+DI7HohDdOPn4j7muo0byo0hjR1BFWhtqF8G6KPjTLaknPysIJ+ad9D/EiLAZDtt4WLs+zO5TP4DtEYnTewBYY0bEehi9ojGs+Qdi8Yx9CD8zDr3PnY+eWjahRvepXfT8MowiwwGGYIuhpU6dOHdy5cwdVatWH/kcVnHtnCG1bO1S0MRR1NrUdcjfZex4ajSUnvXFQtDMDlL3qXMVWdEblRRgx+aeBkwWOj2mI5Wd8sPa8nxCVV3xDMb6lC/rXtRPDSCXQZ3Bg7hBc8wvDkFl/4p2OIYI/AjVrVEP95h44d+wA++cwxZri3bvJMEWIp0+fom3btqJDqn7zNtA2MEGAVV2Y9l4IW9fKmNe1Ag6OrJ+ruHkd/hFTdnsKr5YD91PFDUV8ToxtiEU9KrG4kTM0iX1iK1ccGd1AeA9FxydhzqHH6LTyMh4GRmTZnz7P25tmY93BSzB1qgaoqOJeaAq6ztuLH0aOE+7UDFMcKRSBs2rVKtjb20NbWxvVqlXDxYsXs913wIABIi+d+Va+fHnpPps2bZK5T2xsbGG8HYZRKMLDwzF+/Hj88MMPOHLkCHoMGo4DydVgMXgNTNwbY0gTJ5yd0Bg9a5TOEAHIzNuoOMw++AiNF5zD9psBSEpOQRMXC2HQt6pPNThaGhTq+yruOFsZ4L8f6+C3zu4w0FbHg8AIdFhxCb8eeozotLEXEuhzHdysAp4eXIsxq/fDsvG3OLFxIdavWIJKDVphxdr/icgewxQn5J6iojqAMWPGCJFTr149rF27Fm3atMHjx49ltjf++eef+OOPP6T3ExMTUalSJXTv3j3DfoaGhvDy8sqwjQQUwxQXkpKSEBcXh+HDh+Pff/+FnWsFGJVviKS6faCqpok2FUthmoebdJJ1dkTEJGDthWfYeNkfHxOSxLZa9qaY2MoF1fPQMs7IDyos7lOrjKjFmXPwMQ55BuF/l54L76E5HcujmZtVhv2p82rJj+0wJiwaP6oE4czfQXgdFomRQ77H9gPHcWb/dmiyNxFTTJB7BGfx4sUYPHgwvvvuO7i5uWHp0qUoVaoUVq9eLXN/anW0traW3m7duoX3799j4MCBGfajiE36/ejGMMWFq1evomrVqpgyZQpa9h0O/RJlEVOxO4zbTUKFcq749/taWNuveo7ihqIAK874oP78M1h17pkQN5VsjbBlcE3RvsziRnGwNNDGit5VsXFgDeEQHRj+EYM338Kwf27jTWTWyDV97scWj8XZa7dgX70JVDR14WtSA24dh6JZuy4ICAj4Ku+DYYpMBIc6OW7fvi1Owulp2bKlMBzLCxs2bEDz5s1RpkyZDNs/fPggttFVbOXKlfHLL7+gSpUqMp+DrnLpJiEyMtWYLCEhQdyY/CFZM167wl9Pf39/MV6B0lKenp7wfh6AfWoNYNp3Kcz0NTG2mRO6VyspUhbZPV9cQhL+vfkKay744V106j7OlvoY08wRzd0sxMUDRU6VlaJ8fNZ3MMGRkXWw/KwfNl55gSMPgnHBOxQTWjiiV41SWVKQNcsY4+72Rdhy8XusOOeHB//Nhl98DGp7BKNzq0b4bcYU4WhdXNfza8Dr+WXk5zhUSUk/7a2Aef36NUqWLInLly+jbt260u2///47Nm/enCXFlJmgoCAR7aHwO5mUSbh27Zpofa1QoYIQK5TWotoDmqfj5OSU5XlmzZqF2bNnZ9lOz5vbHzfDKMp4Bfp7oCvvJUuWQENbFwZ1voGue1No6hmioXUKWtkmQyeHS5akZOD6WxUcf6WK8PjUL0JzrRS0KZWMquYpokuKUR4Co4Edfmp48SH1gyujn4KeDkkomU3Q7mMi8O+15zh7/DBint1Cckw4anT5AQM96sLGzIjHPjBKQUxMDHr37o2IiAhRqvLVBQ5Fa6htVcJvv/2GLVu2iK6PnJg7dy4WLVoknkdTM3sjMSqeo3B9w4YNhalZXiI4JJxCQ0NzXSBGtoI+efIkWrRoAQ0Ndq6V53rSsX3mzBmR4qXIyrQNRzBz1mxolGsKTQs7NHO1wJTWzrDLIRVFxcKHHgRj2RlfvHyXOgjT2lALI5qURZcqJaBRxAZhFqfjkz7bf28EYNEpH0THJUFdVQWD6pXBiMZlsx2X4R8ajR/nrMTVo7tg0XUG3v4zAbaWZtj1z0aUc3Up1utZGPB6fhn0/W1ubp4ngSPXFBW9CDU1NQQHB2fYHhISAiurjMVxmSHd9ddff6Ffv345ihtCVVUVNWrUgI+Pj8zf0+RdWdN36Y+V/2A/H14/+a4nRSnp+CchnpAMfIhPxq87L0Gv0SA4W+ljRrtywjclp7+h44/eYPFJL3i/+SC2melpCufiPrVKi3bkokxxOD7p3Q1qUBYeFUti1oFHOPYoGOsu+uPoozf4tVMFNHLOenw42RjjzOrpuOz7Iyas3otXYa/hF/EWDTy6omGdGti4cgksLCyK5XoWJryen0d+jkG5XrqRMKG2cFL/6aH76VNWsjh//rw4wVOBcm7QifzevXuwsbH54tfMMF8bilhScb6pqSmeennBPyAQKTV6w2LgSliVccQvHcvjyKgG2Yob+nu44P1W+KaQGy6JG2ozpq6oC5OaiLlRRV3cFDesjbSxpl81rP+2OmyMtBHw7iP6/3UDo7bdFe3/sqjnaI6L8wdj6e6zsKrXDVGvvHF4z04M3XQZG3fsEzWUDKPMyL1NfNy4ceIqtHr16iJNtW7dOrx8+RJDhgwRv586dSoCAwPx999/ZykurlWrFtzd3bM8J9XT1K5dW9TbULiK0lIkcFauXCnvt8MwcoN8nPbs2YNhw4bh3bt3uBupC51W46FnWgraRuZiSvfoZk4w0s3+CuaW/zssOO4lJlQTuppqGFjPDj80KJvj45iiAbWT1ylrhsUnvLHpynNh1HjOKwRTPdzQs3qpLMNU1dVUMapDHfRvUR0TljfD/nM3cPnRC+yeNBYTS5TBnv92oHb1yl/t/TCMQgucnj17IiwsTMzEoSJJEixUECzpiqJtJHjSQ7m13bt3i+JhWVAHCZmaUeqL2sqpe+rChQuoWbOmvN8Ow8itiJi6AZ8/fw7XqrURqxuOU/5x0CpTGU1dLYWfDU2fzg5yuKUJ3+e83or7mmqq6FO7NIY1doSFQdb0LFN0oVlVNI6jc5WSmLrXEw8DIzF1zwPsufMKv3euACerrIaNRjoaWD+pL6YO7IShczcgRM8IH6CNxo3qo3zV2pgyZthXeS8M8yXItchYUaGoDwmjvBQpMbKL5Eikenh4cE7+C6FWb7JR6Nq1K3bsPyqEumHzIdApW1N8EVGdjaw6Cgm+IVFYfNJbtAsT1Cbco7otRjZ1yjCgsTjBx+cnEpOSxXRyOkZi4pPEoNQhjcpieBPHHNOUx+76YcycpfA6uBratuVQ/ptpcA69hJXz58DExKRQPseiCh+fhff9zcM2GeYrQIXD1CV49uxZ3L17Fw/exAL1v4el4zcwMTLA2OZO6FO7TLYdTgHvYrD0lA/23n2F5BQyvgQ6VCqBMc2dYW+es3MxU3ygFNR3DRzQpoINZu5/iFNPQrD8jC8O3n+N3zpXEHU4smhdxQEPdy3F0j1d8dclP/ic2YnbN/bg5OkzmDR2FMYOHQx1df76YBQbPkIZppCv3qjdu2PHjsI+waFiLei71kdKzb7Q0jNCv9plxLRuY13ZnYPkWkuTpnfcDEBCUmrwtWU5K4xr6QxXa45GMrIh92MqQD72MBgzDzyCf1gM+vzvOrpULYnpHm4w09eSKY4mdG+EAa1qoPeUVzjndwsfkjUwadQQrN34N25fOsd1XYxCU7QMMBhGgTl9+jQqVqyIVatXo0HPIdCxdkB0+U4w6zgF7qUtcHB4HczqUF6muHkXHY/fjzxBw/lnsfXaSyFuGjiZY9/welj3bXUWN0yukLCmSM6p8Y3wbZ0yIuq3506gmBr/360A0X0nC6rPGd6qIu7duYNqTdtBVVsf4dY1UHnQHFRv3BrePr68+oxCwhEchpEz3t7eYto3OXCSueW0XxfCrP9yWHy7FGUtDTC1tTOifW7CSUYRcWRsAv538Tn+uvQcH9ImSFcrY4IJLV1EtwzD5BdDbQ3M6eiOTlVKYtqeB3gaHIWJuzyx+84rkbYqayG7mN3J2hCX/voVB4d8i6XnA3Bhbj/4vw9CnbY94dG0AVYv/BX6+tkXwjNMYcMCh2HkBBXBUQchjQQ5dOgQDCxKwrjhtzCo4gEjXS2Mbu4srqSRnIQjmS6CP8YnYfNVf6w5/wzhMamzV8rZGAovm8YuqfOiGOZLqFraBAdH1seGS8+x9JQ3rvm9Q5ulF0UB8pDGDtDKZup4+5rOaFPNEfNtNmLBvD8Q+doXW9f+iXsvQvHPqgVwL2MlzFcZ5mvDRyHDFDA0AHb//v3Cp6lX795Qq9wBBhWbw6DjDJjU7YH+jcvj3MRUw73MRcTxicn4+6o/Gi44iz+OPhXipqyFHlb2ropDI+ujiaslixumwKDjj7qqTo5tJLr14pOSseSUNzz+vIjrfmHZPo7qc6b1bQW/a8fQbchkaFo7ItyxDeq27ooybpVx6fot/pSYrw5HcBimAKGOKHLfVtfQQOSHaHg+e43XJx/CtM0Y1Hc0F23fLtZZfUioXnjXnUCsOOuHwPDUeVG2JjqiK6pT5RLiC4Vh5EUpU11sGlgDhzyDMPvgYzx7G42e664Jc0CaVJ4dJnpa2PbbCPiM+hbT/7mAXf73EJ0Qhza9vkMlVwf8u24ZSpey5Q+O+SqwwGGYAuDFixfCw4Ycu8nbRkVTByZtxkHHoTrsLQ3xU9tyaOaWNfqSnJwiPGz+uKeGkGuPxDZLAy2MbOqInjVKQ1OdhQ1TONCx2b5SCTR0ssC840/x7/WX2HErACefBKOtjQra5GCZ5mRliJ3j2mF3jeuYMn8VfI+sx2X/++g2vwm+a+iIfu0aQ0enePoyMV8PFjgM8wVER0dj7969+P7775GYmIhuv22HSbuJ0C5dAcYmZhjVzAn969plESrUsXLWKwQLj3vjcVAkfb3AWEcDQxuXFSMZspsEzTDyhkZ6kONxF3JC3vMAPiEfsMVXDX6b7+D3LhVQJofJ9V0bVECHuivxy8ZWWPfPHgSlmGBI364Yb2iMjZv+RlePZpxiZQoNvjxkmC+I2ri6umLo0KGwKO0EjZLlcMHnLQzc6qNv4wo4O7Exvm/okEXcXHkWiq6rr2DQpltC3OhpqaG1bRLOjGuAHxuVZXHDKATV7UxxeFQDjGvuCHWVFFx+FoaWSy5g5VlfUSuWU13PnO864cmh/6FlaRWoaukhLkUdvXr2QNlKteDl96JQ3wdTfGGBwzD55NatW+jRowdMTc2gY2yBeHV9xFXtBYuev6Fhjco4NLIB5napAPNM5mn3AsLR93/X0Xv9ddx5GQ4tdVX82NABZ8c1QJtSKWLiN8MoEiTOhzZywJRKSajrYIq4xGQxzLXd8otisGtOmOhpYsOkPrh97wGqteuH5LhoBAS8wjebPNH9h/EIfhNSaO+DKZ7wGZVh8ggNd6XJ9Vu3bkVAQAA8YwwR02AkrHSMYGdlIgZitipvlSUE/yQoEotOeOPUkzfiPs0D6lWjNEY0dYSVobZwN2YYRcZCB9jUpRoOPwrBL4eewPvNB3RbcxW9a5XG5NauwgwwOyo5WOPKxl+xvV8XLNx/Hf43TuDu8RU4vPtfTJ4xG1OHDYCmpmznbob5EljgMEwexyvUrl1bpKVsy9eEXnlHfCjTAMbm1kKoDKxnl8U35HloNJac9MZBz9eg+kxVFaBLVVuMbuYkulYYRpmgv4HOVWzR2NkSc48+wc5br0Qh8olHbzCzfTm0q2iTbX0Nbf+mWTV0a1wFP6/djaX3jyJFSx+zxv6IZYsX4ubNW3Cw4lEjTMHCKSqGyYGDBw/Czc0NR46fRJW2/aBdwhmJVXvCov149GlSBWcmNBI+IunFDbV5T9ntKSzwD9xPFTdtK9jgxNhGWNi9EosbRqmh1NP8bpWw/YfacLDQQ+iHOIzcdhcDN90UQ2Bzgupz5g7rjpdeD9CmW1+o6hojsWQlNJn+N9xqNcEdz4eF9j6Yog9HcBhGBg8fPhTjFeLj4/Hs2TP0Gz0dxl1mwbJvTdRyMMfP7crBvaRRhse8jYoTBZh0VUuGaURTV0uMa+GcZV+GUXZqO5jh6OgGWH3uGVadfYZzXm/RYsl5jG3ujEEyTCzTY2Gogz0LJ+DWD72w8KQP9i8Yh9gX99CwfU80aVgfm/78A2amJoX6fpiiBwschknHu3fvhBPxH3/8gRMnTsC4jBuM6vWGYc3OKGWmL+ps2rhbZwjFR8QkYO2FZ9h42R8fE5LEtlr2ppjU2gXVypjy+jJFFopckhllu4olMH3vA1x//g5zjz7FvnuvRaF95VLGOT6+urMttjmVxN8OKzFh0hREvgvBoa3rUN7zKXZu/xf1nK2hpsaWCcznwSkqhgGEhw3NjKLxCsNGjoF2nb7QdWsI3dbjUaJpP0xuXxmnxjWCR4VPdQY0/HL5aR/Un38Gq849E+Kmkq0RtgyuKcL3LG6Y4oKjpb445ud3qwhjXQ1RWN951WXM3P8QUbE5F9HT31P/NnURePccho6bAi2LMlCt2g3tB46BlUM57Dt+ttDeB1O0YIHDFHsuXbqEypUrY9369SKCc/DcNZzwi4ZFh0n4pll1nJvQWAwg1NZIvZKMTUjC/y76odH8s1h00htRsYlwsTLAun7VsG94PTRw4mGYTPGDhEqP6qVwelwjYRJItWebr74QtWjHHgYJc8vcWtKXjh+AV8+eol/rWoh+eAphL73Re8g4VGzYBo+8Mk2kZZhc4BQVU2zx9fUVN0pJPXr0CBqGZrDoNhM69lVRw94cP7cvh4q2n0LsCUnJ+O/WKyw77YPgyFixzc5MF2NbpIbo1ahNimGKOWb6Wljcs7LoGJy+7wFehMVgyNY7aO5mhTkdy6OEcc4jG8wNtDGvZ010cr+DYTMXwfPIFjzwv4e2U1ZhZC8P/OBRGwYG+oX2fhjlhQUOU+yIjIzE9u3bMXLkSOjq6aP+tH9g2nKYSEmVtrbAlDauGVpek5JTcPD+azFlmU7WhI2Rtmj37lrNNsdiSoYprtR3MsfxMQ2x4owv1px/JnygyMV7fEsXDKhrl+sFQZ3y9rjz33Ks29cVcxYsRbJDPUwZ0g8/pSRh3rLVGNmvC499YHKEz8xMseL27dtwdnbGnF9+gY6JFWKNyuCuXzAsa7bHxPZVcXp8IzFwkMQNhdSPPQxGmz8vYMyOe0LcmOtrig6qsxMao1fN0ixuGCYHKK07oZULjoxugGplTBATn4RfDj1Gp5WX8TAwIte1o7/DHzs3gd+FPfiusiFIEiUkJmL82DGwcamM8zfu8/oz2cIChykWXL58WcyMKuPgiLgkICQmGTotRsOyxxz0aFRZCBYajEknZBI2F7zfouPKyxiy9bZwbTXUVsfEVi44P7GJaIGV1OMwDJM7zlYG+O/HOmKIJ40keRAYgQ4rLgmxEx2XmKdurZ/7tYCfz1O0/X4ykqJCEfLCBwP/fYTOI37GM/8A/hiYLHCKiinS0EiFv/76CwsWLBCTv8+Em0Gn40wYmtigqr2FiMZUKf3Jb+Om/zsxa+fG89Q5O7qaahhUz14MzczJjp5hmJxRVVURox2al7MU4x4o7bvh0nMcfRCEOR3d0bycVa5LaGtuhP1LJuNcv06YufEwvEJeYN+2X3Bww2KMmrkAv44aAF3dnGt8mOIDCxymSEKFw1RrU758eURFRcHStQZUNE0QY+6KMiVsRJ1Nh7RUFEHh8oUnvIRZmaSjo2+tMhjWpGyWoZkMw3w+lgbaWP5NFXSpWhIz9j3Eq/cf8d3ft4S/1KwO5cV8ttxoXNUF56o4Y+WuU5h2yQ3xKapYMm0E1i76FQdPXUTTSg78ETGcomKKFpRe2rlzJ1xcXPDkeQAc67WFlm05qNX5FiU7jMG4DjVFnU3HyiWFuPF5E4WhW2+j3fJLQtxQ4eM3NUuJ1nDqomJxwzDyoYmLJU6MbYgfGzmIv7ujD4PRbNF5/H3VXxT25wb9/Y7o3gJvfO9j8LBRUNc3QYppGfT/31U4VG+CU5dv8kdXzOEIDlNkuHPnDn766Se8CXkrxiu0HTQBhs2HwapiH3SqUlJMPZa0qL4Mi8HS097YdzcQdC6lQE7HSiWEK6udud7XfisMUyzQ1VTH1DZu6FipJKbufYD7AeH4ef8j7LkTKJyQ3WxyH8Cpo6mBlVN/wIQBXbHkqCc2rViAqLvn0LZzd9Rr0hKbl/6CUja5p7+YogcLHEbpCQkJgYGBAYYNG4br16/DyKESjOr3gUHNzqhsR3425UUHB/EmMhbLz/hg+40AJKZdJbYsZyVaV12sDb7yO2GY4km5EobYM7Qu/rn+AvOPeeFeQLiIqn7XwB5jmjlDRzP3on57GzMsG9QEbcuZ4rvhsQgNC8PZnevheukcNu8/hfYVrKGlpVko74dRDFjgMEoLDcL8+++/xVDM/j+OhH6jQdCN1IZ+4wEoaWsrIjadKpcUxY3vouOx+pwv/r76AnGJqYMwGziZC2GT27wchmHkD6Wpvq1jh5blrDH74CORslp73g+HPYPwayd3NHaxzNPztKpdCS9vncaCTXsxe/pE6NXsjiE/L8LAuwcwf9ESDO3dKeMDol8CcaH5f8Fa5oBe6fw/jik0WOAwSsnx48eFUZ+VTUlRTLxu+35YfTMXJTtPwo+NymJIo7LQ01JHZGwC/nfxOTZc9EN0fOogzOplTIQ3B01DZhhGsbA20sbqvtVw6vEb/Lw/tQh5wMabwp9qRjs3UaScl/qcSQO7YPg37bDp6guM79kScW/9MW7aTCxYtAR/r1uB+tUqpIqbgy5Acqozeb5Q1Qbae7HIUWBY4DBKxZMnT/Dhwwd4PngIHx8f+IcnwKLzNOg41RYnNSNdTZQvYQhVFRXhnkq38JjUYX+0nYRNY2eeFcUwig61jdcpa4bFJ72x8fJz0VZ+3isEU9q4oVeNUiIymxt62poY3sQJLW9fxcCJv+D6yf14/sITHUbOxvDhIzC8ZjKsP0fcEPQ4ivxwFEdhYYHDKCTURUFeNCFRseKKzclYBVv+3oyJEyfC1t4RZr0XwrhRfxhUaQtVLV3p40Kj4sTcGzLmi4xNNRAra6EnUlGty1vn6aTIMIxiQFHYGe3KiVTz1L2eeBgYiWl7H2DPnVf4vUsFYSCYF5xKWePS9pU4cH4ARk76CSm1emDBtNFY+O45tn0PtK9KKTK5vx2mkCmUj3TVqlWwt7eHtrY2qlWrhosXL2a777lz58SVeObb06dPM+y3e/dulCtXDlpaWuLn3r17C+GdMIUBTR6uP+8Mvll/DaO330OnaathU9oei1eth5q2HkJUTBAS/gFGtbtnEDeEpLmUxE1JY20s7F4JJ8Y2gkcFGxY3DKOkVLA1wr5h9YTYIfPNWy/eo+2yi1h43AuxCamp57zQoVEN+F87hrkd3KAS8w7xHz9i/iGg5gzgznO5vgWmKAqcHTt2YMyYMZg+fTru3r2LBg0aoE2bNnj58mWOj/Py8kJQUJD05uTkJP3d1atX0bNnT/Tr1w/3798XP3v06CE6aBjlFzdDt95BUEQsPr64j/DL26BhWhJJ8bF48z4aZj1+Q4luM6Cmk/uV2x9dKqJbNVue8s0wRQB1NVUMrm+Pk+MaobmbJRKSUrDirC9aL72Ay755LxKmC+Y+TSrgjb83xv40AY8DgbsvUq0i/rsOvPyMemOmmAqcxYsXY/Dgwfjuu+/g5uaGpUuXolSpUli9enWOj7O0tIS1tbX0pqb2qU2QnqNFixaYOnUqXF1dxc9mzZqJ7Yxyp6VmH3yM+PBghF/6ByHbpyPi0j9IigkXBcQ2g5bDsGRZTG3jmqfnexcTL/fXzDBM4VLSWAfrv62ONX2rwspQC/5hMejzv+sYt+Mewj7E5fl5DPV0sHBUd3gvBDb/CBjrAv1WAy4TgC0XgZi8PxVTHGtwqI2XpjdPmTIlw/aWLVviypUrOT62SpUqiI2NFeknMm9r0qRJhgjO2LFjM+zfqlWrbAVOXFycuEmgrhsiISFB3Jj8IVmzgl6768/f4aWfL15vHCHua5Vyh6aFHdSNraGmk2r4FZuQjKTkvIWkzXTVleLzldd6Fld4PYvHejZzMUfNkXWx5JQvtt4IwJ67gTjzNASTWzuja5VPY1hyJDERlkZAvwaAdxBQ2xGI/AiM3gJM/w84Mw1wtM7+4TTZHPlcF0VdT2UhP+smV4ETGhoqZgJZWWV0kaT7wcHBMh9jY2ODdevWiVodEiVbtmwR0RmqzWnYsKHYhx6bn+ecO3cuZs+enWX7iRMnoKubsYaDyTsnT54s0OW6HaoCddOS0C5VgWYuwLTlUJGeysxLnycw1lRFuAjQyDqJpcBYE3j7+BqOPEGxXc/iDq9n8VjP6qqAeXlgh58aXsckYOreR9hw+iF6OCTBKpe5m0ZJz9A47d/ONsDZ6cDZx0CrecD7aGDndWBax+wff/nSJUSoBRWp9VR0YmJiFKuLKrOSpnlB2alrmiFENwl16tQRE6EXLlwoFTj5fU5KYY0bNy5DBIfSZBRJMjTM3Qqcyaqg6Y+T0oQaGgU3Yfvjudt4M/t7qKhpwLLnL9l+nq0b1kLtmASM3H5f3E8/tSb1ESr4tUsltCpvVazXs7jC61k81/P7pGThebPszDP4RiZjwQMNDG3ogB8a2kNLPZtqjPd3gVOf7tIpp2l54MYcoM184KZfzv/PevXrAyZViuR6KiqSDMxXFzjm5uaidiZzZIWs9TNHYHKidu3a2Lp1q/Q+1eTk5zmp04pumaGDiw+wz6eg18/JTAtxAQ+hpm8qU9yopJmA1XG0FIXD6upqmHXwMYIjPvlY0O9nti+H1u42UDb4eOT1VGQU/fiklzasiTPaV7LFT/se4rz3Wyw7+wyHHgbj984VZBt7qsv+CkxKBt5EALdz6azSoMd/5poo+noqKvlZM7kWGWtqaopUU+ZQHN2vW7dunp+Huq8odZU+qpP5OSndlJ/nZBQPQwN9lKtcHZrWTlkST5L7JF5I3BAkYk6PayTdZ0P/6rg0ualSihuGYQqGUqa62DSwBpZ/UwXm+lrwexuNXuuuYdKu+wjPY+MBFRx3rAa0rMCfijIj9xQVpYaojbt69epCmFB9DbWIDxkyRJo+CgwMFDOFCCoUtrOzQ/ny5UWRMkVuyPOGbhJGjx4t0lXz5s1Dx44dsX//fpw6dQqXLl2S99th5Eh0dDQe37sFUwsrEYmhVvHcIjOSuVIEzaqRiB+GYYovFAGm0Q4NnSww7/hT/Hv9JXbeeoXTT0LwUzs3YRyYUxFyCRNgVheK0BTqy2YKGLl/fORXExYWhjlz5gg/G3d3dxw5cgRlypQRv6dt6T1xSNRMmDBBiB4dHR0hdA4fPgwPDw/pPhSp2b59u+iumjFjBsqWLSv8dmrVqiXvt8PIERMTE/Tu3RvGxsZYNrkpGi84i4D3HzHNw034X8gSL9FxqW7F2hqqLG4YhsmAka6GSE91qVISU/c8gE/IB4zdcR+7bweKAZ522eQwyBunxgyglBnwchkvqrJSKPp02LBh4iaLTZs2Zbg/adIkccuNbt26iRtTdDAzMxPHCeVYSczoaKZ6H7mXMMxWvHxIEzj6WnypxTCMbKrbmeLwqAZYf9EPf572wSXfULRaegFzGsajp4z9NdRSozhWRryiygxP32AUBl9fX9SvXx+dOnWSGv8ROc2PkkRwaGYNwzBMdmiqq2J4E0ecGNMQ9RzNRHr7z4vvEJuctWiVslckcuiW4zRxLXNecAWGvxUYhYE63WgkB7lYE2n6RkwGzw5JBEdPkw9lhmFyx85cD1sH1xIDOyfs8kRTr7UwUc/Yehwe+BIvQhfjbaIJklqelB1BJnHDk8QVGv5WYBSGxMREvH//XnTfEckpqQonpym/0XGprsacomIYJq9QgXEJY13yE8XrBEtxS0+yvi2sv7WFipo6bkTYoU5ZGS3mjMLDAodRGMj1mtyvJT4HkhRVTt0On1JUOcWSGYZhMhISFZv9uSj6PSKv74aqjiFCorrw0ikpLHAYhYEGpz569AjqaeZbaQEcqOUkcOK5BodhmPxjaaCd7e+SYz8gxusy1AwsctyPUWxY4DAKA1kDkEWAqamp8D+SRHBy8raRRnC4BodhmHxQ094UNkbawgk9/bgXQlVDG1ol3aBrZCb2Y5QT7qJiFIaoqCgcPXoUZ86cyVCDk9NQ4A9pNTjcRcUwTH6gCycyDxXnmEy/S0mIRVzgE6i/82N/LSWGIziMwkCT3cntmmaYZSwyzj2Co881OAzD5BGKDt94/k60io9p7iScjt9ExUl/b2luAus6DeFmb8trqsSwwGEUhpiYGFy9elU6dywvbeLsg8MwTH449jAIsw8+zjAKxsrw0zBmmmNVo5QBAl/VkdYDMsoJf3qMwmBkZCRmi9HIhgxGf3nxwWGjP4Zh8iBuhm69k6XmJiTyU/SmvqM57t29gxo1aqBUqVIZRgkxygULHEZhsLa2xqxZs6RXTZIUVU7zMyVdVOyDwzBMTtAFE0VuMosbIv02sqVQU1MTF1yGhoa8qEoMCxxGYfDy8kKVKlVEiur169dIzlMXFRcZMwyTO1Rzkz4tldN+BhoaKFGiBKysrHhplRgWOIzCoJHppJKvGpy0wZwMwzD5NfbLvJ9GYiyePHmCDx8+8GIqMSxwGIWBQsMkciQpqiRJioqHbTIM84Xk1bCP9nOxdMHZs2fFfDxGeWGBwygM8fHxePHihfhJpOShBoeLjBmG+VJjPwnkmk77vXzhj02bNgnTUbKuYJQTNvpjFAZnZ2fcuHEDR44cEfelTsbZpKhIAEXH87BNhmG+zNhPgqWBltgvLCwMmzdvxq5du3hplRgWOIzC8ObNG8yfPx+rVq3KUIOT3bBNMumSiCAetskwTG60drfB6r5VYZnO94Yw1dNM/amf+lNPTw9169ZF9erVeVGVGE5RMQpDRESEuGKiLipJB1VOXVSS9BShx7OoGIbJo8hxsTJEk0XnoKGmgr8H1UJ4TDyG/nMHWuqp1/zR0dG4cuWK8MFhlBcWOIzCILlqolENEg+cnGpwYtJaxHU11XIsRGYYhklPVFyC+Gmur4U6Zc1w4P5rcV9LXS3DuYjbxJUbFjiMwiC5aqIIjqSDilDNJYKjy9EbhmHyQXhMqsAx1k1NScUlpF4saaZFcMqWLYudO3cKwz9GeWGBwygMBgYGaNGihehcSKdvsvXB+eRizCchhmHyTvjHNIGjoyF+xicli5+SFJWnpyePaigCsMBhFIaSJUtixYoV4qpJUjycUxcVt4gzDPM5RMSkWlEY66YKnLiENIGjoZbBk4tujPLCAodRGJ4+fSod1eDl90K6PTsjY54kzjDMl6Wo0gROYsYIjo6ODipVqsQ1OEoOCxxGYUg/4C45Od32bGpwJAKHB20yDPM5KSojndQanPg0gSOpwYmJicGtW7e4i0rJYYHDKAw0ooFmUVlYWGQsMs42RcWDNhmGKYgITlKGCI6joyMOHz4Mbe28jXdgFBMWOIzCEBcXJwbchYeH57FNnIuMGYbJPxEf4zMUGX9KUalJPbkOHToEExMTNG3alJdYSWGBwygMdNVEA+40NTWlRn8UvMnOyfhDWhcVt4kzDJMf3mcTwZGkqN6+fYvVq1eLFNVvv/3Gi6uksMBhFIZ3796JAXfGxsaYUr5qjh1UBBcZMwzzOZBzsawaHEmKioz+aMgmG/0pNzyLilEogUMD7shgS1KDk5NDcXRaDQ774DAMkx8iPubcRUWmo1evXsXt27d5YZUYjuAwijmqIS1FldMEBvbBYRgmv6SkpGQtMk6Q3SZubW3NC6zEsMBhFHJUg6TIOC8pKm4TZxgmz+eZ+CQkpl1AGUtSVEkZi4xdXFxw7ty5bOv/GOWgUFJUq1atgr29vWi5q1atGi5evJjtvnv27BF2/dQqTH4olAc9fvx4hn2oToMOvMy32NjYQng3jLwjODVq1IDEyDi7FvEMNTg8i4phmHzW31BBsbaGasY28bT79+7dEx1UFSpU4HVVYuQucHbs2IExY8Zg+vTpuHv3Lho0aIA2bdrg5cuXMve/cOGCEDhHjhwR+c8mTZqgffv24rHpIfETFBSU4caeBcpNmTJlRP3NmjVrpKMacqzBiWcfHIZh8oc0PaWjIY3QSFJUmmpcllqUkHuKavHixRg8eDC+++47cX/p0qUiIkMteHPnzs2yP/0+Pb///jv279+PgwcPCht/CXRgcn60aPH48WPpqIbzd71yrcH51EXFwzYZhvm8AuMMKaq0CI4kmsxdVMqNXAVOfHy8iMJMmTIlw/aWLVuKWou8kJycjKioKDFhOj0fPnwQV/xJSUmoXLkyfvnllwwCKLOBHN0kREZGip8JCQnixuQPyZoV9NolJiZKB9zFxSdIU1TZ/X8kRcZaqgX/WorCehZXeD15PXMiLOqj+GmorS49VmLTosFqSBHbyOiPvqPIB6eg/y75+Pwy8vN5yFXghIaGCgGSWQXT/eDg4Dw9x6JFi0TxaY8ePaTbXF1dRR0O5UdJrPz555+oV68e7t+/DycnpyzPQZGi2bNnZ9l+4sQJ6OrqftZ7Y4CTJ08W6DIEBgYK0UrpxwuiTksdCfFxIl2ZGapB/hBLkRsVXL98Hk9TawWVmoJez+IOryevpywuv6GwsBriIt9Jzy3vIlPPJbdvXke4F8RF9cSJE4XpqKzzDx+fXw+aE6ZQXVSZK9GpTS8v1enbtm3DrFmzRIrK0tJSur127driJoHETdWqVbF8+XIsW7Ysy/NMnToV48aNk94nUUTKnCJJ9GXK5F9B05cH1UpRtKWgoMI+X19fkaKqW68+4HkNOjra8PBolGXfj/FJSLl2Wvy7fZuWSt1JJa/1LK7wevJ65sTL836Any+c7W3h4eEutv3x+AIQG4vGDeqhQkkjBAQEiLpPEjgeHh58fCoQkgxMXpDrtwL5mdCE6MzRmpCQkFxzm1ScTLU7//33H5o3b57jvqqqqqLzxsfHR+bvtbS0xC0zknQI83kU9PpRayYNuKPPSlVNTdomLuv/ER77ady4ka52jsXIygIfj7yeikxROT6j4lPPHaZ6WtL3I3Ey1tNO3Uamo1Q/ShfCP//8s1xeR1FZz8ImP2sm15JxUr/UFp45VEz3qYArp8jNgAED8O+//6Jt27a5/n8oIkRX/3TlzygvkgF3p06dyrWL6lOLuFqREDcMwxRum7ix7qe8tsTJWDKLSlJkXL16df5YlBi5x/UpNdSvXz9xoJCnzbp160SL+JAhQ6TpI6q9+Pvvv6Xi5ttvvxV1NZSGkkR/yFnSyMhI/Jvqaeh3VG9D4SpKS5HAWblypbzfDiNHwsLCRHcdCdVuQybl6IMTnTZoU0+JU1MMw3y9NnGjtEnismZRSUxHKYLDKC9y/3bo2bOn+OKaM2eO8Kpxd3cXRVtUTErQtvSeOGvXrhXdNMOHDxc3Cf379xeFxUR4eDh++OEHIX5I9FD3FPnn1KxZU95vh5EjVPBNIliMapA4GWcbwWEPHIZh8k94pjZxGgvzyck4VeBQmpwuoDkroNwUyuXvsGHDxE0WEtEigeyxc2PJkiXixhS96ngacCdGNaSlqLKrRWcPHIZhPocIqdFfxjENhJZGau1f+fLl8fDhQ15gJYdtGxmFQTLgjk4uSbnMopIO2uQxDQzD5IPwj/EyJ4mndzK+c+eOiOI4Ojry2ioxXMDAKAx0MpEMuHv49pPRnyx40CbDMAVRgyOZQ0WnGg01blgoSnAEh1EYHjx4IAbcubm55dpFJY3gcJExwzB5JDYhSRqxkUZwEj7V30j82aiLiuoBqQuYUV44gsMoJJIUVXYd4DE8aJNhmM+M3qirqkjNQaUt4ukGbVIXFdUDcheVcsMCh1EYJN4T1EVF3kY5d1GlRnD0edAmwzCfUX8jidZIW8TTCowJOzs7bNy4UZyTGOWFBQ6jMEi8J6iLStLYkN1ID0mKSpeLjBmG+QIPHEkNjqRFXDIo2t/fX+q9xignLHAYhYGumnbu3Altbe1PPji5tIkr8wwqhmG+jsDJycWYeP36tTCUpRTV2LFj+WNSUvjbgVGoCM6NGzfEANTqdtVzTFF9YKM/hmHySYQkRSXTxVgtS7o8t5mJjGLDAodRGN6+fYuFCxeKFNX6jt/lmKJioz+GYfLLe0mKKq2DKn0EJ32Kikc1FA1Y4DAKWWScm9GfZBYVp6gYhsl3iirNxTh9DU76FBVNrC5RogRHcJQcFjiMQhYZD5K0iWfj1PQpgsOHMMMw+UxR6WY/aJMgR3UaAs0oN2z0xygMkgF3Dg4On4z+sk1RpV51cQSHYZj8FxnLSlF9qsG5ffu2iOLQuYhRXvjyl1EYXFxcpAPuDjwIydOoBl3NTyclhmGYfLeJJ2RtEycfrsTERHFjlBeO4DAKg6enp4jiULu4ZJq4rC4qOvlwDQ7DMPkl/GPWNnHJNPH0AkdSD1i9emo3J6OccASHUUgkPjiyusQ/JiQhTf9wDQ7DMHkmIiZrm7h0FpUGd1EVNVjgMAqDrq6uGHCXvotKVopK4mJMv+IUFcMw+Y/g5FyDU7p0aaxcuRL6+vq8uEoMCxxGYYiJiRED7qiLqntahEaWwJEUGOtpqmfrk8MwDJMeageXDOlN3yYuSVGlbxNPTk7Ghw8f+Pyi5LDAYRQGumqiAXc6OjqIzaEGh03+GIbJLxFp0Ru6JjLQVs+xyPjVq1eYPHmyGNUwdOhQXmwlhQUOozBIBtwZGBjA1Lam2CYrQMMeOAzD5JeIdB1UqukunGQ5GVORMaXLeVSDcsMCh1EYgoODxYA7SlHN2dE5+whOmouxHk8SZxgmv/U36QqMsxu2SaajlC6nCA6jvLDAYRRyVENajXE2RcZpNTha7IHDMEw+PXDStYhnN2xTTU0NRkZGYvAvo7ywwGEUclRDmxy6qCQpKnYxZhgmr4TLaBFPP4sqfYqqSpUqCA8P58VVctjoj1EYJAPurK2tc/TB4RochmE+t8g4fYt4dimqu3fvwtjYGO7u7rzQSgwLHEZhKF++vBhwd+fOnRydjCU+ODxok2GY/E8Sly1w0qeokpKSEBERgcjISF5gJYYFDqMw3Lt3T0RxqF1c4lQsy+dG4mXBKSqGYfJKeNok8cw1ONl1UfGoBuWHa3AYhUIy4E4yTVxNhgSXRnC4i4phmC+N4KT54GTuoqJ6QO6iUm5Y4DAK2UX1qQaHjf4YhpFfDY6sYZu2traYP3++8ORilBcWOIxCdlHVz5PA4cOXYZh8RnAyFxlLh21+qsFRVVUVs/HIVZ1RXvgbglEY6KqJBtzRiSU49ZyT47BNFjgMw+S7BifdHKoMXVTp8uEvX77EiBEjRIqqf//+vMhKCgscRmGgzgUacJeSkiJu2dXgSIZt6rPRH8MwXxjBiZf44GhkLTLmUQ3KDQscRmEICgoSA+4oRTVsXROxLf3MmCwpKi4yZhgmDyQmJSMqNvW8YZKHLiouMi4aFEqb+KpVq2Bvbw9tbW1Uq1YNFy9ezHH/8+fPi/1ofwcHB6xZsybLPrt370a5cuWgpaUlfu7du1eO74ApDCg1RQPuqlevLm0Tl1mDI5lFxTU4DMPko8CYMEw3SZwixbKM/siegiwr6MYoL3IXODt27MCYMWMwffp04Q7ZoEEDtGnTRuQ4ZfH8+XN4eHiI/Wj/adOmYdSoUULQSKAhaD179kS/fv1w//598bNHjx64fv26vN8OI0diYmLEZ3vr1i1pF5WazCJj9sFhGCb/gzYNtNWhni7vnZCUdiWVyeiPLrDj4+Px7NkzXmYlRu4CZ/HixRg8eDC+++47uLm5YenSpaJwa/Xq1TL3p2gNGb3RfrQ/PW7QoEFYuHChdB/6XYsWLTB16lS4urqKn82aNRPbGeUl/YC77EY10BWXJIKjyzU4DMN8SQdVWv1N5hQVXTjT2BiKJjPKi1xrcEgB3759G1OmTMmwvWXLlqIdWBZ0BU+/T0+rVq2wYcMGJCQkiJAh7TN27Ngs+2QncOLi4sRNgsR+m56Pbkz+kKxZQa8didW3b9+Kf/984HHqxpSUDP8fqr+RTBrXUs34O2VFXutZXOH15PXMTFjUR/HTSFsj4/nk46fvBdWUJCSktYx//PhR1ASqq6sX+N8lH59fRn4+D7kKnNDQUNEZk7kSne4HBwfLfAxtl7U/udvS81EBanb7ZPecc+fOxezZs7NsP3HihKj7YD6PkydPFujS+fv7i5SkiYkJGo6jCJ8qfH29cSTWS7pPhOj0VIcKUnD25AnIyGApLQW9nsUdXk9eTwk339KJQg0J0eE4cuSIdPt7oW/UoaaSgqNHj0q3BwQEwMXFRQzcTL8/H5+KUcqgUF1UmecJUZpB1oyhnPbPvD0/z0kprHHjxmWI4FCajCJFlA5h8q+g6cuD0oQFWYRHs6jo4KU0VclSpYCQQLi6uMCjkYN0n+eh0cDty9DT0kDbtq2KxEcnr/UsrvB68npmJuTqC8DXC46lSsDDo6J0u39YNHDnMnQ0NeDh8el8QpkHLy8v8T1BNaF8fCoO+RmAKleBQ5b7VFeRObISEhKSrb+AtbW1zP0pVGhmZpbjPtk9J3Va0S0zXCX/ZRT0+tHVkmRUA5AqVtXV1TL8P+KSVKSDNouaGODjkddTkVHm4zMqLjX1ZKKvmeE9JKWVoWprqGbYTnWgM2fOFBdb8nrPyryeX5P8rJlci4w1NTVFNXrmUDHdpy8yWVCbcOb9KZVExV6SN5bdPtk9J6McSLwnbt68iaRsuqg+tYh/6nhgGIbJiYiYVBdj48wuxglZXYzFfU1N2NnZCXd1RnmRe4qKUkPUxk0ChYTJunXrRIv4kCFDpOmjwMBA/P333+I+bV+xYoV43Pfffy8KiqnAeNu2bdLnHD16NBo2bIh58+ahY8eO2L9/P06dOoVLly7J++0wcoS6FmjAHbmI+mTjgyMx+aMIDsMwTH7axLMdtJluDpWkHnDgwIEiRdW9e3deZCVF7t8S5FcTFhaGOXPmiKp0d3d3UbRVpkwZ8Xvalt4ThwwB6ffUJUVziehLb9myZejatat0H4rUbN++HT/99BNmzJiBsmXLCr+dWrVqyfvtMHKEaqgkA+6S0pz+MjsZS+ZQ6bKLMcMw+WwTN9LJZtBmuhZxgkc1FA0K5TJ42LBh4iaLTZs2ZdnWqFEj3LlzJ8fn7Natm7gxRQeK5NGAO+qU67LwkEwfHInJH7sYMwyT/whO5jENSVlcjNOnyymCwygvHOdnFAbJVRMVGUudjFWzS1FxDQ7DMPmswckyaFN2BIcpGrDAYRQGyVUTRXBs0oJzmVv/JSkqjuAwDJPvCE7mFJVU4GS8YKpataowqmWUG5atjMKQfsBdtl1UXGTMMEw+SE5OkQ7bNMpmVEPmFNWjR4/EEGcaAcQoLyxwGIWhUqVK4qrpxYsXUnPHLDU48VyDwzBM3omK/TTeJXORcXYpKhrt4+vrK7qpGOWFBQ6jMDx8+FB0zVWpUiXbLipJBIdTVAzD5IXwj6mpJl1NtSypqE8pKtldVDxsU7nhGhxGYaB5Y2QbQCTn4oOjp8lFxgzD5GOSeKboTXqBw11URRMWOIzCQB44ZAaZsYsq4z5cZMwwzOcUGBtlahHPqciYxv5MnDhRjI9hlBcWOIzCQIM2ybmauqgat5fU4Mge1cBOxgzD5IVw6ZgGWRGcpGxTVJSeoosuRnlhgcMoDDRElQbcGRgY4Ka0yDhzioqLjBmGyTuSDqrMHjgZZlFlEjh+fn7ChZ+M/tq1a8fLraSwwGEUBsmAOxrVcP09ZBr9fUpRcQ0OwzD5qMGRIXCks6gypagogkPpckpVMcoLCxxGYaCZZDTgjlJUdWfsktkmHsM+OAzDfNYcKhk1OJJZVBpZRzVQupxHNSg3LHAYhRzVkCQjRUWGXeyDwzDM57SJy0xRZVODwxQNWOAwCjmqoUZLZBE4MQmpJyNCj6eJMwwjpzbxypUr4/3791lGxTDKBctWRiGhaE3mGhyJBw5t0s4UUmYYhsmxi0pWDU42beJeXl5o3LixKDRmlBf+lmAUclSDxAcn/QVUeg8cvrJiGCZfPjiyanCySVF9/PgR9+/fx+PHj3mRlRgWOIzC8OTJEzHgjq6cknKI4LAHDsMweSUiLUVlopf3FBWPaigacA0OozBQ9IYG3FEtTgUZoxrYxZhhmPxAQ3slERxjGRGc7IZtSuoBuYtKuWGBwyhkF1W0jC6qGDb5YxgmH9BFkSQaLLuLSnYNjoWFBYYNGwYTExNebyWGBQ6jkF1Urg0lAkfWmAY2+WMYJu8dVBSh0dZQy7YGJ3OKysjICB4eHtDW1uZlVmJY4DAKg6WlpRhwZ2hoiFMyanAkKSpdbhFnGOYLxzTklKKiVDmNaKAUFRmQMsoJCxxGYaARDTTgjq6aTjxN3Za+W4qLjBmG+TwPnKz1N+lTVJltJ2jIZo0aNXhUg5LDAodRGKg9nHwnKEVlP2qrjAiOZNAmp6gYhsm7i7FRNhEc6bBNNbUsbeJ3797lImMlhwUOozDQVRMNuKMi41BpkTGyRHDIB4dhGOZLXIwzDNvMFMGh7qvExERxY5QX/qZgFIaYmBgx4I4iOKVqy+iikhQZcw0OwzBfWIOTmJQs7bDKXINTsWJFBAYGQi1TZIdRLtjoj1FIkmX64EhSVKzLGYbJz5gGWS7GqdEbWV1Uz549Q7du3TBkyBBeZiWGBQ6jMFSoUEEMuCNH45xmUbGTMcMw+UlRGclIUUk6qAhNtaxGfxRNvn37Ni+0EsOXwozCQK2ZVGRMJlvJDSZmqcGRtolzkTHDMHlA6mKcg8mfuqoK1DMJHDIdpXpAKysrXmclhgUOozBIBtxRDY55vbQaHBkRHE5RMQyTnzlUstrEszP5Sx/B4VENyg0LHEYhRzU8k1GDwykqhmE+p01cVgQnO5M/wszMDP3794epqSkvuBLDAodRyFENxpXSanBkFRlzFxXDMF9Yg5PdHCqCLrIGDBgALS0tXmclhgUOozDQSYUG3NEcmN1pRcbp9M2nNnHuomIYJj+TxGXW4CTJ9MAhvLy80KRJEx7VoOTItYuKOmL69esnvrDoRv8ODw/Pdv+EhARMnjxZdNNQuqJEiRL49ttv8fr16wz7NW7cWFj4p7/16tVLnm+FKQQMDAzEgLumTZsizedP2kVFXVUx8exkzDBM3ohNSJamoWS2iUtdjLN+DdK4mHLlysHJyYmXW4mRq8Dp3bs37t27h2PHjokb/ZtETk5Gb3fu3MGMGTPEzz179sDb2xsdOnTIsu/333+PoKAg6W3t2rXyfCtMIfD8+XMx4I5EbVKawpEIHMkkcYKLjBmGyWv9DXVJ6WnKmCSejYux5GKbjP7evHnDC63EyC1FRV4mJGquXbuGWrVqiW3r168XrXcU/nNxccnyGIrynDx5MsO25cuXo2bNmmKia+nSpTPY+ltbW8vr5TNfAcmAO0pVPUnJmKKKTqu/IcEjqyiQYRhG5pgGXY0MQ3szR3Bk1eAkJSUhIiIChoaGvKhKjNwEDrXYkWCRiBuidu3aYhsVksoSOLKgg4wOTmNj4wzb//nnH2zdulX4FLRp0wYzZ84UKQ5ZxMXFiZuEyMhIqUqnG5M/JGtW0GtHnzUNuKMiY9WKqduSk5LE/yc8OlbcpyuxojYfRl7rWVzh9eT1JMKiPoqfhtoaMv+2YuJSIzwUwMn8e/p+evToEdTV1Qv875KPzy8jP5+H3AROcHAwLC0ts2ynbfS7vBAbG4spU6aIVFd6Jd2nTx/Y29uLCM7Dhw8xdepU4Z+SOfojYe7cuZg9e3aW7SdOnBBRA+bzyG69Pxc/Pz8hXj58+ADJp33m1CnoaQAvPtA9dagmJ+DIkSMoihT0ehZ3eD2L93reD6OojRpS4j7IPGfcCkn9fcT7sCy/p/TUunXrxPfO+PHj5fL6lG09FQUqZZGbwJk1a5ZMsZCemzdvip+ywoJU2S5ruyyVRoXDycnJWLVqVZb6Gwnu7u6iEKx69eqibqdq1apZnosE0Lhx4zJEcMjAqWXLlhyC/Azos6E/zhYtWkBDQ6NAn7dLly6g1HjTtQ/FtlYtW8BQRwNX/cKAB7dhbqQPD496KErIaz2LK7yevJ5E1K1XgPdj2JWwgIdH1u+F8BsBwLMnsLWxhodH5Qy/oxENdNFM3xPU+MDHp+IgycDIReCMGDEi144lOzs7eHp6yizQevv2ba7213SC6tGjhyg6PXPmTK4ihEQNfTH4+PjIFDjkZSDLz4Aew18on09Brx993uQ9YWpqBrinDrnT1Ez9f8QlpYpife2i+5nx8cjrqcgo2/EZFZdaY2OipyXzdSempJ5TtDXVs/yeSinIdJS+q+T1npVtPRWF/KxZvgUOFYDSLTeomJhqKm7cuCGKhInr16+LbXTg5CZuSKycPXtWOErmBuVK6XFUu8Eod+iRaresbWyg5Y6MXVQ8aJNhmM9xMZYxpiE3J2OJ6SiPalBu5FaD4+bmhtatW4t0kqSF+4cffhBtwOkLjF1dXUWNTOfOnUX9BY2op1TToUOHRCW7pF6HLLM1NTXFGHsqMKawIQmtx48fixxplSpVUK9e0UpdFDeoHoqEsbGpGR4j46gGyaBNPR60yTBMfuZQyTD5y2D0J0PgmJiYoHv37nm6wGaKqZMxCZFRo0aJWheC/GxWrFiRYR9qGaeoDvHq1SscOHBA/Lty5Yw5UYrmkMEfiZzTp0/jzz//FMWopLDbtm0ruqjU1LK2+zFKGMGx/hTBkQgc6aBNHtPAMEw+28RlIRnVIGvYJjWwTJgwQXzfMMqLXAUORV2olTsnqOg4fe1O+vuyIEFz/vz5AnuNjOJAxwsNuNPRM8DRtG2SYeI8SZxhmM8btJlbikpNpo8bWZzQ9w15sDHKCc+iYhRK4FCR8cdE4Oip6Aw1ONJBmzyHimGY/ERwZAzazC1FRZGbMmXKiHFBjPLCAodRGHx9fcWAOytrG2j3Xy+2SSwFPhUZcxqSYZjcichh0GaGWVQyBA5lEuLj48WNUV5Y4DAKg2TAnZGpGV6ni96kn0XFERyGYfIXwckmRSWZRSVD4FBXLs04JCdjRnnhT49RGCQD7uISUsVMOn3DNTgMw+SZ2IQkfExITUEZ5RLB0dLIGhWmCy0aG8M+NcoNTy1kFAbJgLuoqMgMHVTph23qcw0OwzC5EJmWnqKLJINszhnSGhy1rF+Dr1+/Fq79S5Ys4bVWYjiCwygM5J3k7e2NN1Hx6LvTP4PAkfjg6GpyDQ7DMDkTniZwjHQ0oJo+FCyjTVyLpm1mfnx4OPbv389Gf0oOCxxGYQgICMCwYcOgpWcIuA6SWYPDERyGYfLugZO9j01OTsZ6enrSUQ2M8sICh1EYyLjx1KlTsLSyho7rIKSfyco+OAzD5JXwmHhpBCc74nLwweFRDUUDFjiMwiC5atIxMIZvpi4qSYqKIzgMw+Q1RZVdi3j6GhxZbeI0bJPGAVlYWPBiKzEscBiFQXLVZGFlDd3Kn4qMk5JTEJvW8cBt4gzD5HkOVQ4RnJxSVLa2tpg/fz53USk5LHAYhcHY2FgMuIOWHm6kn0OVVn9D8LBNhmG+dExDbimqR48eoUaNGjyqQclhgcMoDBQOpgF3L97H4sbZKEi6NyX1NxpqKjJPRgzDMLKKjPNSgyMrRUUGf3Q+4mniyg0LHEZh8PHxEQPuLCytoTvwf1kmievyJHGGYQqoBienFJWampqowzEwMOD1VmJY4DAKg2TAnb6xGT6kS1FJBm1ygTHDMPmqwclDkbEsH5y4uDgxG49+MsoLCxxGYUhOThbD7RIT0lxIM6WouP6GYZh81eBkM4cqOTkFCUkp4t+aMpyMXV1dcfnyZXHRxSgvLHAYhSExMVEMuEtISoYehYmlERwetMkwzGfU4GQTwZEM2sxuFlVISAhWrFghanCqV6/OS6+k8CwqRmFwcXERA+6Wbd4t7lOKilrEPV+Fi/uJScniPsMwzJe0iUsGbWZXg/Pu3Tts27ZNjGtglBeO4DAKQ3BwsBhwl6CuCzj2QUx8EurPO4OgiFjx+weBkeL+zPbl0Nrd5mu/XIZhFBCKAEelRX2zaxOX1N9QkFhdxqwqMh2tU6cOj2pQcjiCwygMNEmcrpiuXzwj7gdHxkrFjYTgiFgM3XoHxx4GfaVXyTCMIhOR1kFFGGpnN0n8UweVSvqZMOlMR69evYrbt2/L8ZUy8oYjOIzCjWpQ0TbAq2z2oQQVnY5mH3yMFuWsM4xzYBiGkdTfGGirQ11GAXFuJn/isQYGaNSoESwtLXlBlRgWOIzCjWowMrOEcY3s9yORQ5GdG8/foU5Zs8J8iQzDKDgRUhfjz5tDRdjZ2WHTpk3CD4dRXljgMAqDoaGhGHD3+qMa3udh/5CojOkrhmEYSQQnuxbx3Ez+iAcPHvCohiIA1+AwCoONjY0YcNe+//A87W9poC3318QwjJIKnBwjODkLHFVVVejq6kJHR0dOr5IpDFjgMAqDl5cX3N3dsXJi/zztf+B+IKJiPxUUMgzDSMY05G0OlewUlJaWFpycnGBvb88LqsSwwGEUBsmAO0MjE+m2zCXE6e9vuxGAVksu4KxXSKG9RoZhFJuImDzU4CQk5RjB+fjxI+7fv4/Hjx/L6VUyhQHX4DAKA4WFacCdhr4+YgA4WugjOj4xQ6u4tZG28MEx0tHE5N2eePkuBgM33kTXqraY0c4tW98LhmGK2aDNnGpw0pyMsxM4FL05ceIEtLU5Da7MsMBhFAaaQ0UD7kzMrWAouqVSsLBbJRG2Cf0QJ2puatqbSlvDj41pgEUnvPHX5efYfecVLvi8xa+d3NGqvPXXfisMwyhyDU6ak7GsMQ3E+/fvsXPnTpiamqJBgwZyeqWMvOEUFaMw0FXToi37Yd55urj/7G00+my4jgn/3RdXWtQSnt73RldTHTPalcOuIXVR1kIPb6Pi8OOW2xjx7x2EfeApwAxTrCM4OURzpTU42fjkhIaG4n//+58Y18AoLyxwGIVh16VH+Pn3xQi6dSxf7sXVypjg8KgGGNa4rBBAhzyD0GLJBRy4/xopKTy7imGKZQ1ODkXG8Wk+OFoaqjmajvKgTeWGBQ6jENAQzUWH7iL6yXl89LmW4XcSiULuxdkN29TWUMOk1q7YP7weXK0N8C46HqO23cUPW27jTST75TBM8YvgfH6buMR09NatW3J6lUxhwAKHUQjIlTgsTgVaJVyhae2Yo3txTriXNMKBEfUxtrkzNNRUcPLxG7RYfB7/3QrgaA7DFAMKwgeHPHBq1KiBSpUqyelVMkovcKhQq1+/fqIzhm707/Dw8BwfM2DAADH8LP2tdu3aGfaJi4vDyJEjYW5uLkKJHTp0wKtX2U0vYpQBciVOSYhD3OuniA/2zXG/3CD79dHNnXBoZANUsjVCZGwiJu7yRP+NN/HqPfVnMQxTFKEIb2SaNxZ1WubuZKyWbT3goUOHsHHjRjm9UkbpBU7v3r1x7949HDt2TNzo3yRycqN169YICgqS3o4cOZLh92PGjMHevXuxfft2XLp0CR8+fEC7du2QlJSaV2WUD+qQUtHUhlYpdxHFyY781NS4WBtg99C6mNrGVYieC95vhW/OlmsvkJxNqothGOWFjD8lp4jsjP5IBD0P/SD+Tc0IstLe5IFjZWWFqlWryvcFM8rZJv7kyRMhaq5du4ZatWqJbevXr0edOnWEY62Li0u2jyUXSWtr2a2+ERER2LBhA7Zs2YLmzZuLbVu3bkWpUqVw6tQptGrVSk7viJEn1P5ta1sKKm3HAirZ626KxPiFxoiCYqq7yQ2aJvxjo7JoXs4Kk3d54taL95ix7yEO3X+NeV0rws5cr4DfCcMwX4uwD6kFxlpqKrj94n0GWwmCGhWolk/irXXQM0icE8hbq7W7zVd73YySCZyrV6+KtJRE3BCUaqJtVLyVk8A5d+6cGFNvbGwsRtb/9ttv0rH1t2/fRkJCAlq2bCndv0SJEsLin55XlsChlBbdJERGRoqf9Dx0Y/KHZM0Keu36OKtg9M+DoaKuiVLjdov0JEH/pWssFyt9eL35gGWnfbDndgB+8nBFU1cL6X45UdpYC/8Mqo6tNwKw8IQ3rj9/h9Z/XsC45k74tnbpDCfBorKexRVez+K5nscfvcHMg0/Ev+OSUvDN+muwNtQS54lW5a3E70duvy9tWsjcpbm8VyWxn8STq1y5cnB0dCzw960s66mo5Gfd5CZwgoODpaIkPbSNfpcdbdq0Qffu3VGmTBk8f/4cM2bMQNOmTYWwocgOPVZTUxMmJp/s/AkKJ2b3vHPnzsXs2bOzbCenSiomYz6PkydPFujSqYcFip/a5qXwZutEmDT/AVo2zjDSTEEXu2RUNA3HfSMV7PVXxavwWAz59x7KGSejq30yzPNoOGpOUSB3YNszVfhEAr8f9cK/F5/im7JJsNYtWutZ3OH1LD7reT9MBX95SyK/ny5WgiNjMWL7PQx0ThbnjVRxk/FiJnVbCn7acw9xfgnw9fEWY2Mo00A1o5lLJIrDeioyMTEx8hM4s2bNkikW0nPz5k3xU9aVNdVQ5HTF3bNnT+m/KSpDPgQkdg4fPowuXbpk+7icnnfq1KkYN25chggOpbQoCmRoSJ65TH4VNP1xtmjRAhoa2Xcq5Bf6DOn5Fi5cCG/vp2is+wrtGtVAXbdSsC1ZQuzTFsDouESsPv8cf13xx+NwVfg+UMf39e3wYwN76GjmnrYi+qakYOftQMw95gX/D0lY+FADo5qUxeD6dtDIxvxL2dazuMLrWbzWk2po5i66QHEbGb9VEXJmf6A2wuNzuvJXwds3QZgxayW8Hnnixo0bmDdvHipXroyGDRsWq/VUdCQZGLkInBEjRqBXr1457mNnZwdPT0+8efMmy+/evn0roi15xcbGRggcHx8fcZ9qcyh8SB1a6aM4ISEhwphJFhT5oVtm6ODiA+zzkcf6ff/992jbti3mz5+P336bI6J3Qx4/xooVK8RxR5+jsYYGprYthx41S2PWgUe46BOKlef8sP9+EH5uVw4tylnlKW3Vt449mrpZY/reBzjr9RaLTvni2OMQzO9WEeVLGKGw4eOR11ORUdTj89azMARHZu9cThGad9HZi5vkhFiEHVsu/LcMXMuLzlzKHowfPx7FcT0VnfysWb4vVak129XVNccbDSijYmIqCCYlLOH69etiW3ZCRBZhYWEICAgQQoeoVq2aeIPpw3vUafXw4cN8PS+juFBN1dKlS0UoUmIVsHLlSpETv3z5snS/shb6+HtQTazqUxU2Rtp49f6jMPYbtOkm/EOj8/b/MtbBXwNqYHGPSqLr4tHrSHRccRmLT3ghLs3tlGEYxcXnTdRnPY6ETcT1XUByMpKiwpCSEI+GzVqKi+lu3boV+OtkCh+5xeLd3NxEuzddkVMnFd3o39TOnb7AmAQRtXwT1O49YcIEUaDs7+8vio3bt28vRFXnzp3FPlSkPHjwYKGuT58+jbt376Jv376oUKGCtKuKKRpYWFiIwvGDBw8iMDAQL1++FMfCmjVr8ODBA7EPiR+PCjY4Pb4RhjYuK8z9KBrTcskFIVI+xucuUug5ulS1xclxDdG6vDUSk1Ow7Iwv2i+/hHsBOfs2MQzzdXgd/lF0RM459DjfqfAYn+sI+W8Wws9tQsS1/2DWchjch6/E6kW/5yvDwCg2ci02+Oeff4TwoFoXulWsWFG0d6eHCrkoqkOoqamJL66OHTvC2dkZ/fv3Fz9J8BgYGEgfs2TJEnTq1Ak9evRAvXr1RKEwfQnS45mihaqqquiko+PkwIEDSExMFGnSKlWqiLQVRfgkgzcnt3bFsTEN0cDJHPFJyUKkNF98HiceBefJP4e8eNb0qyYiQmZ6mvB+8wFdVl3G3CNPEJvA0RyGUQTIrJPSyo0WnBWeVnRBoqmWc0qaHIt/ausm6nFC983F2z2/QEVVDWoGFtCyKgtN89JYMKTTV+2mZAoelZRiOI2QipQoEkTCiouMP69IjjoLPDw8Cj2HTFEcKhintCWlJakmh1Kf5DwqgQ7pYw+D8cuhx3id5nfR2MUCs9qXz7PvDc2ymnPwEfbdey3u25vridqcGnamRWo9iyK8nkVzPQPexWDVOV/suv0KCUmpX1u1HUwxupkzIj7Gi1bv7L7M/mhTGqe3LkPJSvXx95n78D+4CsYN+sKgaluUMDMsVB8cRVnP4vD9Lbc2cYaRB6VLl8auXbuEbcCgQYNEZI/qs8j/iNKb1JlAKac2FWzQyMUCK874Yv1FP5yjtJXvBfzYyAHDGjvm2m1lqqeJpb2qoF3FEpi+7wGeh0ajx9qr6F/HDhNbuUBPi/90GKYweBkWg5VnfbH7zisRrSHqljXD6GZOqOVgJi5oTj0JgaWhFt5kKja20lNDxcir2Pj7MmEL4uh4Do8ePMSJWwMBfVMRtc1sBsgUHfgszSglVGx+584dhIaGioJkOnk9fvxY1HhRhIciOpS2ognjXavZSrutlp/xxZ47gfi5fTm0zEO3FTkg17A3xe+Hn2DHrQBsuuKPU0/e4I8uFVHfiVx1GIaRB9QoQMJmz91A6TiF+o7mYs6cJJJ69+V7zD3yFDf8U4fwGuuoo2PlkqhY0ghmOio48vdyLFg4XzQoUMRk2rRp0NHWQsf6FfhDKwawwGGUFqq5ooLA4cOHi1ocSllRAfL+/ftFkTqFgqnlU9JtdfxRMOYcfIzA8I/4ccvtPKetqLtqXreKaFfJBlN2PxDdWn03XEevGqUwra0bDLU5zMwwBQVFS5ef8cH+e6+lwqahswVGN3NEtTKmUvGz4LgXDj8IktbYDK5vjyGNyyLk1QuMGPG9SGP8+eef2Ld3r4juUk0n1fQxxQcWOIzSQ35IVHj+9OlTjB07VnTe0VwySmGRnw4NeKVIDeXY6URJV4XrLuQ/bdXAyQInxjbE/GNPsfnqC2y/GSCe4/cu7mjqyp0XDPMlPHv7QaSU998LhGT+JV2EjGrmhKqlTaTDMWlUyz/XX4p0FQVgu1W1xbiWztBOjsXPUybi0aNHOHPmjHAjpr9/Oi+wsCmesMBhigxkOSCxVadOPDKa3LhxI5YvXy4EUP369UXaamIrV3StaouZn5G2otqb2R3dRWv65N2e8A+LwaBNt9C5SklhMmiip1mI75hhlB/fkCjxN3jw/mupsGnqaimETeVSxuI+2T38dfk5Vp97hg9xiVLxM6WNKxzNdcVwZzJ7pYgNRXbJvX7AgAHCJJYpvrDAYYoUEnFChciUrtqxYwdu3bqFtWvXomzZsiJtRYXKDl+YtqLixqOjG2LJKW/876If9t4NxEWft/ilo7socGYYJme830SJaAylmSS9vM3drETxcAXbVCdxSlHtuh2AxSe9pQXE7iUNMbWNG+o5mou0dNWWvUV3JRn0kT8aNRxQswHDsMBhiiQ0kHXUqFFivAPNT5s+fbrIw5Op5KJFi0Q+nvyTviRtRb+b5uGGNu7WmLTLEz4hHzD0nzvwqGCN2R3cYWGQdTwIwxR3ngZHYvlpXxx5+EnYUOSUIjbuJVOFDXVGnfUKwR9Hnwo/KqKksQ4mtXZB+4olEBYWiq5duwoDWeqkpBQUiR2aY8cwErjiiinS0PT6VatWCQdkckOOjY3Fzp07M6SzJGmr4+lMApenmQRShCc3q6gqpU1waFR9jGzqKNpNjzwIRosl57HvbmCeDAYZpjjwJCgSQ7feRuulF6VRG3IOPzKqAdZ9W10qbjxfheOb9ddE6pfEDRX5k0nfmQmN0MzRCCtWLBfNA9RFSWnoH374QURvmjRp8rXfIqNgcASHKRaQIeDZs2fF6AfKzZNRIOXqKaJDKStqO0+ftvrl0BNp2qqRswVmdSgvzP6yfX51NYxv6YJW5VOjOY+DIjFmxz1RV/Bb5wqwNtIu1PfLMIrCw8AIkYo68Th1+DJlkT3cbTCymSNcrQ0z+N0sOOEl/mYITXVVDKxnh2GNHGGgrSasIKiJgIqGyeBt8+bNMDU1hbu7+1d7b4xiwwKHKVb1OTTagybdk7CpXr268MsJDw/HH3/8IdJW1HaePm21/sJznPd+i1ZLLuCHhg4Y3iTntBVdhe4fUQ9rzz/DstO+OP00BDcWn8f0tm7oWaNUnqacM0xR4MGrCPx52kf4RhF06LetYCNSUc5Wn0bvvI+OFxHTLdf8hUMx7UdF+3TBQGmppKQkMa7l0qVLaNq0qaijs7a2RsOGDb/iu2OUAU5RMcUOHR0d9O7dW5w4yfzL0dERy5YtE2JHMq1ckrY6NqaBEDuUtlpxNjVtRWMgcko9aaipYkRTJxweVV90gUTFJWLKngfot+GGsJtnmKLM/YBwDNp0E+1XXBLihkyCO1YugRNjGmJF76pScUPz3agrquGCs6JDisQNpYgPjayPxT0qA9HvxEUHGXjWqlUL+vr6oqaO2sBpkDPD5AZHcJhiXZ+zdetWcQKlEylFVypXriym0/fs2VO4IlPaavPAGjj+6I2YbUVpqyFbbwvRMzuXtJWTlQF2D62LjZefC1OyS76haLX0ghgK2q92GaiyPTxThCBXYYrYUKE+kSpsSmJEU0dhtimBOqOo63DRCS8Epc2Kc7OhzihX8XdFdXJz587FsWPHcOHCBQQFBYmuSOqQorEsDJNXWOAwxR6ycaeBnS9evBBT6f/55x/hhvzdd9/h+++/F79v7W4tanEk3VYX0tJW3ze0F2krivjIgoqOv2vggGZuVsI358bzd8J/55Dna8zrWlEIKIZRZm6/SBU29DchOeYpYjOiiWOG45uinhd8QjH3yBM8DY4S20oYaWNCKxd0qlxSpKYo/TRw4EBs374ddevWFYXDv/32m6i54cHITH5hgcMwdLWpqgp7e3vRbTVlyhTcu3dPzLgiLx3y2KCTq46WljgZ02wrEil0Ql959hn23X2NGe3KoVX57E0CKdKz/fva+Of6C8w9+hQ3/d+jzZ8XMb6lM76tVYo/A0bpuOn/Dn+e8hGRSYmw6VKlpBD8mX2kqNCYWr4l+xpoqwsB1L+uHbQ11PDgwQNh60BmnGTnQKliGsFCqWSG+VxY4DBMOshTg8Ljfn5+IiROTqg00LNChQr4+eefMXToUCFWPidtRSmpfnXs0NjFEtP2PhAuyr8feSqiOR5m/DEwysF1vzARsbnyLEzcV1dVEc7gJGxKm+lm2PfV+xgsPO6FfffSOqPUVPFtnTJiX3L9fvv2LcbMmIHnz5/j3LlzuH//PiZPnoxnz55BQ4NnvDFfBgschpGBg4OD6LSiQmS6spSMfaCanV9++QXNmjX77LRVKVNd0Y7+361X+OXwY3i+isSjQDV8NHuGEc2cRZEywyga15+/w4pzfrjm904qbLpXtxWGmHRMpyciJgErz/li02V/UaBPUNpqQksXsS+lokjE0Mw4chmndm+KnA4ZMkQUEzNMQcACh2FygLxyaL5N+fLlsW/fPpw8eVK4IteuXVuIHhJCkrTVrAOPREv5p7SVm/DFkZW2om09apQSUZ9pezxxxust/jzzDCeevMWCbhWlpmcM8zWhuhmK1Cx7qIZnV2+JbRpqKuhRvRSGNi4LW5OMwoY6o/6+6i+GZkbGps6MquNgJhy/JeMXqN6NvKjo4oEiNhcvXhRmfdz2zRQ0KinF0Go1MjISRkZGiIiI4MK1z4CuvsgFmFqsi1MY+d27d2LsQ58+fUSHx++//y4KIH/88UeR2qI/JTIzk8y2IkjAzGpfLsdi4vj4ePy65RgOBmrjfUyCqGUY0sgBI5s6ifoEJn8U1+OzIKFjmeplqMbm1ov3UmHTq0ZpIWxKGOtk2D85OQX77wdi4XFv6bHvYmWAKR6uaOxsIQQ9FfEPGzZM/KTUFP0/Tp8+LdK/xQk+Pgvv+5tj4QyTRyiMTn45NWvWFEXIJEyOHj0KZ2dnMdSTTuIUsTk1rpEooKR6A0pbkTX9guNPEROfekWbGXpcNfMUHB1ZVxihURstRYHaLb+EOy9Tv1wYpjAg0UFRyK6rrwjfJhI35CjcwDoZp8c2wC+d3LOIm0s+ocLzZuyO+0LcWBtqiyjkkdEN0MTFUnwRbdmyRYjN8+fPw9vbG7NnzxbjFYqbuGEKF05RMUw+IUGyZ88ecbKeOnUqgoODERISIsZA0O/q1KkjM221904gfm5fLtu0lZm+Flb2qYr2D4Pw075H8A35IL5oBtezF66uOTkoM8yXChvyr6Hi4XsB4WKblroqetcqjcF1S+P2pTOwyTRu5PHrSPxx7Km0PdxASx1Dm5TFwLr24lilFNSJE6fx7bffinQupaI2bdqEihUriosChpE3LHAY5jMggdK4cWPR+UEnbTqJ0zyrJ0+eiG4r8s+xt7XFpoE1MqSthmy9I9xaqdsqu7QVjYqo7WCGOYceY8+dQPzv0nOcfPJG+ObQdoYpSGFz5mmImBV1/1WE2KatoYo+tcrgx4YOsDTUFimV9LwO/4hFJ7yx5+4rMTCTUld9a5cRKVVTPU1pGoHGK9CEbxL8JiYmYnu3bt34w2MKDU5RMcwXDvGkGhz6EiBjMnJHpu4rFxcXUZScPm1F08YpbUXt4eRoPP9Y9mkrY11NYVe/cUANceX8IiwGvdZdw4x9D/EhTvZjGCY/wubk4zfosOIyBm++JcSNjoYavm9gj4uTmgpfJxI36Yn8mCC8bBovPIfdd1LFTbuKNuLYntm+vBA3/v7+wqiP0relSpUSU7+pgJjmv5HHDcMUJhzBYZgCgIrd/ve//4kTPI16oC8QiuhQ6yv9pCtXSjN1qWqL2QcfiXTAqnPUbRWIqW1cxJeFLJq4WuL42IaYe+Qptt14iS3XXogr7rldKogCZobJD1QMTBFFitjQxHtCV1MN/eqUwfcNHGCur5XlMXGJyTgXpIKZSy4h/GNqNKemvanojKJZa0R0dDRWrFiBbdu2ic4oXV1drFq1SlwAWFjwccp8HVjgMEwBYmdnJ2oNKFVFM3TmzZsnto8ZM0aksapUqSKiMnT1PDstbTVy+324GqmiXK1oONukfmGkx1BbQwgaulqesscTAe8+4tu/bqBHdVtMb1sORjrcKcTkLmyOPwoWNTaSMQl6mmr4tq4dvqtvL+q/ZD3m0IMgEWl89Z7qvxLgaKkvZkY1dbUU0UkS8nRr2bKlqEGjad80XoGimra2tvyxMF8VFjgMU8DQiZ/mV8XExIi2crKdX7lypfDToatb+l3L8tZo4GSB1ed8sfr8MzyNUEXbFVfEVTQNJ5RlEljP0RzHRjcUgzs3X/XHzluvRCTot84V0KKcFX+OjEyRcvRhsIjYeL1JFTb6WuroX7cMvqvvINyEZXHlWahIR3mm1eUYaqRgskd59KxZBuppRpQ3b97E6NGjhTkf/Xz9+rUwxezQoUO2I0sYpjBhgcMwcoLC9DNnzkRgYKCYr0PtsiVLlhT1OeQHMnLkSIxr6YIOFa0xctMFPAlXlaatqAaCnJIzf1HoaaljVofyaFvRBpN3ecIvNBrf/30LHSqVENslRZ5M8YasBg4/CMLy0z7wCfkg7XIaWM8Og+rbixovWXgFR2HesaciDSqJ8lBdTomop+hc3VaIG4pMUvdgQEAArl69irCwMBGxJGGjrZ2xbodhviYscBhGzpCoodoEKrxcsmSJsKhfvXq18NAZO3YsWrRogR9dk6HlUBW/HvESaauh/+TcbVXDzlT4jCw55Y31F/xw4P5rXPYNxeyO5YWXDl9BF19hQ7PNlp/xFTYDksGWg+rZi5uRrux0ZnBELJac9MZ/twOQnJI6hoFaxEc1c4KRliqOHHmK2NhYIWYWL16MzZs3w9HRURy/JN5pWC2LG0bRYIHDMIWEpqam+DKgokvqtDp06BAeP34sbm/eBGNQW0s0drUWaas1F/yk3VbZpa3I5XhqGzd4uNtg0i5PkYIY8e9dHCz/Gr90dM/SBcMUXRKTknEwTdj4vY0W2wy11TG4vgMG1LPLtk4rMjYBa88/w4ZLzxGbkDozqo27NSa2cpEKaxLmNO173LhxQtSQqSUds2SHUKtWrUJ8lwyTP1jgMEwhz7YaNGiQ6Kr69ddfUbVqVRw4cADDhw8XE5VnzJgh0laSbquz6bqtsktbVSpljIMj64uhn3SjKec0EPHnduXQpWpJjuYUcWFDk7rpc38emipsjHU1ROFw/7p2MNCWLWziE5Px7/UXWHbGF++i48W26mVMMNXDDdXKpHrWEE+fPhXHJtXXvHr1SkRxPn78iMOHDxfSO2SYz4cFDsN8pbby+fPni39Txwm5vt64cUN0Yc2dOxffffcd/hpQA6eehAih8+r9p7QV1dqUzZS2Ijv9sS2chefOpN338TAwEuP/uy+u6n/vXCGLvT6j3CQkJWPv3UAhbMgjiTAhYdPAQQgbKiSWBXU8HXkQjPnHn0of52ChhymtXUWhukQ8UyqKhmKam5vjzJkzIvpIbd89e/bkad+M0sACh2G+MuQfQnU69IUSGhoqOq1oTg8VJ7do3FiIGorirDn/TKStWi+9IL7IRspIW5UrYYh9w+ph3UU/LD3lI7qsWi65gKkervimRmmoqnJ3i7ILmz13XmHFWV9hF0BQYTmlMcnLJjthQ9x4/g6/H3kiHcVAnjdjWzihZ/VS0s4oMqyk2jCa9k2RGi8vLzF/jfxsyOaAh5cyyoRcnYzfv3+Pfv36icmfdKN/h4en/nFlB11ByLotWLBAug9Z5Gf+fa9eveT5VhhGrlSqVAn//fefEDvUWk7Fm+QnQjU7Qa9eYlwLZ5wY0xBNXCyQkJSC1eeeodmi8zjyIEhclaeHvqyGNXbEkVH1UbW0sXA+nr73Ifr87zpepl21M8pFakrpJRovOIfJux8IcWOur4lpHq64NLmJmPCdnbjxDYnCd5tvocfaq0LckLHfmOZOOD+xsRjJIBE3NOWbjkNKoZILMc2LonM4tYFbWbENAaN8yFXg9O7dW0xdPnbsmLjRv0nk5AS1IKa//fXXX0LAdO3aNcN+NOsn/X5r166V51thGLlDV8dU70DpK0pVUVcKRXXc3Nywfv162JnribTV+m+rw9ZEB0ERsRj2zx0x9VnSMZMeR0sD/DekrqjdoflCV/3CRNHyX5eei24bRvGJS0zC1msv0GThOUzb+0B02FHk5ae2bmKkwg8Ny8r0TCJCImMxdc8DEcE79eQN1FRV0KdWaZyb2BhjmjsLywGCooV03NGYEeqGIkhk37lzRwgehlFW5JaiIl8EEjXXrl2TVtrTSZoGr1HYk7xAZGFtbZ3h/v79+8WVrIODQxaPkcz7MkxREToUyZk+fTr69++PuLg4uLu7ixoIAwMD9OnTBw2cGknTVpd8Q9HmzwuiY4bSVpIvLoK+1AbXt0dzN0tM3u0pio9piCd5pNDwTnKmZRSP2IQk/HcrQHzGJGQJCwMtDGlUFr1rls5xsjxF7NZd8BP2AR8TksS2luWsMKm1a4bPm3yZ6AJy0aJFIh3q5OSEnTt3wsbGRjock2GUGbkJHDKAorRU+jbC2rVri21k6Z2dwEnPmzdvRLU+eS5k5p9//sHWrVtF6LRNmzbCUI1O/rKgLwi6SaBJt5J8c+ZJuUzuSNaM106+60kFntRKfvv2bfGl07x5c+GOfOnSJVEPMbJxTXSoaIVfDj/Fee9QIXb23X2FaW1c0Lr8p4JRooShJjb3r4Ydt19h3nFv3H7xHh7LLmJUk7IYXO9TmqIooMzHZ1xCEnbcDsS6i8/xJjL1nGVloIUfGtqjR7WSwhoASEZCWkt35vqcnbcDsfzMM4SldUZVLmWEya2cRYeU2CchQaQ0ExMTRQefn5+fuICsUKECmjVrJkSOZL+isJ6KCK/nl5Gf41AlJXMCv4D4/fffsWnTJnh7e2fYTnldmjZLTpi5QV0mf/zxh2hRTG8iRZEge3t7EcF5+PCheC7yZzh58qTM56Fw6+zZs7Ns//fff0UkiGEUHfIiOXjwoOi0opQsiXRqKafZVioqqnj4XgV7/FXxLi5V1DgbJaObfTKsZDRPvYsDdjxTFeMhiFJ6KehdNgkl9Ar7XTES4pOAKyEqOB2oisiE1M/QWDMFzUsmo7ZlCjRy0J90Bvd8p4JDL1UREpv6WHPtFLQvnYxKpilI7ypA/jUbNmwQKX4qZj9//rzo4uNUFKMs0EUelb9QBJLS+QUqcLITC+mhGSUnTpwQkRdKR6WHrhAGDx4spiznhqurq3B5Xb58eY770RVu9erVxU+6KslLBIeK6KhjJbcFYmQraBKT9NlwV0Xhricds5MnTxZCh2405LBt27ZimCfUNLD2wnOsu+QvilI11FQwsG4ZDGvkkCFtRdCf/d57r/HbES9ExiaKfYc2dMCPDe1Fy7kyo0zH58f4JGy7GYD1l/wR+iE16mJjpC0+h25VS0Irl8/izstwEZGjn4SpngZGNimLntVtoZEuKvfy5UshiOknzUajaM2+fftElI9awIvKeioDvJ5fBn1/U3Q7LwIn3ymqESNG5NqxRAWSnp6eIsWUmbdv3+apIp8mMpM4ItfM3CBRQ394VCwnS+BQiyPdMkOP4T/Yz4fXr/DXk1JVf//9N6Kjo0VEh4qQ6bi/deuWSFuN79IF3WuUFpPKaZ7Quov+OOgZjJ/aloNHhYwmgT1r2qGJqzWm73soppsvO/sMJ56EYH63iqhom3WqubKhyMdnTHyiKB6mWhmJsClprIPhTRzRtRoJm+xrbAi/tx8w/5gXjj0KFvepiJxaxX9o6JDB3I+udikdNXHiROGeXaNGDVFQTC7E+vr6RWY9lRFez88jP8dgvgUOKSe65QYVE5PCoqvMmjVrim10MqZtdevWzfXxFEatVq1ankKnjx49EqqYTv4MUxzQ09MTpmtkEEhX4rt27RJRU6qpUE9OTjUJfPwGs9JMAof/ewf1HVNNAtMXmtI4h3X9quGQZxBmHniEp8FR6LTysujOoVbi1JoPpqCIjkvElmsvRAGwpE6GOuJGNHEU7tW5Rc/eRsXhz9Pe2HYjQHTCka1Rj+qlhMmjVbrRHBSho/pFGuraqVMnYR5JV75kt0FpTYYpDsityJhaW1u3bi1yvZIW7h9++AHt2rXLUGBMaSj64+vcubN0G/0hkicIVfdnhgYVUoGxh4eHEFqUUx4/frz4o61Xr5683g7DKBwUjaGOKvoCo783qlMjO326gKC/OxqKeGpcI+GZszqHbit6nvaVSqBuWTPMOvgYB++/FgXLJx4HY0G3iqhWxvRrv1WlhzqbNl/xx/8u+uF9TGqRZBkzXRGx6VylZIZ0UnbC6H8Xn2PdhWeIpoIdILUzrrUrnKwyNlfcvXtXTKqn44Emfh85ckTUM546dUqO75BhipmTMQmRUaNGiToBokOHDqL9NT2UhqKoTnq2b98urkC++eabLM9J+eLTp0/jzz//xIcPH0QtDdUgUBcVzflhmOIYzZkzZ474N/2klASliClVTAX4o0eOFDOp5hx8jNNPQ4R42X8vMEvaykxfC8u/qYL2FW3w076HYmhjtzVXMaCunRi+mJ3fCpM9UbEJqcLm0nOEpwkbOzNdjGjqhE6VS+TavUazpnbeeiWmxlP0hqhkayRmRtV2MMuwb3BwsJhn5u/vL+psaKgrNXr06NGDJ30zxRK5nrFMTU1FK3dOyKpxpkgP3WRBgoYq/xmGyQrVVpDT9+rVq8WXHKUpyEvn6a3r2DCglUhbzT70SDjhUtqqnqMZZndwz5C2alneGrXszfDr4cf47/YrbLzsL4zi5nWpiLqOuaenmdQp3Zsu+4sp3REfU4WNg7memArfoVLuwobOizSH7I+jT/AsbTp4aVNdTGrtgrYVbDLUUlEDxZ49e0Q3FBVd0oBMagQhR2JbW1v+OJhiC1+SMUwRo2HDhiJd27RpU5GuoqjOkiVLMHToUDEC4uTYT2mry75hIm01qL49RjV1kqatjHQ1sKB7JbSrVAJTd3sKQdT7f9fRu1ZpTG3jmu2U6uIOiRlyiv7r8nNExSaKbWUt9DCqmRPaVSwhjBdz4+7L95h75Clu+L+TDtGkx9NYhcw1OtTqTZPpCfIYozpEauQgocswxR0WOAxTBKF0LdXhSNK66urqIh1cvnx5cXVPqav0aau15/2w/+5r/NTOLUOEoJGzBY6PbYh5x55i67WXYh7S2ach+L1LBTRxsfzK71JxCI+JF8KGol1RcanCxslSHyObOYn1zIuw8Q+NxoLjXsJlmqAWcXKhHtK4LAwzCUqqPSTXYSoBePfunfiMyZaDzCAl4xYYprjDAodhijhklkmtwZMmTRLdhjS5nFrMycZh3cCBOOsVKk1bjfj3LrY5vsRs0W2VWrxK0ZpfO1VA2wolxLiHl+9iMHDjTSGQfm5XDsa6OfuoFGXeR8eLNNSmK/6ikJhwttIXERcPd5s8TW8P+xCH5Wd8Rdt4YnKqMV+3qrYY19IZNkYZnRpJzJBBKX2WNO2bOlTJ7ZpGeWTn5M4wxRUWOAxTDKDatW3btonJ0OSFQl1XZPpGYx+oVoPSVlR8TLOPKG3VeulFDG6QMW1Vp6wZjo1pgEUnvEUKZs+dQFzwDsWvndzR2r14zYV7Fx2P9Rf98PcVf2lXk6u1gRA2rctb50nYkMkfrSOlCyXiqLGLheiMcrPJamBG4oac4MPCwkSdlbGxsehCzTynj2GYVDiWyTDFiEaNGgmjLKrFoat+GoZL27Zu3ojRzZxwamwj0X5MkQRKWzVbdB6HPF9LmwGok4qmk+8aUlfUloR+iMOQrbdFwTL9u6hD0Za5R5+g/rwzQpiQuClnY4g1favhyKgG8KiQe9SG/Gt23gxA44VnRUqKxI17SUP8810tbBpYM4u4oa5RmuNHpn00d48+NxqFQ8Z9LG4YJns4gsMwxQwSODTa4bvvvsOECRPE+Afy0iFbBvKo+nPSJFx7WVqYBGaXtqpWxgSHRzXA8jM+WHPeD4c9g3DFN1QYCVKXUPoun6IAtWhTxGbL1RfSCd3lSxgKUdiiXMbBptlBIvGc11shkLzffJC6F1NnVPuKJbIII/L8+uWXX3Dnzh08ePBAFIuTzQbZAlBNFcMwOcN/JQxTTCGr/jVr1mDevHmitZjGotAX55MnT4TgOTGmB9Ze8BORCllpK3I5ntjKFW3cbTDhv/vCBXn09ns4eD8Iv3V2z+Csq6yERMWKSNY/118gNm2Cd0VbI7EGzdws8yzkPF+Fi86oq35h4r6RjoYwW+xXp0yWsQxRUVHCt4a6o+7duyemfZNxH80BpE4phmHyBgschinm0JcmFavSuIfdu3cLo01yEr/l4oK+lcugSxVbzDn0SPiyyOq2ci9phAMj6osaHorokGfO9edhmNG2HLpXt1XKaM6byFjxfqhrLC4xVdhUKmWMMc2cRJ1MXt/Ty7AYLDjhJdyhCWrzHljXDsMaO4pW/PQkJyeLdafo2rRp00RxOLlR04263xiGyR8scBiGEV/YXbt2FSNQ6MuUHHHJU4VqPCiSsHTpUtysmTrEk7qoKG31b9mXmNMxNW1FX9xUYNuqvDUm7bqP+68iMGm3Jw56vsbcLhVga6KrFKscHJEmbG68FBPZiSqljUUqilrm8ypsqLuKOqO2XPNHQlJqZxSNZBjf0kWkpTJz9epV0fJNRn3kSEwu8LSNXOCVUSAyjCLAAodhGCk6OjqYPn26+Pf69etFuoQG5tLgW6rZOTRiFDZefYVV53xx5Vla2qq+vfB70ddSh4u1AXYPrStapxed9MZFn1C0WnIBU9q4CqO6vHQXfQ1eh38Uwmb7jQDEJyVL64xI2DRwMs+zyIhNSBJeOLQ+EqM/ejy9//IlsqaXaHYYjZw5fvy4mAhfrlw5rFu3TkyGZ2HDMF8GCxyGYWRCRoEVKlQQQuevv/7Chg0bRJTBNckPJ8c2xJxDT0Q6iup09qXNtmpX0UaMIfixUVk0L2eFybs8cevFe8zY/0hMLJ/XtSLszPUUZsUDwz9i1Vlf/HfrlVTY1LQzxejmTmL4aF5FBnVG7b0biEUnvBAUESu2UTcUuT43dLbIsj952JCfDa1nxYoVxSBMEpPkb2NlZVXA75JhiicscBiGyRZqT6b6HGolt7a2Fm65NPKBBjhSq3LvWqUw60Bq2mrktrvYdiO124omXJe10MfOH+tgy7UXwgn5+vN3aP3nBUxo6YKB9ezz5O4rLwLexQjPn123A0QKiahlnyps6jjkXdhQZ9QFn1DMPfJEFFkTJYy0MaGVCzpVLikzYkWz9Gg+GAkZch0m1+nY2FgsWrSogN8lwxRvWOAwDJMj9CVMKRNi/vz50NLSEmZzZDpHfjon5s4TxceStFWbPy+mzrZKS1v1r2uHpq6WwgWZfv/r4ScimrOgW0UhhAoTKvpdedYXu++8El4/BAkaEjaZp3PnxsPACPxx9Cku+YaK+wba6hjRxFG8X+owywx1RJHvEM0Ko9QUFRUfPXpUzA3jdBTDFDxs9McwTJ6hFAq1lJMTMn1BUzvz4wf3ofrkGI6OrIvmblZCOKy7QCaB53DgfqpJYClTXWFkRwXHJHruBYSj7bJLQmwkpKWG5MmLsGhM/O8+miw6hx23AsRrrO9oLiJM236onS9x8+p9DMbuuId2yy8JcaOpporv6tvjwsQmIjWXWdyEhISIIu2qVauKdm8SjDQqg9axfv36LG4YRk5wBIdhmHxhZ2eH5cuXi64rKj4md93Lly/j4sWLom6nd63q0rTVKEpbXU/ttqJozTc1S4s262l7HuCs11vh5HvkQRDmd6soswj3S3keGo0VZ3xFjRDVyUiKfsc0d0K1Mqb5eq6ImASsPOeLTZf9pfU6HSuXECk3EnCyIBFTq1YtIQZpRAZ1pZH/EK0ZwzDyhQUOwzCfBc1Doi/uvn37IiAgQMy2orZmMqQ7Me0nadqKzO3Sp61ogORfA2oI0UFC6NHrSHRccRnDGpfF8KaOWYzvPoc3H4EJux7goGcQ0nSNEFb0/69a2iRfz0WdUeRgvOKsLyI+JkjTWtM83FDBVrYoO3LkiPCv2b9/v4h20ZDTtWvXonLlyl/83hiGyRsscBiG+Wwo3UIDPGlg5+TJk0WdSZcuXfDLrJ+RmJiIvUPHYvH5V6LbitJW++8FYnrbcmhf0Qadq9iinqM5ft73CMceBWPZGV/xc363SqhcyvizXo9vSBT+POWNQ55qSEGQ2Eb1PyRs8vucyckpIsVGUSbqtiJcrAwwxcMVjbPxxCEX6CVLlojaGqqzoUgXiRwLCwuxVgzDFB4scBiG+WI0NTXFFzs58FJH0IIFC0TUwtfXF23btkWvb9tgzqGnMtNWa/pVE2mqn/c/FDOauqy6jO8aOGBcC2dpPQull248fydGJ1gaaKOmvWmGLiyfN1FCIKUOBqUtKmjqYoExLZxR0Tb/Yumybyh+P/JERJcIa0NtjGvpjK5VbWV2f5GXDc34oqgW1dxQJItmew0bNkyY9zEMU/iwwGEYpsCgSAUVFdOka2opJwO7PXv2iDEQuwY3xfZ7oaKwOHPaiqZwU6HvnIOPsO/eaxHtOfn4jfDNeRcdJxyUJf4yhI2RNma2Lyc8dZaf9sWRh0FpwgZo4WaJyhqv8UP3KkJ05IcnQZGiM+q891tx30BLHUOblMXAuvbQ0cyaOqP3Sh5BU6ZMwdatW/HTTz8JT5uFCxfCycnpS5eTYZgvgAUOwzAFCqVuKGrTokULkaKhehQqqnV3dxcFtztm/4GV194KAZM5bbW0VxW0q1gC0/c9EAXCPdZelfn/ILEzZOudDNtal7fGyGaOcLbQxZEjqbOf8uNkvOiEN/bcfSWEkoaaCvrWLoORTZ1gqqcp8zFnz57FuHHjYGxsLCaykwMxCTnqlGIY5uvDAodhGLmlrcaPHy9EAIkcf39/kb7y6d0V7du3x+qe/fHHqed4EfYpbTW7Y3nhgFzD3hS/HnqM/26/yvX/08bdCqOaOQvnYIJSY3mFioZpWvrGy8+lQzXbVrTBpFYuKGMm23GZ5nTRSAtyeKaaowYNGmDlypWig4z9bBhGcWCBwzBMoUR0bt++LdI5K1asgI+PD3yGDsWMasDjZCesPPdMpK08/ryIgfXsMLq5M7pUtc2TwPm2jr1U3OSVuMQkbL32Ukw/D49JFURU10OdUdkVI1OdDb1+8gIiJ+d58+bB3NwcP//8s/jJMIxiwQKHYZhCoUqVKli2bJmIeERHR+PBgwdo0bw5mjRpgs1LV+Ove1E48fgN1l98LrqXaDJ5XqDC4/x0Rh16EIQFx58i4F1qZ5Sjpb6YGUXdVtlFYA4cOCC6xVxdXREXF4fXr1+LUQv0fhiGUUxY4DAMU2iQgKDoB0HDOynVQ146Daq5i/lM/xs7B78c8RJpq7+vvsjTc1JXVV64+iwMc48+geeriLTHaYlOrW7VbMWAUFlcv35diBkaTxEUFAQjIyMxboHmc3E6imEUGxY4DMN8FQYPHizaqVetWiUGUEZGRsLdTBUeuA3Nxm2w5tILxCWmtUbJgGIt1kapLeM54f0mSnRGnXkaIu7raaphSKOyGNzAHrqask+BgYGB2LRpk+iKMjMzg7e3N7Zt24bOnTsLscMwjOLDAodhmK9GqVKlMHfuXFF0XLJkSUyfPl3UubRtew5T+w7CnrcW8AxM9aJJjySRNKNtuSz+OBKCI2Ox4uwT/Hc7QLgZq6uqoHet0qIt3Vw/e5FCYotqhkxNTUWrNw3DpChTr1695LIGDMPIBxY4DMN8derWrSt+0qRt6riKiYnBwG+64ocffkCDLmOx5ryfdJYUQQKlS9US+OVwVn+ccc0dcfKlKiYvvYTYhNTOqDbu1pjYygUOFvoy//8S757du3eLkQrU+m1rayuiODQ1nWEY5YMFDsMwCgPV4VCk5PfffxcDPD08PBAWeAWtIu9CrVp3HPP9gISkFIR+iMPaC8+zPJ7EzsTdD2mIBJUUo3oZE0z1cEO1MtnPn/L09BQt3xQ5ImHVsWNH8f8uXbo019kwjBLDAodhGIWCalxmz54tfGVMTExgb2+Pt2/fol07fwxo2Ax+ZrVw2vtdjs+hihQs61kZbSuXzFakUDdUeHi4mPIdHx8vzAirV68u0lN6erI9cBiGUR54+hvDMAoJpYhIaPz7779o166daCv/adJY6PqcQHJCzq3hyVCBsZ6GTHFD6ailS5eKCA0JHBJS3bt3x+rVqzFnzhwWNwxTROAIDsMwCk3z5s3FEEuqjSERUq9dL6z/phtUdQxg2uwHqBtayHxcSFRclm3Hjh3DzJkzRdEwDcWkDq4///wTampZ50wxDKPcyDWC89tvv4niQV1dXVG0lxfo6mrWrFkoUaKE8MigE9ujR4+yhJZp3gu5h9IVXocOHfDqVe6OpwzDKCfq6uoYPny4qJeJexeEuFePEet3G6GHF+P9+U1IjovJ8hjyuZFAbd5UXzNjxgzcuHFDnDtodtTixYtZ3DBMEUWuAofy2hT6HTp0aJ4fM3/+fHHSITv3mzdvwtraWgzti4qKku4zZswY0fGwfft2XLp0SVioUwg7KSlJTu+EYRhFQFVVFX3bNkLFUetgULUd4l4+QOT1PUgMD0ZcsC9SUig5BRhrpogC44iICNGGXr58eXFeIedhmo9F5w5KTXHkhmGKLnJNUVGhIEGtlnlBkhsnL4wuXbqIbZs3bxaW6JSH//HHH8UJixxQt2zZIkLXxNatW4WfxqlTp9CqVSs5viOGYb42aqoqmPedB4ZoWUHbtjziQvygoqmNoA3DoGlpD4v2k9ClliW2b/tXzI2qUKECEhMTRQ0PnVvq1Knztd8CwzDFrQaHpvQGBwcLd9P0HRWNGjXClStXhMChgX00LTj9PpTOcnd3F/vIEjiU0qKbBHJMJeh58jN5mElFsma8dgUDr2f+aeZijhXfVMavR7QRHFkLMd5XoaKmAXWVFLz/ZxxeJHaHlrOz6L6iSPLZs2eFYR8JHYaPz68J/71/Gfn53lEogUPihqCITXro/osXL6T7aGpqivbRzPtIHp8ZClFLoknpOXHihKgPYj6PkydP8tIVILye+WdyOeBZpAoinWoCjVfD59px7Ni+XRrNpbodGuZJkV8yEGT4+FQU+O/986BaOrkJHCoAliUW0kO1M+Qn8blkbu2k1FVug+1y2mfq1KkYN25chggOpbQoCmRoaPjZr7M4K2j646TaKA0Nja/9cpQeXs8CZGhvjBg+XKSxqUGBZkcxXwYfnwULr+eXIcnAyEXgjBgxIteZLHZ2dvgcqKCYoEiMjY2NdDu1c0qiOrQPhZzfv3+fIYpD+0js3jNDaS5ZA/Loy5m/oD8fXr+ChdezYKBJ36GhoeL8wH/fBQcfnwULr+fnkZ+/6XwLHGqvpJs8IMdSEjAUHahSpYrYRmKGht/NmzdP3K9WrZp4g7RPjx49xLagoCA8fPhQdGAxDMMwDMPItQbn5cuXePfunfhJLdz37t0T2x0dHaGvnzr0ztXVVdTIUCiZUkzUAk5zaGiKL93o31Qn07t3b7G/kZERBg8eLFo9zczMxMTfCRMmiE4JSVcVwzAMwzDFG7kKnJ9//lm0eUuQRGWoo4Hy44SXl5coAJRAbZ0fP37EsGHDRBqqVq1aohjYwMBAus+SJUuE8RdFcGjfZs2aiVZ09rRgGIZhGEbuAodER24eOFQcnB6K4lAhM92yQ1tbG8uXLxc3hmEYhmGYzPCwTYZhGIZhihwscBiGYRiGKXKwwGEYhmEYpsjBAodhGIZhmCIHCxyGYRiGYYocLHAYhmEYhilysMBhGIZhGKbIwQKHYRiGYZgiBwschmEYhmGKHHJ1MlZUJO7J+Rm7znwiISEBMTExYv14WvOXw+tZsPB68noqMnx8fhmS7+3MUxBkUSwFTlRUlPhZqlSpr/1SGIZhGIb5jO9xGr6dEyopeZFBRYzk5GS8fv1aDPCk2VdM/hU0icOAgAAYGhry8n0hvJ4FC68nr6ciw8fnl0GShcRNiRIloKqac5VNsYzg0KLY2tp+7Zeh9JC4YYHD66mo8PHJ66nI8PH5+eQWuZHARcYMwzAMwxQ5WOAwDMMwDFPkYIHD5BstLS3MnDlT/GS+HF7PgoXXk9dTkeHjs/AolkXGDMMwDMMUbTiCwzAMwzBMkYMFDsMwDMMwRQ4WOAzDMAzDFDlY4DAMwzAMU+RggcPkid9++w1169aFrq4ujI2N8/QYql+fNWuWcJzU0dFB48aN8ejRI15xAO/fv0e/fv2EYRXd6N/h4eE5rs2AAQOE83b6W+3atYvleq5atQr29vbQ1tZGtWrVcPHixRz3P3/+vNiP9ndwcMCaNWsK7bUWtfU8d+5cluOQbk+fPi3U16yoXLhwAe3btxfnPVqXffv25foYPj7lAwscJk/Ex8eje/fuGDp0aJ5XbP78+Vi8eDFWrFiBmzdvwtraGi1atJDOAivO9O7dG/fu3cOxY8fEjf5NIic3WrdujaCgIOntyJEjKG7s2LEDY8aMwfTp03H37l00aNAAbdq0wcuXL2Xu//z5c3h4eIj9aP9p06Zh1KhR2L17d6G/9qKwnhK8vLwyHItOTk6F9poVmejoaFSqVEmc9/ICH59yhNrEGSavbNy4McXIyCjX/ZKTk1Osra1T/vjjD+m22NhY8dg1a9YU6wV//PgxWTOkXLt2Tbrt6tWrYtvTp0+zfVz//v1TOnbsmFLcqVmzZsqQIUMybHN1dU2ZMmWKzP0nTZokfp+eH3/8MaV27dpyfZ1FdT3Pnj0rjtX3798X0itUXmid9u7dm+M+fHzKD47gMHKBrkqCg4PRsmXLDAZXjRo1wpUrV4r1ql+9elWkpWrVqiXdRqkm2pbb2lB6wNLSEs7Ozvj+++8REhKC4hZJvH37dobjiqD72a0drXfm/Vu1aoVbt24hISEBxZnPWU8JVapUgY2NDZo1a4azZ8/K+ZUWXfj4lB8scBi5QOKGsLKyyrCd7kt+V1yh908iJTO0Lae1obTBP//8gzNnzmDRokUi7de0aVPExcWhuBAaGoqkpKR8HVe0Xdb+iYmJ4vmKM5+zniRq1q1bJ1J8e/bsgYuLixA5VHvC5B8+PuVHsZwmzqRCBcCzZ8/OcTnoS7R69eqfvWRUZJceitpm3lbc1pOQtQa5rU3Pnj2l/3Z3dxefS5kyZXD48GF06dIFxYn8Hley9pe1vbiSn/UkQUM3CXXq1EFAQAAWLlyIhg0byv21FkX4+JQPLHCKMSNGjECvXr1y3MfOzu6znpsKiiVXJ3TFJ4FSKpmvFovbenp6euLNmzdZfvf27dt8rQ2tKwkcHx8fFBfMzc2hpqaWJbqQ03FFx6Ks/dXV1WFmZobizOespywoxbp161Y5vMKiDx+f8oMFTjE/udFNHlDLKf3hnjx5UuTqJfl+aoecN28eivN60hVvREQEbty4gZo1a4pt169fF9uoFT+vhIWFiSvn9AKyqKOpqSnamOm46ty5s3Q73e/YsWO2633w4MEM206cOCEiYBoaGijOfM56yoK6r4rTcViQ8PEpR+RYwMwUIV68eJFy9+7dlNmzZ6fo6+uLf9MtKipKuo+Li0vKnj17pPepg4q6pmjbgwcPUr755psUGxublMjIyJTiTuvWrVMqVqwouqfoVqFChZR27dpl2Cf9etI6jx8/PuXKlSspz58/F50sderUSSlZsmSxW8/t27enaGhopGzYsEF0pI0ZMyZFT08vxd/fX/yeun/69esn3d/Pzy9FV1c3ZezYsWJ/ehw9fteuXV/xXSjvei5ZskR0Bnl7e6c8fPhQ/J6+Snbv3v0V34XiQH+rkvMjrcvixYvFv+kcSvDxWXiwwGHyBLUo0x9r5ht90UoPJkC0kadvFZ85c6ZoF9fS0kpp2LChEDpMSkpYWFhKnz59UgwMDMSN/p257Tb9esbExKS0bNkyxcLCQnwZlS5dWnwmL1++LJbLuXLlypQyZcqkaGpqplStWjXl/Pnz0t/RujRq1CjD/ufOnUupUqWK2N/Ozi5l9erVX+FVF431nDdvXkrZsmVTtLW1U0xMTFLq16+fcvjw4a/0yhUPSRt95hutI8HHZ+GhQv+RZ4SIYRiGYRimsOE2cYZhGIZhihwscBiGYRiGKXKwwGEYhmEYpsjBAodhGIZhmCIHCxyGYRiGYYocLHAYhmEYhilysMBhGIZhGKbIwQKHYRiGYZgiBwschmEYhmGKHCxwGIZhGIYpcrDAYRiGYRimyMECh2EYhmEYFDX+D0AR1z2FVIUKAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "plt.plot(uNum.real, uNum.imag, 'o-')\n", "plt.axis(\"equal\")\n", "\n", "times = np.linspace(0, tEnd, nSteps+1)\n", "uExact = u0 * np.exp(lam*times)\n", "plt.plot(uExact[0].real, uExact[0].imag, 's', ms=10, c=\"orange\")\n", "plt.plot(uExact.real, uExact.imag, ':', c=\"k\")\n", "plt.grid()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, we can observe some numerical error of the time-scheme, since we do not retrieve exactly the hexagon showed by the exact solution.\n", "One accuracy indicator is the $L_\\infty$ error in time (computed over all time-steps) :" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "L_inf error : 0.11925\n" ] } ], "source": [ "print(\"L_inf error : {:1.5f}\".format(np.linalg.norm(uNum-uExact, ord=np.inf)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let say we want to use a collocation method instead of RK4, using $4$ Legendre Radau-II nodes (_i.e_ Radau nodes including the right time interval bound). The only line we have to change in the previous code is :" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# replace : \"rk = Q_GENERATORS[\"RK4\"]()\" by\n", "coll = Q_GENERATORS[\"coll\"](nNodes=4, nodeType=\"LEGENDRE\", quadType=\"RADAU-RIGHT\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the **exact same code as before** can be used to build the time-stepper, compute the numerical solution and the error. In practice, you don't have to, as the $Q$-generators in `qmat` already provides dedicated methods for that :" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "L_inf error : 0.00001\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT3RJREFUeJzt3QlYVGXbB/C/7KKACwqYKLjiLoIiGO6gkOab5pJlWWppman19opLbilpalYuZWmmZvqVmpWo4I6Cu+K+o7gh7qgk63zX89gMDIOsMwyc+f+u6zjMPWfGMw+HM/c8axmVSqUCERERkYKYGfsAiIiIiPSNCQ4REREpDhMcIiIiUhwmOERERKQ4THCIiIhIcZjgEBERkeIwwSEiIiLFYYJDREREimMBE5SRkYEbN27Azs4OZcqUMfbhEBERUT6IuYkfPXqEatWqwcws9zoak0xwRHLj6upq7MMgIiKiQrh69SqqV6+e6z4mmeCImht1Adnb2xv7cEqd1NRUhIeHIzAwEJaWlsY+nFKP5cnyLMl4frI8S5LExERZQaH+HM+NSSY46mYpkdwwwSncBc/W1laWHROcomN56hfLk+VZkvH81I/8dC9hJ2MiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiUoyUtAz8FHUZv8eayVtxn4hMk0ETnF27dqF79+5yUSwx6+Aff/yR53N27twJLy8v2NjYoFatWvjuu+909lmzZg0aNmwIa2trebtu3ToDvQMiKi1Cw07BY8JGTN94DpHxZvJW3BdxIjI9Bk1wnjx5gmbNmmHevHn52j82NhbBwcHw9/fHkSNHMHbsWIwYMUImNGrR0dHo27cvBgwYgJiYGHnbp08f7Nu3z4DvhIhKMpHEfL8rFhkq7bi4L+JMcohMj0HXogoKCpJbfonamho1amDu3LnyfoMGDXDw4EHMmjULvXr1kjHxWEBAAEJCQuR9cStqfUT8119/NdA7IaKSSjRD/RAZm+s+4vGPAz1gZcFWeSJTUaIW2xS1M2KF6qy6dOmCxYsXywXKxMKOYp9Ro0bp7KNOinKSnJwst6yrkQriNcVGBaMuM5adfrA8i2Zp1GWdmpvsxONL91zE235uRfzfTA/PT5ZnSVKQz50SleDEx8fDyclJKybup6Wl4c6dO3BxcXnuPiL+PKGhoZg8ebJOPDw8XK6KTYUTERHBotMjlmfhRMaKWpm8a2Yij5yB0wP2xyksnp/6xfIsnKSkpNKZ4OS0BLpKpdKJ57RPbkuni2as0aNHa9XguLq6ytoie3t7PR696WTQ4o9TNBWKWjVieRrTpe0XERl/Mc/9/Jt7ILgNa3AKin/v+sXyLBp1C0ypS3CcnZ11amISEhJgYWGBypUr57pP9lqdrMRoK7FlJz6c+QFdeCw//WJ5FkxGhgqDlx3EtjMJ+dp/+qZzqGxXFr28qhfq92PqeH6yPEuCgnxml6ged76+vjrVdqIZydvbW/OmnrePn59fsR4rERnP8r1XUGtsWL6TG7WPf4uB25gNOH/rkcGOjYhKBoMmOI8fP8bRo0flph4GLn6Oi4vTNB29+eabmv2HDh2KK1euyOak06dPY8mSJbKD8SeffKLZ56OPPpIJzYwZM3DmzBl5u2XLFowcOdKQb4WISoAT1x/KBGXCHyc0MVsrc8RMDMR7bd1hlq2lWtwX8d+G+mrFA77aBb/QrUhKSSuuQyeiYmbQJioxxLtDhw6a++p+MG+99RaWLl2KmzdvapIdwd3dHWFhYXKU1Pz58+UEgd98841miLggampWrVqF8ePHY8KECahduzZWr14NHx8fQ74VIjKiR09T4T9zOx4kaY+g+HN4GzStXkH+HBLcUA4FF6OlRIdif08PDGxTWzM0/PIXL2H+9gv4cvNZef/Gw6do+NlmvOVbE5N7NDbCuyKiUpvgtG/fXtNJOCciycmuXbt2OHz4cK6v++qrr8qNiJRNXD/++/sx/H7omlZ8YveGeLuNu87+IpkRQ8HFaKlgPzdYZpv35oMOdfBu21rot2gvDl25L2M/R1+R2w9veiOg4fP78hFR6VKiOhkTEan9FXMDH/56RKtA/Os6YunbrWCevS2qACzNzbBmmB9uPvwHvqHbNPEhyw7K293/64DqFTl9BFFpxwSHiEqU2DtP0GHWDp34/nGdUNXORm//j4tDWdlstf1MAt5eekATf3HGdjR+wR5rh7XhzMdEpViJGkVFRKbraWo6Os7aoZPcrBzsIxMRfSY3WXXwqCpff/CLmU1eJ64not74jfh6y3mD/J9EZHhMcIjI6EI3nobHhE24dOeJJjaiU12ZePjVcSyWYxjfrSHOTO0Kd8dymthXW87JUVt7L90tlmMgIv1hExURGc3Oc7fx1pL9WrFG1eyx7n3jNA/ZWJpj+yftcfH2Y3SavVMTF52ShQPjOqOKne6koURU8jDBIaJil72Dr1rkpx3gWsn4HXxrVykva4/WH72Oj1Y9m8dLaDltC9rVq4IlA1sWqaMzERkem6iIqNikpWeg93dROsnNd294yYSiJCQ3WfVo/gJiQ4PRq0V1rVqn2mPDsCz6slGPjYhyxwSHiIrFdzsvos64jThw+dn8M8KA1jVlAtG1sXOJ/S2IhXxn92mGY5MCYWeTWen92fqTsn/O8WsPjXp8RJQzNlERkUEdunIPvRZGa8Wc7W2w7ZN2sLUqPZcgextLHJ/URSY03eft1sTFz3bWFtgT0lHuQ0QlQ+m5uhBRqXL/SQq8Po9ARrbJzDePbIv6znYorZpUd5DNaaKJStTiCI+S09B0Ujh6tngBs3s3k7U+RGRcbKIiIr3KyFDhveUH4TlVO7mZ+WpTmRiU5uQmqzd93XBxejDa16+iia09fB3uIWGyczIRGRcTHCLSm1X741BrbBg2n7ylib3UxEX2s+nj7aq4khYjqcTSEWL4eFZi5JXonyOGmxORcbCJioiK7Ex8IrrOjdSKWZmbyeUVKthaKb6Exdw4onZKTAionjNHEHPp1KxsK5vlxBw7RFR8WINDRIX2JDkN3p9H6CQ3a9/3w7lpQSaR3GTVulZlmeiM6lxPE7tyN0nO0vz536eMemxEpoYJDhEVmEqlQsja42g0cTPuPE7RxMcFN5Af8C1qVDTpUv2oc12c+zwITV5w0MR+3B0rm63E4p5EZHhsoiKiAtl4/CaG/XJYK+ZbqzKWD2oFC3N+Z1ITS0389eGLuHY/Sa5QrqZeuTxqTEdUq1CWZx+RgTDBIaJ8ibubhLZfZn5Qq+0b2wlO9oZZ6VsJqle0lbVaEaduYciyg5q43xfb4FmjAv7vPV9YMjEk0jt+3SKiXCWnpaPLV7t0kptl77SSH9xMbvInoKGTLK+3fGtqYkfiHqDuuI2Yv/0Cz0IiPWOCQ0TPNWvzWdQfvwlnbz3SxIa1ry0/qNvWy5z/hfJvco/GOD2lK17I0jz15eazsn/Ogcv3WJREesImKiLSsfv8HbyxeJ9WrJ5Tefw5/EUOd9aDslbm2DOmI87feoSAr3Zp4r2/e7akxeEJAahUzrRGoBHpGxMcItJISHyKVtO36pTIjk/aw82xHEtKz+o62cnasDWHruHj32I08RZTI9C5gRMWDfCCmRmXfSAqDDZRERHSM1To/8NeneRmfv8W8gOYyY1h9fKqLmd7frlZNU1sy+lbclbolfvieIYSFQITHCIT92PkJdQeG4aoi3c1sX4tXeUH7ktNXYx6bKZELND5zWueiPksEDaWmZfmseuOy/45p24kGvX4iEobNlERmaijVx/gP/P3aMUql7PCzk87oLw1Lw3G4mBriTNTgxBz9QF6ZPn9BH8TKfvl7OLvhyhfeBUjMjEPk1LRcvoWpKRlaMXDRvijYTV7ox0XaWvmWkE2Dy7eHYup/y7zcO9JChpP3Iw+3tUxo1dTWetDRDljExWRCS2v8MHKw2g2JVwruZn+ShP5QcrkpmQa9KI7Lk4PRps6lTWx/zt4De4hYfj72A2jHhtRScYEh8gE/H7o2QfihmM3tSaeuzQ9GP19ahj12Chv5mZl8Mvg1tg/tpNWfPjKI7J/zuU7T1iMRNmwiYpIwbLPs6LGeVZKp6r2NrK2bc+FO3j9x8x5itrP2oE6Vcvj7w85TxGRGmtwiBQoKSUNvqFbdZKb34f6yg9ITiJXurWp4yh/jx92rKOJXUh4DI8JmzBj0xmjHhtRScEEh0hhJq4/gYafbcbNh081sU+71pcfiN5ulYx6bKRfHwfWx9nPu8LD2U4TW7jjomy22nXuNoubTBqbqIgUIvxkPN5dfkgr5l2zIn59tzVXq1YwawtzbBrZVme19zeX7Je3XO2dTFWx1OAsWLAA7u7usLGxgZeXFyIjI5+778CBA+XQx+xbo0aNNPssXbo0x32ePs38xkpkKq7eS5Lf2LMnN9EhHfH7MD8mNyaiRmVbWUu38PUWWnGf6VvR5/topKVrTwtApHQGT3BWr16NkSNHYty4cThy5Aj8/f0RFBSEuLicpx//+uuvcfPmTc129epVVKpUCb1799baz97eXms/sYkEishUiKHe3b6NhP/MzG/twk8DW8oPOheHzNWqyXQENXGRs1C/nmV03P7Ye6gzbiN+2HXJqMdGpKgEZ86cORg0aBAGDx6MBg0aYO7cuXB1dcXChQtz3N/BwQHOzs6a7eDBg7h//z7efvttrf1EjU3W/cRGZCq+3nIe9cZvxInrmdP3D/F3l4lNB4+qRj02Mj5xfZz2ShOcnNwFVe2sNfFpYadlbd/huPtGPT6iUt8HJyUlBYcOHcKYMWO04oGBgYiKisrXayxevBidO3dGzZo1teKPHz+WsfT0dDRv3hxTp06Fp6dnjq+RnJwsN7XExGcfCqmpqXKjglGXGcuu+MtzX+w9vLHkoFbMrbIt/vrAFzaW5vyd8PzUYmUG7Pm0Hc7GP0K3+dGaeM8FUbAwK4Po/7VHBVtLvZ2flDeWZ9EU5DwsoxLTmxrIjRs38MILL2DPnj3w8/PTxKdPn46ff/4ZZ8+ezfX5otlJ1PasXLkSffr00cT37t2LCxcuoEmTJjJZEc1aYWFhiImJQd26dXVeZ9KkSZg8ebJOXLyura1tkd8nkaElpgATDul+HxnbPA1ObImifIq+VQarLplrxZpXysDAehngqg9UGiQlJaF///54+PCh7Kpi9ARH1Nb4+vpq4tOmTcPy5ctx5kzu8zWEhoZi9uzZ8nWsrKyeu19GRgZatGiBtm3b4ptvvslXDY5InO7cuZNnAVHOGXRERAQCAgJgaZn7tz8qWnmmZ6gwZPlhRF7IXOlbmNO7Cbpzpe8ClyeJ66UKI1bHYPOpBK3imP6fRujt9QLL08B4fhaN+Px2dHTMV4Jj0CYqcRDm5uaIj4/XiickJMDJySnX54q8a8mSJRgwYECuyY1gZmaGli1b4vz58zk+bm1tLbfsxMWPF8DCY/kZtjx/jrqMiX+e1NqnV4vqmNWbiyzy/Cya799siftPUuA9bYtMooWxf5yU2+aRbVE/y7w6zzs/qWhYnoVTkHPQoJ2MRWIihoWLb1NZiftZm6xysnPnTtkMJToo50UkQ0ePHoWLi0uRj5nI2I5feyg7gmZNbuxsLHBsUiBm92nGFaRJLyqWs5KLeK4Zpn0t7jJ3F3ymb5GzYROVZgaf6G/06NGyFsbb21s2Uy1atEgOER86dKh8PCQkBNevX8eyZct0Ohf7+PigcePGOq8p+tO0bt1a9rcR1VWiWUokOPPnzzf02yEymMR/UtHu8214nKz9wSLWF2r8ggNLngzCq2ZFOfru+50XEbrxWbeBW4nJcjZsMdR84kv1WfJUKhk8wenbty/u3r2LKVOmyE7DImERHYLVo6JELPucOKJtbc2aNbLzcE4ePHiAd999VzZ9iWHlYvTUrl270KpVK0O/HSK9EzWQy8+b4aNo7flspvZohAG+bixxKhbvtauNQS+6o/8P+7D/8j0Z+2VfnNzeqVcGwfw9UClj0E7GJZWo9RGJUX46KVHOneREkhocHMw2+SL648h1jFx9VCvWoX4VLH6rJczMyvD0KwSen0UX//ApWodu1YlHftoBrpU48rQoeH4W3+c316IiMoKLtx+j0+ydOvGD4zvDsbxuh3ii4uTsYCObrXaeu423/l3TShCzZjdwsccfH/jJNbCISjKuJk5UjJ6mpqPtzO06yc2HDdNwfmogkxsqUdrVqyLPy07VMtexOn0zEfXHb8Kc8NznMSMyNiY4RMVk6t+n4DFhE+LuJWliowPqyQ+QOuxDTCXYyzUzcOKzTqhVpZwm9s22C3K0X9SFO0Y9NqLnYRMVkYFtO3ML7yzVXl6hWXUH/DbUD1YWZpwCn0oFa0tzbPu4PWLvPEGHWTs08f4/7pO3+8d1QlU7LnhMJQcTHCIDufHgH/h9sU0nvvt/HVC9IjtqUunk7lhO9s/5K+YGPvz1iCbeatpWvFjHET+/0wrm7CBPJQCbqIj0LDU9A/+Zv0cnufnhTW/5wcDkhpSge7NqiA0NRh/v6prY7gt3UHtsGH7aE2vUYyMSmOAQ6dH87RdQd9xGHL36QBMb6OcmE5uAhrkvT0JU2pQpUwYzX22G45MCUTHLquST/zol++fEZPk7ICpubKIi0oP9sffQ5/tordgLFcpiy+h2KGvF4bSkbHY2ljjyWSBO3niIl77ZrYn3mL8HZS3NsXdsJziU5TpWVLyY4BAVwd3HyfD6fItOfMvotqhTVXfBQiIla1TNQdZWrth7BeP/OCFj/6Smo9nkcLzcrBq+7teca6lRsWETFVEhZGSoMPjnAzrJzezezeQFnskNmbI3WtfEpenB6Nygqib2Z8wNuIeEYe3ha0Y9NjIdTHCICuiXfVdQa2wYtpxO0MRe/rfDZS+vzA6XRKZMLDXy41stcWh8Z6346P+Lkf1zLiQ8MtqxkWlgExVRPmXvXyDI/gUhneCQpYMlEWWqXN5a1mpm76fWec4uVK9YFhGj2E+NDIM1OER5ePQ0Fc2nhOskN+s/aIPTU7syuSHKh1bulWSi898u9TWxa/f/QYPPNmHSnydZhqR3THCInkOlUuG/v8WgyaRwPEhK1cQ/69ZQXqibuVZg2REV0Acd6uD8tCA0z/L3szTqsmy22nLqFsuT9IZNVEQ5+PvYDQxfmTlLq8BZWon0w9LcDH980EZntu/By54tabJnTEc5zQJRUTDBIcoi+zo7alxnh0j/qlUoK2tDt56+hUE/Z67X1uaLbWha3QFrhvnJZIioMHjmEAF4mpqOjrN36CQ3Kwf7yAswFxEkMpxODZzk39k7bdw1sWPXHspZwb/dep5FT4XCBIdM3hcbz8BjwiZcuv1EUxYjOtaRF1y/Oo4mXz5ExeWz7g1xZmpX1KycuRjt7Ihzsn/O3kt3+YugAmETFZmsnedu460l+7ViDVzs8ccHfrC24PIKRMZgY2mOnf/tgAsJj9F5zk5NvN+ivfL24PjOcCxvzV8O5YkJDpmc+IdP0Tp0q0488tMOcK2U+c2RiIynTtXyshZ13ZFrGLU6RhP3/nwLOtSvgsVvtZSTCRI9D5uoyGSkpWegz3fROsnNd2+0kBdSJjdEJc8rntXlLOH/aV5NE9t+9racTXz53itGPTYq2ZjgkEn4fudF1Bm3Efsv39PE3mhdQ144uzZ2MeqxEVHuypQpg7n9PBEzMRDlrTMbHib8cUL2zzlx/SGLkHSwiYoU7dCV++i1MEorVtXOGts/aY9yWS6URFTyOZS1xInJXXDs2gO8PG+PJt7t292wt7GQ8+fY2XDZFHqGV3hSpPtPUuA9bQvSM1Ra8c0j26K+s53RjouIiq5p9QqyWXnpnlhM+uuUjCU+TZOzjvdqUR2zejeVtT5k2thERYqSkaHC0OWH4Dk1Qiu5mdmrqbwgMrkhUo6BbdxxcXow2taroomtOXwN7iFh+DPmhlGPjYyPCQ4pxuoDcbLj4aaT8ZpYUGNnXJoejD4tXY16bERkGOZmZbDsnVY4MK6zVnzEr0dk/5xLtx+z6E0Um6io1DsTn4iucyO1YhZmZeQFr2I5K6MdFxEVnyp21rKWNuriHfT/YZ8m3nH2Trg7lsPGj/zlHDtkOpjgUKn1JDkN7WftwO1HyVpxsX6NV82KRjsuIjIev9qOMtH5KuIcvv53mQexxpyYrfy9trUQEtxA+wlP4oDkOwX/j6wdgXI19HTUZAhMcKjUUalUGPfHCazcF6cVHxvsgXfb1jbacRFRyTEqoB4+6FAHPRfuwYnriTL2/a5Lcvvp7ZboUL/qs+Tmr/pAxtOC/wdmNkD3s0xySjAmOFSqbDpxE0NXHNaKtXKvJBfFtOCqw0SUhZWFGf7+0B9X7yXBf+Z2Tfztnw7I2+hhVeFSmORGEM8TNT+sxSmxmOBQiZGSloHl0Zdx5V4SalayxQBfN3mBEuLuJqHtl5kXKLW9IZ3g7GBjhKMlotJCzFIumq02n4zHe8sPaeK+CxPgZTsTq2uPgUWZDE38n3QrTL/5Ni4nV4Ob9Q2MdfkJZc1TjHT0VKJHUS1YsADu7u6wsbGBl5cXIiO1O4RmtWPHDjl/QfbtzJkzWvutWbMGDRs2hLW1tbxdt25dMbwTMpTQsFPwmLARUzecxrLoK/JW3v/7JLrO3aWT3IhRE+KCxeSGiPKrSyNned1407emJnYoqSHqHP8TCxN6yftDLo9Fg5NrsPxed0Q+8ZK34r6IU+li8ARn9erVGDlyJMaNG4cjR47A398fQUFBiIvT7j+R3dmzZ3Hz5k3NVrduXc1j0dHR6Nu3LwYMGICYmBh526dPH+zbl9lznkpXcvP9rlhkm5NP3l+8+zLOxD/SxIa2qy0vUFnnvSAiKogpPRrj1JQuqGaX+RE4I/5tuB37GxGJvjk+R8SZ5JQuBk9w5syZg0GDBmHw4MFo0KAB5s6dC1dXVyxcuDDX51WtWhXOzs6azdw8c3ifeI2AgACEhITAw8ND3nbq1EnGqfQ1S/0QGZvnfrUdbXFmaleMCfIoluMiImWztbJA1DAnhNd7P9sjYgbk7LMgl9EkOaL5ikoHg/bBSUlJwaFDhzBmzBiteGBgIKKitNcHys7T0xNPnz6VzU/jx49Hhw4dtGpwRo0apbV/ly5dnpvgJCcny00tMfFZj/rU1FS5UcGoy0wfZbc06rJOzU1O+rasDnNkIDU1s51cKfRZnsTy1DdFn59paahnE4fLTbvhtQufIzqpeS47P0tyRN+cqdW/lz+npqWJginQf6no8iwGBSk3gyY4d+7cQXp6OpycnLTi4n58fOZss1m5uLhg0aJFsq+OSEqWL18ua2dE35y2bdvKfcRzC/KaoaGhmDx5sk48PDwctra2RXiHpi0iIqLIrxEZa5avisTII2fg9ODZmjNKpY/yJJanoSjx/HRIv4j2//5sUSY9X88RHY/V9uzejYfmNwv1fyuxPItDUlJSyRpFlX3RMzGPyfMWQqtfv77c1Hx9fXH16lXMmjVLk+AU9DVFE9bo0aO1anBEM5moSbK3ty/0+zLlDFr8cYpmQkvLoq3ceyvqMiI3nstzP39PDwT7uUGJ9FmexPLUN0Wfn/ePAFue/ShGS4lOxXkR+6m1efFFoKJngf5LRZdnMVC3wBg9wXF0dJR9Z7LXrCQkJOjUwOSmdevWWLFihea+6JNTkNcUI63Elp04uXiCFZ4+ym9gm9r4YtO5XJupzMo828/y3yHjSsXzkeVZkiny/LTI/AgUQ8GX3+v2772cviyrNPupWYrnF7JMFFmexaAgZWbQTwwrKyvZ1JS9Kk7c9/Pzy/friNFXoukqa61O9tcUzU0FeU0qGcQ8N0P83XPdRzyung+HiMgQxDw3AfbR/97L/o3r2X3xOOfDKT0M3kQlmobEMG5vb2+ZmIj+NWKI+NChQzXNR9evX8eyZcvkfdFR2M3NDY0aNZKdlEXNjZjzRmxqH330kWyumjFjBnr06IH169djy5Yt2L17t6HfDhlASHBDebtoV6zWZUXU3IjkRv04EZEh/eA2XQ4Fz2mouEhuxONUehg8wRHz1dy9exdTpkyR89k0btwYYWFhqFnz2URLIpZ1ThyR1HzyyScy6SlbtqxMdDZs2IDg4GDNPqKmZtWqVXJ01YQJE1C7dm05346Pj4+h3w4ZiEhi3BzLIWTtCXl/wksNtGYyJiIqDiKJWZLQDVPin30JH1DpL85kXEoVSyfj999/X245Wbp0qdb9Tz/9VG55efXVV+VGymFhlpnMDPKvZdRjISLTZWWeOaJKPSScSh9+PSYiItNk7fhsVfDCEM8Tz6cSi4ttEhGRaRIrgXc/+2xV8KyOPAGu/zscuWvm4pxaRHLDlcRLNCY4RERkukSSkj1RKXcFwLP+gKjUwiiHRUXHJioiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4hRLgrNgwQK4u7vDxsYGXl5eiIyMfO6+a9euRUBAAKpUqQJ7e3v4+vpi8+bNWvssXboUZcqU0dmePn1aDO+GiIiIYOoJzurVqzFy5EiMGzcOR44cgb+/P4KCghAXF5fj/rt27ZIJTlhYGA4dOoQOHTqge/fu8rlZieTn5s2bWptIoIiIiIgsDF0Ec+bMwaBBgzB48GB5f+7cubJGZuHChQgNDdXZXzye1fTp07F+/Xr89ddf8PT01MRFjY2zszN/g0RERFS8CU5KSoqshRkzZoxWPDAwEFFRUfl6jYyMDDx69AiVKlXSij9+/Bg1a9ZEeno6mjdvjqlTp2olQFklJyfLTS0xMVHepqamyo0KRl1m+i478bvM/n+YAkOVp6liebI8S/K1iOdn0RTk92HQBOfOnTvyRHFyctKKi/vx8fH5eo3Zs2fjyZMn6NOnjybm4eEh++E0adJEJitff/012rRpg5iYGNStW1fnNURN0eTJk3Xi4eHhsLW1LdR7IyAiIkKvxXAsoQwAc/mzaKI0NfouT1PH8mR5FtaJeMNfi3h+Fk5SUlLJaaJSNydlpVKpdGI5+fXXXzFp0iTZRFW1alVNvHXr1nJTE8lNixYt8O233+Kbb77ReZ2QkBCMHj1ac18kRa6urrImSfTloYJn0OKPU/SVsrS01Fvx/XP4OlZePCl/Dg4ONplfi6HK01SxPFmeRfVg/1X8FnvaINcinp9Fo26BMXqC4+joCHNzc53amoSEBJ1anZw6J4u+O7/99hs6d+6c675mZmZo2bIlzp8/n+Pj1tbWcstOfJjwA6Xw9F1+4lzJ+tqmhucjy7MkM6XzsziuRaZUnvpUkDIz6CgqKysrOSw8e1WcuO/n55drzc3AgQOxcuVKvPTSS3n+P6JG6OjRo3BxcdHLcRMREVHpZvAmKtE0NGDAAHh7e8s5bRYtWiSHiA8dOlTTfHT9+nUsW7ZMk9y8+eabsl+NaIZS1/6ULVsWDg4O8mfRn0Y8JvrbiOoq0SwlEpz58+cb+u0QERFRKWDwBKdv3764e/cupkyZIueqady4sey0JUZACSKWdU6c77//Hmlpafjggw/kpvbWW2/JjsXCgwcP8O6778rkRyQ9YvSUmD+nVatWhn47REREVAoUSyfj999/X245USctajt27Mjz9b766iu5EREREeWEa1ERERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDjFkuAsWLAA7u7usLGxgZeXFyIjI3Pdf+fOnXI/sX+tWrXw3Xff6eyzZs0aNGzYENbW1vJ23bp1BnwHREREVJoYPMFZvXo1Ro4ciXHjxuHIkSPw9/dHUFAQ4uLictw/NjYWwcHBcj+x/9ixYzFixAiZ0KhFR0ejb9++GDBgAGJiYuRtnz59sG/fPkO/HSIiIioFDJ7gzJkzB4MGDcLgwYPRoEEDzJ07F66urli4cGGO+4vamho1asj9xP7iee+88w5mzZql2Uc8FhAQgJCQEHh4eMjbTp06yTgRERGRhSGLICUlBYcOHcKYMWO04oGBgYiKisrxOaJ2RjyeVZcuXbB48WKkpqbC0tJS7jNq1CidfZ6X4CQnJ8tNLTExUd6K1xMbFYy6zPRddunp6Tr/hykwVHmaKpYny7MkX4t4fhZNQX4fBk1w7ty5I08UJycnrbi4Hx8fn+NzRDyn/dPS0uTrubi4PHef571maGgoJk+erBMPDw+Hra1tId4ZCREREXotiGMJZQCYy5/DwsJMrpD1XZ6mjuXJ8iysE/GGvxbx/CycpKSkkpHgqJUpI06WTCqVSieW1/7Z4wV5TdGENXr0aK0aHNFMJmqK7O3tC/huSGTQ4o9TNBOKGjV9+efwday8eFL+LPphmQpDlaepYnmyPIvqwf6r+C32tEGuRTw/i0bdAmP0BMfR0RHm5uY6NSsJCQk6NTBqzs7OOe5vYWGBypUr57rP815TjLQSW3biw4QfKIWn7/IT50rW1zY1PB9ZniWZKZ2fxXEtMqXy1KeClJlBOxlbWVnJ4d7Zq+LEfT8/vxyf4+vrq7O/aEry9vbWvLHn7fO81yQiIiLTYvAmKtE0JIZxiwRFJCaLFi2SQ8SHDh2qaT66fv06li1bJu+L+Lx58+TzhgwZIjsUiw7Gv/76q+Y1P/roI7Rt2xYzZsxAjx49sH79emzZsgW7d+829NshIiKiUsDgCY6Yr+bu3buYMmUKbt68icaNG8tOWzVr1pSPi1jWOXHEhIDicTFKav78+ahWrRq++eYb9OrVS7OPqKlZtWoVxo8fjwkTJqB27dpyvh0fHx9Dvx0iIiIqBYqlk/H7778vt5wsXbpUJ9auXTscPnw419d89dVX5UZERESUHdeiIiIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxmOAQERGR4jDBISIiIsVhgkNERESKwwSHiIiIFIcJDhERESkOExwiIiJSHCY4REREpDhMcIiIiEhxDJrg3L9/HwMGDICDg4PcxM8PHjx47v6pqan43//+hyZNmqBcuXKoVq0a3nzzTdy4cUNrv/bt26NMmTJaW79+/Qz5VoiIiKgUMWiC079/fxw9ehSbNm2Sm/hZJDnPk5SUhMOHD2PChAnydu3atTh37hxefvllnX2HDBmCmzdvarbvv//ekG+FiIiIShELQ73w6dOnZVKzd+9e+Pj4yNgPP/wAX19fnD17FvXr19d5jqjliYiI0Ip9++23aNWqFeLi4lCjRg1N3NbWFs7OzoY6fCIiIirFDJbgREdHy4RFndwIrVu3lrGoqKgcE5ycPHz4UDZBVahQQSv+yy+/YMWKFXByckJQUBAmTpwIOzu7HF8jOTlZbmqJiYmaJjGxUcGoy0zfZZeenq7zf5gCQ5WnqWJ5sjxL8rWI52fRFOT3YbAEJz4+HlWrVtWJi5h4LD+ePn2KMWPGyKYue3t7Tfz111+Hu7u7rME5ceIEQkJCEBMTo1P7oxYaGorJkyfrxMPDw2VNEBXO88q7sI4llAFgLn8OCwuDqdF3eZo6lifLs7BOxBv+WsTzs3BEVxaDJTiTJk3KMVnI6sCBA/JW1Lxkp1KpcoznlKWJjsMZGRlYsGCBTv8btcaNG6Nu3brw9vaW/XZatGih81oiARo9erRWDY6rqysCAwO1EifKH/G7EX+cAQEBsLS01Fux/XP4OlZePCl/Dg4ONplfh6HK01SxPFmeRfVg/1X8FnvaINcinp9Fo26BMUiCM3z48DxHLLm5ueHYsWO4deuWzmO3b9+WzUp5nQB9+vRBbGwstm3blmcSIpIa8cFw/vz5HBMca2truWUnnsMPlMLTd/mZm5trvbap4fnI8izJTOn8LI5rkSmVpz4VpMwKnOA4OjrKLS+iM7HoP7N//37ZSVjYt2+fjPn5+eWZ3IhkZfv27ahcuXKe/9fJkyfl81xcXAr4boiIiEiJDDZMvEGDBujatatsThIjqcQmfu7WrZtWB2MPDw+sW7dO/pyWloZXX30VBw8elJ2IRUcv0V9HbCkpKXKfixcvYsqUKXKfy5cvy/bR3r17w9PTE23atDHU2yEiIqJSxKDz4IgkRUzaJ/q6iK1p06ZYvny51j5iyLio1RGuXbuGP//8U942b95c1sioNzHySrCyssLWrVvRpUsXmSiNGDFCvvaWLVu0qhWJiIjIdBlsFJVQqVIlOZQ7N6LTcda+O1nv50R0Dt65c6fejpGIiIiUh2tRERERkeIwwSEiIiLFYYJDREREisMEh4iIiBSHCQ4REREpDhMcIiIiUhwmOERERKQ4THCIiIhIcZjgEBERkeIwwSEiIiLFYYJDREREisMEh4iIiBSHCQ4REREpDhMcIiIiUhwmOERERKQ4THCIiIhIcZjgEBERkeIwwSEiIiLFYYJDREREisMEh4iIiBSHCQ4REREpDhMcIiIiUhwmOERERKQ4THCIiIhIcZjgEBERkeIwwSEiIiLFYYJDREREisMEh4iIiBSHCQ4REREpDhMcIiIiUhwmOFRipGVkaH5eHHkJKWmZ94mIeC2iEpPg3L9/HwMGDICDg4PcxM8PHjzI9TkDBw5EmTJltLbWrVtr7ZOcnIwPP/wQjo6OKFeuHF5++WVcu3bNkG+FDCw07BTGrj2huT91w2l4TNgo40RExUVccyb/mXnd4bWo9DJogtO/f38cPXoUmzZtkpv4WSQ5eenatStu3ryp2cLCwrQeHzlyJNatW4dVq1Zh9+7dePz4Mbp164b09HQDvhsy5AXl+12xUGWLZ6gg40xyiKg48FqkLBaGeuHTp0/LpGbv3r3w8fGRsR9++AG+vr44e/Ys6tev/9znWltbw9nZOcfHHj58iMWLF2P58uXo3LmzjK1YsQKurq7YsmULunTpYqB3RIYgmqF+iIzNdR/x+MeBHrCyYIsqERkGr0XKY7AEJzo6WjZLqZMbQTQ1iVhUVFSuCc6OHTtQtWpVVKhQAe3atcO0adPkfeHQoUNITU1FYGCgZv9q1aqhcePG8nVzSnBEk5bY1BITE+WteB2xUcGoy0wfZbc06rKsqcmNeHzpnot4288NSqTP8iSWp76ZyvlZXNciUylPQylIuRkswYmPj9ckJVmJmHjseYKCgtC7d2/UrFkTsbGxmDBhAjp27CgTG1GzI55rZWWFihUraj3Pycnpua8bGhqKyZMn68TDw8Nha2tbqPdHQERERJGLITJW1MrkXTMTeeQMnB4ouz+OPsqTWJ6GovTzs7ivRUovT0NJSkoyXIIzadKkHJOFrA4cOCBvRQfh7FQqVY5xtb59+2p+FrUy3t7eMtnZsGEDevbs+dzn5fa6ISEhGD16tFYNjmjSErVA9vb2ub4XyjmDFn+cAQEBsLS0LFIR3Yq6jMiN5/Lcz9/TA8EKrsHRV3kSy1PfTOX8/GPFYSD+jsGvRaZSnoaiboExSIIzfPhw9OvXL9d93NzccOzYMdy6dUvnsdu3b8valvxycXGRCc758+flfdE3JyUlRY7QylqLk5CQAD8/vxxfQ9T8iC07cXLxBCs8fZTfwDa18cWmc3lWDX+99QL6+bjB3ka5FwSejyzPkkyp5+ehK/fRa2FUvvY1K/PsmmWph/6ASi1PQytImRU4wRFDs8WWF9GZWHQI3r9/P1q1aiVj+/btk7HnJSI5uXv3Lq5evSoTHcHLy0u+QZEB9+nTR8bESKsTJ05g5syZBX07ZGSi4/AQf3c5Wio3T1Iy0HRSOHq2eAGzezfLtRaQiCgv95+kwHvaFqTn9e0qC3Gt4mCH0sNgw1IaNGggh3sPGTJEjqQSm/hZDOfO2sHYw8NDDvkWxHDvTz75RHZQvnz5suxs3L17d5lQvfLKK3If0Ul50KBB+Pjjj7F161YcOXIEb7zxBpo0aaIZVUWlS0hwQ7zX1l1+O8pK3BcXlLb1qmhiaw9fh3tIGNYfvV78B0pEpV5GhgrvLT8Iz6kRWsnNzF5NcfmLl557LRJxca2i0sNgnYyFX375BSNGjNCMeBIT8s2bN09rHzFkXNTqCObm5jh+/DiWLVsmJwQUtTYdOnTA6tWrYWdnp3nOV199BQsLC1mD888//6BTp05YunSpfD6VTuLCIYaCL4++jCv3klCzki0G+Lppvi3dfpSMltO2aPb/aNVRuW37uB1qVSlvxCMnotJi9YE4/G/Nca1YUGNnzO/fAmb/ZjV5XYuo9DBoglOpUiU5R01uROdgtbJly2Lz5s15vq6NjQ2+/fZbuZFyiAvIIP9aOT5Wxc5afruKunAH/X/cp4l3nL0TbpVtsWlkW9hYMsElIl1n4hPRdW6kVszCrAwOjOuMiuWsCnQtotKDKSmVKn51HGWi81GnuprY5btJ8JiwCdM2KHsYOREVzJPkNFnzmz25WTPMDxemB+eY3JByMMGhUmlUQD2c+zwIjarZa8147DZmA7afSTDqsRGRcYmWgbHrjqPRxM2yeVstJMhDfkHyqqk9jxopk0GbqIgMSVQjbxjhj6v3kuA/c7sm/vbSZ/MwRY3piGoVyvKXQGRCNp24iaFiTpssWrlVwsohPrAw53d6U8IEh0o910q28lvZ5pPxeG/5IU3c74ttaFGjAla/5wtLXtiIFC3ubhLafpn5RUdtb0gnODvYGOWYyLiYzpJidGnkLBOdN31ramKH4x6g7riNWLDjglGPjYgMIzktHV3n7tJJbn5+p5W8HjC5MV1McEhxpvRojFNTuqBalm9tMzedlf1zDly+Z9RjIyL9mR1+FvXHb8KZ+Eea2HvtasnEpl2W+bPINLGJihTJ1soCUSGdcO7WIwR+tUsT7/1dtLw9PCEAlTiCgqhU2n3+Dt5YnDldhFC7SjnZJ4/TRZAaExxStHpOdvLb3G8Hr+K/vx/TxFtMjUBAQyd8/4aXZoIvIirZEhKfotX0rTrx7Z+0h7tjOaMcE5VcbKIik9Db2xWxocHo1vTZmmZCxKlbqDU2DL/ujzPqsRFR7sSSCq//uFcnufn2NU/5BYbJDeWECQ6ZDLFA57z+LXD0swCtaddD1h6X/XNO3Ug06vERka4lu2NRe2wY9ly4q4n1/fcLS/dm1Vhk9FxsoiKTU8HWSk4SeCTuPl5ZEKWJB38TicrlrLDz0w4ob80/DSJjOnr1Af4zf49WrKKtJXZ92gF2NpZGOy4qPXgVJ5PlWaOirN7+MfISPt9wWsbuPklB44mb5TfEL3o1kbU+RFR8HialotX0LUhOy9CKbxjxIhpVc+CvgvKNTVRk8gb718KFaUHwrVVZUxarD16Fe0gYNhy7afLlQ1RcyysMX3kYzaaEayU3015pLL+IMLmhgmINDpH4QzA3w6/vtsatxKfwydKR8YOVh/HBSmDnf9ujZmWO0iAyhDWHruHj32K0Yp0bVMWiAd4c5UiFxgSHKAsnexv5bTHy/G0MWLxfE2/35Q7UrVoef334IufZINKT87ceISDLPFVqh8Z3RuXy1ixnKhI2URHlwL9uFZnofNChtiZ2PuExPCZswsxNZ1hmREWQlJIGv9CtOsnN/73nK//umNyQPjDBIcrFf7t44OznXVHfyU4TW7DjohxWvuvcbZYdUQFN+vMkGn62GTcePs3yd1ZfJjat3CuxPElv2ERFlAdrC3NsHtUWV+4+kU1Vam8uedaEtW9sJ9m0RUTPJybWHLLsoFbMs0YFWWtjac7v2qR/THCI8kl0MhbfMsOO38T7vxzWxEWnZB/3SvhlsI/srExEma7dT8KLM7RX+haixnREtQplWVRkMLwaExVQcBMXOYvqa61cNbF9sfdQZ9xG/LDrEsuTCEBKWga6f7tbJ7lZMtBbflFgckOGxgSHqBDEBIChPZvixOQucMwy2mNa2GnZP+dw3H2WK5msb7aeR73xG3H8+kNN7J027jKx6ejhZNRjI9PBJiqiIhBLOhwc3xmnbyYi6OtITbzngihYmpfBgXGd5dIQRKYg+uJdvPbDXq1Yzcq22DyyLadXoGLHBIdIDxq42Mtvpyv3xWHsuuMylpquQvMpEXipiQvm9ffksg+kWHceJ8P78y068S2j26FO1fJGOSYiNlER6VF/nxq4ND0YXRplVsNvOH5TLvvwfwevsqxJUTIyVBj4036d5GZu3+Yy4WdyQ8bEBIdI339UZmXw/QBvHJ4QgKxrdX76+zHZP+ds/COWOZV6y6Ivo9bYMOw4mzkfVE/PF2QH/P94vmDUYyMS2ERFZCCVylkhNvQlHLx8D69+F62Jd5m7C872Ntj2STvYWvFPkEqXE9cfotu3u3X6ou0Z0xEOZS2NdlxE2fHqSmRg3m6VZHX9wh0XMePfZR7iE5/K2VzfaF0DnwXX5++ASrxHT1PR7vNteJScphX/c3gbNK1ewWjHRfQ8bKIiKibD2tfG+WlB8K5ZURNbsTcO9T6LQMzdLG1ZRCWISqXCivNmaDFtu1ZyM/nlRjJxZ3JDJRUTHKJiJKak/32YH6JDOmrFl5wzR90J4bh6L4m/Dyox1h+9LhPwA3cyPyra1quCi9OD8Zafm1GPjSgvbKIiMgIXh7Ly2+/2swl4+6cDmrj/zO1yyPn6D9rAyoLfP8g4Lt5+jE6zd+rExbxOVewyJ7YkKsl4BSUyog71q+L81EB0dMnQxMSkgWIW2DkR5/i7oWL1NDUd7b7crpPcDG+YLs9TJjdUmhg0wbl//z4GDBgABwcHuYmfHzx4kOcU+DltX375pWaf9u3b6zzer18/Q74VIoPq4ZaB4591grtjOa3p7sWw8qgLd1j6ZHCf/30KHhM24crdzGbSjzrVlYlNXQcVfwNU6hi0iap///64du0aNm3aJO+/++67Msn566+/nvucmzdvat3fuHEjBg0ahF69emnFhwwZgilTpmjuly3LVWmpdLOxNMf2T9rj0u3H6JjlG3T/H/fJ2/3jOqGqnY0Rj5CUaPuZBLy9NLOZVGjyggPWDPOTzaSpqalGOzaiEpngnD59WiY2e/fuhY+Pj4z98MMP8PX1xdmzZ1G/fs5DY52dnbXur1+/Hh06dECtWrW04ra2tjr7EilBrSrlZf8c0cHzo1VHNfFW07bCv64jlr7dCuZmHHVFRXPjwT/w+2KbTjzy0w5wrWTL4qVSz2AJTnR0tGyWUic3QuvWrWUsKirquQlOVrdu3cKGDRvw888/6zz2yy+/YMWKFXByckJQUBAmTpwIOzu7HF8nOTlZbmqJiYnyVnwz4beTglOXGcvOsOUZ3KgqgqYEYMy6k1h75IaMRZ6/g9pjwzA+uD7e8q2ppyNQFp6feZRPegZeX3wAR65mrvQtLOzfHJ0bVNUqQ5Ynz8+SpiCfOwZLcOLj41G16rM/lqxETDyWHyKxEUlLz549teKvv/463N3dZQ3OiRMnEBISgpiYGEREROT4OqGhoZg8ebJOPDw8XNYEUeE8r7xJv+XZzgbwaQlMPmyOpPRnNTefh52V28dN0lCDaxny/MzvOXa9DP6OM9eK+Ttn4FX3DKTEHkRYbMHPTyoclmfhJCUlGS7BmTRpUo7JQlYHDjxrzxWdf3OaNCqneE6WLFkikxkbGxud/jdqjRs3Rt26deHt7Y3Dhw+jRYsWOq8jEqDRo0dr1eC4uroiMDAQ9vb2+ToW0s6gxR9nQEAALC05NXtxlWfPl4GTNxLxn4V7NbHZxy1ga2WOyE/awp7T5PP8fI6DV+7jtR+1+9m4ONhg0wi/PJcL4d+7frE8i0bdAmOQBGf48OF5jlhyc3PDsWPHZBNTdrdv35bNSnmJjIyUfXVWr16d574iqREfDOfPn88xwbG2tpZbduI5/IAuPJZf8Zdn85qVZf+c5dGXMWH9SRlLSkmH1/TteLlZNXzdr3m+v0AoHc9P4N6TFLSYqlvzEj6qLeo55dykz/IsHjw/C6cgn9kFTnAcHR3llhfRmfjhw4fYv38/WrVqJWP79u2TMT8/vzyfv3jxYnh5eaFZs2Z57nvy5EmZFbu4uOTzXRCVbgN83fC6T00MXnYQ284kyNifMTfk9lXfZnjFs7qxD5GMKCNDhfdWHELEKe0vmbN6N8OrXjw3yDQYbB6cBg0aoGvXrrI5SYykEpv4uVu3blodjD08PLBu3TqdKqjffvsNgwcP1nndixcvyuHhBw8exOXLlxEWFobevXvD09MTbdq0MdTbISpxzMzKYMnAljg4vrNWfNTqGDl/zoWER0Y7NjKelfviUGtsmFZy062pC2JDg5nckEkx6Dw4YqTTiBEjZF8X4eWXX8a8efO09hHNUKJWJ6tVq1bJvjqvvfaazmtaWVlh69at+Prrr/H48WPZl+all16So6jMzbU7zxGZAsfy1rLZat+lu+i7KLN/Tuc5u+BaqSzCR7ZDWSv+bSjdqRuJCP4mUitmbWGG/WM7w8GWfeXI9Bg0walUqZIcyp0bkchkJyYEFFtOREKzc6fuGilEps6n1rP+Od9uPY/Z/y7zcPXeP2jw2Sa83cYNE7s3MvYhkgE8Tk5D25nbZX+brP74oA2au1ZgmZPJ4lpURArzYae6OPd5EJpVd9DEftpzWTZbbT2t2/GfSifx5fB/vx9D44mbtZKb8S81kIkukxsydVxNnEiBxBT764e/iGv3k/DijO2a+KCfD8rbPWM64oUKXN6ktNpw7CY+WHlYK+ZXuzKWvdMKFub83kokMMEhUrDqFW3lt/ktp27JEVdqbb7YJmt4fh/mB0t+IJYal+88QftZO3Ti+8d2QlV7rlNGlBVTfSIT0Lmhk0x0Bvq5aWIx1x6i7riNmLftvFGPjfL2NDUdnefs1EluVgzykb9XJjdEupjgEJmQSS83wukpXbWap2aFn5P9c8QoLCp5Zmw6A48Jm3Ah4bEmNrxDHZnYvFg37znJiEwVm6iITIwYMi764Ih5csRQcjX1EPND4zujcnndmb+peO06dxtvLtmvFfNwtsP64W1gbcFh/0R5YYJDZKLqVLWTtQBrD1/D6P+L0cS9Pt+Cjh5V8eOb3nIyQSpetxKfwmf6Vp34zv+2R83K5fjrIMonNlERmbieLarLWW7FWlZqYvkHMRvu8r1XjHpspiQtPQN9v4/WSW4Wvt5CJqJMbogKhgkOEckFOr95zRMxnwWirGVm88eEP07I/jknrmvPNk76tWjXRdQZtxH7Yu9pYv19asjEM6gJ19gjKgw2URGRhpjS//TUroi5+gA95u/RxLt9uxsOZS2x+38dYGfDaf/15XDcffRcEKWz9IZojipnzcszUVHwL4iIdDRzrSCbRZbsjsWUv0/J2MN/UtFkUjh6taiOWb2bylofKpwHSSloOW0LUtO1l6rZ+JE/GrjYs1iJ9IBNVET0XO+86I6L04Phn2U48prD1+AeEoY/Y26w5AqxvMIHvxxG8ykRWsnNFz2byISSyQ2R/jDBIaJcmZuVwfJBPtg/rpNWfMSvR2T/nNg7T1iC+fB/B67KxHDD8ZuaWJdGTrg0PRj9WtVgGRLpGZuoiChfqtrZyFqGPRfu4PUf92niHWbtgLtjOdm8YpOlgzI9czb+EbrMzZxvSBCj7w+OD0ClclYsJiIDYQ0OERVImzqOMtEZ0bGOJiZqccRsu6Fhp1ma/3qSnAaf6Vt0kps1w3xxKfQlJjdEBsYEh4gKZXRgfZz9vKtWv5Hvd12SzVY7ziaYdD+b8X8cR6OJm3ErMVkT/19XD5kYetWsZNTjIzIVbKIiokITSwaIpqmr95LgP3O7Jj7wpwPydm9IJzg7mM4q15tOxGPoikNasZZuFfHrkNaw4KrtRMWKCQ4RFZlrJVtZO5H9A7516FZ416yIVe8q+wM+e4KnFh3SES4OmQubElHxUe4Vh4iKXdfGznL23TdaZ44KOnjlvpyl97udFxX3G0lJy0DQ15E6yc1Pb7eUCR+TGyLjYYJDRHolJgD8/D9NcHJyFzjZZ65K/sXGM7J/zsHLmcsRlGZzIs6h3viNOH0zURN7r20tmdh0qF/VqMdGRGyiIiIDEUsN7BvbWWeY9KvfRUNMgnx4fAAqlsJh0lEX7qB/lmHyQq0q5RA2gsPkiUoS9sEhIoOq72wnazXERHefrjkmYyoV4Dk1Qk50t/B1L5iJiWFKuIRHT9FqmvZK38K2j9uhVpXyRjkmIno+NlERUbHo09L12erYjZ01sc0nb6HW2DD8uj+uxP4W0jNUeOPHfTrJzbevecrEjckNUcnEBIeIirV/zsI3vHD0swBYmmfW2oSsPS7752Ttz1IS/LQnFrXHhmH3hTuaWG+v6jJR696smlGPjYhyxyYqIip2FWytcH5aMA5duY9eC6M0cTEiybG8FXb8twPKWxvv8hRz9QF6zN+jFXMoa4nd/+sAOxtLox0XEeUfExwiMhqvmhVlM8+iXRcxPeyMjN15nILGEzfjtVaumP5KE1nrU1we/pMK39CtSEpJ14r//eGLaPyCQ7EdBxEVHZuoiMjo3m1bGxemBcHHPXMZg1/3P1t9e2OW1bcNubyCWB292eRwreRm6n8aywSMyQ1R6cMaHCIqEcRMx6vf80X8w6dyBmS1Yb8clrc7/9seNSuX0/v/u/bwNYz+vxitWEePqvjxTe9SMbqLiHLGBIeIShSxdpWoNdl57jbeWrJfE2/35Q7UcyqPvz58Ua6BVVQXEh6h8xztlb6FQ+M7o3L5zAkKiah0YhMVEZVI7epVkYnOsPa1NbFztx6j/vhN+HLzs/46hfFPSjrafLFNJ7lZ/W5r+f8xuSFSBiY4RFSi/a+rB85M7Yo6VTMn05u//aIcVr77fObwbfXaUD9FXcbvsWbyVtzPavJfJ9Hgs024/uAfTezjgHoysfGpVbkY3g0RKSLBmTZtGvz8/GBra4sKFSrku7PfpEmTUK1aNZQtWxbt27fHyZMntfZJTk7Ghx9+CEdHR5QrVw4vv/wyrl27ZqB3QUTGZmNpji2j22HHJ+214m8s3icTnYTEpwgNOwWPCRsxfeM5RMabyVtxX8S3nr4l9/tpz2XNc5u5VsD5aUH4sFNdI7wjIirVCU5KSgp69+6NYcOG5fs5M2fOxJw5czBv3jwcOHAAzs7OCAgIwKNHjzT7jBw5EuvWrcOqVauwe/duPH78GN26dUN6uvbQTiJSFjfHcrK2ZV5/T614q+lb8f2uWGSotPcX90V80M8HteJ7xnTE+g/awNKcldhESmXQv+7Jkydj1KhRaNKkSb5rb+bOnYtx48ahZ8+eaNy4MX7++WckJSVh5cqVcp+HDx9i8eLFmD17Njp37gxPT0+sWLECx48fx5YtWwz5doiohOjWtJqcTbivt2uBnidGRokE6YUKZQ12bERUMpSoUVSxsbGIj49HYGCgJmZtbY127dohKioK7733Hg4dOoTU1FStfURzlkiGxD5dunTReV3RpCU2tcTEZ9PBi9cRGxWMusxYdvrB8iy8z3s0gGsFa8zaciHPfccG1UO7upV43hYQz0/9YnkWTUE+d0pUgiOSG8HJyUkrLu5fuXJFs4+VlRUqVqyos4/6+dmFhobK2qTswsPDZf8gKpyIiAgWnR6xPAtnX6xZviqjI4+cgdODU4X8X4jnp36xPAtHtOgYLMERHYBzShayEn1nvL29UVjZp2YXTVd5Tdee2z4hISEYPXq0Vg2Oq6urrAWyt7cv9HGacgYt/jhF3yhLS67Lw/I0rltRlxG58Vye+/l7eiDYz61YjklJ+PfO8ixJ1C0wBklwhg8fjn79+uW6j5tb4S4iokOxIGpiXFxcNPGEhARNrY7YR3Revn//vlYtjthHjNjKiWjmElt24sOZH9CFx/LTL5Zn4QxsUxtfbDqn08E4KzEhsdjP0oKdiguL56d+sTwLpyCf2QX+axdDsz08PHLdbGxsUBju7u4ygcladSeSmZ07d2qSFy8vL/kGs+5z8+ZNnDhx4rkJDhEpl5WFGYb4u+e6j3hc7EdEpsOgfXDi4uJw7949eSuGcB89elTG69Spg/Lln03aJRIi0UfmlVdekU1MYgj49OnTUbduXbmJn0U/mf79+8v9HRwcMGjQIHz88ceoXLkyKlWqhE8++USO1BKjqojI9IQEN5S3P0RqDxUXNTciuVE/TkSmw6AJzmeffSaHeauJId3C9u3b5QR+wtmzZ+XQb7VPP/0U//zzD95//33ZDOXj4yM7A9vZ2Wn2+eqrr2BhYYE+ffrIfTt16oSlS5fC3Lzo69MQUekkkpiPAz2wdM9F2aFY9LkRzVKsuSEyTQZNcETSIba8OgdnJWpxREdmsT2PaAL79ttv5UZEpCaSmbf93ORoKdGhmH1uiEwXG6WJiIhIcZjgEBERkeIwwSEiIiLFYYJDREREisMEh4iIiBSHCQ4REREpDhMcIiIiUhwmOERERKQ4THCIiIhIcQw6k3FJpZ49uSDLrlOm1NRUJCUlyfLjauxFx/LUL5Yny7Mk4/lZNOrP7eyrIOTEJBOcR48eyVtXV1djHwoREREV4nNcLL6dmzKq/KRBCpORkYEbN27IBTzF2ldU8AxaJIdXr16Fvb09i6+IWJ76xfJkeZZkPD+LRqQsIrmpVq0azMxy72VjkjU4olCqV69u7MMo9URywwSH5VlS8fxkeZZkPD8LL6+aGzV2MiYiIiLFYYJDREREisMEhwrM2toaEydOlLdUdCxP/WJ5sjxLMp6fxcckOxkTERGRsrEGh4iIiBSHCQ4REREpDhMcIiIiUhwmOERERKQ4THAoX6ZNmwY/Pz/Y2tqiQoUK+XqO6L8+adIkOeNk2bJl0b59e5w8eZIlDuD+/fsYMGCAnLBKbOLnBw8e5Fo2AwcOlDNvZ91at25tkuW5YMECuLu7w8bGBl5eXoiMjMx1/507d8r9xP61atXCd999V2zHqrTy3LFjh855KLYzZ84U6zGXVLt27UL37t3ldU+Uyx9//JHnc3h+GgYTHMqXlJQU9O7dG8OGDct3ic2cORNz5szBvHnzcODAATg7OyMgIECzFpgp69+/P44ePYpNmzbJTfwskpy8dO3aFTdv3tRsYWFhMDWrV6/GyJEjMW7cOBw5cgT+/v4ICgpCXFxcjvvHxsYiODhY7if2Hzt2LEaMGIE1a9YU+7EroTzVzp49q3Uu1q1bt9iOuSR78uQJmjVrJq97+cHz04DEMHGi/Prpp59UDg4Oee6XkZGhcnZ2Vn3xxRea2NOnT+Vzv/vuO5Mu8FOnTompGVR79+7VxKKjo2XszJkzz33eW2+9perRo4fK1LVq1Uo1dOhQrZiHh4dqzJgxOe7/6aefysezeu+991StW7c26HEqtTy3b98uz9X79+8X0xGWXqKc1q1bl+s+PD8NhzU4ZBDiW0l8fDwCAwO1Jrhq164doqKiTLrUo6OjZbOUj4+PJiaamkQsr7IRzQNVq1ZFvXr1MGTIECQkJMDUahIPHTqkdV4J4v7zyk6Ud/b9u3TpgoMHDyI1NRWmrDDlqebp6QkXFxd06tQJ27dvN/CRKhfPT8NhgkMGIZIbwcnJSSsu7qsfM1Xi/YskJTsRy61sRLPBL7/8gm3btmH27Nmy2a9jx45ITk6Gqbhz5w7S09MLdF6JeE77p6WlydczZYUpT5HULFq0SDbxrV27FvXr15dJjuh7QgXH89NwTHI1cXpGdACePHlyrsUhPkS9vb0LXWSik11WotY2e8zUylPIqQzyKpu+fftqfm7cuLH8vdSsWRMbNmxAz549YUoKel7ltH9OcVNVkPIUCY3Y1Hx9fXH16lXMmjULbdu2NfixKhHPT8NggmPChg8fjn79+uW6j5ubW6FeW3QoVn87Ed/41ESTSvZvi6ZWnseOHcOtW7d0Hrt9+3aBykaUq0hwzp8/D1Ph6OgIc3NzndqF3M4rcS7mtL+FhQUqV64MU1aY8syJaGJdsWKFAY5Q+Xh+Gg4THBO/uInNEMSQU/GHGxERIdvq1e39YjjkjBkzYMrlKb7xPnz4EPv370erVq1kbN++fTImhuLn1927d+U356wJpNJZWVnJYczivHrllVc0cXG/R48ezy3vv/76SysWHh4ua8AsLS1hygpTnjkRo69M6TzUJ56fBmTADsykIFeuXFEdOXJENXnyZFX58uXlz2J79OiRZp/69eur1q5dq7kvRlCJUVMidvz4cdVrr72mcnFxUSUmJqpMXdeuXVVNmzaVo6fE1qRJE1W3bt209slanqKcP/74Y1VUVJQqNjZWjmTx9fVVvfDCCyZXnqtWrVJZWlqqFi9eLEekjRw5UlWuXDnV5cuX5eNi9M+AAQM0+1+6dElla2urGjVqlNxfPE88//fffzfiuyi95fnVV1/JkUHnzp1TnThxQj4uPkrWrFljxHdRcoi/VfX1UZTLnDlz5M/iGirw/Cw+THAoX8QQZfHHmn0TH7SakwmQw8izDhWfOHGiHC5ubW2tatu2rUx0SKW6e/eu6vXXX1fZ2dnJTfycfdht1vJMSkpSBQYGqqpUqSI/jGrUqCF/J3FxcSZZnPPnz1fVrFlTZWVlpWrRooVq586dmsdEubRr105r/x07dqg8PT3l/m5ubqqFCxca4aiVUZ4zZsxQ1a5dW2VjY6OqWLGi6sUXX1Rt2LDBSEde8qiH0WffRDkKPD+LTxnxjyFriIiIiIiKG4eJExERkeIwwSEiIiLFYYJDREREisMEh4iIiBSHCQ4REREpDhMcIiIiUhwmOERERKQ4THCIiIhIcZjgEBERkeIwwSEiIiLFYYJDREREisMEh4iIiKA0/w9swGSt7gUiRgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "uNum = coll.solveDahlquist(lam, u0, tEnd, nSteps)\n", "plt.plot(uNum[0].real, uNum[0].imag, 's', ms=10, c=\"orange\")\n", "plt.plot(uNum.real, uNum.imag, 'o-')\n", "plt.axis(\"equal\")\n", "plt.grid()\n", "print(\"L_inf error : {:1.5f}\".format(coll.errorDahlquist(lam, u0, tEnd, nSteps, uNum=uNum)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the error is way lower, even if we used $Q$-coefficients of the same size ... but there is a major difference here between the two $Q$ matrices :" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Q for RK4 :\n", "[[0. 0. 0. 0. ]\n", " [0.5 0. 0. 0. ]\n", " [0. 0.5 0. 0. ]\n", " [0. 0. 1. 0. ]]\n", "Q for Collocation :\n", "[[ 0.11299948 -0.04030922 0.02580238 -0.00990468]\n", " [ 0.234384 0.20689257 -0.04785713 0.01604742]\n", " [ 0.21668178 0.40612326 0.18903652 -0.0241821 ]\n", " [ 0.22046221 0.38819347 0.32884432 0.0625 ]]\n" ] } ], "source": [ "print(\"Q for RK4 :\")\n", "print(rk.Q)\n", "\n", "print(\"Q for Collocation :\")\n", "print(coll.Q)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While $Q$ is **dense for the collocation method**, it is **lower triangular** for RK4, \n", "which allows to solve for the first node (or stage) first, then the second one, etc ...\n", "Hence we need to solve instead $M$ equations sequentially instead of solving a system of $M$ equations **all-at-once**,\n", "reducing thus the computation cost for one time-step.\n", "This is one reason why RK methods have been generally favored in scientific computing against collocation methods.\n", "\n", "> 🔍 In this case (Dahlquist), solving the _all-at-once system_ is easy and cheap as showed above. \n", "> But for large scale non-linear problems, this can quickly become unfeasible, as solution at each time node\n", "> may represent thousands or millions of degrees of freedom ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, the high accuracy of collocation methods motivates to solve the _all-at-once system_ in a cheaper way, \n", "which is the main idea of **Spectral Deferred Correction (SDC)**, a.k.a **Iterated Runge-Kutta methods**.\n", "Those use fixed-point preconditioned iterations to solve the all-at-once system, where the preconditoner is built using a **lower triangular** approximation of the $Q$ matrix, named the $Q_\\Delta$ **matrix**.\n", "\n", "The second main feature of `qmat` is then to [generate those approximations ...](./03_qDelta.ipynb)" ] } ], "metadata": { "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 2 }