Skip to content

Instantly share code, notes, and snippets.

@mscroggs
Created October 28, 2022 10:31
Show Gist options
  • Save mscroggs/45ab606d6e69b811122b2697821267b1 to your computer and use it in GitHub Desktop.
Save mscroggs/45ab606d6e69b811122b2697821267b1 to your computer and use it in GitHub Desktop.
lecture4.ipynb
Display the source blob
Display the rendered blob
Raw
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"provenance": [],
"authorship_tag": "ABX9TyPfiLuyW+WXLNMB4yNiJDhx",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/mscroggs/45ab606d6e69b811122b2697821267b1/lecture4.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {
"id": "FY3gDECWsGp6"
},
"outputs": [],
"source": [
"import numpy as np\n",
"import math\n",
"\n",
"from scipy.sparse import coo_matrix, linalg"
]
},
{
"cell_type": "code",
"source": [
"def make_matrix(N):\n",
" A = np.zeros(((N+1)**2, (N+1)**2))\n",
" b = np.zeros((N+1)**2)\n",
"\n",
" h = 1/N\n",
"\n",
" for i in range(N+1):\n",
" A[i, i] = 1\n",
" b[i] = 0\n",
" for i in range(N**2+N, (N+1)**2):\n",
" A[i,i] = 1\n",
" b[i] = 0\n",
" for i in range(N + 1, N**2+N, N+1):\n",
" A[i,i] = 1\n",
" b[i] = 0\n",
" for i in range(2* N + 1, (N+1)**2-1, N+1):\n",
" A[i,i] = 1\n",
" b[i] = 0\n",
" for i in range(1, N):\n",
" for j in range(1, N):\n",
" index = j * (N+1) + i\n",
" A[index,index] = 4/h**2\n",
" A[index,j * (N+1) + i-1] = -1/h**2\n",
" A[index,(j-1) * (N+1) + i] = -1/h**2\n",
" A[index,j * (N+1) + i+1] = -1/h**2\n",
" A[index,(j+1) * (N+1) + i] = -1/h**2\n",
" b[index] = 1\n",
" return A, b"
],
"metadata": {
"id": "CVRdpBQ0sK0L"
},
"execution_count": 70,
"outputs": []
},
{
"cell_type": "code",
"source": [
"N = 4\n",
"A, b = make_matrix(N)"
],
"metadata": {
"id": "Vatz07_i4Uq7"
},
"execution_count": 71,
"outputs": []
},
{
"cell_type": "code",
"source": [
"b"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "1V-yIjEv4Wu1",
"outputId": "018ced50-fe16-4b29-c6b9-b6acfb79923a"
},
"execution_count": 72,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 1., 1., 1., 0., 0., 1.,\n",
" 1., 1., 0., 0., 0., 0., 0., 0.])"
]
},
"metadata": {},
"execution_count": 72
}
]
},
{
"cell_type": "code",
"source": [
"A"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "atU0RYXh7izD",
"outputId": "c3a5b0b9-9ae0-4bea-8f40-5d02c75cd289"
},
"execution_count": 73,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([[ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., -16., 0., 0., 0., -16., 64., -16., 0., 0., 0.,\n",
" -16., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., -16., 0., 0., 0., -16., 64., -16., 0., 0.,\n",
" 0., -16., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., -16., 0., 0., 0., -16., 64., -16., 0.,\n",
" 0., 0., -16., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., -16., 0., 0., 0., -16.,\n",
" 64., -16., 0., 0., 0., -16., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., -16., 0., 0., 0.,\n",
" -16., 64., -16., 0., 0., 0., -16., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., -16., 0., 0.,\n",
" 0., -16., 64., -16., 0., 0., 0., -16., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" -16., 0., 0., 0., -16., 64., -16., 0., 0., 0., -16.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., -16., 0., 0., 0., -16., 64., -16., 0., 0., 0.,\n",
" -16., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., -16., 0., 0., 0., -16., 64., -16., 0., 0.,\n",
" 0., -16., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,\n",
" 0., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 1., 0., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 1., 0.],\n",
" [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
" 0., 0., 1.]])"
]
},
"metadata": {},
"execution_count": 73
}
]
},
{
"cell_type": "code",
"source": [
"sol = np.linalg.solve(A, b)"
],
"metadata": {
"id": "gCyWtfV67l1K"
},
"execution_count": 74,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sol"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "v6Un2-aR8RjU",
"outputId": "ab5a6baf-db81-490c-cb45-e3db3593f5cb"
},
"execution_count": 75,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"array([ 0. , -0. , -0. , -0. , 0. ,\n",
" 0. , 0.04296875, 0.0546875 , 0.04296875, 0. ,\n",
" -0. , 0.0546875 , 0.0703125 , 0.0546875 , -0. ,\n",
" -0. , 0.04296875, 0.0546875 , 0.04296875, 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ])"
]
},
"metadata": {},
"execution_count": 75
}
]
},
{
"cell_type": "code",
"source": [
"from matplotlib import pyplot as plt\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from matplotlib import cm\n",
"\n",
"u = sol.reshape((N+1, N+1))\n",
"\n",
"fig = plt.figure(figsize=(8, 8))\n",
"ax = fig.gca(projection='3d')\n",
"ticks= np.linspace(0, 1, N+1)\n",
"X, Y = np.meshgrid(ticks, ticks)\n",
"surf = ax.plot_surface(X, Y, u, antialiased=False, cmap=cm.coolwarm)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 466
},
"id": "UZB1KLnK8SY1",
"outputId": "e50f2980-ebed-42ff-ba89-8dc65090fbce"
},
"execution_count": 76,
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 576x576 with 1 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAHBCAYAAADkRYtYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9e2wb95n3+x2SongTb5Ily6Rs2VZ8k235IsqbdNPd19muT727arcNtulp0gJJgD0HCRC0B2i6f2wQBAtst7vvP4ucPUXx5o9gsam8SbE1utg3b9PinDbFNrbjNJYoyxIlUZZE3UVS4v0ynPOH9zchKV5mhjPkiPx9gACxRA1/vMx853l+z/N9GI7jQKFQKBQKpTKaRi+AQqFQKJT9ABVMCoVCoVAEQAWTQqFQKBQBUMGkUCgUCkUAVDApFAqFQhEAFUwKhUKhUASgq/J72nNCoVAolFaDKfVDGmFSKBQKhSIAKpgUCoVCoQiACiaFQqFQKAKggkmhUCgUigCoYFIoFAqFIgAqmBQKhUKhCIAKJoVCoVAoAqCCSaFQKBSKAKhgUigUCoUiACqYFAqFQqEIgAomhUKhUCgCoIJJoVAoFIoAqGBSKBQKhSIAKpgUCoVCoQiACiaFQqFQKAKggkmhUCgUigCoYFIoFAqFIgAqmBQKhUKhCIAKJoVCoVAoAqCCSaFQKBSKAKhgUigUCoUiACqYFAqFQqEIgAomhUKhUCgCoIJJoVAoFIoAqGBSKBQKhSIAKpgUCoVCoQiACialZeE4DhzHNXoZFApln6Br9AIolEaQy+WQTqeRSCSg0+nQ1tYGrVYLnU4HhmEavTwKhaJCmCp32PT2m9JUcBwHlmWRyWTAcRzS6TQYhuGjTYZheOHU6XTQarXQaGgihkJpMUreNVPBpLQMRCBzuRwvkul0ukAQiXDmnxcajQZtbW1UQCmU1oEKJqV1yWazyGazfBTJMAyflq0kgMUCyrIsdnd30dPTw6dxqYBSKE1HScGke5iUpobjOF4sGYYRLW5EXAmZTAbLy8twOBxIp9MAHkWg+fugxX9DoVCaAyqYlKaFRJD5USWBZVnMzMwgFArB4XDAbrfDZrNBp6t8ShDB1Wq1AMBHnul0mhdQhmH2pHCpgFIo+x8qmJSmg0SVS0tLSKfTOHr0aMHvI5EIvF4vent74Xa7EYlEEAwG4ff7AQB2u50X0La2torPRYSwlICmUikApfdAqYBSKPsPKpiUpqK4sKf4d4uLi1hZWcHZs2dhNpuRTqfR1dWFrq4uAI/2Ond2dhAKhfDw4UNwHAebzcaLKCkWKkexgOavKT8CJVW4Op2OCiiFsk+ggklpCvLbRYDP9h6JuKVSKXi9XphMJoyMjECr1SKXy+05jk6nQ2dnJzo7OwE8EtDd3V2EQiEsLi7yz7GxsQG73Q69Xl91baRVJX+tmUymYK1UQCkU9UOrZCn7HiJALMsW7FWura0hFovBarViZmYGJ06cwIEDB/i/E1IlW0wqlcKnn36KAwcOIBwOg2VZWK1W2O12OBwOQQJaav35lbhEYEkalwoohVJ3aJUspfmoVNjDcRw2NjYQDocxPDyM9vb2PX8vVoiIkB07dgzAZ20m4XAYKysryGQyBQJa6jlLraF43SzLFrTB5EegpBKXQqHUFyqYlH1JtXaRSCQCn88Hk8mES5cuySYwxXuYWq0WDocDDocDwCMBJwI6NTWFdDqNjo4OXkANBoOg5ygloJlMhv85FVAKpf5QwaTsO3K5HDKZDF/YUywupLCnv78fyWSyrmKi0Wj4AiGy1kgkgnA4jOnpaaRSKVgsFr6VxWAwVF0feY3kpqBYQNPpNLLZLJxOJxVQCkVBqGBS9g3lCnsIxYU9oVAIiURC1jVUq5ItRqPRwGazwWaz4ciRI8jlcohGowiHw/D5fEgmkzCbzbyAGo1G0QIaj8cRDAYL/jZ/D5QKKIUiD1QwKfuCcoU9hI2NDfh8voLCHiHiRvYI64VGo4HVaoXVasXhw4fBcRwvoLOzs0gkErBYLHyUajKZBAsoMV3gOA65XA7JZJJ/DBFQOpGFQpEOFUyK6qnm2DM9PY1EIgGPx7OnSrWaYIoVDrERppDjdXR0oKOjA319feA4DrFYDOFwGPPz84jH4zCbzbyAms1mwQJKyBdQsvbiPVDqh0uhVIcKJkW1CCns8Xq9cLlcOH369B4hESIsaoNhGFgsFlgsFrjdbnAch3g8jlAohIWFBcRiMZhMJl5ALRYL/3eVjllKQFOpVEU3IgqFUggVTIoqqVbY8/DhQ6yuruLcuXO8aBQjdzRIjllPGIaB2WyG2WzmBTSRSPBGCtFoFFqtFlqtFru7u+jo6JAUgXIct0dAqaE8hVIIFUyKqhBT2HPlypWKkZASgtloGIaByWSCyWSCy+UCx3FYXV3F5uYmlpeXEY1G0d7ezrexWCyWqtFiOQHNt/OjAkqhUMGkqIhiH9hyhT0nT57kvV8r0YyCWQzDMGhvb4fFYsHx48cBAIlEAuFwGIFAAJFIBHq9nhfQjo4OSQIK0IksFAoVTIoqIFFlpcKeZDJZsrCnHGqskq0HRqMRRqMRvb29AIBkMolwOIzV1VVMT0+jra2NF1Cr1SpIQAE6kYVCoYJJaSgcxyESifApv+KL9+7uLrxeL9xud8nCnkoIEcxWuKgbDAYcPHgQBw8eBPAorR0Oh7G+vg6fz8e7Fdntdlit1gKj+FIInciSTqdhtVqpHy6laaCCSWkYpLDnwYMHOHr0KKxWK/+7/MKe8+fPly3sqUQrpGSl0N7ejp6eHvT09AB4FCmGw2Fsbm5idnYWWq22YCZoNQEFSk9kGR8fx+XLl/nf04kslP0OFUxK3Sku7Cm+ICeTSXi9XlgslqqFPUKei1IZvV6P7u5udHd3A3gkoDs7O9ja2sL8/DwYhikQUGKQUAmSVs9P42az2YJirnwzBepGRNkPUMGk1JVShT0Mw/CzKcUW9lSC9hJKQ6/X48CBA7xjUiaTwc7ODoLBIPx+PwDwAmq320UJKCF/IguBGspT1A4VTErdICYExYU9DMMgm81icnISqVRKVGFPNWiEWTttbW3o6urib2Cy2SzC4TDC4TAePnwIjuNgs9l4AW1ra6t6TCETWagfLkVtUMGkKE41x55MJoOpqSn09/fD7XYrNoqLIg86nW6PgO7s7CAcDmNxcRG5XA42mw2ZTAbpdFrQzU+piSy5XA6JRIIKKEU1UMGkKEq1Ac8LCwsIBoN47LHH4Ha7ZX1uKpj1QafTobOzE52dnQAetQHt7OxgY2MDXq8XLMvCarXylbhUQCn7FSqYFEWoFlWSwp6Ojg4cOnQI7e3tsq+hlQRTTWKh1WrhdDrR3t6OS5cugWVZfqh2IBBANpstGKot5LMvJ6DlJrJQP1yKElDBpMhONcee9fV1zM7O4tSpU+js7ITP51NE2FpFMNX+Gkmfp8PhAPAo60AE9P79+8hkMujo6CgYql2NShNZCNRQniI3VDApslHNB5ZlWTx48ADpdLqgsEcpYWsVwdxvaDQavkCov78fuVwOkUgEoVCI/35YLJaCodrVkGIoTwWUIhYqmBRZKE7BFkeVxLGnr69vT2GPRqOhgtnCaDQa2Gw22Gw2AI8iUDJUe2ZmBqlUqmCottFolDyRJZ1O43e/+x0GBwepoTxFNFQwKTUjpLBnfX29rGNPfh+mnFDB3J9oNBpYrVZYrVYcPnwYHMchGo0iFAphdnYWiUSiQEBNJpMoAU0mk3x0SSeyUMRABZMiGSGFPRMTE7BarRgZGSmbAqMpWUolGIZBR0cHOjo6eAGNxWIIhUKYn59HPB6H2WzmBdRsNgseHl7OUJ6IJTWUp+RDBZMiiUoDnoG9hT2VaISwcRyH5eVl+P1+fr/M4XAISvdRGgvDMLBYLLBYLOjr6+MFNBwOY2FhAbFYDCaTia/CLRbQUp+vUAEt3gOl35XWggomRRTVCnuy2SwePHiATCYj2LFHo9HUNSWbyWQwOTkJrVaLy5cvI5VK7Un35Rec0IuieOp5A5QvoG63GxzHIR6P805E0WgURqORF1Ahays3kSWTyRR896mhfGtBBZMiGHLBYFm2ZFS5s7ODyclJHD58GC6XS/DFQ8k9zGJIK8PRo0dx8OBBpNNptLW1FUQr+ftlyWSS3y8jEShF3TAMA7PZDLPZDJfLBY7j+KHaS0tLiMfjuHfvHn9TZLFYBFXMlprIQgW0taCCSRFEtcIev9+PjY0NDA0NwWw2izp2PVKyZI2bm5u4cOECTCZTyecstV8WiUQKKjaTySRWV1fhcDgE9Qy2ImoazM0wDEwmE0wmE3p7exGLxXDixAmEw2EsLy8jGo1Cr9fzN0UdHR2SBTR/IguAghQudSPa/1DBpFSEXAR8Ph+sVis/wYIgtLCnEkq1lRBSqRQmJibQ0dEBj8cjao0MwxRUbOZyOdy+fRvpdJrvGSRN90Jda1oBNQlmPuR7ZjQaYTQa0dvbC+DR9zgcDmNlZQWRSARtbW0FQ7WFCiidyNLcUMGklCW/sIc4qeQjprCnEkqlZIFHe6off/xx2XFhYi9YGo0GWq0WR44cwZEjRwqa7olrjdVqFWX7RqkfHMeVFD+DwYCDBw/i4MGDAB7dZIXDYaytrcHn80Gn0/FVuFarVfBQbTqRpbmggknZQ6nCHq1Wy4ualMKeSiiRks3lcvD5fEilUvj85z+vmHDlN90T15rd3V2EQiGsrq4ik8nwo68cDodsY8uKUduFVs0RppB1tbe3o6enBz09PQAeVcuGw2FsbGxgdnYWWq22YKi2GAEtZyifSqXAsiycTicVUJVCBZNSQLnCHhIFSi3sqYTcghmPxzE+Po6enh6YTKa6Rnn5tm8A+PeMGI9LmdxRDTX2mqpxTYB0Idfr9eju7kZ3dzeAzwR0a2sLc3NzBZ+7zWYTNVSbCGg0GkU0GoXJZOIfQyNQdUEFk8JTqbCHYRhsbGxgaWlJUmFPJeRsK1ldXYXf78eZM2dgt9uxuroqy3GlotFo+P3No0eP8pM7QqEQlpeXwbIsbDYbL6BChi/vF9R4cc/lcrJ4yBYLaCaTQTgcRjAYhN/vB8MwBUO1hQgox3F8cRD5N53Ioi6oYFIEOfYsLi7CaDRKLuyphBwRZjabxdTUFFiWhcfjESw89U4dFk/uILMjQ6EQFhcXwXFcwYV2vwrofk/JiqWtrQ0HDhzgi+Ky2SzC4TDfCyrkcy0WczqRRX1QwWxxqjn2rK2tYW5uDt3d3Whvb1fkhKxVMImxu5xp4mrIdeElsyOdTieARwJKLrT5AkoiUCGRihpoNcEsRqfToauriy80y2azfGqefK75qfm2traq0S+dyNJ49sfZR5Gd4sKe4hMrP2IbGRnB1tZWwZ2tnEgVTI7jsLi4iJWVlbLG7vsNrVaLzs5OvuqYXGhDoRAWFhYAgC8gErpXRvmMRgm5Tqcr+FxJZoGYKbAsC41GA7PZjHQ6LWhvu9JEFmoorwz0bGtBqg14LlXYo2Trh5Q9zHQ6Da/XC4PBgJGREUFVinJBBL5ekUqxgBbvlen1emi1WrAsW9f3oRJqjTDJd77RlMos+Hw+ZDIZeL1evjiMpHCFFK6VElCg9EQW8h8VUHFQwWwxSFRZzrFnfn4em5ubewp7lPJ7BcRHmMFgEFNTUxgYGODL/qU+L6Dei3spilN9mUwGi4uLCIVC+OSTT/hqTRKBNkpA1fqeluvDbDRarRYGgwFOpxPd3d0F7UkrKyvIZDK8QYbdbhfkMFXNUB4ovQeqxs9NLVDBbBGqFfYkEglMTEzAbreXLOxR0o1HqGByHIe5uTkEg0FcunRJFl9XKelgNY0Na2tr451ojh49yldrknYH0i/ocDgEN9w3M2oVcqCw6KdUexIxyMh3mMofql2NcobydCKLcKhgtgBCC3tOnz7Np4iKUTIlK+TYyWQS4+PjcDqdGB4eVmWUoAaKqzVLNdyTKl2hlm9SUKswqXVdQOWWl3yDDPJYMiSAeBznD9UWOmWHGsqLgwpmE1NsxVWtsKdSC4OSKdlq0evGxgZ8Pl9FQaeUplTDfSgUwvr6OmZmZiR5pgpBrcIkVx+mEohZm0aj4T2Ojxw5smfKDhlTRwTUZDLVJKDxeBzLy8s4duxYSwsoFcwmpVphDxlzdeTIERw6dKjql74RKVmWZTE9PY1kMlmTBZ/cF281pWQJQl+fXq8vsHwjs0DX1tYwMzMjaWrHfkKtQg7UJualpuxEo1GEw2HMz88jHo/DbDbzn61YAU2lUryFXytPZKGC2YSQvcpKhT1bW1v8mCsh1LvoJxqNYmJiAocOHcLp06eb+iRsJO3t7QWm48VTO/R6PZ/CFTo3ElCvMKl1XYC80W++gJI5r7FYDOFwGH6/H7FYDCaTiRdQs9lc8X3J5XIl21JabSILFcwmQmhhj8PhkDTmqh5tJRzHIRAIYHFxEWfPnoXVaq3p2Eq0gKgtwpRzLcVTO5LJJG/jF41G0d7eXiCg5d5XtQqTWtpKSqFkBS/DMLBYLLBYLHC73eA4DvF4nHciisViMBgMBUO1898nIpiljlttIsuPf/xjfOlLX8KhQ4cUeW31hApmk8CyLILBIP9FL74orK6uYn5+XvI+YD1SstlsFpOTk9BoNBgZGZGlKV9t4rbfMBgM6O3t5edGJhIJhEIhLC0tIRKJwGg0lr3IqhG1tpUA9RVzhmFgNpthNpvhcrnAcRwSiQTvRBSLxfibI7vdjmw2K2omaP5Elvfffx9f+MIXlH5JdYEK5j4nf8q71+vFE088UfD7bDaL+/fvg+O4qoU9lVA6JZvJZHD79m309/fLeieqhGC2sgiTwcuHDh3ac5ElkzbU7IGr1sgXaGxBEsMwMJlMMJlM/GdL0vPLy8sIh8MAwO9xC93fZhgG8Xhc8NaP2qGCuY/JL+wp9eUlhT39/f3o7e2t6UKhVEqW2NvF43F87nOfk/3EUkowKaUvsiQCXV1dRSQSQTab5VO4QgtNlETNKVk1VfAyDMPfHPX29mJjYwO7u7vQ6/X8/rbQCut4PN4UtpUAFcx9SakBz3IU9lRCiZRsKpXCxMQELBYLf+GVG6WiwVaNMCuRL6AWiwWrq6vo6+tDKBTiC03MZjN/kW2EgNIIUxosy0Kv1xek5/MrrH0+H3Q6Hd/Gkm+SkUgkKhorvP/++3jllVfAsixefPFFfO973yv4fSqVwje/+U3cvXsXnZ2duHHjBvr7+/Ev//Iv+Pu//3v+cePj4/jkk09w4cIFBd6BR1DB3GeUG/BMqKWwpxJyR5hbW1uYnp7GiRMncODAAfznf/6nbMfOp5XTp42E7BWSfTJSaEIqNUmrA+kVdDgcgpvt5ViXGlGzYJYq+imusE6lUgUmGWNjY3xbSjqdLpmmZ1kWL730Ej744AO43W54PB6Mjo7izJkz/GPeeustOBwO/pivvvoqbty4gW984xv4xje+AQCYmJjAl7/8ZUXFEqCCua+oNOAZeOQr+sknn+DMmTP8vEW5kGsPM5fLwefzYXd3F8PDw4JMpWuB7mE2hlKRXKlKzVgsxjfbJ5NJPgJ1OBwwGAyyC6iaBRNQb7o/l8tVLcJrb28v6PF1u934xS9+gV/84he4evUq2tvb8eSTT+LP//zPMTw8DAC4ffs2BgYGcOzYMQDAM888g5s3bxYI5s2bN/H6668DAJ5++mm8/PLLe75fP/7xj/HMM8/I+ZJLQgVzH1CtXYQU9mQyGTzxxBOKFFzIIRLxeBwTExM4cOAAhoeHVXFxoOKnDELe03wBJb2CxK3G5/MhmUzCYrHwKVw5vIPVnJJVMyzLir657e3txbPPPosf/ehH+OijjxAKhfDhhx8iFArxjwkEAujr6+P/7Xa7cevWrYLj5D9Gp9PBZrNhe3ubH0AAADdu3MDNmzelvDRRUMFUOdV8YMPhMCYnJ3H06FFEo1HFqhNrvciQtpbBwUHeULoe0AizcYj9zpRyq4lEIgiHw7xfKjEcJxGoWNRc9KNmak0XMwwDp9OJL33pSzKu6hG3bt2CyWTC2bNnZT92MVQwVUq1wp5cLof5+Xlsb2/j4sWLMJlMWFhYUN0ddDabxYMHD5DNZmtqa5FKNXELhULQ6XRVnU4o4pDje8gwDO+Xevjw4QLD8fyJHSSFKyQCUntKVq1InbVaTWhdLheWlpb4fy8vL8PlcpV8jNvt5geqk/mwADA2Noavf/3rotcmBSqYKqRaYQ9JbXZ2dhYU9pBKVrVc+CORCCYmJgoGUdebSj61U1NTSCaT0Gg0iMVifPqv2v6ZGiNMtXzmSlJsOJ4/8opsSVitVj6FW0pA1XR+7CekRpjVKmQ9Hg98Ph/8fj9cLhfGxsbwzjvvFDxmdHQUb7/9Nh5//HG89957uHr1Kv8Z5nI5/Ou//is+/PBD0WuTAhVMlVGtsGdlZQV+v79kYY9GowHLsg2/g+Y4DktLSwgEAjh//nxDe7Aq+dS6XC6cOHGC/z0pQCH7Z2Kjl0ahNvEG6iNM+SOv+vv7Sw5dttlsvIDq9XqakpWI1AgzHo9XFEydToc333wT165dA8uyeP755zE4OIjXXnsNw8PDGB0dxQsvvIDnnnsOAwMDcDqdGBsb4//+17/+Nfr6+viiIaWhgqkSqhX2ZDIZTE1NAQCuXLlSsmJNSTceoaTTaUxOTqK9vR0jIyMNH1hcLJgrKytYWFjgfWozmQx/cc8vQCHpv2AwyEcv5OLLcZwqRUpNNCKSKx66zLIsdnd3ebcalmX5dVmtVsnTb5RA7d8nqREm6b2txPXr13H9+vWCn73xxhv8/xsMBrz77rsl//YP//AP8dFHH4lel1SoYKqAaoU9JOV09OjRirZxjRZMss6BgQG+tLzREMEkKViyl1qtRD4//Ueil52dHYRCIYTDYXi9XjidTj56kcP3liIv+cOyjx49CpZlcf/+fSSTSXi9XrAsWxCBNtLOT+2pYqkRZiKRaBpbPIAKZkORUthTiUYJJsdxmJubQzAYxKVLl2Qp/5cL4mU5OTkJt9sNt9st6cKk0Wj4i28ymYTL5UI2m0UoFMLCwgIYhuF/n+9y0qqoUQC0Wi3a29vR3d0Nu90OlmX5m6DFxUVwHFcgoPW8CVKzaQEAyVs9zeQjC1DBbBjVBjyXK+ypRCMEM5lMYnx8HA6HA8PDw6o76ePxOHw+H4aGhmoeFZaPVquFzWbjq/UymQzC4TA2NzcxOzsLnU7HC2gzDmKuhhoFEyisktVqtXA6nfz0HlKBSUZecRzHp3iVFlC1C2a58V7VoIJJqRkSVUop7KmE0oJZPFtyY2MDPp8Pp06dKijzloqcF9n8KtjBwUFZxbJUIVFbWxsOHDiAAwcOAPjMZ5MYVQudI9ksqFUwKxX96HQ6dHZ28t9lIqAkiwCA7wG12WyyCqja212kFksJ2cPcT1DBrCNCCnvu378PhmHKFvZUoh6CSU6c6elpJBIJeDweWYon5Bz0TKpgSfq1EReiYp9NMsUjfwyWmqZ4tApivmOlBDQcDiMYDMLv94NhmAIBrSUNvx+qd6WsLx6PU8GkiIdlWb7nr1Jhz7Fjx/hpAGJRWjA1Gg0ikQimpqbQ29uLU6dOyXaSyzUNJRAI4OHDhzh37hw6OjowNTWlCqef4jmS8XgcoVCINyHPb2GR4mCjNtQaYdayLp1Oh66uLt6SjaTht7e3MT8/z+9j2+120QKq9pSs1PeMFv1QREEKe4LBIBYWFjA0NFTw+1wuh7m5OYRCoZoLZpQUTLLnOjk5iXPnzsma4gQ+i16l3qWTCkgyKJtE59XErRH7vvnT7okJOXGwmZ6eRiqV4hvwHQ6HoAhebeKk1jYJOSO54jQ8EdCtrS3Mzc1Bq9XyEWi1QjC1C6bUz5MYgjQLVDAVJL+wR6fT7bkwk8Kerq4ueDyemk9kpS7+2WwWk5OTyGQy8Hg8ipwAtbjnkBRsX1/fHkehaseV8pxyi1Oxh2p+A34gEKja/qBWcVKbiAPK7hUWC2g6nS4Yd0XmRRIBzV+H2gVTKolEQjUtZnJABVMhyF4lSQFptVpezDiOw8rKCh4+fIgzZ87IZkae/xxysbOzg8nJSb4XUamTWqpgFqdg5TpuNZQUqfwGfNI/WNz+YLfb4XQ6YbPZFFtHLTRjSlYser0e3d3d6O7uBvBIQMnA5ZmZGbS1tfFZBDU4dJWjlveMVslSKlKusEer1fLVsaSwR0gDvRgYhgHLsrIci+M4LCwsYH19HUNDQzCbzdjc3FRMKMTuYWazWX5/stL7qNS0knpSqv2B7J3Nzc2BZVkYDAbeaEENF161CmYji2v0en3BvEhSSb26uopgMAgAvIiqqRWplhtlWvRDKUslxx6NRoN0Oo3bt2/XVNhTCblSsqlUCl6vF2azGSMjI/zJQvYZlUDMscWYugu5OEq5uDcyDVpcfLK8vIzd3V2sr69jZmYGer2+oAe0UQKhRsFUU/tGfiX15uYmwuEwDAZDQSsSSeFaLJaGrVuqyw9ABZNSgmrtIqSwJ5FI4Pd///cVc8LRarW8a5BUtra2MD09jRMnTvB7MQQlC2SERoLLy8tYXFwsm4ItdVy516y2aSU6nQ4WiwWHDx8G8MhMIhQKYXl5GdFoFAaDgRfQeo0xU2uEqeZ1tbW17WlFIj640Wi0Yb28tUSYtEqWUoBQx56uri6YTCZFbeNqEbRcLgefz4fd3V1cvny5ZGuDXK0fpagmQtlsFvfv3wcA2VPZzYbBYEBvby96e3vBcRzfA7qwsMA3khMf3EpjzJoRtQpmKVEirUjkcyQ3QktLSwU3Qna7XVEBpRHmZ9CrjkSq+cCWKuxZX19XdE1SBTNf1IeHhyvOgVQqwqwkxvkpWLfbLeq4Su1hNiLC/LVzmP//zwc/FvQ3DMPAZDLBZDLB5XKB47i6jDFTqzCpdV3VojiGYfb08pIINN8Mg6Rw5cwk1BJh0rYSyp4UbPEXM5PJYHJyElqttq7RkBTBXFtbw9zcnCAbPqVTssXH5jgOgUBAVAq21HHVlD6VQr5Q5v9MqGjmI2aMmcPhkDzBQ63CBKhzb1WsKOXfCOWbYYTDYT6TIJebVC0RJk3JtoyMi/8AACAASURBVDjVBjwHg0FMTU3h+PHj/F5EvRAjaMXjroRcGOuZkiUp2Fqrife7YJYSy/zfSRHNfCqNMVtaWpI8wUPNgqlGSK+2VPLNMEgmodhNymw285+jGAGlVbKfQQVTIEIKe2ZnZxEOhxs24kqoYJIUZ19fn6hxV0pXyRJhI+s7cuQIXC6XbMeVi3qIcCWhLH7cifv/Ltvz5o8xA/YakNMxZsogd49zKTcpkoonAmqxWPgUrtFoLHsdqCXCzGazqhrUXStUMAVQbcBzLBbDxMQEuru7Kzr2EMFRqjy8mmByHIelpSUEAgGcP39e9N6CkilZjUYDlmWxvLyMpaUlSesrxX6MMIWKJWHmzJ/i6L1/U2QtxQbkQseY0QhTHEo7/RSn4okdYzgcxuzsLBKJBCwWS4GfMfn8mtWFSApUMCsgpLCH7LENDg5WdV0hgtMIwcxkMvB6vdDr9RgZGZF0x6hkSpbjOMzPz8NoNEpeXymqCaaUi7qSIixWLAn+oT9HX43pWSEIHWOWSqWaqthDaeotSvl2jPkCGgqFMDMzw39+DocD6XRa0vm4325UhUAFswwcxyGTyYBl2YqFPTqdTvAeG4milCoCKieYZBJKrfuqSqVkI5EIVldXcejQIZw6dUrWYyuZRpYTqUJZfIxa9zTFUm6M2c7ODoLBINbW1ugYMwE0OiIv9jPmOA6RSIR3ImJZFpFIhE/hipmo00yfORXMEihV2KOE12s+xYJJoratrS1Z9lXlTslyHIfl5WUsLy+jp6eHt36TG7XvYcohlsXHqrdwEkjrQyKRgM1mg9FobOoxZnKhtrQnwzB8MRjw6MbIaDQiHA7jwYMHSKfTVduRaITZ5Agt7NnZ2ZEkQCTCVIp8QU4mk5iYmIDdbofH45HlZNRoNMhmszUfB/hsAgppvVlYWFDkZkJNF6FSyCmWxcdtlGgCn0VMco8xa1bUJpj5kLF7NpsNNpsNR44cQS6X4yNQ0o5U/FmmUqmqN0Xvv/8+XnnlFbAsixdffBHf+973Cn6fSqXwzW9+E3fv3kVnZydu3LiB/v5+AMD4+Dj+8i//Eru7u9BoNLhz547iN2FUMP8LMYU9lZr7K1GPCJNlWWxsbMDn8+HUqVN8sYYcyJXe3N3dhdfrRX9/Pw4dOsQfW8n9UTmRY61KCWXxczRSNIvPkVrHmDUzahbMUlWyGo2GF1DSjkQ+y5WVFXz3u9/l+0M3Nzf32GyS47700kv44IMP4Ha74fF4MDo6ijNnzvCPeeutt+BwODA7O4uxsTG8+uqruHHjBrLZLJ599ln88z//M4aGhrC9vV2X70vLCyaJKgOBAHp6emou7KkEmViiJLFYDEtLS/B4PLLfsdeaks1PwRZXwSolmGqskq2HWOY/VyNEU8ienNAxZkRAay0EU9v3IB81C6aQteV/lgDw05/+FDdv3sQ//uM/4i/+4i8QjUbx+OOP4+WXX8aJEycAALdv38bAwACOHTsGAHjmmWdw8+bNAsG8efMmXn/9dQDA008/jZdffhkcx+HnP/85zp8/j6GhIQCQNTCoREsLJinsyWaz8Pv9eyaIpNNpTE5Ooq2tTRbHHiXbMmKxGMbHx8EwDC5duqTIRnstVbLZbBZer5cvkip1x6rEeyPEozadTosuYpDyPtRTKIuft96iKeX9KTfGLBgMYn5+HlqttqAHVKzAUMGUhpQ+TLPZjEuXLuHs2bO4ceMGEokEfvvb3xZsYwUCAfT19fH/drvduHXrVsFx8h+j0+lgs9mwvb2NmZkZMAyDa9euYXNzE8888wy++93v1vAqhdGygplf2FPqi7q9vY0HDx7I6tijVIRJhigPDg5icnJSsao0qSnZUinYUseud4QZDocxOTnJ3wjIGc0U0yixzH/+eotmrd/D4jFm6XQa4XBY8hgzNY32KkbNgil1bfnDo41GI65evSrbmrLZLH7zm9/gzp07MJlMeOqpp3D58mU89dRTsj1HKVpOMIUU9lSb2iEVuaOoek/wELv+fKMEMoS6HPUUTI7jsLi4iNXVVQwNDUGn0yGXy5WMZpxO555hvmLX2mixJNRTNJVok9Dr9eju7kZ3dzcA8WPMGt26UQk1C6ZUpx/iZ1sOl8uFpaUl/t/Ly8t7nL3IY9xuN+861dnZCbfbjc9//vP8zdT169fxySefUMGUE6GFPT09PZILeyohZ4S5s7ODycnJilGb3IhJyZIULElnVzvhlOqXLOVRS9bl8XgAPOqpLRXNBINBviGfXIydTqfg90AtQplPvUSzHuIkdoyZmgVTzWuTI8Ishcfjgc/ng9/vh8vlwtjYGN55552Cx4yOjuLtt9/G448/jvfeew9Xr17lU7E/+MEPEI/Hodfr8atf/Qrf/va3Ra9RLC0hmEIce9LpNO7du1dzYU8l5IgwOY7DwsIC1tfXq0ZtciNU1EgK9ujRo3v2hcuh0WhqHn5dilIetfk3GeVuYPR6Pd+Qn38xnp+fx87ODoxGIziOg9PpLNmDpkaxJDS6elYJhIwxM5lMSKfTSKVSso0xkws1p4ulRpiJRKLi9Umn0+HNN9/EtWvXwLIsnn/+eQwODuK1117D8PAwRkdH8cILL+C5557DwMAAnE4nxsbGAAAOhwPf+c53eCvS69ev40/+5E8kv0ahNL1gVhvwTAp7WJaFx+NRtDS51j7GdDqNiYkJmEwmjIyM1P0EE+NVK1bMlUzJAsDKygoWFhYkVecWX4wXFxf5YiHSg0b2P6dOflH216AESotmoyOmUmPMtre3MTc3J+sYs1ZAaoRJovxKXL9+HdevXy/42RtvvMH/v8FgwLvvvlvyb5999lk8++yzotdVC00vmOSCWEosSWHPwMAA0um04ie4VqtFOp2W9LdkrSdOnCjZ01QPKqVkiVWg0BRsMUoJZi6Xw9bWFhKJRNl9XrGfu0ajgcFgwKFDh3DkyBG+HcI78MdyLbsuKOkK1GjBLEaj0cBiscBsNuPcuXOyjTFrBaR+ls02CxNoAcEE9l7oSxX2LC0tKerzCkjbw8x3FxJahKTUxapcSpbsp4pJwRajRFtJPB7H9PQ09Ho9hoaGFLuAa7XafSeW+TRjirYU+ecFHWOmPLFYTDG7y0bREoKZTzQahdfr3VPYUw9TAbGiEI/HMTExga6uLsFFSOQ5lDjBS3nVLi4uYmVlpeb9VLkjTOJ2dOTIEUQiEVnFMn+tat6rFIPcoqm2CBMAvy1TCqljzFoBqZ8j2TduJlpCMElkRGYtnj17ljcVJtRDMMU8x9raGubm5nDmzBn+LlgISgsmEQoyLqy9vV2WcVxyCSbJHkSjUXg8HsTjcezu7tZ83FI0i1gS5BRNNQqmmMIaoWPMHA4HLBaL6l6rnEg9L4XsYe43WkIwSbGMXq/HlStXSl7clfZ5BYRFmCzL8tMARkZGRBcjKPk6yI0HScEeO3ZMNlMHOQQzlUrh3r176Ozs5N2OEomE7Hujc+e+JOvx1EQzp2drEfFyY8wWFxcRjUZhMpmacoxZLecO3cPcp8zOzuLQoUPo6ekp+xg1RJiRSARerxculwt9fX2STjolJ6IwDIN4PI6pqSnZW1pq3cMkI9eKDefVPIpLrcghmmqNMOVaExljRgzG4/G45DFmap7XWku7SzwepxHmfmRwcLCqiDRyDzO/HePs2bPo6OiQ/BxKGQBkMhlMTEyAZVlFWlqkChvHcfD7/dja2ipZFCWnYLaCWBJqFU01CqZSbjq1jjFTew+m1LVVc/rZj7SEYAqhXoJZ/BxkL1Cv18uyF6jVamVPQe7s7MDr9eLYsWNIJpOKXXTErpuIuNlsxvDwcMl17ZdRXGqk2dKz9RJxsWPMAPXOba2lHiKRSBT0PDcDVDD/C6WHOwN79xfJ8FU5Dd7lfB0cx+Hhw4dYW1vDxYsXYTKZ4Pf7ZTl2MWInoRARHxgYqJhqB2rbh2lVsSRIFU01RpiNWlO1MWbEsnN7e1sR4/9akOryA9CUbFNTzz1MjuMwPz+Pra0tXojkQq5+RhL5GgyGurgKCU0l56evhbx3QtZd6iLa6kKZjxTRpIJZnuIxZpFIBNPT07KNMZOTWtLYVDD3KUJOEq1Wi1Qqpfg6WJbFxx9/DJvNBo/HI/vJIIdgkrFXcka+1RCSOiXTWTQajaj0tdgIk4rlXpohPVupD7ORaDQaGI1GPPbYYwBqH2MmJ7VEmOl0WvYh9o2mJQRTCPWIMLe2thCPx3H69GnFJoTXIpilUrD1oppgRqNRTExM4PDhw3tGANVy3GKoWJZHjGiqJZrLR63FNcVRXKUxZpFIBEajseIYMyXXJhY1vt+1QAXzv1BSMHO5HKanp/lxN0qJJSBdMEkBjdFobJixezlhW11dhd/vx7lz50RXEAsVTCqUwhAqmmoVTLWtCagu5ELGmJHRcwaDQdbXKLVKVglfaDXQEoIpNCWrhGCSGZsHDx7EqVOn8Nvf/lb258hHimCSFKyQAhqlKLWHmcvlCkwcpPj8ChHMDzs9oo/byuzX9KxaU7JiREnIGLP8HtBax5jV6hqmxve7FlpCMIWghGCSkVJKztgsRoxgktmaGxsbdU/BFlMsbIlEAuPj4+jp6cHp06cln3iVBPMXlvOSjkmpLppqjObUmpKt5b0qNcYsGo0iGAzKMsaslj3MZoQK5n8hp2CS4hQAJSMjJS8mQgUznU7D6/XCaDQqUnwklnxh29zcxMzMjGgf3WrHzYeKZe382jmMS4v/X8l9NCqYwpHTUEGj0cBqtcJqtaK/v3/PGLNcLsfPbhUyxkzq2pqx4Aeggskjl2ASn9UjR46ULE7Jn8+pBBqNBplMpuJjak3BKrF+IvQ+nw87OzvweDyynHClBJOKpXx8cvgP0f4//29YLJaCfTQ1okYRB5RzIAJqH2PGsqyk4drxeBxGo1GW16AmWkIw67GHSSpMV1dXK/qsEmMBpU6QSu0xcqRgSXGO3BeedDqNaDSKrq4uXL58WbbjFx+HiqX8pL74Es48/H8LrOBSqRQ2NzfhdDolXXCVQK17mEoKZjFix5hJ3cNsxkklQIsIphBqEUwyDcVkMuHKlSsVv/xKT0Up5/RD1mg2m2tKwZJIUM4TnDge6fV6DAwMyHZcoDDCpGKpHL878t/w+eDHvBXcrVu3EIvFsLy8DI7jYLfb4XQ6YbPZGrYn1gopWbFUG2PGsizsdjuMRqOoMWbNOKkEoILJI9aajbC9vY0HDx7gscce4/umKqF0v2epalMiSELXKPb4UsmPeC9duoTf/e53shy31PNQlIcUAmk0Guh0Ohw7dgzAozRgOBzG1tYW5ubm+CjG6XTWtRG/FVOyYikeYzY1NYW2tjbRY8xohEkpIJfLYW5uDuFwuOSUjHLIZV1XjvwItliQ5NhTUMJ6T8miIxJh0uiyPpSqntXpdOjq6kJXVxeAz6IY0ohPLsJOpxNGo1ExUVOzYEppmaoHDMOgu7sbVqtV1Bgzuoe5j5H7JCEtD11dXRgeHhZ1fKUjTCJocqVgSx2/1ohtd3eXn36itPWeGi+Qzc6vncMw/q//p+zv86OY/Ivw3NwcP+HC6XTK0keYD93DFE/+2oSMMTMajbhz5w6sVmvFCPP999/HK6+8ApZl8eKLL+J73/tewe9TqRS++c1v4u7du+js7MSNGzfQ39+PhYUFnD59GidPngQA/N7v/R5++MMfKvcGFNESggnINxdxbW0Nc3NzklselI4wNRoNEokE7ty5I0sKtphaUrIcxyEQCGBpaUn2AdSViP/x/1GX56F8RuLa/wkIMDcodRGORCJ8H2E2m4XNZoPT6RTUBlEJte5hqjXyBSqbKpQaY7aysgKfz4df/epX2N7eRjabxdWrV/EHf/AH/CgzlmXx0ksv4YMPPoDb7YbH48Ho6CjOnDnDH/utt96Cw+HA7OwsxsbG8Oqrr+LGjRsAgOPHj+PTTz9V/sWXoGUEUyjlvrwsyxa4zkit/FMywuQ4DisrK9jd3cXjjz+uSEpEquCzLFvQm0qboZsfKY5ADMMU9BGSUVjBYLCgDcLpdIqe5KFWYVJ7hCn0XNVoNHC73fjbv/1bjI2NYW1tDU888QR++ctf4gc/+AH+4z/+AzabDbdv38bAwAC/x/3MM8/g5s2bBYJ58+ZNvP766wCAp59+Gi+//LIqahGoYOZRrkcyEonA6/XC5XKhr6+vppNOqQiTpGDb29v5qjYlkJKSjcViGB8fR19fH9xutyLrKsePf9uOA3V9Rko+tdroFY/CymQyCIVC/CSP9vZ2XkCrGZFTwRSP1Ba4eDwOu92Op556Ck899VTB7wKBAPr6+vh/u91u3Lp1q+xjdDodbDYbtre3AQB+vx8XL16E1WrF3/zN3+DJJ58UvT6pUMHMQ6fTFXxBOI7D8vIylpaWJBl/l0KJQdX5VbBWqxWTk5OyHj8fsSnZtbU1zM/P4+zZs7BarYqtqxwHvnCy7s9JKURO79m2traCSR7EiPzhw4eIRqMwm838/mfxTaNahUmt6wKke8nG43FFfKl7e3uxuLiIzs5O3L17F1/+8pcxOTlZt2tLywimkD1MImZtbW3IZDKYnJyETqfDlStXZEshytmHyXEc/H4/Njc3+SrYdDqt+B6pkOOTCS3JZBIej6chzes0ulQPShm2G41GGI1GHDp0iDciDwaDmJmZQSqVgtVq5QWURpjikVooVWl4tMvlwtLSEv/v5eXlPa5o5DFut5t3J+rs7ATDMHwh2OXLl3H8+HHMzMxgeLg+04ZaRjCFQPYXScSmxABluSLMdDqN8fFxdHR0FFTB1qOoqNqNB6ki7u7uxqlTpxp2kaLRpbpQespJvhE5KULZ3d1FMBjE8vIyYrEY3yZht9tVs4+uZsEEpFWaVzIu8Hg88Pl88Pv9cLlcGBsbwzvvvFPwmNHRUbz99tt4/PHH8d577+Hq1atgGIZ3j9JqtZifn4fP5+P3QusBFcw8NBoNHj58iEgkotj0Dq1WW9XrtRrBYBBTU1M4ceIE79BBUFowq6Vkt7a2MD09LYtxei3Q6FKd1HM0mEajgd1u56sz7927B6vVimAwiPn5+T02cI0SLbULphQqRZg6nQ5vvvkmrl27BpZl8fzzz2NwcBCvvfYahoeHMTo6ihdeeAHPPfccBgYG4HQ6MTY2BgD49a9/jddeew1tbW3QaDT44Q9/yO9v14OWEcxqd0nJZBLBYBBOp1PRRnqtVotkMinpbzmOw/z8PLa3t8uaJcjRJ1mJcsfnOA5zc3MIhUIYHh6WtX9OCjS6VC+NmqfJMAw6Ozv59F+xDZzBYODTt5VcbORGre0utRCPxysGHNevX8f169cLfvbGG2/w/28wGPDuu+/u+buvfvWr+OpXvyrfQkXSMoJZCTJOym6349ChQ4p+eaVGgKlUChMTE7BarRgeHm7YCVZq/SQ9bLPZRBs5FCPHPhONLtVPI0Sz+LtVbKBACoiKXWycTqeiN4BqNVQApBt/VIow9zMtLZi5XA4zMzOIRqMYHh7G0tKSoi48gLQ+zEop2HpTnJIl+71yrE2uSSg0utwf1Fs0K323GIaByWSCyWSCy+XiDRTI95sYKJAUrpxWdmpNydaSqaKC2WTEYjFMTEygp6cHJ0+eBMMwitvWAeIiTCEp2HpD1k/Gma2vr8vmU0vEWOrFg8zTpOwf6imaYiK5fAOFI0eO8AYKpIUlf46kzWarSfDUKpi1jCFMJBJUMPcz+SfKysoKFhYWMDg4CJvNxv+8XoIp5DnUkoItRqPRIJvN4t69e9Dr9bLu99ZiX5hKpTA+Po6d//aCLGuh1I96iWYte4WlDBTC4TA2Njbg8/mg1+v5/U8xY7DIutSYkpXagwk8CkgsFovMK2o8LSOYwKMxQ/fv3wfHcRgZGdmTVpGjgrUaQvowyciwkydP8hMe1EIqlcLy8jJOnjyJ3t5eWY8tVTDD4TAmJydx4sQJ3JN1RZR6UQ/RlFOYiudIJpNJhEIhfgyW2WwumMBSDTUKZi0RZiqVUkVGTG5aRjB3d3fx6aef4vDhw3C5XCW/oLVUsAqlUoRJKk2DwaBqUrD5BAIBBAIBHDx4UHaxBMQXRHEch6WlJaysrODixYv4z+7fk31NlPqhtGgqGckZDAb09vait7eXN1AIhUIFBgokhavX6xVZg9zUEmE2Y+Uv0EKCqdFocP78+YppAjldeMQ+B0kp2u32mlOwte4FFsOyLKampsCyLAYGBpBKpWQ5bjFiIkyWZTE5OQmGYeDxeFTThE6pjXqYGyhNvoFCX18fcrkcP4ElEAggl8vBbrfD4XDwPaJqROo1RA0m6UrRMoLZ0dFRNd2q1WqRzWYVXUepfVK5U7AkUpNDMOPxOO7du8cbz29sbCCRSNR83FIIFUyyJrfbDbfbDYZh6IDoJqJRfZpKodFoYLPZYLPZcPToUd7qjRgoxONx+P1+OJ3OhhooFMOybE03ompMM9dKywimEOoRYeanHfOb/eVMwcrl9rO+vo7Z2VmcPXuWL45S0klIiLE76ZkdHBxU9d05pTaaTTTz0el06OzsRGdnJwDg1q1bMJlMBQYKZP+zngYKxdSyh9msUMHMo15VshzH7UnBynlS1CpqpD81Fovtmf2ppJNQpWPn31x4PJ6CfSAaXTYnzSya+Wg0GvT09KCnp6fAQMHv9yMWi/EGCg6Ho651DVL3MLPZbNNukbSMYAoRJCVGb5Uim83i448/VqwKthbBTCaTGB8fR1dXF9+fmo/Y8V5iKJeSzWQyGB8fh8ViweXLl+ldbwvR7KJZ/H0vZaAQjUYRDAbx4MEDZDIZfgKL3W5XdApQLbMwm7EHE2ghwRQCmYepFBzHYXZ2FqlUCk8++aRid4tSBZPspZ4+fbqsobHSKdniC8ju7i68Xi+OHz9ecr4ejS6bn2YWzWq1BgzDoKOjAx0dHThy5AhyuRy//7m4uAgAsNvtcDqdNRsolFqb1FmYSg2wbzRUMPNQMsJMJpOYmJiA3W6HyWRSNLUipT1DqKOQ0inZ/HWvrKzg4cOHGBoaato7VoowmlU0xRbnaTQaPj0LfGagsLm5idnZWeh0Ot5gQayBQq1rI9AIs0VQag+TjLw6deoUOjs7sbm5Kftz5COmeCmdTmNiYgIWi0VQO0s9UrK5XA4PHjxAOp2Gx+Mp69tJo8vWohlFs9Zq9mIDhVQqhWAwiKWlJUQiEZhMJt6ByGg0ihJQlmUlpXxjsRiNMPc7Qr4oclejkRRsOBzeM/JKySZqoaJGHHIee+wxdHd3Czq20inZVCqFO3fuoKenB6dPn27K0vRWZ/ODackG+c0mmnL7yLa3txcYKMTjcQSDQczOziKZTKKjo4MX0GoGClL3MJvVRxZoIcGsN6R4xul07qmCJZGUUmJQLbXMcRwWFxexuroqelC2kinZZDKJ9fV1nDt3rupQWBpd7m82P5jm/1+seP7aOYzND6bx1csRuZdVd5S+cTabzTCbzSUNFFiWLdj/LM7kSN3DjMVioq4p+4mWEsxazL3FUJyCLYakfpUcUl3udWazWXi9XrS1tUlyyFEiJctxHBYWFhAOh/kJ65TmJF8oi38mRjgPfOEkfiJCNNXqPlPPSSXFBgosyyIcDvMtLPn7o1arle5hlqClBFNpcrkcZmdnsbu7uycFmw8RTKVKwsulTSORCCYmJtDf349Dhw7JemypEAHX6/VwuVyC3hMaXTYnYoVTjGiqdUhzI0d7abXaAgOFdDqNUCiEtbU1zMzMIJPJgGEY6HQ6mM1mwe9fPB6nEWYrISVNQlKwnZ2duHz5csW/V3IfkBy/OCVLRpqdO3cOHR0dNR1brrv1aDSK8fFxXsDn5uZUGwlQaqdUdFnpcUKEU6hoqtUMXE2zMPV6PW+gAAD37t2DTqfDwsICP66L7H9WqqSne5hNgpCULIn+xExUJ3ZtlfoXSz2HUuQLMsuyfMNzqZFmYpErJbu2tob5+fkCARcixjS6bB2ECqcQ0VTzzEm1CGYxDMOgt7cXBoOBN1AIhUJ8BbvNZuNTuPmZoVgsVnIrqhlQ5yfVQMSIGbGQe/jwIYaHhwXvvdUjwszlcojH47hz5w4sFguGhoZqFsv8Y0sll8thenoaKysr8Hg8BdGuki0rlMYiNLos97fkv3Ic+MJJ/ORu+cwJFUzx5NdZEAOFw4cP48KFCxgeHkZ3dzefJfr444/x0Ucf4Wc/+xnfzlKO999/HydPnsTAwAC+//3v7/l9KpXC1772NQwMDODKlStYWFgo+P3i4iIsFgv+4R/+QdbXK4SWijCFIFQwxaRgi1Hagk+j0SAUCmF5eVl2k/JaLjrEP9fpdOLixYslbfcqRZg0uqRUijorRZpqFSa1rguoXCVbbKCQzWYxPj6On/3sZ/jtb38Lm82GxcVF/NEf/REuXrzIH4dlWbz00kv44IMP4Ha74fF4MDo6ijNnzvDHfuutt+BwODA7O4uxsTG8+uqruHHjBv/773znO/jiF7+o4Csvjzo/qQYiRDA3Nzdx9+5dDAwM4Pjx46JFRMmpKLlcDmtra7xJuVomeoTDYXz88cfo7+8v+57Vq4qZUl9qiS4rHbPUcctFmjTCFI+Ytel0Oly6dAk/+tGP8Kd/+qf4q7/6K7jdbrz55pv4yle+wj/u9u3bGBgYwLFjx6DX6/HMM8/g5s2bBce6efMmvvWtbwEAnn76afzyl7/krws//elPcfToUQwODsr0KsXRUhGmkBOmkpjlcjn4fD5EIpE9EzPEoFRKNpVK4d69e2hvb8fBgwdVMdmd4zgsLS1hZWWlas9npfeFRpeUUpSKOEtFmmoVTLWuC5C+tng8jsOHD+OJJ57As88+W/C7QCCAvr4+/t9utxu3bt0q+xidTgebzYbt7W0YDAb83d/9HT744IOGpGOBFhNMIZQbIl1LCrbUc8idkg0Gg5iamsKpU6fAcRy2t7dlQ9LBwQAAIABJREFUPb4UWJbF5OQkGIYR1PO5n/Ywa2m8byWUiC4rPQ/5LIho/vGJFVgsFtUKk5ojTKnvl1JtJa+//jq+/e1vw2KxyH5soVDBLKJUhCm2CrYacu5hchwHv9+Pra0t3jg9GAw2XHji8Tju3bsHt9sNt9st2JqwVEpWTdFlucZ7KprqoPhGZvHDtxGNRtHe3o5sNotkMlnXmZLVyOVyshTjKYHU7ZFEIlFW1FwuF5aWlvh/Ly8vw+VylXyM2+1GNpvFzs4OOjs7cevWLbz33nv47ne/i3A4DI1GA4PBgJdfflnSOqWgzk9KIYSmZImYkRRsNBqtKQVb6jnkELRMJoOJiQmYTKYC43Slq3CrQW4wxBYcKWm7VyvVIiUqmnupV3RZ6fk/PxgBx3FYX1/H8vIypqamkM1mYbPZ+JmSjRQsNUeYUonFYmX7MD0eD3w+H/x+P1wuF8bGxvDOO+8UPGZ0dBRvv/02Hn/8cbz33nu4evUqGIbBhx9+yD/m9ddfh8ViqatYAi0mmEIg0V8ikcD4+DgOHDiAS5cuyZrO0Wq1SKVSNR1jZ2cHXq8XAwMDe+ZEKllUVAmO4zA3N8cXHIm9wSgVYTY6uhRz0aeiqT5+crcDX70cgdFoREdHB06ePAmWZfmZkgsLC3zFp9PphNVqrWvqVq2CWcuNayXjAp1OhzfffBPXrl0Dy7J4/vnnMTg4iNdeew3Dw8MYHR3FCy+8gOeee463yRwbG5O8FrmhglmETqdDOBzGwsICzpw5w5dNy0ktESApogkEAmWLaBiGUbRtpRSZTAbj4+OwWCy4fPmypIuAWvYwa+0ZpKLZ+OiymHxrPK1Wy8+MBB5ZwgWDQaysrODBgwf8SCyn06n4mCo1C6bUdSUSiYp7mNevX8f169cLfvbGG2/w/28wGPDuu+9WfI7XX39d0tpqhQpmHrlcDuvr64jH47hy5YpiVaZS9zCz2Szu378PjUaDkZGRskU0lczX5SK/iGJ3dxderxfHjx/fE+2KoTjCrHd0KddFnoqm+qgkAHq9HgcPHsTBgwcLRmLNzMwglUrBarXylnBy+z+rVTBZlpU0qQSQPuVkP9BSglkp1UJSsEajUfGWDCkp02g0iomJCRw+fHjPJnkxSu9h5o8nCwQCWFxcxNDQUM3+kY3aw1SqT5CKpjr4yd0O/GH/tuDCs+KRWLu7u/xQZgB8w77NZqtZ7KjH7f6ipQSzHBsbG/D5fDhz5gxYllW8JUNshLm6ugq/3y/YOL0eTkLZbBazs7NIp9PweDyyFE7kp2TrEV0qnTZsVdFUWzoWkN5TqNFoYLfb+eK1TCaDUCjEXzPa29v59K3JZBL9HEqO+asFqRGmWov25KKlBZN4wcZiMb5IJRQKKb7/JzTCzOVyePDgAVKplCjjdKUjTI7jcPfuXfT29uL06dOyFUnUy+mnnhf0VhVNtSHXeK+2tjZ0d3eju7sbwKPMVDAYxPz8POLxODo6OngBFZKlatb+UDW+JjloKcHM/xBJCra7uxsnT54sKAhQWjCFRIBkfT09PaJFSUnB3N7eRiQSwfnz52varywFEUwlostGRj2tJJpqjC4B4MOl4xju8sp+XKPRCJfLBZfLBY7jEIlEEAwG4fV6kcvlYLfb4XQ6YbPZSkZsak19So0w1Tp3VC5aSjAJ+SnY4irYeghmtQiT9DFKrdJV4gvLcRwWFhawsbEBu90Oq9Uq+3MosYeplgt4K4mmWlH6Qs4wDKxWK6xWK/r7+5HNZhEOh7G1tYW5uTnodDo++rRYLPwWhBoFU+q6qlXI7ndaSjBJijMej5ftE6yXYJZ6Do7jMDs7i52dHVmNEmolm83C6/VCr9fD4/FgfHxckQiWYRis/d7/Lsux1CKU+TS7aKrxPc+n3pGPTqdDV1cXurq6ADyy1wyFQlhcXEQ0GoXFYkEikUAmk1GV+xAgPcKMxWJUMJuFXC4Hs9lckIItpl4p2WLBIaOvHA5HzV61ckLm3fX39+PQoUMAlKtmleM1q/2i3eyiqWZub5xBX1/5IdNKYzAY0Nvbi97eXn4g88TEBHw+H1iWVY37ECA9wqSC2US0tbUVOOWXohF7mKFQCPfv38fJkyf5u1E1sLa2hvn5+T3VuUoZDHx08AlJf6d2kSymGUVzv30GjYYMZNbr9bhw4QI4jlON+xAgPcKs5PLTDLSUYAqhHr2A5M4tf1/w0qVLiruKCIV46JLq4eJmbbV4vu7ni3QziiZFPKRKVqPRqMZ9CJAeYcbjcdVcx5Sg5QRTLUOKOY7Dp59+CoPBAI/Ho5qNf5IadjqduHjxYsk7WyWqcMVUxu5nocynWURzP30exFtWTZQ6xxrpPgQ8ijCl1FDE43EaYVLkZXd3F7FYDMePH8fBgwcVex6xPV7hcBiTk5M4ceIEDhw4UPZxjfJ83U8XZqE0i2hSlKPe7kMA3cMsBxXMOsJxHAKBAJaWlngLPqXIt68Tsq6lpSWsrKyUNXTPR+4Is1p02YxCmc9+Fs1m/2zUiNLuQwDdwyxHywmm0JSs3A4cLMvi/v374DgOIyMjuHXrlmzHLgXp9ax2l8iyLCYnJ8EwDDwej6CTpB57mK12Id7PornfUGNathbkdh8CatvDpBFmiyEmOhNCLBbD+Pg4+vr64HK5+OMqaYslJAqMx+O4d+8e3G433G634LXImZItji5bTSjz2W+i2cqflZqp1X0IqK0PU4jf9X6FCmYJdDqdbKbIpDXj7NmzBe44RNCUGoNTzX6PuAkNDg7yqR0xx5Z7D5NefB+x30STom6kuA8BtTn99Pb2yv0yVAMVzBIQsaml+iyXy2F6ehrJZLJkawbp91RSMEuJGsdxmJubQygUkuwmJFel8VNPf4QXqVDuYT+I5n6/wVFDWrYR1fpC3IccDgfS6bSkaxNtK2kyhKQdazUvyDd2P3XqVN1aM6odP5PJYHx8HB0dHRgeHpacDlZ67ZT9IZqU2lCDj2wp96FgMIhIJMI7j4lxH4rH47BYLHVYeWNoOcEUQi2CubW1henp6arG6fWYWZkvaru7u5iYmMDAwEDNU0bkWHssFsOL/9cf1HSMZketornfo0u1oLbRXsR9qKOjA9vb2xgcHEQsFhPlPtTsRT/q6JZXGVIEkxin+/1+DA8PV50yInQmplTyBTMQCMDr9eLChQuyjOSqtUp2Y2MDo9+aqHkdrQAVJ+X4yd3GFqeoIcIsRy6XQ1tbG5xOJwYGBjA8PIyzZ8/CaDRiZWUFt2/fxsTEBAKBABKJBP931dpK3n//fZw8eRIDAwP4/ve/v+f3qVQKX/va1zAwMIArV65gYWEBAHD79m1cuHABFy5cwNDQEP7t3/5N9tcshJaLMJVIyabTaYyPj8NmswlOddYjwsxms7h//z7S6bSoAdTVkFoly3Ec5ufnEQwGZVlHq6CmSJMKuHyoXTCL11bJfej999/H1NQUYrEYstlsyWOyLIuXXnoJH3zwAdxuNzweD0ZHR3HmzBn+MW+99RYcDgdmZ2cxNjaGV199FTdu3MDZs2fx8ccfQ6fTYXV1FUNDQ/izP/uzupvUq/PTajBiBDMUCuHOnTs4cuQIHnvsMcEpFqUjTFJ0ZDKZMDQ0JOsXS8oeZjabxaeffopsNovZ7O/LtpZWgQpV86FmwawGcR/q6+vD0NAQXnnlFXz961/H5uYmXn75ZXzuc5/DX//1X+Pu3bv839y+fRsDAwM4duwY9Ho9nnnmGdy8ebPguDdv3sS3vvUtAMDTTz+NX/7yl+A4DiaTib+GJZPJhqWy9+enpTBCBJMYp09PT+PSpUsVreRKoWSEub29jbW1NfT29qK/v1/2L5fYlGwsFsPt27dx8OBBnDx5Ev/jv/9K1vW0ClQ05aeRadn9LJjFGAwGXLt2DSaTCT//+c/x7//+77h48SL8fj//mEAgUDAtyu12IxAIFBwn/zE6nQ42mw3b29sAgFu3bmFwcBDnzp3DD3/4w4aMQGu5lKwQqglm/kDlkZERSV96JSLM/Oknbrdbsc13MSnZjY0NzM7O7ulDpUijkelZKtjy0kyCSSBFP3q9Hl/5yldkPfaVK1cwOTmJqakpfOtb38IXv/jFug/ebq5PSwC17mFGIhHcvn0bPT09OHPmjOQvvNxzN7PZLO7du4dEIsH3VyqV8hWSkiVFUIuLixgeHubF8qmnPwIAWiFbA1S4moNmFMxK/esul4s3iweA5eVluFyuso/JZrPY2dlBZ2dnwWNOnz4Ni8UCr9cr8+qr01yflkyUE7NAIICJiQmcP3++ZjcLOXsZo9Eobt++je7ubl7EleyVrJaSzWaz+N3vfgeWZXHp0iVJ5giUytRbNJtZpBuVllWrYCrV7uLxeODz+eD3+5FOpzE2NobR0dGCx4yOjuLtt98GALz33nu4evUqGIaB3+/ni4kePnyIBw8eoL+/X/Y1VoOmZEtQLJgsy2Jqagosy8pWbSpXhEms986dO1fg4ajVapHJZGo+fikqpWSj0SjGx8dx9OjRPTcVP/5tuyLraVXUVD273wkEAnUbzkxQWx8mQaotaLW6Bp1OhzfffBPXrl0Dy7J4/vnnMTg4iNdeew3Dw8MYHR3FCy+8gOeeew4DAwNwOp0YGxsDAPzmN7/B97//fbS1tUGj0eCf/umfeLeiekIFswT5YkYMyl0uF/r6+mT7gms0mpoELZfLwefzIRqNlrTeU3KiSLnolexXFos3gRb7yE89RLOZo0tCLpfjhzMTc3Kh7ja1PKcaI8xaPa4rXSOvX7+O69evF/zsjTfe4P/fYDDg3Xff3fN3zz33HJ577jnJa5KLlhNMoXuYuVwO6+vrfMGKzWaTdR21VMmmUimMj4/D6XTi0qVLZa33lKrCLfaSJf604XAYw8PDNAVbZ2ikWTu3N87gq5cjyOVyCIfDvLuNVqstaU4uB2oVTKke143wxq03LSeYQmAYBsFgEMlkEiMjIzWZsJdDapVsOBzG5OQkTpw4UbGVRek9THLsbDaL8fFxmM1mXL58uewFhRT7ALTgRwmUEs1WiC7z0Wg0vEACj25Og8FggTm50+lEZ2dnzTeGahVMqetKJpN1r1qtN1Qwi0gmk/xA5XLRmxyIjQA5jsPS0hJWVlZw8eLFqi0j9RDMSvuVlPpDI035aW9v32NOvr29vWe2pN1uFy0yahVMqXuYsVisqX1kgRYUzEoCuL29jQcPHmBgYABLS0uKbsiLiTBZluVF3OPxCEqXKC2YyWQS4+PjZfcr86HFPvVDTtFsteiy2sivfHPy4tmSs7OzaG9v56NTk8lU9fqRy+Ua0nxfDal7mPF4vKKPbDOgvk+rARCP0+3tbVy+fBk6nY43/VUKoVWypOjI7XbD7XYLFnGlBJPjOPj9fiSTSTz55JOC0lK02Ke+0EizPhTPlkwkEggGg5ifn0cikYDVaoXT6YTD4Si5rdNsEWazTyoBWlQw84tW0uk0JiYmYLFYMDw8zFeXKmmMDggTtM3NTczMzEgqOlJCMMl+pclk4t08KOqkVtFstehSDoxGI1wuF1wuF3K5HHZ3dxEMBvlGfIfDgc7OTnR0dPDnpxoFU2qEmUgkqGA2M6SA5rHHHkN3dzf/83r0RlWKMEnVaSgU4l17xCK3YJL9ymPHjqGnpwcfffRR9T9CYbEPpb7QSFM81dKyQtFoNLDb7bDb7QAeDW8PhUJYWVlBJBKB0WhELpdDe7v6tiukVslGo1EqmM0Ix3F4+PAhVldXBRXQKEE5QctkMhgfH0dHR4fgUWFiji+FUv2VUkvIaYVsfZEimjS6lJ+2tjZ0d3eju7ubH401PT2NQCCAQCAAm83Gp29r6YGUA6mRb7VZmM1ASwqm1+sVVUCjBKWqZHd3dzExMYGBgYGaBz3LIZjED3ZnZ0dSfyUt9lEHNNJUF2Q0lsViQXd3Nzo6OrCzs4NgMAi/3w+dTqdY76cQpEaYdA+zSRkYGGh4v1BxlWwgEMDDhw9x4cIFWe7Sap2Gkslk+L3dSv2VlaDFPupBqGjS6FK+tGw1SCSXb44A7O397Ojo4H9fj7qBXC4n6XlolWyTYjKZBImJkl6P5Li5XA4PHjxAOp2WzacWqC3CzN+vPHjwoCzroTQeGmmqi3KpTyG9n52dnbDZbIoUDdVSJUsFs0UhRTlKe0neuXMHPT09OH36tKziLGZmZT7r6+uYm5sT1F9ZCVrso04qiSaNLuuLkL3Ccr2fGxsb8Pl8fO9nZ2cnjEajLNeQWvowi0dxNRstKZhiZmIqJZjb29uIx+MYHh7mUzFyIvbEIfuVu7u7Jc3cKc0DjTSrU4+0rJTimnK9n3Nzc4J6P4UgdQ+TtpW0MHIPeCZwHIeFhQVsbGzAZDIpIpZiya/MlcMOsFyxD62QVQ/Fokmjy/rDcVzNKdVqvZ9k79NqtQo+r6VWycZiMZqSbVWUEMxsNguv14v29nZ4PB7BvYxKosR+JS322R/QSLOx5HI5WbdhSvV+BoNBrKysYHp6mr9BdzqdFYsepe5h0raSFqbWKtNiiDD19/fj0KFD/M8bOUS21v1KtQ7ApQjnb/+3HwEAXmzwOtSI0mlZpZ1+2tra0NPTg56eHr73k/hlZzIZ2Gw2dHZ2wm63F6Rga9nDpCnZJkToHmY2m5Xl+dbW1jA/P79HmEgvZr0NmOXYryT2gsXvJS32oVCEUU9rPNL7aTabcfjwYbAsi52dHWxvb2N+fr6g9zObzVLz9TK0pGAKQY4IM5fLwefzIRqNlhQmuaNYIci1X6lmL0yKeP7Hf/8V3WNuAI3K0FTq/dzd3cXMzAw6OztF9X62gmDSq10Zat3DTKVSuHv3LnQ6HS5dulQyihM7E1MK+RZ20WgUd+7cgcvlwokTJ2o6WYlJfT6VnH3oxZiyH/nJXemtVfsJ0vs5ODgIk8mEvr4+JJNJeP9/9t48rqkz7f//JIQ9AZKQILtABATZhKi1y9iOrRYlPs/zc1xmvl2my8xYu810cTrt9LHtdOo8XWap3aa1HTttxbaAaFut1mnrdKZVUVmVTZBdE8IWIGQ55/z+cM4pSCAhOYFA7vfr5avSnJxzI+G8z3Xf131d1dUoKyvjaltP9oBPhDlHcUQUrsisr68PZWVlmD9/PpKSkia8nrsjzNHFCy5dusT1r+QjucfWPk+S7EMgzA1CQkIwf/58LF68GNnZ2QgJCYFWq0VZWRkqKirQ3t6O4eHhMQ/NRqMRYrF4wnMeOnQIKSkpUKlU2LFjx7jXTSYTNm7cCJVKhaVLl3ItFo8cOYLc3FxkZGQgNzcX//jHP3j/fh2FTMlOgEgkmrIwGYZBW1sbOjs7HSrq7s4mz+z5KYpCU1MT7/sr3T12wvRDpmUJthCJRFAoFFAoFAAui1Gv16OxsREjIyOoqKhAYGAgKIqa8P5CURS2bt2KI0eOICYmBmq1GhqNBmlpadwxu3btglQqRWNjIwoLC7Ft2zbs3bsX4eHhOHDgAKKiolBdXY1Vq1aho6NjWr73K/HKCNMRphphUhSFqqoq9Pf3Q61WO5Qt5q69nqOpqKgAwzATTgs7y5VTsiTZZ/ZBZgQcw1umZR0lMDAQMTExyMzMRF5eHhYuXIgzZ86gvb0d119/PZ555hkcP358zL3txIkTUKlUSExMhJ+fHzZt2oTS0tIx5y0tLcVtt90GAFi/fj2OHj0KhmGQk5PD7SxIT0+H0WiEyWSavm94FESYEzCVCHN4eBgnTpyAVCrFokWLHM4wc+ca5uDgIAYGBhAREeHyeqUtnC29RyAQ5g5CoRDXXHMNXnjhBURGRqKoqAjJycl44403sGfPHu64jo4OxMbGcl/HxMSMixJHHyMSiRAaGgq9Xj/mmKKiIixevHjG+oh65ZQsn2uYOp0O9fX1WLRoEUJDQ6c0DnetYbLbWNgize7AVtIPYfZDpmWnB0/93XF1XAqFAhs3bsTGjRt5GtH31NTUYNu2bTh8+DDv53YUEmFOgL3pUnYv44ULF6BWq6csS0euMVUYhkF9fT06OjqgVqvh5+fntihwdIRpbzqW3IAJsx2+p2U9teiHs1vF7Ik2OjqaK9cHAO3t7YiOjp7wGKvViv7+fu6Bv729Hf/93/+Nd999F0lJSVMeH18QYU7AZDKzWCw4ffo0aJp2qrEyC5+JM+yYBAIBt17pzixckvRDIDiPp+5hdrbKj9lsnvQ+qFar0dDQgObmZpjNZhQWFkKj0Yw5RqPRYPfu3QCAjz/+GDfccAMEAgH6+vqwZs0a7NixA1dfffWUx8YnnvcTmwam0q3kSgYGBnDixAnExMTwspeRjwjTYDDg5MmTiImJwYIFC7gxuVNq7LlJss/sZLKEH5IM5H48WZjO9sIMDAyc8HWRSISdO3di1apVWLhwITZs2ID09HQ8+eST2L9/PwDgzjvvhF6vh0qlwksvvcRtPdm5cycaGxvx9NNPIzs7G9nZ2dBqtc59gy7ilWuYjmBLmB0dHWhpaUF2djYvG3T5mJJl1yszMzPH7YFyZ1IRWxqPQPAW+Kwt66nCdLa1lyOdSvLz85Gfnz/m/z399NPc3wMCAvDRRx+Ne98TTzyBJ554YspjcgdeK0x7N/zRMqNpGrW1tTCbzViyZAlvtV+FQiHMZrNT72UYBg0NDTAYDBPur5yOCJNA8Cb4Wnv0ZGGSTiUT43k/MQ+BFebIyAhOnjyJoKAgZGVl8Voo3dk1RlvrlbZwtzA339Nk9ziS8DM7IdOytjl58iTq6uqg0+lcas7Ad2svvnB2DXNoaGjOdyoBvDjCtIdQKITFYsGpU6ewcOFCtzR6dmbK1GAwoKqqCklJSYiIiLB7fncKk0DwNtp8fogFylb09PSgpaWFK2Iul8sRHBzssAT5aB7tDpyNML2htRdAhGkThmFw4cIFmEwmXHvttZM2W3WFqUaYk61X2sKdwrzrkUtuOS/B/TgaPZI9mbaRSqWQSqUAvu/yceHCBQwNDSEkJARyuRxSqXTSylqeOiVLemFOjtcKc6I1TKvViqqqKgQEBCAoKMhtsgQcT/pxZL3SFu4SZmdnJ+/nJBBmI2yXj8jISDAMg4GBAfT09HD7CdnoUyKRjIk+PVmYzkaY3rCG6bXCtMXg4CAqKysxf/58REVF4d///rdbr+eI0Nj+lSEhIVPuX8lOK/MFWxhheHiYt3MSCLONibJlBQIBQkNDERoaioSEBFgsFvT09KC9vR0GgwFisZjrMempwnQ2S5ZEmF4GO92ZkZEBieT7qh7urMhhbw1zKuuVE52fz8IIrLgffnZmCh8Tph8yLes8vr6+iIiIQEREBBiGweDgIPR6Paqrq2EymeDn54e+vj6EhIR4jDydXcN0ZFvJXMBrhclKkKZpNDQ0YHBwcNx0Jztt6y5hTraGOdX1yqmefyoMDQ2hoqICiYmJ/+mlqbf7HoBkyBIILAKBABKJBBKJBPPnz0dnZyf6+vpw8eJF1NfXIzAwkIs+3bkMZA+app2qXGY0GhEZGemGEXkWXitM4PKCfWVlJWQymc3pTrZjibue/mxFmM6uV9qCj44ibHH5jIwMhISEkMo+BAJcL2LACjQ2NhYMw2B4eBh6vR61tbWwWCyQSqWQyWQICwub1uiTZMlOjtcKs7e3F1VVVUhOTuYao14JKzQ++0iO5soI0GKxoKKiAqGhoVNer7SFK1OybKawTqfjCrkTZj/O7K8k07L8M3oNUyAQIDg4GMHBwYiLiwNFUejr60N3dzcaGxvh7+/PRZ/ulpKzWbLeUrjAa4VptVqRk5Mz6QfQ3Q2eRwvN1fVKWzg7JUtRFGpqaiASiZCXl+cx6ysEwlxhsn2YPj4+kMvlXKcOo9EIvV6PxsZGjIyMICwsDDKZDFKp1Cm5TYYrWbIkwpzDKJVKu5U63C1MNoLkY73SFs5EmCMjIygvL0d0dPSYhq+A/TZeBII34cq0LE3TDlcNCwwMRExMDGJiYkDTNPr6+tDT04Pm5maIRCJOrkFBQS7PSjmbJUsiTILbhckwDEZGRtDZ2enyeqUtplpJqK+vDzU1NW6rbESYvZBpWX5xNpITCoWQyWTc76fJZIJer0dzczOGh4cREhLCve5MGU9XsmRJhOnluFOYZrMZlZWVEAgEyMnJcUsm7lQizI6ODrS2tmLx4sU22/RMNbokN1cCYWL42ofp7++PqKgoREVFgaZpGAwG6PV6tLa2cnKVy+UQi8UO3WNcqfTD5+yYp+K1wnSlJ6arsOuVKpUKjY2Nbt3naU+YNE2jrq4OJpMJarWa1+LyBM+CFFTnH2enZd1RuEAoFHKFExITE2E2m9HT04PW1lYMDg5CIpFw0edESXzOjstoNJII09txhzC7urrQ3NzMrVc2Njbyev7R2BMmm5UrlUqRmprqkd0TCJ4DmZblj+mo9OPn54d58+Zh3rx5YBgGBoMBPT09qK6uBk3TXPQZEhLC/e670g+TRJheDp/CZMvKDQ0Njeup6a5fnsmEyZYBVKlUUCqVk56HJPsQCPwy3aXxBAIBQkJCEBISgvnz58NisaC3txednZ2ora1FcHAwZDIZrFarU+MymUwzWnBhuiDCnAQfHx9earGy65VhYWHj1ivZrR/TKUytVovGxkbes3IJBG/EmWnZma4l6+vrC6VSCaVSCYZhMDQ0hJ6eHhiNRpw6dQpSqRRyuRyhoaEOj9Mbtp/N/e9wAqZrDdNgMKCsrAxxcXFQqVTjrutMT0xHubIjC8MwOH/+PFpaWpCXl+eQLJ2JLsm03dyFrIPygztLbk4VgUAAsViMuLg4BAYGIicnB2FhYdBqtSgrK0NlZSU6OjpgNBptvt9W16e5itcK0xFcFWZXVxeqqqqQmZk54bQnX/VebTH6F5KiKFRUVMBsNiM3N5dU7vEyiOg8i5mOMCdDJBIhPDwcKSkpUKvVUKlUoGka9fX1OHHiBOrr66HX68fdGyd6ADh06BA+U+xWAAAgAElEQVRSUlKgUqmwY8eOca+bTCZs3LgRKpUKS5cuxYULFwAAer0e119/PcRiMe69917ev09nIFOyk+CsMCdbr7wSd0aYLEajERUVFdzmZwKBwC9TnZb1ZGGORiAQICgoCEFBQYiNjQVFUejv74der0dTUxPef/99KBQK+Pj42IyaKYrC1q1bceTIEcTExECtVkOj0SAtLY07ZteuXZBKpWhsbERhYSG2bduGvXv3IiAgAM888wyqq6tRXV093d+6TTz/JzaDOCNMs9mMU6dOwcfHBzk5OXa3abgzwgQulwA8ffo0UlNTpyxLkuxDsAWJVl1ntgjzSnx8fCCTybBgwQKo1Wrce++9CAwMRFdXF3JycrBlyxaUlpZiaGgIAHDixAmoVCokJibCz88PmzZtQmlp6ZhzlpaW4rbbbgMArF+/HkePHgXDMAgODsY111zjUclEs+8nxhOOrmFORWb21isnuoa7Isy2tjaYTCbk5uYiLCzMLdcgEAhTxxOFSdP0lNdVVSoV7rjjDqhUKpw8eRKbN2/G8ePHue1yHR0dY0psxsTEoKOjY8w5Rh8jEokQGhoKvd6xFoLTjVdPyV6ZFHMlU5HZlfsrHYXPJs8sNE2jtrYWVqsVQUFBTj2hkeiSQJgaU5mW9VRhOlvlJygoCL6+vrjuuutw3XXXuWF0noFn/cQ8DEeEyVbK6erqwpIlS6a8TYPvNUx2SjgwMBAZGRnTnolHMmS9AzIt6zqekiXL4o5OJdHR0Whra+O+bm9vR3R09ITHWK1W9Pf3c51aPA0izEmwJ0yz2YzTp087vF450TX4ijDZKeH58+cjISEBAoGAlybShNkNkRvBEZyt8jOZMNVqNRoaGtDc3Ayz2YzCwkJoNJoxx2g0GuzevRsA8PHHH+OGG27wuIcJFq+ekrXHZNHfwMAAqqqqsGDBAruVcpy9xlRgW4RlZWWNabPDTvlO5cmRTMcSCM5x849P4+AHi2d6GE7hjghTJBJh586dWLVqFSiKwh133IH09HQ8+eSTyMvLg0ajwZ133olbbrkFKpUKMpkMhYWF3Pvnz5+PgYEBmM1m7Nu3D4cPHx6TYTvdeLUw7a1hTvSUw65XZmdnu9wDztUIk2EYNDY2YmBgwGaLMHeskRIIAKkta4vL/x7O9cicaVypIzvZfTA/Px/5+flj/t/TTz/N/T0gIAAfffSRzfeyezI9Ba8W5lRhN+8ajUa7+ysdRSgU2m1kPRFWqxVVVVUICgrC4sWLbQp+qsIk0SWB4J24I8Kca5A1TAdhk2l8fX2RnZ3NWxssZyPM4eFhnDx5EhEREUhJSZkwGp7OCJNEGwTC5WnZ2Ygra5iuzrTNFrxamI4uLA8MDODkyZOIj49HUlISrwvSzuzD1Ov1OHPmDNLS0hAVFTXpsWRK1rtxd8IPSSgaj70HR0+tvUp6YdqHTMnawWKxoLq6mpf1SltMRWgMw6C1tRUXL15EXl4e/P39eT0/mY4lEPihtrYWcrkcUql03GyUJxVeH40ra5jh4eFuGJHnQYQ5Aex6pdVqxdKlS91WnsnRLFmapnH27FkwDAO1Wj2lljskwiQQppdfPj2Mnc/4oqWlBT4+PpDL5ZDL5QgKCvLIogWAaxGmt0zJEmHawGw2o6KiAjKZbEw3cnfgyBqmyWRCeXk55s2bh7i4uCmNx1FhkuiS4CwkW3Y8dz30AyQlXc6WNZlMXLFyo9EIiUQCi8XidETnLiiKcio3w5uSfrxamLbEw+6vTE5OhkKhwMDAgFu7idiLMPv7+1FdXY3U1FSnql+QCJNAmFn8/f0RFRWFqKgo0DSN7u5u9Pb24vTp0/D19eWiz8DAwBmdqqUoyqm2fyTC9FI6Oztx4cKFMeuV7iyOzp5/IqF1dnaipaUFOTk5Tj/BTZcwSYTheZCEnJnl5h+fxqfvZY+Z5hQKhZBIJJBIJMjIyMDIyAj0ej0aGxsxMjKCsLAwyGQySKXSaY8+na0lOzQ0RCJMb2Ky/ZXubr9lK8Ic3U9TrVa7tIXFkfGT6ViCq5Bp2fHc9dAPQFE93O+3UCjkSlWyEg0ICEB0dDSio6NB0zT6+vqg1+vR3Nw8JvqcDiG5sg+TRJhegEAgGLNeaWs/o4+Pj9OFBRzhSqFZLBZUVlYiJCQEOTk5Lk/RkClZAmHm8PPzA8MwoGmaE6fJZAIwXlBCoRAymQwymQzA5anOK6NPuVyOsLAwt0Sfzq6pkilZL8FgMODUqVPceqUt3B1hjp7yHRoaQkVFBRITEzFv3jxezm9PmCS6JBDcx+rNp3B4rxo+Pj7w9fXFyMgImpqaEB8fD4qixkWfowUaGBiImJgYxMTEgKIo9Pf3Q6/X4/z58/D39x+z9skHJMK0j1cL08/Pz+7+SnevYbL1bHU6Herr65GRkYGQkBDezs93+zACYSLItOx4Lv97DAO4HFlWVFRApVJBLpeDpmkwDAOKoriHWoqiOHGOlpePj4/N6LO+vh5ms3lM9OnslhVS6cc+Xi1Mf39/ux+u6RCOyWRCc3Mz1Gq1U1lqkyEUCmGxWHg9J4FAmBpGoxEVFRVISUmBVCoFAO7ew0qKnbZlJcr+3cfHx2702dfXh+7ubjQ2NiIgIAByuRwymWxK0Sep9GMfrxamI+uDIpHIbcKkKAo1NTWgaRp5eXlu2cw82ZQsX9OxJKrwPGYqQ5ZEmeO5aePJ//zNB4f3Sic8bnRUOdXok52eBS5HfHq9HnV1dbBYLJBKpZDJZHajT2cjTFbq3oBXC9MR3BWhjYyMoLy8HFFRURgcHHRb5Q+S9EMgzCx3PfQD7gHme3kCh/eqJ3yPK9FnUFAQgoKCEBsbC4qi0NvbC51ONyb6lMvl46qXeWoFIk+CCNMOIpEIw8PDvJ6zr68PNTU1WLhwIWQyGdrb23k9/2gmEiZJ9iEQZpbR8gTsC9TZ6DM8PBzh4eFgGIaLPs+dOwer1QqpVAq5XI7Q0FBQFDVlYTIM47G1cd0BEaYd+F7D7OjoQGtrKxYvXjxmfcFdHzoSYRKmGzIt6xyuRp/sfyfKvBUIBAgODkZwcDDi4uJgtVrR19cHrVaLhoYGGI1GdHV1QS6XO9TYwRvxamE6Iii+smRpmkZdXR1MJtO4YgRspux0CVOv1/N+HQKBMDGjp2Udgc/o02q1cseMjiBFItGY6PP48eOwWq04e/YsKIrios+QkJAJI09notLZjFcLE/heVhPBhzAtFgsqKioglUqRmpo6Toys1NzxwbtyH2lrayt++qtO3s5PIgnPg5TEm3tMR/QpFAoRFxfHRZ+9vb24ePEi6urqEBwczGXejo4+vSlDFvDyBtKO4KowBwcHcfLkScTFxU3YfNqdez3ZUlxse7C+vj63XIdAGA2Rtvu4aeNJ7o89hEIhfH194e/vDz8/P/j6+nIP6BRFwWw2j4lGWUQiERQKBVJTU7FkyRLMnz8fZrMZNTU1KCsrQ1NTE86cOYP+/v5JhXno0CGkpKRApVJhx44d4143mUzYuHEjVCoVli5digsXLnCvPffcc1CpVEhJScHnn3/u+D+QG/H6CNMershMq9WisbERGRkZkEgkbrmGPYRCIaxWK06fPg2ZTIa7H9W65ToEAmFypjot6whTnboFbEefIyMjAC7PhtmKPsViMcRiMeLj42G1WtHT04MdO3bgu+++g1AoxN///nesXr16TMU0iqKwdetWHDlyBDExMVCr1dBoNEhLS+OO2bVrF6RSKRobG1FYWIht27Zh7969OHv2LAoLC1FTU4POzk6sXLkS9fX1M759hUSYdnBGZgzD4Pz582hpaUFeXt6ksgTcm5hjMpnQ29uL2NhYJCYmuuUaBALBM3Am+mQYBmfPnoVKpeLuRVarlevZaSv6VCqV+Otf/4r33nsPKSkpaG9vx/r167F8+XJcunQJAHDixAmoVCokJibCz88PmzZtQmlp6ZhzlZaW4rbbbgMArF+/HkePHgXDMCgtLcWmTZvg7++PhIQEqFQqnDhxgqd/Jefx+giT7zVMiqJQVVUFf39/5ObmOrQu6a5qQt3d3aitrUVwcDAiIiIm/T4JBL4h2bIziyPRJ7sffHQFotFbVdjlHPb+xEZ47H1tZGQEUVFReOyxx/DYY4+ht7cXoaGhAC7vCIiNjeWuFRMTg+PHj4+5/uhjRCIRQkNDodfr0dHRgWXLlo15b0dHh/P/GDxBIkw7CIVCh0VjNBpx4sQJhIeHY+HChQ4n8bijwHtrayvOnz+PxYsXc9/Dyh8dt/9GwqyGrB16NjP5AHFl5Gk0GlFeXo7U1FROlsD3Wbej1z5FIhH3YE9RFBd9XtkLUyqVzumsWa+PMPmit7cXZ8+eRVpa2pgPnyPwGWHSNI3a2lpYrVau3J6taRU+INEDgTA7YWvbLly4kIsIJ8LWthU28jx27Biamppsvi86OhptbW3c1+3t7YiOjrZ5TExMDKxWK/r7+yGXyx1670wwdx8FppG2tjbU1dUhNzd3yrIE+IswLRYLTp06hYCAAGRkZHDTJxaLxaE1DQKBb0jE61kc3qvG8PCww7K8EqFQyLUq+/zzz/HPf/4TO3futHmsWq1GQ0MDmpubYTabUVhYCI1GM+YYjUaD3bt3AwA+/vhj3HDDDRAIBNBoNCgsLOQaUzQ0NGDJkiXOfdM84vURpivFAthozmKxQK1WO53BxUeEOTg4iMrKSqhUKiiVSgCXk4+sVisyMjIA1Lt0fgKBwA/uyJZ1hAPvZnKyTE9Pd6mN4P79+/HnP/8ZBw8e5FqOXYlIJMLOnTuxatUqUBSFO+64A+np6XjyySeRl5cHjUaDO++8E7fccgtUKhVkMhkKCwsBAOnp6diwYQPS0tIgEonwyiuvzHiGLECE6TBXVuIxm82oqKjg1itdEa+rEWZ3dzfq6uqQmZkJiUQyZspEIBBM+SmSQCDMLf76BwUqKipgMBgQGRnJTa06c98qLS3Fyy+/jE8//XRCWbLk5+cjPz9/zP97+umnub8HBATgo48+svnexx9/HI8//viUx+dOiDAdgE2aYT9cBoMBVVVVWLBgwZh9R66c39kIs6WlBRcvXkReXh78/f3HyVIgEJDpWMKMQrJlZw42M3ZwcBBdXV3Izc2FyWRCR0cHzp07B7FYjPDwcMjlcod68e7btw+vvPIKPvnkE7uynIt4vTCnUk9WKBTi4sWLaGpqQmZmJsRiMS9j8PHxmXILMZqmce7cOVAUBbVazUn9Slm6C3ID9DzIeuHsYTqmZUfLsqqqChkZGdw9i91mZjAY0N3djYqKCgCAXC5HeHg4JBLJuPtHSUkJXn31VXzyySdO5WrMBbxemI7AVsu5cOECBgYGoFar4evry+v5pxJhjp4Onj9/PreXlO2RN1qWJLokELwPVpYGgwHV1dXIzMxEcHDwmGMEAgFCQkIQEhKCxMREmM1m6PV6tLS0YHBwECEhITh37hxWrFiBY8eO4fXXX/dqWQIkS9YhBAIBqqurQdM0Fi9ezKssgamtYQ4ODqKsrAzx8fFISEjgZGm1WsEwDFfWikDwJEj0O32wshwYGJhQlrbw8/NDZGQkMjIysGzZMkRGRuLf//43Vq5ciQcffBA33XQTOjs7vboAChGmHYaHh9Hb2wuZTIaUlBS3teByJMJkp04yMjKgVCq5qNJqtY7ruk4gEDwbdyxrjJbl2bNnkZWV5ZAsr0QgEEAmk2H58uWIiIjAN998g9jYWGzfvh05OTkYGhrie+izAq+fkp1MgHq9HrW1tZDJZG5d4LYXYTIMg5aWFmi1WqjVavj5+Tm0XkmmYwkE74GVZX9/P86dO4esrKwxTeqnykcffYS33noLBw4cQFhYGFJTU/HTn/4UVqt1TD9fb4KEJDZgBdXY2Ijc3FwEBQW5rZsIMHmESdM0ampqYDAYkJeX57AsCd7FbJjynA1jnK2wsuzr6+NFlh9++CF27dqFTz75BGFhYWNe81ZZAkSY42AFxSb3BAQEuLX9FjBxhGk2m3Hq1CkEBwdj0aJFXCYsW+puMlm6M7okGbIEAj/w8bvEyrK3txe1tbXIzs52SZZ79+7FO++8gwMHDpA93FfgvY8K/2G0cEwmE8rLyxEREYH4+HjuNXcL01Z7r8kq97DvIRAI3g0ry56eHtTX1yM7OxsBAQFOn6+wsBC7d+/GgQMHXKoENFchd93/0N/fj7KyMqhUKm6rBst0RJijz6/T6UhyD2FOQqZl+WO0LBsaGpCTk+OSLPfs2YPdu3fjk08+IbKcAHLnBdDV1YWzZ88iJycHcrl83OvTJUyGYXDhwgU0NzdDrVbbLHPnyHolSfYhEGYPzkzLsrLU6/VoaGhAdnY2/P39nR7DBx98gHfffReffPKJ3Yb33ozXT8mazWbodDqo1eoJF7N9fHwwMjLitjGwST81NTUAwLXlIsk9BEcgUZt3wcqyu7sb58+fR05OjkNl7Sbi/fffx/vvv09k6QBeH2H6+/sjMzNz0swvd0eYFosFQ0NDEIvFSE9Pn1Jyz5WQ6JLg6RDBOw8rS51Oh6amJl5k+cEHHxBZOojXC9MR3CnMwcFBnDp1Cn5+fmPK3DlTuWc6ZEkyZAkE/nHk94qVpVarRXNzs8uy/Pvf/449e/bgwIEDvNXFnut4/ZQsAE5SE+EuYep0OjQ0NCAzMxOVlZVOT8GSqJJAmNuwsrx06RJaW1uRk5PjUonOd999Fx9++CEOHDjgVCUgb4VEmA7gar/KKxmd3JOXlwexWExkSfAqyLSs41wpy+zsbKdlyTAM3n33XXz00UdElk5AhOkAfEaYbGGEwcHBMZV7GIaBVqvlhOkIRJYEwtzB1rQsK8uuri60tbW5FFmOluX+/fuJLJ2ACBP2e2LyJUyz2YyysjKbyT1paWno7e3FyZMnUVFRga6urkl7ZBJZEgASqc1l9r2TDuCyLDs7O5Gdne10WTqGYbB7924UFRXxIss77rgDSqUSixYtmvB6999/P1QqFTIzM3H69GmXrucpEGE6AB/CNBgMOHnyJBISEmwm94SGhiIlJQXLli1DUlISjEYjzpw5g1OnTqG1tRVGo5E710zJkiT8EPiEyH5idr0wD3V1dfjnP/+JxsZGxMfHO12whJVlSUkJb5Hl7bffjkOHDk34+sGDB9HQ0ICGhgb89a9/xZYtW1y+pidAkn4cYKoNnq9Eq9WisbERWVlZXDbaRFtGBAIBxGIxxGIxEhMTMTIyAp1Oh3PnzuHXf3Df1hYCgTDz3PXQD7BhyTCAy/cCmqYRExPDJQgGBwcjPDwc4eHhDmXIMgyDd955B/v370dpaSmCgoJ4Ged1112HCxcuTPh6aWkpbr31VggEAixbtgx9fX3o6upCZGQkL9efKYgwHcDZggFsck93d7fTnUYCAgIQGxuLOx++6NQYCATC7IGVZVtbG7q7u5GdnQ0fHx9ERESAYRgMDg5yfXEBQC6XQ6FQQCwWj7uXMAyDt99+GwcOHMC+fft4k6UjdHR0IDY2lvs6JiYGHR0dRJhzAXdU0GGTewQCAXJzc12q3EPWKwlzlbde/JpM9f8HVpatra3Q6/XIysoaMw0rEAggkUggkUiQkJAAs9kMvV6P5uZmDA0NITQ0FBKJBGFhYRCLxdi1axc+/fRTlJaWutS9hPA9RJhuwGw2o7y8HEqlkut6wib3MAxDZEngBbIGOHdgZdnS0oLe3t5xsrSFn58fIiMjERkZCZqm0d/fj2PHjmH79u0IDAyE2WxGSUnJjMgyOjoabW1t3Nft7e2Ijo6e9nHwDUn64Rl7yT2eVrmHQJhpvF38rCwvXLiAvr4+ZGZmTjnBRygUQiqVQqPR4J577oFCocBdd92FBx98EGq1Gi+99JI7hj4hGo0G7777LhiGwXfffYfQ0NBZPx0LkAgTgONTsmx0OBFTSe6xh6fJkkybEQj8w8qyubkZBoMBGRkZLmXDvvnmmzh8+DAOHDiAgIAA/OpXv8LQ0BCampr4HDY2b96Mr776Ct3d3YiJicFTTz3FbYP7xS9+gfz8fHz22WdQqVQICgrCO++8w+v1ZwoiTAdht5bY2gfFR3IPi6eJkkAguAdWlk1NTRgaGsKiRYtckuUbb7yBL774AsXFxWP6YgYHByMjI4OXMbPs2bNn0tcFAgFeeeUVXq/pCZApWQeZaC8mTdOorq7G8PAwcnNziSwJBCfwtmnZDUuGwTAMzp8/j+HhYV5kefTo0XGyJPALiTAdxJYwTSYTysvLMW/ePMTFxZHkHsK04W2CmUuMlqXJZEJ6erpLW9def/11fPnllygqKiKydDNEmHBsDfNKYRoMBlRWViIlJQXh4eEAwCX3AJjS0yKRJYHgHbCybGxshMViQVpamkuyfO211/D111/j448/JrKcBsiUrIOM7lii1WpRVVWFrKwsTpYURcFqtUIgEMw5WZKEH8J0MNejZlaW9fX1sFqtWLhwoUuyfPXVV3Hs2DEiy2mERJgO4uPjA4vFgqamJuj1epLcQyAQHIaVZV1dHQAgNTXVJVnu3LkT//rXv/DRRx/B39+fz6ESJoEIE45NyQqFQjQ1NUEsFpPKPQQCwWFYWdbW1kIoFCI5OdklWb788sv49ttviSxnADIl6wAmkwldXV0ICgpCWlramLZcRJYEAn/MtWlZVpbnzp2Dj4+Py7L8y1/+gu+++47IcoYgEaYdBgYGUFVVhfDwcISFhQEgyT2EmWWuSWWuwsry7Nmz8PPzg0qlckmWf/7zn3Hy5El8+OGHDnUqIfAPEeYkXLp0CefPn0d2djZ6e3tBUdScqdxDIBDcByvLmpoaBAQEICkpySVZ/ulPf8KpU6ewd+9eIssZhEzJYvwaJrtHqrW1FWq1GsHBwRAKhbBYLE5Nwc5mWZIMWcJ0M9sj6A1LhrmCJoGBgS5Hln/84x9x+vRpFBYWElnOMESYV0BRFKqqqmAymZCbmwtfX18wDIOAgAB0dHSgsbER/f39YBhmpodKIBA8jNGyFIvFSEpKcvpcDMPgxRdfRHl5Ofbs2eOyLA8dOoSUlBSoVCrs2LFj3Outra24/vrrkZOTg8zMTHz22WcuXW8uQoT5HwQCAUwmE8rKyhAWFjYuuUcikWDp0qWQSqXo6OjAd999h7Nnz6K7u5vbn2mLw3vV0/hdEAiEmWKxrBJ9fX2orKxESEgIEhISnD4XwzB44YUXUFlZyYssKYrC1q1bcfDgQZw9exZ79uzB2bNnxxzzu9/9Dhs2bMCZM2dQWFiIe+65x6VrzkXIGuZ/GBgYQGVlJVJTUyGXywGMT+4RCoVQKBRQKBSgaRp9fX3Q6XRoaGiAWCyGUqmEXC4fV6D90/eyUVVVBalUip//unvavzfC3GG2T1c6ymxrLP0/iweg0wWjqqoKNE3D19cXWq0WcrkcPj4+UzoXwzB4/vnnUVNTgz179sDX19fl8Z04cQIqlQqJiYkAgE2bNqG0tBRpaWncMQKBAAMDAwCA/v5+REVFuXzduQYRJi5/QJuampCdnY3g4GAA9ttyCYVCyGQyyGQyMAwDg8EArVaL5uZm+Pv7Q6lUQqFQwGq1oqqqCvPnz0dERAQO7x371Dmb1zcJBAI7DSvEpUuXEB8fj9jYWPT390On06G5uRm+vr7cg7a9ijwMw+D//u//cO7cOXzwwQe8yBIAOjo6EBsby30dExOD48ePjzlm+/btuOmmm/Dyyy9jaGgIX3zxBS/XnksQYeLyk1VOTg5omnaqGIFAIEBISAhCQkKgUqkwNDQErVaLU6dOwWg0IiYmBqGhoTbfO3rKlsiTQJhdbFgyDIqiUFlZifDwcE5KYWFh3Da04eFhdHd3o6amBlarFXK5HAqFAiEhIWPuLwzD4A9/+APq6+vx/vvv8yZLR9mzZw9uv/12PPTQQ/j2229xyy23oLq62ukuKnMRIsxRsLKkKApCodDpzLbg4GD4+/vDx8cHubm5GBgYQE1NDSiK4p402QbTo/E0ec6mKTHC3MPTp2VZWVZUVECpVCImJsbmcUFBQYiLi0NcXBysViv0ej3a2tpgMBjg6+uL2tpaFBQU4NVXX0VjYyP+/ve/8y7L6OhotLW1cV+3t7cjOjp6zDG7du3CoUOHAABXXXUVRkZG0N3dDaVSyetYZjNEmP9hdFsuV2TJbkkZGhpCbm4ufHx8EBoaitjYWFgsFuh0Opw/fx5GoxFyuRxKpXLckybgefIkEAjfM1qWERER4+QzESKRCBEREYiIiOAazx8/fhw7duyAyWTCY489hosXL46ZPuUDtVqNhoYGNDc3Izo6GoWFhfjggw/GHBMXF4ejR4/i9ttvx7lz5zAyMgKFQsHrOGY7JNbG5dJ3P/nJT/DBBx+gv7/faVmyUzMMwyAzM3PcYr+vry+ioqKQlZUFtVqN0NBQtLW14bvvvkNtbS30er3NjNvDe9XcH4L34i0JP54OK0u2F66jsrwSgUDA5TYsW7YMX375JWiaxp133okVK1ZMmn0/VUQiEXbu3IlVq1Zh4cKF2LBhA9LT0/Hkk09i//79AIAXX3wRb775JrKysrB582b87W9/c/peOFcR2NlP6BWbDRmGQVVVFUpKSvDpp58iNDQUa9euRUFBASIiIhz60IyMjKCyshIxMTFTzi6jaRq9vb3QarXo6+uDRCLhMm4ny7Bzd+TpydNh3oi3CtOTPocblgzDarWivLwc0dHRiIyMdPpcDMPgd7/7HVpbW7F79+4x2fUjIyOkZdfMYvOmT4R5BWzGbElJCfbv3w+aprFmzRoUFBQgISHBpjzZNcrU1FRIpVKXrz8wMACtVgu9Xo/AwEAolUqEh4fbXdfgW6CedKMiEGHONKwsz5w5g9jYWMybN8/pc7GybGtrw9/+9rdxW9EIMw4R5lRhGAYXL17Evn37sG/fPvT09GD16tXQaDRYuHAhhEIhiouLIZVKsXTpUgQFBfF+fTbjtnIAPyYAACAASURBVLu7GyKRyOH0dFfl6Sk3KcL3EGHOHBuWDMNisaC8vBxxcXGIiIhw+lw0TeN3v/sdOjo68M477xBZeiZEmK7S29uLAwcOoKSkBE1NTYiKisKlS5dQVFTk0i+QoxiNRuh0Ouh0OtA0zcmT3Ts6Ec7I0xNuUoSxeKswgZn9PLKyPHPmDObPn+9S1ihN03j66adx8eJFvP3220SWngsRJl+YzWbcfffd6OrqgkKhQEVFBa6++mpoNBpcc80107J/ymw2Q6fTQavVwmw2cxm3Eolk0jVXR+VJhOlZeLMsgZn7PG5YMgyz2Yzy8nIkJCS4lDVK0zSeeuopXLp0Ce+8886UKwARphWbN1HyeDNFKIrCmjVrsHbtWtx///0QCAQwm8348ssvUVxcjG3btiE7OxsajQY//OEPERgY6JZx+Pn5ITo6GtHR0dzerpaWFgwODkImk0GpVCIsLGzC7SoMw2DVpjK3jI1AmAu4Q5ZarZbIchZDIkwnaG5unrCwMkVR+Pe//42SkhIcPXoUSUlJKCgowOrVqyes9sMnNE2jp6cHWq0W/f39CAkJgVKphEwm435JaZpGfX09KIri1mKvjDxJhOlZeHuECUzvZ5KV5ZkzZ5CUlITw8HCnz0XTNLZv3w69Xo+33nqLyHJ2QKZkpxuaplFRUYHi4mIcOnQIMpkMGo0Ga9asgUKhcPseJ4Zh0N/fD61Wi56eHgQFBSE8PBxdXV2QSqUTZv0CwIcn+E1gIrgGEeb0CHPDkmEAl/dml5eXQ6VScc0YnIGmafzv//4vent78eabbxJZzh6IMGcShmHQ0NCA4uJi7N+/HyKRCGvXroVGo0FsbOy0yLOnpwc1NTUQCoUICgqCUqmEUqm02zqIyHPmIcJ0vzBZWY6MjKC8vBzJycmQyWROn4+maTz55JPo7+/HX//6VyLL2QURpqfAMAw6OjpQUlKCffv2YXBwkNuukpqa6hZ5Dg4Oorq6mrsJDA8Pc0lDADh52ltzJfKcGYgwL+MuaV4py5SUFJf2VNM0jSeeeAKDg4N44403iCxnH0SYnkp3dzf279+Pffv2oa2tDStXroRGo0FOTg4vnQJ6enpQX1+PRYsW2Sz6bjKZOHlaLBaEh4dDqVRCLBZPKm8iz+mByPJ73CFMVpZGoxEVFRVITU3lOo04A03TePzxxzE8PIzXX3+dyHJ2QoQ5GzAYDPjss89QUlKCmpoaXHvttdBoNFi+fLlTe7Y6OzvR0dGBzMxM+Pv72z3eYrGgu7sbOp0Ow8PDXMZtaGgokecMQYT5PXwL80pZLly40KXkPJqm8Zvf/AYjIyN47bXXiCxnL0SYsw2TyYSjR4+iuLgY3377LdRqNTQaDVasWOFQI9qmpiYMDg5i0aJFTv3iUhTFZdwODAwgNDSUy7idLPIl8uQXIsyx8CVNVpbDw8OorKzkRZaPPfYYLBYLXn31VZdnhw4dOoQHHngAFEXhrrvuwq9//etxx3z44YfYvn07BAIBsrKyxnUgITgNEeZsxmq14ptvvkFxcTG++uorpKSkQKPR4KabboJEIhlzLE3TOHv2LHx9fZGcnMzLmihN0+jr64NOp0NPTw/EYjFXIN5e5EsE6hpEmGPhQ5h54dVQKBQQCoWorKxEeno6QkJCnD4fTdP49a9/DavVyossKYpCcnIyjhw5gpiYGKjVauzZswdpaWncMQ0NDdiwYQP+8Y9/QCqVQqvVkt6V/EGEOVegaRqnT59GUVERPv/8c0RGRqKgoAD5+flgGAYPPfQQnnrqqQn3iroKwzAwGAxcjVt/f38olUooFAqScesGiDDH44o0/zunH93d3ejq6kJvby+USiWioqIglUqdEh1N09i2bRtomsYrr7zCS97Bt99+i+3bt+Pzzz8HADz33HMAgMcee4w75tFHH0VycjLuuusul69HGAep9DNXEAqFyMvLQ15eHn7/+9/j3LlzKC4uxrp163Dp0iWsW7cOfn5+YBjGLRm3AoEAISEhCAkJgUql4grEV1RUQCgUQqFQQKlU2pw2ZqfBACJPRyCy5JfLnz9fSCQSXLhwAXl5ebBardDpdKivr4dYLIZCoYBcLneoxCVN03j00UcBgDdZAkBHR8eYJtIxMTE4fvz4mGPq6+sBAFdffTUoisL27duxevVqXq5PsA0R5ixHIBAgLS0Nw8PD+Pjjj/H666+jqakJd999N0wmE/Lz81FQUIAFCxa4ba9ncHAwEhISkJCQgJGREeh0OtTU1ICiKK5AvK3sXCJPwnTCft4GBwdRVVWFjIwM7nMpl8vBMAwGBweh1WrR2toKHx8f7vNra7sVTdN45JFHIBQK8fLLL/MmS0exWq1oaGjAV199hfb2dlx33XWoqqpyKcOXMDlEmHOA7u5u3Hfffdi3bx/mz58PAHjggQeg0+lQWlqK3/zmN7h06RJuvPFGrFu3DhkZGW775Q4ICEBsbCxiY2NhsVig0+lw/vx5GI1GrkB8SEjIOHlvWDIMhmFw4cIFnNSlu2VshLnDWy9+PaVpWVaWBoMB1dXVyMzMHNflRyAQQCKRQCKRICkpiXv4O3fuHCwWC+RyOfz9/REVFQWBQICHH34YIpEIf/nLX3j/fYqOjkZbWxv3dXt7O6Kjo8ccExMTg6VLl8LX1xcJCQlITk5GQ0MD1Go1r2MhfA9Zw5wjWK3WSZNv+vv78emnn6KkpAR1dXVYsWIFNBoNli5dOi2p7xRFQa/XQ6vVwmAwQCqVQqFQcOtGDMOgrq4ONE0jNTWVuwF5e+RJpmQnxlFhOiJLe7ANDt58800UFhYiPDwcERERKCwsnPK5HL1ecnIyjh49iujoaKjVanzwwQdIT//+YfLQoUPYs2cPdu/eje7ubuTk5KC8vNylUn4EDpL0Q7iM0WjEkSNHUFRUhFOnTmHZsmUoKCjAD37wA7tJO3xA0zR6e3uh1WrR19cHsVgMo9EIqVQKlUpF9nuOgghzYhwRJivLgYEB1NTUICsry6VG7zRN48EHH8TQ0BCio6Pxj3/8A/Hx8fj5z3/O+/rhZ599hgcffBAUReGOO+7A448/jieffBJ5eXnQaDRcgt+hQ4fg4+ODxx9/HJs2beJ1DF4MEaY97O17MplMuPXWW3Hq1CnI5XLs3buXmwKdrVgsFhw7dgzFxcU4duwYFi1ahIKCAtx4441ueXK+ErPZjNOnT8PPzw9msxmBgYFQKpUIDw+3m3Qx1+VJZGmfyaTJyrK/vx/nzp1DVlaWS+32KIrCL3/5S4jFYrz00kvcLEhdXR0GBgbIVOjcgghzMhzZ9/Tqq6+isrISr7/+OgoLC1FSUoK9e/fO4Kj5haZpnDhxAsXFxTh8+DDi4+Oxdu1a5Ofnu1RXcyJGRkZQUVGBxMREKBQKMAzDZdx2d3dDJBJxSRf2CjXMRXkSYdpnImGysuzr60NtbS1vspRIJHjxxRenPcGHMO0QYU6GI/ueVq1ahe3bt+Oqq66C1WrFvHnzoNPp3N5pZCZgGAY1NTUoKirCZ599hpCQEKxduxYFBQWIiIhw+Xtmi8FPVrfTaDRCp9NBp9OBpmlOnvYi37kiTyJM+9gSJivL3t5e1NXVITs72+4D12RQFIUHH3wQoaGheOGFF4gsvQOyD3MyHNn3NPoYkUiE0NBQ6PV6l5rLeioCgQCLFi3CokWL8OSTT6KpqQklJSW47bbbwDAM1qxZg4KCgkl7ak4EeyPLyMiYVH6BgYGIi4tDXFwczGYzt1fObDZzGbcSicRmxi3LXJEnwTZXZsu6Q5YPPPAApFIpnn/+eSJLL4cIk2AXgUCApKQkPPzww3jooYdw8eJF7Nu3Dw899BB6enqwevVqFBQUIC0tze4N5dKlS2hpaUFOTo5DxeBZ/Pz8EB0djejoaC5jsaWlBYODg1yB+LCwMCJPL4b9WbPdeab6GbsSiqJw//33Izw8HH/4wx+ILAlEmCyO7Htij4mJiYHVakV/f7/XpXALBAJERkZiy5Yt2LJlC3p7e3HgwAHs2LEDTU1NuP7666HRaJCXlzduu0prayu6u7uxePFipzqvsIhEIkRERCAiIgI0TaOnpwddXV2ora1FSEgIVyD+yuuPlidABDqXYH+2er0ejY2NvMjyvvvug1KpxI4dO4gsCQDIGiaHI/ueXnnlFVRVVXFJP8XFxfjwww9ncNSexdDQED7//HOUlJTgzJkzWL58OdatW4errroKjzzyCJYtW4bNmze77ebDMAz6+/uh1WrR09ODoKAgLuN2thWIJ+uXjnN47+Xs1O7ubpw/fx45OTkubY+iKAr33nsv5s2bh+eee47I0jshST/2sLfvaWRkBLfccgvOnDkDmUyGwsJCJCYmzvSwPRKz2Ywvv/wSH3/8MQ4cOIDk5GTcc889uPHGG13KVnSU0WXOuru74evrC6VSCaVSOSsKxBNhOgYrS51Oh+bmZmRnZ7ssy61btyIqKgq///3viSy9FyJMwvRiMBjwox/9CDfeeCOWLl2K4uJifPHFF1CpVCgoKMDq1atd6j84FYaHh6HT6aDVagGAk6c9ec+UPIkw7cPKUqvV4sKFC8jJyXGoYPpEUBSFe+65BzExMXj22WeJLL0bIkzC9LJ161Zce+21Y6qP0DSNyspKFBUV4eDBg5DL5SgoKMDatWuhUCimZYuOyWTi5GmxWBAeHg6lUgmxWGzz+jRNo6qqCnWmq9w+NhYizMkZLcuWlhZkZ2e7JEur1Yp77rkH8fHxeOaZZ4gsCUSYhOmFpulJbzwMw6ChoQHFxcU4cOAAfHx8sGbNGqxbtw6xsbHTIk+LxYLu7m7odDoMDw9zGbehoaEQCASwWq2orKyEQqEYs+3I3ZEnEebEsLK8dOkSWltbeZHlli1bkJCQgGeeeWZO7qsmTBkiTILnwjAMOjo6UFJSgn379mFwcBCrV6+GRqNBamrqtNzEKIpCT08PtFotBgYGIJFIMDAwgPj4+HEZ01fCp0CJLCeGlWVXVxc6OjqQnZ3tUsa11WrFz3/+cyQlJRFZEkZDhDkbsVff9qWXXsJbb73FlZF7++23ER8fP0Oj5Q+9Xo/9+/ejpKQEbW1tWLlyJTQaDXJycqZlusxoNOLUqVNcYXixWAylUgm5XO72jFsiTNuMlmVnZyeysrJ4kaVKpcLTTz9NZEkYDRHmbMOR+rZffvklli5diqCgILz22mv46quv5lR9W+By8tDBgwdRXFyMmpoaXHvttdBoNFi+fLlLN8yJGBoaQlVVFVe2j2EYGAwGLuPW398fSqUSCoXCLRm3RJjjYWXZ2dmJrq4uZGdnu9SWzmq14mc/+xmSk5Px1FNP8SJLew+3LEVFRVi/fj1OnjyJvLw8l69LcAtEmLMNR+rbjubMmTO499578a9//WvaxjjdmEwmHD16FMXFxfjuu++Ql5eHgoICXH/99S6VQGNh20BlZGRALBbbPGZ0gXihUAiFQgGlUslbgXgizLE89+jlf2OKotDb24usrCyXZXn33XcjNTUV27dv50WWjjzcApcf/tasWQOz2YydO3cSYXouNj8UJBXMg7FV37ajo2PC43ft2oWbb755OoY2Y/j7+yM/Px9vvfUWysvLcfvtt+Prr7/GihUrcNttt6GoqAgGg8Gpc/f09ODs2bPIzs6eUJYAEBwcjISEBKjVaqSnp0MgEKCmpgYnTpxAc3MzBgcHbb5vw5Jh7g/BMQ7vvfxv3NfXh5aWFlgsFrS0tMBgMMDOw75NLBYL7rrrLqSlpfEmSwA4ceIEVCoVEhMT4efnh02bNqG0tHTccb/97W+xbds2Xh7uCNMPKY03R3jvvfdQVlaGr7/2nuhEJBJhxYoVWLFiBWiaxunTp1FcXIw//vGPmDdvHjQaDfLz8x0qjj96L99USqoFBAQgNjYWsbGxsFgs0Ol0OH/+PIxGI1cgPiQkxOEatyS6/J7RRQkoisK1114LmqbR3d2N5uZmDA0NQSaTQaFQICwszO7aNivLjIwM/Pa3v+V1zdKR5g2nT59GW1sb1qxZg+eff563axOmDyJMD8aR+rYA8MUXX+DZZ5/F119/7VL9zNmMUChEXl4e8vLy8Oyzz6K2thbFxcXYuHEjAgICsHbtWmg0GkRFRY27UXZ0dKCrq8vlje++vr6IiopCVFQUKIqCXq9HW1sbDAYDpFIpFAoFpFLpuBv7aHlu2KvGTRtPOj2GuQIry5aWFm4aVigUwsfHB5GRkYiMjOTqCF+8eBF1dXWQSCRcYtaVU7YWiwV33nknsrKy8MQTT0x7gg9N0/jVr36Fv/3tb9N6XQK/kDVMD8aR+rZnzpzB+vXrcejQISxYsGAGR+uZMAyD1tZWbruKyWRCfn4+CgoKsGDBAjz11FNYuHAh1q9f79K62GTQNI3e3l5otVr09fVNemO/Em+UJyvLCxcuoL+/HxkZGXajR7aOsE6ng16vR0BAAEZGRhAXFweFQoE77rgDOTk5ePzxx90iS3v5Bv39/UhKSuKm+i9evAiZTIb9+/eTdUzPhCT9zEbs1bdduXIlqqqqEBkZCQCIi4vD/v37Z3jUngnDMNDpdCgtLUVxcTHq6+u5bhS5ubnTsl2FYRgMDAxAq9VCr9cjMDCQKxBvL7r1BnmysmxubobBYMCiRYuc+rkMDQ3ho48+whtvvAGDwYDU1FT85S9/gUql4nvIABx7uB3NihUr8MILLxBZei5EmAQCcPnm9otf/AJ+fn5Yvnw5SktLUVdXhxUrVqCgoADLli1zW7Q5GoZhxmTcsntpFQqF3aSQuShPVpZNTU0YHBx0WpYsZrMZd9xxB1JTU5GQkIDS0lJcvHgRd999N37+85/zNWwOew+3oyHC9HiIMAkEAPjlL38JhUKBxx57jJueMxqNOHLkCIqLi1FWVoZly5ahoKAA11133bStCxuNRuh0Ouh0OtA0zckzODh4wvcYDAb8f3fVTsv43MnhvWowDIOmpiYYjUaHmpFPhtlsxk9/+lMsXboU27Zt437Og4ODaG9vR2pqKl9DJ8xNiDAJBODyTXOybSMWiwX//Oc/UVRUhGPHjiE9PR0ajQY33njjpPLiE7PZzBWIN5vNXMatRCLhbv59fX2ora1FZmYmgoIuZ9rOxsiTleX58+dhMpmQlpbm0jqj2WzG7bffjquuugqPPvooqeBDcAYiTAJhqtA0jZMnT6KoqAhHjhxBbGwsCgoKkJ+fD6lUOi1jsFqt0Ov10Gq1GBwchEwmQ0BAADo7O5GTkzPp9K2nC5SVZWNjIywWCxYuXOiyLG+77TZcffXVeOSRR4gsCc5ChEngB28tAcYwDGpqalBUVITPPvsMEokEBQUFKCgoQERExLTcnGmaRnNzM9ra2uDn54fQ0FAolUrIZLJZl3HLyrKhoQEURblcZN9kMuG2227Dtddei4cffpjIkuAKRJgE1yElwC7DMAyam5tRXFyM/fv3g6IorFmzBgUFBUhMTHTbzXp0SyuRSIT+/n5otVr09PQgKCiIy7i1V2N3puXJyrK+vh4MwyAlJYUXWV533XV46KGHiCwJrkKESXAdR+vbPvjgg7jxxhvx/PPPz/lsQIZhcPHiRezbtw/79u1DT08PVq1aBY1G43Lyymg6OzvR2dlps6UVwzAYHBzkMm59fX2hVCqhVCrtFoifbnmysqytrYVQKERycrLLsrz11luxYsUK/OpXvyKyJPCBzQ8RqfRDmBKkBNh4BAIBIiMjsWXLFmzZsgW9vb04cOAAduzYgaamJlx//fXQaDTIy8tzertKW1sbdDodcnJybJ5DIBBAIpFAIpEgKSkJw8PD0Ol0qKioAABOnoGBgePey27nANwvT1aW586dg0gkwoIFC1yW5S233IIbbrgBv/zlL4ksCW6FCJPAK6QEGCCVSnHrrbfi1ltvxdDQEA4fPoy3334b9913H5YvX45169bhmmuucbgMH1vxJjs72+FoNSgoCPHx8YiPj4fJZIJOp8O5c+dgsVgQHh4OpVIJsVg8TjDulCcry7Nnz8LPzw8qlcolwY2MjOCWW27BypUr8eCDDxJZEtwOmZIlTAlSAsx5zGYzvvrqKxQXF+Obb75BdnY2CgoK8MMf/pDbFjIadquF0WhEeno6L1O7FosF3d3d0Ol0GB4ehkwmg1KpRGhoqF3huCJQVpY1NTUICAhAUlKSy7L8f//v/+Gmm27CAw88QGRJ4BuyhklwHVICjB8oisK3336L4uJiHD16FElJSVi7di1uvvlmhIaGgqZpvPPOO7jqqqtc3mox2Rh6enqg1WoxMDAwJuPWnpynIs/De9WgaRo1NTUICgpCUlKSS+MeGRnBT37yE6xevRr3338/kSXBHZA1TILriEQi7Ny5E6tWreJKgKWnp09YAoxgGx8fH1xzzTW45pprQNM0KisrUVRUBI1GA6lUCqvVCqVSidtvv91tQvDx8eGqCdE0jb6+Puh0OjQ0NEAsFnMF4m1l3Do6dcvKsrq6GmKxGImJiS6NmZXlzTffjPvuu4/IkjCtkAiTQPAgLBYL1q9fj+HhYQwNDUEkEmHNmjVYt24dYmNjp0UQDMPAYDBwGbf+/v5QKpVQKBRTyrhlZVlVVYWQkBAkJCS4NC6j0Yif/OQnWLt2LbZu3UpkSXAnZEqWQPBkTCYTNm/ejKuuugqPPPIIGIZBZ2cn15psYGAAN998MwoKCtw2TWuL0QXihUIhFAoFlEql3QLxrCzDwsIQHx/v0hiMRiN+/OMfQ6PR4J577iGyJLgbIkwCwZOxWCw4cuQI8vPzbb6u1+uxf/9+lJSUoK2tDStXrkRBQQEWL148La3JgMtTomyNW4qiuCndK2vzstPMMpkMcXFxLl3TaDRi8+bN+K//+i9s2bKFF1naq1b10ksv4a233uI6yLz99tsuS58wqyDCJBDmCgaDAQcPHkRxcTFqampw7bXXQqPRYPny5Xar/PCFxWLhuqsYjUauQHxwcDCqqqoQHh4+Zs+uM7hDlo5Uq/ryyy+xdOlSBAUF4bXXXsNXX32FvXv3unxtwqzB5gdteh5LCYRp4NChQ0hJSYFKpcKOHTtsHvPhhx8iLS0N6enp+PGPfzzNI+QPiUSCDRs2oLCwEGVlZVi7di2KioqwfPlybNmyBQcPHsTIyIhbx+Dr64uoqChkZWVBrVYjNDQUra2tOHbsGCiKQlBQEGiadvr8w8PD2LRpE/7nf/6HN1kCwIkTJ6BSqZCYmAg/Pz9s2rQJpaWlY465/vrrua0+y5YtQ3t7Oy/XJsxuSJYsYU5AURS2bt06JmpgS9OxNDQ04LnnnsO//vUvSKVSaLXaGRwxf/j7+yM/Px/5+fmwWq345ptvUFJSgqeeegrJyclYt24dbrrpJkgkEreNwcfHB3K5HO3t7UhOTkZgYCC0Wi3q6+shkUi4jFtHKx2xsvzRj36En/3sZ7yuWTpSrWo0u3btws0338zb9QmzFyJMwpxgdNQAgIsaRgvzzTffxNatW7m2XEqlckbG6k5EIhFWrFiBFStWgKZpnDlzBkVFRfjTn/6EiIgIaDQa5OfnIzw8nNfrUhSF8vJyREZGIioqCgAgl8vBMAwGBgag1WrR1NSEwMBArkD8RJWOhoaGsGnTJmzcuBF33333jCb4vPfeeygrK8PXX389Y2MgeA5EmIQ5gSNRQ319PQDg6quvBkVR2L59O1avXj2t45xOhEIhcnNzkZubi2effRa1tbUoLi7Gxo0bERAQgLVr10Kj0SAqKsolKVmtVpSXlyM6OhqRkZFjXhMIBAgNDUVoaChUKhWXcXvmzBkuoUahUHAZt0NDQ9i4cSM2b96Mu+66yy2yjI6ORltbG/d1e3s7oqOjxx33xRdf4Nlnn8XXX38Nf39/3sdBmH0QYRK8BqvVioaGBnz11Vdob2/Hddddx217mOsIBAIs/P/bu7uQpv4wDuDfQ3VRGdKkXDSVYkRCWLsI9SKHvcJq566jhs61hRlBGhG9DGIF0VXdeBOElCHNi+aaF6eZFb1clFDNIhZhYEgR5Yxsi5ZznP/Fn0aSttN0L7nv59JzmI9XX3/6PM+vtBQOhwMnT57E8PAwPB4PmpqaEIlEYDKZYDab/3oZ+s+w1Ol00Gq1CWvIy8uLLzD4/v07RkZG0N/fj+PHj6O6uhpPnz5FY2Mj7HZ7yk6WGzduxODgIIaGhrBy5Up0dXXh2rVrk97x+/3Yv38/fD7fnPxLBCWHTT80J6g5Neh0OoiiiAULFmDVqlVYs2YNBgcH011qxgmCgJKSErS2tuLu3bvwer3QarVwOByoqqrCmTNnMDAwkLBhJxqNwu/3o6ioKGFYTmXhwoUoLi6G0WhEZ2cnXrx4gVAohIsXL+LUqVPw+/1I0MWflF+3VZWWlkKSpPi2qp6eHgDA0aNHEQ6HsXv3bmzYsIEbrAgAx0pojlCz49bn88HlcqGjowPBYBAGgwEDAwMoKCjIYOXZZWxsDLIso7u7G69fv4bRaIQoiqioqJjUsDM+Po7nz5+jpKRkxiewcDiM2tpa1NfXw2azIRQKwefzwePx4OzZszPeEESUBM5h0twmyzJaW1vjO24dDsekHbeKouDIkSPw+XyYN28eHA4HamtrM1121opEIujr64Pb7caTJ09QXl4OURSxdu1a1NXV4cqVK9Dr9TP6HuFwGDU1NbBYLNi7d+8sVU40YwxMIkpONBrFw4cP0dnZCa/Xi02bNkGSJGzbtg2LFy9O6jPD4TAkSYLVaoXVap3dgolmhoFJRMn79OkTzGYznE4nNBoN3G43+vr6UFRUhF27dsFkMkGj0aj6rFAoBEmSYLPZ0NjYmOLKif4aA5OIkrdnzx7YbDZs3bo1/rWfl0K73W7IsowlS5bEx1UKCwun7HT9GZZ2ux0WiyWdPwKRWgxMIkpeLBb746YeRVEwNDSE7u5u9PT0IBaLYefOnTCbzVi9ejUEQYiH5b59iW17yQAAA1pJREFU+9DQ0JDG6on+CgOTiNJDURR8/PgRN27cgMfjwefPn2E0GnH//n20tLSgvr4+0yUS/QmXrxNlSqLF8MPDw6iurobBYEBZWRlkWc5AlbNHEARotVo0Nzejt7cXt27dQmFhIbZs2cKwpH8WT5hEKabmOqmmpiYYDAYcOHAAgUAAJpMJb9++zVzRRLmNJ0yiTFBznZQgCPj69SuA/5cH/FxgTkTZg7tkiVJMzWJ4p9OJ7du3o62tDd++fcPt27fTXSYRJcATJlEWcLlcsFqtePfuHWRZRkNDw4wuXyai2cfAJEoxNYvh29vbIUkSAKCyshKRSATBYDCtdRLRnzEwiVLs1+ukxsfH0dXV9dvtF8XFxbhz5w4A4NWrV4hEIli2bFkmys06iTqMf/z4gZqaGuj1epSXl7NZilKGgUmUYmqukzp//jwuXbqE9evXxxebp+o+yH9JLBbDwYMHcfPmTQQCAbhcLgQCgUnvtLe3Y+nSpXjz5g0OHz6MY8eOZahamus4VkJEWevRo0dwOp3o7e0FAJw7dw4AcOLEifg7O3bsgNPpRGVlJSYmJqDVajEyMsJfOGgmOFZCRP+WqTqM379/P+078+fPR35+PkZHR9NaJ+UGBiYREZEKDEwiylpqOox/fWdiYgJjY2MoKChIa52UGxiYRJS11HQYi6KIjo4OAMD169exefNm/v+SUoKBSZTDbDYbli9fjnXr1k35XFEUHDp0CHq9HmVlZXj27Fla61PTYWy32zE6Ogq9Xo8LFy5MOXpCNBvYJUuUwx48eIC8vDxYLBa8fPnyt+eyLKOtrQ2yLKO/vx8tLS2/rfUjmoPYJUtEk1VVVUGj0Uz73Ov1wmKxQBAEVFRU4MuXL/jw4UMaKyTKHgxMIpqWmrEOolzBwCQiIlKBgUlE01Iz1kGUKxiYRDQtURRx9epVKIqCx48fIz8/HytWrMh0WUQZwQukiXJYXV0d7t27h2AwCJ1Oh9OnTyMajQIAmpubYTKZIMsy9Ho9Fi1ahMuXL2e4YqLM4VgJERHRZBwrISIiShYDk4iISAUGJhERkQoMTCIiIhUYmERERCokGivhHTlERETgCZOIiEgVBiYREZEKDEwiIiIVGJhEREQqMDCJiIhUYGASERGp8B+sXSIMpWHlFAAAAABJRU5ErkJggg==\n"
},
"metadata": {
"needs_background": "light"
}
}
]
},
{
"cell_type": "code",
"source": [
"def make_matrix_sparse(N):\n",
" rows = []\n",
" cols = []\n",
" data = []\n",
" b = np.zeros((N+1)**2)\n",
"\n",
" h = 1/N\n",
"\n",
" for i in range(N+1):\n",
" rows.append(i)\n",
" cols.append(i)\n",
" data.append(1.0)\n",
" b[i] = 0\n",
" for i in range(N**2+N, (N+1)**2):\n",
" rows.append(i)\n",
" cols.append(i)\n",
" data.append(1.0)\n",
" b[i] = 0\n",
" for i in range(N + 1, N**2+N, N+1):\n",
" rows.append(i)\n",
" cols.append(i)\n",
" data.append(1.0)\n",
" b[i] = 0\n",
" for i in range(2* N + 1, (N+1)**2-1, N+1):\n",
" rows.append(i)\n",
" cols.append(i)\n",
" data.append(1.0)\n",
" b[i] = 0\n",
" for i in range(1, N):\n",
" for j in range(1, N):\n",
" index = j * (N+1) + i\n",
" rows += [index, index, index, index, index]\n",
" cols += [index, j * (N+1) + i-1, (j-1) * (N+1) + i, j * (N+1) + i+1, (j+1) * (N+1) + i]\n",
" data += [4/h**2, -1/h**2, -1/h**2, -1/h**2, -1/h**2]\n",
" b[index] = 1\n",
"\n",
" rows = np.array(rows)\n",
" cols = np.array(cols)\n",
" data = np.array(data)\n",
" # Note: The error we saw in lectures was in the next line: data, rows, and cols in the wrong order\n",
" A = coo_matrix((data, (rows, cols)), ((N+1)**2, (N+1)**2))\n",
" return A, b"
],
"metadata": {
"id": "F-Otnjz78nIM"
},
"execution_count": 77,
"outputs": []
},
{
"cell_type": "code",
"source": [
"A2, b2 = make_matrix_sparse(4)"
],
"metadata": {
"id": "IFCPfu2e_gWB"
},
"execution_count": 78,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# Do a sparse solve\n",
"sol2 = linalg.spsolve(A2, b2)\n",
"\n",
"# Check that the two solutions are the same\n",
"assert np.allclose(sol, sol2)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Tf_GwfPe_ibm",
"outputId": "043e08f5-3a96-4e0f-b529-8588f8995bad"
},
"execution_count": 81,
"outputs": [
{
"output_type": "stream",
"name": "stderr",
"text": [
"/usr/local/lib/python3.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:145: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format\n",
" SparseEfficiencyWarning)\n"
]
}
]
},
{
"cell_type": "code",
"source": [],
"metadata": {
"id": "rN1I9JgKJEsp"
},
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment