{ "cells": [ { "cell_type": "markdown", "source": [ "# networkinference tutorial\n", "\n", "This tutorial explains some of the core functionality of the `networkinference` package and how it is intended to be used in practice. The package obviously heavily features the econometric work of the package author, who is in need of advertising due to his refusal to tweet. \n", "\n", "We will load data from a single network, estimate a linear-in-means model via TSLS, and construct dependence-robust confidence intervals with several different methods, each with different advantages and disadvantages. Inference with OLS or IPW estimators involves the same commands, other than the estimator class.\n", "\n", "## Working with Python\n", "\n", "For users new to Python, there are several options for executing Python commands. First, Stata17 users can call Python from [within Stata](https://www.stata.com/new-in-stata/pystata/). Second, one can install a Python IDE such as [PyCharm](https://www.jetbrains.com/pycharm/). Third, one can work in command line by typing \n", "\n", " python3 \n", "\n", "or execute Python commands saved in a `.py` file by typing \n", "\n", " python3 my_python_script.py.\n", "\n", "There are many tutorials available online for learning Python. A crash course put together for a Ph. D.-level economics class on networks can be found [here](https://github.com/mpleung/Python_tutorials).\n", "\n", "## Load and format data\n", "\n", "This part of the tutorial partly serves to illustrate some of the utilities in the `networkinference` package but mostly serves to explain how to get data into Python and into the right format for analysis.\n", "\n", "We will load data from two files. The first is a csv containing outcome and covariate data for 1000 nodes in a single network. The second is an adjlist file that stores the network.\n", "\n", "We use the `pandas` package to load csv data into an `R`-style dataframe object." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "import pandas as pd\n", "node_data = pd.read_csv('node_data.csv') # file assumed to be in this \n", " # notebook's working directory\n", "\n", "print(node_data.head(5)) # show first 5 rows of dataframe\n", "print(node_data.describe()) # print summary statistics" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " Y X\n", "0 8.304543 1.726436\n", "1 -0.029856 -0.445788\n", "2 1.312531 0.815567\n", "3 2.271223 -0.544314\n", "4 -1.367210 -1.678043\n", " Y X\n", "count 1000.000000 1000.000000\n", "mean 1.894746 0.002677\n", "std 3.269896 0.983470\n", "min -9.054192 -3.158637\n", "25% -0.315210 -0.630493\n", "50% 1.672337 -0.008020\n", "75% 3.815592 0.619925\n", "max 18.022742 3.138272\n" ] } ], "execution_count": 1, "metadata": { "collapsed": false, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:48:59.232Z", "iopub.execute_input": "2022-04-10T00:48:59.246Z", "iopub.status.idle": "2022-04-10T00:48:59.647Z", "shell.execute_reply": "2022-04-10T00:48:59.702Z" } } }, { "cell_type": "markdown", "source": [ "The first column above is the outcome, the second a scalar covariate. \n", "\n", "Next, we load the network data using the `networkx` package." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "import networkx as nx\n", "A = nx.read_adjlist('network.adjlist', nodetype=int) # load data\n", "\n", "nx.draw(A, pos=nx.spring_layout(A, k=1/20, seed=0), node_size=40, \\\n", " edgecolors='black', node_color='white') # plot network" ], "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACr10lEQVR4nOy9e3xT9f0/fsqlSU5OTi5N0yRt0jQpiJRrS0vT9EqLqAiCIFcRaYtQgQilbuJlm7rJ1K+bm246b/O6ObfppoLgnKLFG2IFtdQar8W1KGKFoCiFPn9/+Hu/PyfJOSdpm9K0nOfjsYcjTc55n/c55/16v16v5+v5SgIARoECBQoUKDhNMGygB6BAgQIFChScSiiGT4ECBQoUnFZQDJ8CBQoUKDitoBg+BQoUKFBwWkExfAoUKFCg4LSCYvgUKFCgQMFpBcXwKVCgQIGC0wqK4VOgQIECBacVFMOnQIECBQpOKyiGT4ECBQoUnFZQDJ8CBQoUKDitoBg+BQoUKFBwWkExfAoUKFCg4LSCYvgUKFCgQMFpBcXwKVCgQIGC0wqK4VOgQIECBacVFMOnQIECBQpOKyiGT4ECBUMSR48eZT744APm6NGjAz0UBQkGxfApSCgoi5WCvuLEiRPMFVdcwTgcDmbmzJmMw+FgrrjiCubEiRMDPTQFCQLF8ClICCiLlYJ4YdOmTczevXuZffv2MYFAgNm3bx+zd+9eZtOmTQM9NAUJgiQAGOhBKFBwxRVXMHv37mUefPBBxmazMR0dHczy5cuZiRMnMrfccstAD0/BIMHRo0cZh8PB7Nu3j7HZbPTzjo4OJicnh2lra2M4jhvAESpIBCiGT8GAQ1msFMQLH3zwATNz5kwmEAhE/G3UqFHMli1bmNGjRw/AyBQkEpRQp4IBR3t7O2M2m0OMHsMwjM1mY1JSUpj29vYBGpmCwQa73c589dVXTEdHR8jnHR0dzKFDhxi73T5AI1OQSFAMn4IBh7JYKYgXOI5jamtrmeXLl9PniYTNa2pqlMiBAoZhFMOnIAGgLFYK4onNmzczEydOZHJycphRo0YxOTk5zMSJE5nNmzcP9NAUJAiUHJ+ChMCJEyeYTZs2Mffddx+TkpLCHDp0iKmpqWE2b97MjBgxYqCHp2AQ4ujRo0x7eztjt9uVzZOCECiGT0FCQVmsFChQ0N9QDN8ghWIgEgex3AvlfsUXynwq6AuUHN8gg1LonTiI5V4o9yu+UOZTQTygGL5BBjFVij179jCrV69WZL76CVIyanIKIeQ3GzduPOUqIuTcBw4cGHLyb4oqi4K4AAoGDYLBIAwGA9rb20M+b29vh0ajgV6vR0NDA7q6ugZohEMLXV1daGhogMFgQHZ2NgwGA51fuXvBcRz0ej3cbjfUarXod4xGI4LBYL+N1+FwQK1Ww2KxDJnnQm7O+2M+FQxdKB7fIIJcoXd6ejrzzDPPKLvfOELOu5C7FzzPM8888wzz7LPPMunp6aesMF843ra2Nubjjz9mJk6cyCxZsiTiuRiMYuDRhA4CgcCguyYFA4SBtrwKYkcsO15l9xsfRJvrjo4Oyb8bDAYEg8FT6qFEO1cgEIDRaERnZ2eEF1tTU4POzs64jaW/EIuXTa7J7/ejublZeQ8UiELx+AYRYin0VmS+4gM578JoNDIdHR3M1KlTmQULFoTciwULFjC1tbUMx3GS9+vCCy9kli1bFlc2YjRvqLu7m0lJSRHNOba2tjKZmZkJTxLhOI65+OKLI+Z88eLFjNVqZVpaWpiWlhZm8eLFzL333sucffbZCvlFgTgG2vIq6BlIHsdoNMJqtYbknQAl30EQDAbR2tra63mIlk/V6XTQ6/W44IILoNFokJ2dDaPRCJ1Oh7a2Nvr9rq4u1NXVQa1Ww2azwWAwwGKxoLm5uU/X15PxCj0+vV4v6aVWVFSgoaEhruOKNy6//HJMmDABRqORzjnP83TOGxoaMH36dHqN7e3tmD59esJfl4JTC8XwDVIEg0FUV1ejoqLitH7Jww2cHCGlJ+js7MSkSZNQWFgYMr9VVVVoaGhAe3s7SktLYTKZoNFo0NjYiGAwCL/fD5/PF3FPamtrwfM89uzZ028bE6/Xi5KSkpBzl5SUIC8vD9OnT0dNTQ2ys7NFf5udnY3GxsaE3DSRe/ziiy8iLS0NBw8epJ81NTXRa1LILwpihWL4BjGE3h/Z/Q4F9l4sCDdwer0e1dXVWLVqFSorK3u9GSDH1Wq1SE9Ph1arBc/zcLvdYFkWXq83xLvW6/Xw+/3Uy+jq6kJeXh5YloXH44HRaERdXR0qKytRV1cHn88Xcox4IRgMQq/Xo66uLuR5qKurg0qlgt/vR2dnZ1TDkJ2djdbW1riOrbcIv8cajQazZs0KmTuhsWttbZU17IlyXQoGHorhGwLoa1hvMIKEtNra2tDQ0ACe52EymWTLB5qamqLOUUNDA6qqqiK8vJqaGhouFB4jOzsbzc3NIRsQlmVhMpmQnJwMi8VCywoMBgM1gvH2ysmiHwwG0dTUFHKtwkW/oaEhwiskG4NE84xiDVuS7wUCAcXjUxATFMOnYNBBuMsni15dXR2Kiorg8XhEf2O1WuF0OmVDn7GEyoRGJHxBJRuQlStXwufzYc+ePdDpdNizZ0/IxqQ/FuLOzk7wPB/CbGxoaEBbW1vIud5++20YjUZoNBrY7faQ7yVSmLwnYUth5MNisYiGmhPluhQkBhTDF2ecjt7XqYbQuzEYDHSnL7fjJyUGZCFcu3YtXn75Zdx5551Yu3YtKioqYDKZYLVaRc9JcmDhxwlfUAcq9NbQ0ICysrIITzU7OztkjDU1NaisrERxcTFWrFhBw4iJFibvzdwFg0E0NzfD7/cPSPhfefcHDxTDFyfEi1ShIDqIcWlsbKSLIFkkxcJjhJBC0N7eDrVajfz8fNTU1OA3v/kNtm/fjtbWVlnDWVhYCI7jZBdU4VhOFdkiWn1bZ2cngsEgXn31VXAcB7PZjP3799PfJuJi3de5O5XXpbz7gw+K4YsThjqNOtEWyIaGBlRUVECv14d4esKwV1ZWFjQajegiJOU1iN1Hn88HlmWh1+uh1WpRXV0tWfAdvmCfiucimndUXV0dImM2Z86cQbEoD5Z3arCMU8H/QTF8ccBQplEn6m6WjIvjOPh8PtTV1YUsPoFAAHl5eVCpVD26L8LjpqenR+TAyHnkFjXhQkjq+FiWhcPh6JfQWzSPr6KiAoFAAK2trQgEAoNmUR4MrOWh/O4PZSiGLw4YyjTqRN/NlpeXg2VZaDQamM1mqNVqZGRkgGVZZGZm4pJLLunx+ElpAKnNIwgvBpda1MQW7P6W0JIK8ep0OtTV1YVsXMi/B8uinGjRBiGG8rs/lKEYvjhgqO76xK4rGAyGkDwGEsFgEG63G3a7HR6PB3v37sX1118Pg8GAX/3qVzhx4kSvvIZYFrNYFrWByDMJr7Ompgapqamiht9isSiLchwwVN/9oQ7F8MUJie4Z9QZCAyBWTFxTUzMgYSfhWKxWK9RqNc466yyqXPLBBx9E/CbcCMkZpVjlvxJxURNeV0dHh2Rdo1qtRkdHxwCNcmhhKL77Qx2K4YsTBkM+oqcQq5cTvtyVlZXw+/2nfFxSBJSzzjoLJ0+ejLiG3kiaSS1mseT4EgWtra2w2+2if3M4HIrHFycMxXd/qEMxfHFGIucjegMhe1LKc/D7/afsJY81tNTV1QW/3w+e5+F0Omkz1vr6+ph258LFTNjUNVHIPbFg27Ztkh4fKXNQED8MtXd/KCMJAAa6Q4SCxMWJEyeY/Px85sCBA7QVjBDZ2dlMWloaU1RUxNxyyy39Pp4PPviAmTlzJhMIBCL+NmrUKGbLli2M2+1mSkpKmAMHDjCdnZ1Mamoqc/DgQSYlJYX54osvmEAgENK+p6Ojg8nJyWHa2toiWgUdPXqUaW9vZ3ieZ44cOcLY7fa4thPqL3R2djI5OTlMZ2cnk5+fz/ztb39jbDYb09HRwVx88cXMp59+ysyZM+eU3DMFChINSj8+BbL4/vvvmU8++YQ5duxYhOHr6Ohgvv76a+bBBx9k7rvvvlPS+Zr0GhQby6FDhxi73c5s3LiReffdd5ns7GympaWFCQQCzO7duxme5xmtVtujjugcxzGjR49mrFYrM3r06IQ2eqSrekdHB7NkyRLmxIkTDMdxzN69e5mxY8cyo0aNYnJycphJkyYxL7zwwim7Z4MFYl3phZ8Nxq71CiQw0C6ngsQGIbg0NDTIdj04FdTtjz/+GB6PByzLSrZjCgaD4DiOhviEOT232w21Wo26urqQUOVgZ+CRa9Tr9VQU22q1QqVS4ZxzzkF6erpoGE6h2/8IsbxvfX096uvrYTAY4PF4oNPpoNVq4fF4BlW4W4E4FI9PgSzsdjvz1VdfMX6/n8nJyWHcbjeTnZ3N5OTkMBMnTmQ2b94c4m31F5566ilmwoQJzFdffcW89NJLTF5eHpOTk0O9GDKW9vZ2xmAwMAaDgbHZbMymTZtox/GPPvqI+fjjj5n33nuP8fv9DMP8X0f0MWPGMGq1ut/G358g17hkyRJm4sSJzMcff8x0dHQwn3zyCdPZ2ckcOnSICQaDIR7rqbhngwXCZ4R0pX/rrbeYJ554gnnttdcYi8XCTJgwgQkEAsyHH37I7Nu3j9m7dy+zadOmgR66gt5ioC2vgsSHkOHo9/vh9XoRCAQA/OgtTZs2Le4sR+KhfP3116ivr4dOp0NmZiY+/fTTiO8IvRhSfK7VamVFq0lHdI7j+q1V0KlAuFC32LXyPB/Rasnn8w0IIzfREE6WIs9UIBCARqMBx3HQaDRKnd4Qg2L4FESFkOHo8XjA8zw4jqNsR47jcODAgbgw2oRhJ4/HA41GA5PJhPz8fBw6dCimYzQ0NCA7Oxv5+flwOByi33E4HHj00UdDui0MxoVMWEwvVXTv8XhQU1MTQre3WCxobm4+xaNNPJB5Ewt3EiUgqVZXSqh48EIxfApihtDDCgaDeP/99+H3+6FWq8GybJ+1PIPBIKqrqyPyd6WlpVi/fn3Mxzl27Bjy8/OhVqujFnALFy+PxzPoFrKWlhawLBtTSyahNzMYjXx/gHh84VqvQnFynufR1tYW8jvhvCoYfFAMn4I+YePGjSgtLY1aFydX4xSuxEJq7ojx7Ik31tXVBa/XC5ZlkZGRAaPRKCnZ1djYSI9LjGF1dTXeeOMNdHd3ix5/IGq1xM753Xff4bLLLsPIkSPB8zwqKipEF++ysjJkZ2crqiIy8Pv9YFlWctNQXFwcMYdKqHhwQzF8CkQRbYEPBoNoamqSLGwnBiUWpZRYJJ9iDSv5/X7agTsYDILnedTV1YWE+erq6qDX61FeXo6GhoaQ3T3HcbQzw89+9jO0tLSgtbUVnZ2dEddRU1PTr0XgYnO3ceNG3H777dBqtVCpVPjZz36GAwcOYMaMGVCpVEhJSYFaraadIAg7UVEVkcabb76JjIwM0b+RBsQcx9HwO8uy8Hq9yhwOYiiGLw4YbIoNZLwdHR0R445mqIgiil6vh9vtjtrvTsqobdy4Efv378djjz0GjuNkjSf5/2LjFaKzsxM6nS7kWA0NDaiqqgppy+Pz+cBxHHieh9vtpgtZW1sbfD4fKisrMXXqVGg0GqhUKkqCCe9wXlxcDJ7n+82QSEmzqdVqLF++HF9//TW9VxkZGVCr1Vi8eDH2798vSvqJ5RkdbM9yXyAsA5EKiZNn0Gq1IiMjAzzPn1KlIgX9A8Xw9QGJ2qtOCmL1Xg6HI2Tcct4XCSMSj4r8XazDOTFUUnkn0kaopKQE6enpouMlu+3p06fD6/WKznN3dzdee+011NTUQKvVwmq1il6zRqOB1WqFwWDA2LFjoVKpYLFYoFKpoNPpUF9fj66uLsqCXLduHWVCEqaoVCisoqIi5PrjYTzkpNlIbime4si9fZYHs6EUzp9cnSp5npuamgbldSqIhGL4+oDBpspOxiuWC5o+fTr8fn/UhqZyuRDinZE5kGMaut1uvPvuu1i5cqXkbluj0cBgMMDr9UYsSpWVlZg4cSIMBgNUKhXUajWmTZsW4fEJj7Vs2TIEg0H4/X4UFhbSfnvh9y01NTWkgW2063juueeg0+lw7733Yvny5b3aCJ08eRLNzc247777cMstt2Dp0qWSG4K0tDRZ4k5viCtyz7KYcQs3lHq9ftB4QmJhemEDYrvdHtGAOFHfaQW9g2L4eonB1ocr1novqQU+LS0NI0eOlKR2W61WOJ3OkBxSZ2cntFqtJKsyMzMTFosFdrs9wrBVVVWhpqYmqtd40UUXYceOHTh+/DgA6fAgycnE0nJIpVKFeI5yv9FoNNDr9TCZTHC5XCguLo4gl5x11ln4z3/+g/379+Pw4cPYvXs3HnnkEVx77bVYsGABJkyYAK1WC7VajfT0dGi1WlRUVEgacaPRiMbGxrjR7KWur62tDTzPU9WbaJEB4TwnIsRUfMI3Jp2dnTCZTNDr9UpOdAhDMXy9RE86LydCOCiWei+32y0Z0tPr9cjIyJAtCB8xYgSys7OxaNEi3HXXXaiurkZ2dnZE8TTxOn0+H3JycvDEE0/gvPPOg06ng8PhgFarxRlnnIG8vDzYbLaI8CWB2AIf3iLGYDCEeCLR6t2mTp2KRYsWRVxnTU1NhFEjIV4hOUbOOA4bNgwMw0Cj0cDpdKKsrAxr167FrFmzROfI6/XKemE93XiJPYdHjx7Fr371K9H2RQ0NDRFh7bKyMqxcuRI8z4uem2XZhGU7SnWplwrTD/Q7q6D/oBi+XiKaF3Dffffhhx9+SJgcYCwen9FohN/vlw2DStU7+Xw+LF++HFOmTEFqaiqSkpKgUqnQ1tZGa/08Hk/IDposlNOmTcNFF12EtWvXYuPGjfjDH/6Ap59+Gm+88Qbee++9XnnWUpsNufsmbLEUvki2tLRAq9XSPGL4vST3XWw8mZmZuO+++/Dpp5/i0KFDeOONN/Dggw/iqquuwuzZsyWVQQihhud5uFwu6PV6rFy5EgcPHkR3d3fMoXax/N0ll1yC1atXw2g04pxzzokwZHJ5zeTkZKSlpYnOOxE4SDSDEUvONNFTFQriB8Xw9QFSC8/ixYuRn5+PjIwMTJs2LWFygHI5vpKSEmzYsEG2qSYhAJDyADlq9/vvv4+srCy0traiqakJbrdb1BDFEpaLdy41lqa6YvPg9XpRUFAAq9UqurBbrVY0NjaGfBbNQMt5oE6nE8uXL8fChQtRWlqK3NxcZGVlged5jBgxAqmpqbDZbGBZFunp6eA4DkVFRbjpppvwwAMP4Omnn8Zrr72G2tpaUcmy0aNHo66uDhdffDFGjx6NoqIimvdsbGyUbGJLSEFSmxFyrxMJcvMsFqZXMLShGL4+QM5IHDlyRDZHcyp2xHLdx4WsTq1WC51Oh+nTp9M8WFNTUwSLTXi9brdbktpNSh6Il0e0M3s7F/HucN2T4wnnUHhdcgSgnhjo3uaKf/jhB3R0dODdd9/F1q1b8fvf/x633XYbfvWrX2HDhg1YtmwZzj33XOTl5cmq1xgMBthstpCuDmq1mtYJSkU05s6dGxEGJZuqRM5xS82zwtg8vaAYvjhALKzWkxxgvNHZ2YmamhrJEGt4Hd+RI0dw0003ITk5GaNHj44amo2WsxTzqLKzsyPq4HrqtcU7V9rb4/n9flHqe319fa8MdH+yg2N5DqVyXzqdLsJTrKqqAs/z6OzspAo5JISd6GLfg42FraD/oBi+fkJ/sT5jkf7ieV6UiFFTUyN73tmzZ6OoqKjXCwPxFKUYgkQZZbCz5aJ5jD01qPH2aIWI9hzKsWZZlkVmZmbIRigrKwuFhYV03H6/nwoBJPo97c95VjC4oBi+fkS8C4z9fj90Ol0EtVx4voqKCklSAmEXSnlxvTXUwhCq0+mUZWE2NTUNGbZconig0SD3HDY1NcHpdErmLBcvXkzvKyE3hYszJwJruSc4FeMdbHNyukExfP2IeO0wu7q6UFBQQMNKZAGqqKiA3+9HMBjEzp07wXEctm/fLhvaIkoo4ca3L6FZ4cIqxwZMxNzP6QCx51DYYVxKGJywVMMXcaUdjzTEGLR+vx/Nzc3Ks59AUAzfKUBfdn/E6AmFlxsbG6kaPyEiEGICIaqItVEhhicQCIDneXR0dISMsbdlA+GGjuhjKrmUxILwOZSraSP/P5qGak/PLUaYGmoQzmtXVxfq6urAsmyENKCCgYVi+BIcfr8fGo0GbW1tITtJnueh1WphNpvh9XojShPC26hUVVVR8oVer0dqamqIwHIwGERNTY2kwRITtn799deRnp4eUdMl1MccDLmf0w3RalANBgMmTpwIl8vV5w1MV1cX6uvraf2jRqMBz/NUG3UoIXxeFTJN4kIxfAmMYDBIc3rkJSIGUK/Xw263Q61WIzc3F8eOHaO/I9R6vV4Pq9VK1TTWr19PdRWF/y0sLKTEBZ1OB61WS5l6JCQmFLYmnQBUKhVVIpEqDB7qO/zBiGg1bW63G5WVlTh8+HCfQ/UNDQ0RbN6qqipkZ2cPOQMgnNfBJml4ukExfAmM1tZWKiNGwoliu8iSkhJ4vd6Q39psNmi1Wmg0GtjtdvA8D47jInbwpaWlcLlcIQaV1HJdfPHFuOyyyySL3ouLi1FTU6PsbAcZoqmYFBUVhXS8722oPtp5hloHc+H1DmQ5k4LoUAxfAoO8SIsXL0Z6enpUuS2SsyMhq/Lycvrd7du3SxYykwR8uPEqKioCy7JRZc5Ik1aFJj54sGbNGtECdGEbnr4aJbJxE0N2djacTueAGID+ZFySTWC0d2YoGfzBCMXwJThIqEir1aKxsVFyF2mz2bBt2zbqxYULJjc1NcmWGci182lqaopp96pQuBMfnZ2dqK6uBsuyMBqNId3ahZuVeHglcgzfgfD4TkX/TCGD1mKxSG4uFAwsFMOX4BAWpXu9XsmFRK1Ww2w2g+M4VFdXRxi5YDAoKRtGuqmLwWq1Yvv27crudZCDPEdarRbp6engeR7V1dXYvXs3vF5viE5pPO5rMBjEW2+9BZfLJdq4eCByfKcyJB8MBtHc3Ay/369EQhIQiuEbJCAyZBzHRSwkJSUlmDhxIiWSBAIB0bDmqlWrIn7r8/lk2xER/Ump5rXK7nVwgJSYCHO5pPdfXV0d9Hp9XDoUhLM4VSoVzTUPJKtzIMkmSiQk8aAYvkGGgwcPYvLkyVCr1bDZbFCr1bBYLPj4448BhPZyCyeyVFZWwuv1UpFpo9GIwsJC2komXOZMqD8ZLmyt7F4HD4SLvlT9ntFojEuHAikWp9vtxrJlywaM5auQTRQIoRi+QYqOjg5s27YN+/fvFyWWHDt2TJJwEt5xgHgAZrNZ0rCJ1fEpGFjE6kmQRT8aOWrnzp19Dm8mEotTWDT/+eefx7UgX8HghmL4hgjkGq/Gsjgqhm3woKckjc7OzqjkKIfD0WevJ1FYnCTcyvN8SIhVp9OhpKRECdcrwDBGwZAAx3HM6NGjGY7jYvpc6vdWqzWm7ysYOGzatInZu3cvs2/fPiYQCDD79u1j9u7dy2zatEn0+9deey1jMpmYa665hjl48CDT0dER8veOjg7m6NGjjN1u79O47HY7c+jQIdHjf/XVV8yRI0f6fI5YsGnTJuapp55iCgoKmI8++ojp6OhgPvroI2bChAnM8ePHmZycHGbUqFFMTk4OM3HiRGbz5s39PiYFCYaBtrxDFUpCW0F/IFaSBmEVEj1X0jiY4zjRXG68vJ6BZnGSEgo5ofRYoxpy77Dyfg9uKIYvzpAKQ3V2diovioI+IxpJo7m5OYSMFG6EysrKMHnyZBiNRipqHi+G5dNPP42kpCQqnD4QLM7W1lY4nU7JOXK73VHDrXKh5FNRC6ig/6EYvjgjnDXX1taG7OxsaLVa5UVR0CsIvYvDhw9Lig0QYtK0adOiKofs2bOH5rz66ol1dHTgySefhFarRVJSEn72s58NWDcGom/bl9ZYcvV+ijzf0IBi+OIIsTCU0qJHQW8R7l3o9XpkZWUhMzMTlZWVEc9UbW0tFQyX8ww9Hg+Kiop6LE8WHt47duwYvF5vSGkNqQccKDz//PPQ6/WinSXKyspC3juxcKVcKJnneclNh8IMHVxQDF8cEb7Y9KRoVskZnN4Qu/9i3kVZWRnmzp2L6dOnQ61W0y4aDQ0NaG5uhsfjod6WXOmC3++PWZ5MKrxXWFgYwZIUE0yP5dmOx/P/6quvwmg0Ii0tDSUlJSGsTo7jaLhV6no++ugj/OQnP5GU9rNarTCbzaJ/U2oBBxcUwxdHhBs6uV23y+XCvn37RHf11dXV6OzsPIUjVzBQEN5/t9sNnU4Hv9+Pjo4O8DwvarhYlsWcOXOQkZERYiz2799Pw3ykZ2N4Mfm0adMi5MkMBkNIU+JwCIWXW1tbEQgEUFVVFaEHS46nUqmwc+dO7N+/HzU1NaL5MGKcd+3aBb/f3+ec2Z49e2AwGKDX6/Hwww8DkG5+K7ahKC4uhk6nQ2pqqqSYu1qtxogRI9DY2BhyPMXjG3xQDF+cIXyp5ER61Wo1kpKS4HA4UFFRISolpeQChz4aGhpQWVmJuro6GAwGeDweaDQacBwX0eCXIDs7m3p0gUCAakLyPB9CZiH5ZdKaijQGbmtrA/B/knUWi0XS4JBnmIyPGCfCFhUu9sSIk56NLMtGkGsqKytRWFhIyS9arbZXQs5CD7G1tZUyOV988UXZ33322WeyOVKTyQSTySSqPuNyucCyLOx2O/R6PRoaGtDW1qakLgYhFMMXZwjV2aV23dOnT8eSJUvg8/mgUqlkpaSUF2rogkQIxHRQfT6fpKg4aQXl9XrBsiwcDgdYloVOp6NGTfh9rVYLlUqFuXPnYuPGjTAajfQ3dXV16OrqCjE4X3zxBZ566ilcc801KCoqgtFoFCV0GI1GNDY20nMJ89lSm766ujpq6OQ2hlIeVLiHzHEcOI5DSkoK9u3bF/H9trY2PProo1i9ejVycnKg1Wphs9lE74fVasWwYcOgUqmg0+mooSeh5La2Nqo+Q7xEnueVDeoghGL4+glkJyrXqy4WKSklhDJ0QZROpO69VqtFRUVFzOzC4uJi0Y1SdnY2UlNT8dZbbwH4kYWp0+kQCAQizqnRaKDX6zF9+nRcc801eOihh2RDfyTHF27ExML8PUkFOJ1O3HHHHfjggw9w4sQJ+rmUhzxlyhQcP34czc3NuOuuu3DRRRchMzMTZrMZc+fOxY033oi///3vePzxx2Wv54wzzoBer4fdbhfNOwpzecr7OXihGL5TBDkGmZyUlJI0H7og1HuPxyP6d4/Hg+rq6ohNU2dnp6wmZnj+yWAwQKPR0NyenMFxu91oaWmh/25tbYXD4RD9rsPhQHV1NQwGA0wmUwgpRGwz19raSsk3RB5PjoBDxj18+HCkpaUhNzeXdpMQ85A5joPJZMKZZ56J/Px8TJkyhZYSsSyL7Oxs5OXlISUlRdSDTU1NxRNPPIFdu3bJettShpBct0JSS3wohm+A0dDQgIqKij7VHSkYvPD7/ZIkEXLvwxdTOcOVnp5Ow48kp+ZyubBkyRJ6vG+++SZmweZowtPBYBB+vx9lZWURz3B4KU9LSwu0Wi0l3xgMBtEuIqWlpcjJycG0adMwefJkuFwuGI1GJCcnIzU1VXI8LMti7dq1uOWWW/Doo4/ihRdewPvvv4/Dhw+ju7ubXg/JWQo3FMSDJKoufr8/YlzhqQfhfCmF7YMLiuEbYAgbzfanlFQ8oexq44euri54vd4eETzkjBHHcTRUR1RTdDodOjs7kZ2djZdffhllZWXweDwx15eeddZZEc9mSUkJXC5XSKG8mHhDVlYWtFotnE4nOI5Dbm4uDbESlqnFYompX180D9nlcqG5uTnqnJNwqZClSlp2CQ3XmWeeCbVaTclmWVlZIcSgyspK6kUrhe2DC4rhSxCQRrP91a05Hkoayq62f9DV1RVC6Y/l3osttD6fDyaTCXq9HsuWLcPWrVtRUVFBC9X1ej0sFguuv/56fP/99zG1rdq+fTvS0tKwcuVKGI1GpKenQ61Ww2g0QqPRQKVSUbJIOLFLpVIhLy8vpOWVx+MJeW5IWNPr9dIyATmjMWXKFFqkLwTx+NLT03HjjTfiyy+/lJ3vhQsX0vCn0WiE1+sVzZnOnj2btv8i3rndbqd9MPV6Pb13xJCS90uJ2CQuFMOXYIi3NxXeEVutVkOn04XsqGM9Z3/tagebB9lf4+3JccONDDFCer0earUaVquV1gS2tbWhrKwMKSkp2LFjh+Q5xWpKDQYDXnjhBQDAJZdcgvz8/BCPLT8/P8IQBYNBPProozAYDJg+fTrq6upElWbIc2O1WkPYoeQ74UbjN7/5DcaMGYP8/HxJD/mtt95CdXU1Nf5vvPGG6Pz94Q9/wHnnnYempqaouUZCQGpoaMDKlSuRl5cXMge5ubkwGo2iG0IlR5+YUAzfEIdYR2xSu7V+/fqYi+fjxTyVW2jDBb0TrTdgInq8wWAQW7duRVJSEoYNGwaTyQSz2Ryy0dFqtRg1alTEvQuH2MamvLwcpaWlGDdunCQbUqvVori4OMRjKywshFqtjqoZGggEoNFoRO+x0Gg8/vjjSE9Px6effhqTh/zVV1/hlltuQVZWFqZMmYIHHngAx44do/eQ4zg4HA4YDAYsXbpUksDj8XjQ2tpKc3xiNYC1tbWihpjkERPl+VXwf1AM3xCGnLEidV9VVVUxFc9H6woQbVdLQrlCoyFGbCgrK4NOp4PD4QgJJ8XLwPTFW5OqtaypqRnQxW3BggVgGAaFhYWw2WwRi3B5eTk2btwoewy5Z0Wj0SArK0u0/q2rq4uqnRBDy7IsvF4v0tPTo2qGFhQURCX3vPzyy0hNTcXbb78dMeZo9/LEiRN4+umncfbZZyM1NRVFRUURz1xhYaFk+FRouMhcxEr+YVk2RCVHQeJAMXxDGHIdsbOysugOPpbi+d56fFLkHbLTl2MLCnfOfQ2p9tVbi2YY4mmcY0UwGMS+ffswbNgwDB8+HOnp6b1mBzc1NSEzM1P0by6XCw8++KCohJqwIJ2cq6ysLCTvJRdGnDp1Ks4++2wUFRWFHINIq23ZsgVmsxnPPfdcT6cnAm+//bZkmQLP8zGRfcLDsqREQwwOhyMmso2CUw/F8A1hyClj8DwPq9XaI4MmRajwer2SC75UuYbcgiFWJBwIBPoUNpLLT8biOUTzeBsbG08Zi09oxDMzM6FSqcBxHP7973/3yCvv7u7G66+/jtraWuj1+qiel9/vDzFygUBA9jd+vz9k4yKc+6KiIhqSJZ0diEwby7LgeR48z8Nms0Gn0/V4UyF2T1taWuByuUS/n5GRgezs7JBOE0TVRnhdWq0WXq835jlQwpyJiVNm+AYbgWGoIHyxIgt+bW0t1Gq1bPF8WloaUlJSUFZWhquuugovv/wyCgsLwbIs7QpAiAtiC75cgX5PDC5ZtHtLFIiF/i/mBQqf2c8//1zSQ+U4jnqop2KxEzPixEOKZU6/+uor3HbbbRg3bhw8Hg82b94s6fkLjXlXVxfVq/R4POB5Hg6HQ1LhJLwpLukXaDQaYTabQwxIe/uP9Xt6vT5Cc7QnJCoxz37jxo14+umnMW7cOMl7qFarwfM8li1bhscffxzr1q0T3eQRY81xHJUk1Ol0vR6vgoFBvxu+RCQEnE7o6uoSNVZEZ7CkpETSK9Tr9Vi9ejUKCwuRkpIChmGgUqlipm1Hk2QTC5MVFxejpqYmYhx98fh6UvA9ffp01NfXY+PGjTAYDMjKygLHceB5HizLihY18zxPC82dTieampro8eO94Yu2YSBeVvgivHHjRjz//PNYtGgR9Ho9li5dihdffBEnT56kxwhnioqRRubNm0fzw5mZmWBZNmLj0NbWFnKvhEot0RRbOI7rU887MeNdUlICk8mEyspKmM3miGeuoqKCaqZec801OOOMM5CVlQWfzweDwSCqa1pRUYHq6mr4/X6UlJRQIovH46E5TmWNS1z0u+FTCjsHHseOHcOYMWOQnJyMtLQ0GlYiizbHcTEVz7///vuSOUO3241t27ZJJv7FnoOKigp4vV660BoMBmi1WiqZRXbZOp0OLperVyGv8HEIIcwnCj8T6xiQm5uLtLQ02gXB7XZTw0CkxfR6PaxWKwwGA+rr61FfXx/3DV+0kCvxsjQaDa2ZKysrg8vlwoQJE3D77bfj66+/jjpfUsZ69OjRtHvBkiVLUFZWFtGuKDs7W/b9lrsGu90Oi8Ui+jebzYaVK1fi9ddfDzHYwnHL6Z7ecsst2LNnDxXqDhegPv/88zF27Fi89tpraGpqQkNDAywWC0aOHIk9e/ZEHFOn00Gr1VI2K5m3voblhwISPcLXr4ZPEV9ODIQrVZBch0ajoYvy5MmToxZQi93Prq4u2qKG0MPJb7u6upCZmQmfz0eZowaDgRYAZ2dno62tDcFgEI2NjSgqKkJeXl7EzrmlpQV5eXlYsWJFrzdNsYo6B4NBUdp+S0tLiHdDCpfb2trAcRzKy8tDjp2dnS3alSMeNY/R3qnjx49j5MiRKCkpgcFgwKpVq/Dmm29S2a7e4vDhw7SAe926dVCpVKLtirRaLTo7OyUXv2gbEdKxJPxver0eGzduxJgxY5Ceno41a9bg+eefx/HjxwEAu3btkiToxKKp2d3djcceewxpaWmor6/H+vXrodfrYbPZIjYuXV1dMBgM9DkO//vpWr83WCJ8/Wr4+kqBV9B3iC0y4RqKZFH2+/1Rd2nhBqSuri6iC3dZWRnq6urg9/tRVFQUEgbSaDTIycmBSqVCVlZWyAuSmpoKjUYDi8UCnU6HtWvXUq/J4/FArVajtrY2wkuLBcL6rczMTBiNRvA8H9HGp7GxUbQDd0NDQ4QXKOzRJpxfOVJRPDZ8a9asEWVSVldX48orr0RaWhqSkpJw//334+jRo306lxAvvvgiXC4XVq5cidbWVlgsFtFoTkpKCi644ALJ3CkA1NTUSEYZTCYTpk2bFpFfKygoQHNzM4LBIFpaWnDjjTciLy8POp0ObrcbGo0mLkSTL7/8EuPGjYsYHwlvBoNB0TA9Gf/pvLEfLBE+xeMb4gjffPT1ngjzQG63O8I7It5bcnJyyCIUHgZyuVxwuVzQarUwm83QaDSw2WxU3ioYDIq+RJWVlbBYLL3eNK1duxZXXHEFgsEgFi5cGLF4kcLrWA2ZRqNBampqyOfhHQ2E3oXD4ejzhi8jIwMulytk05CVlUU9opdffhlms7lP5xDDLbfcAofDgSeeeAIdHR0R89TV1YVVq1aJNqANX/wOHjwIjuNEe94lJyfTjgrCvHRRURFtmrt27VrceuutGDt2LDweD+bMmUPD4qWlpbLnjhaGC39HyDNPQtl6vR5arVa096HBYKAycacbBtN6r+T4hjjCH8Z4eeHBYBDbtm2jxyJqGmSXr1arkZKSIhriIIxAUuNlsViwdetWqNVqaDSaqIofarUaHR0dPZkGutht3LgRa9asQUtLC1JSUqDT6ShtnnibS5cuDfE4GhsbYbfbRY9rtVojQnPEKAjDuyQ8yrIs9u/f36OxC8f/0EMP0fMJF3Dh4vL+++9j1KhRPT5HNMybNw8qlQqHDx8WbVfU0NCA4uLimLzdhoYGZGdno7y8PETxpaysjHr+UiSqQCAAn8+HCRMm4KWXXgoJ4XZ0dGDGjBlgWRZWqxUsy2LGjBloa2uLOQwX/o5IlfGIrWFWqxXV1dUJF9o7FRhMEb5TxursL/FlBdHR0NAg2TCUoDe7MmJU29raRDsM+Hw+1NXVRZxHqGhBQph+vx9WqxUmkwlqtRp5eXmyhcGxvkTCxY6owaSkpCA5ORlmsxlWqxWZmZnQ6/VwOp0444wzoNFoaDcB8je5/mxz584NkexqbGyEyWRCdnZ2REjZ5/OhqKgopnxbMBhEc3MzLVPIysqCSqWC0WgUfX/I4vLCCy9g/Pjxcd9h22w25Ofn07EJNyZypSvh4xM+N8K1wWAwQKfTyXZgIMcQe17Du7PrdDqcd955WLRoEa15DM/FVlZWYunSpXj66afx0EMP4bbbbsOVV15J73e0fKRQ9D0RPZtTCcXjE0Gis3yGMvbt2wej0UgXGJ7nI4gXJSUlqK+v7/GxGxoa4Ha7JXMrLMuGCPoKC97J30nejuM4cBxHvSgpbcievERSu3XSCker1WL48OEoKyvDPffcg7fffpuSJYTPrFTkgrA3iTamRqOBTqeDSqWSNJYajQbl5eX44osvEAxGds0QLuAZGRnQaDQhVHqpDYXBYKCs0/T0dNm6xJ7i66+/xsiRI7F582bRuW1tbUVWVpbs4kfKPoRlLuGarE6nEy6XC3q9XpQtKczvpqenY8uWLSFepFR06dChQ6LKM+R+uN1uFBQUYNasWbj00ktRVlaGiooKNDY2Sup4Wq1WSuiSq2c9nTBYInyKcstpgJ/97Ge4/PLL6QLS2dkZ4YWfccYZvTJ8nZ2dUKlUkjt0u90OjuOQnp5OjZxwASdGkHgT4czI8HyNz+fD2LFjZRfwY8eOYc+ePbj//vtlG66SkNno0aOxZs0aXHnllfjlL3+J3/3ud7jvvvvw+OOPY+vWrXj55Zfx5ptvUgMtjFzU19eLyr1NnjwZ6enpouOzWq0YPnw4rVkjvei0Wi02bNiARYsWyTJCxTYU06dPF22tQ4xzX5l2//nPf6BSqfDOO+/Qz6TyvVKbDZZlsWzZMnz++efQ6XSi9X96vR4cx1FFF7KRIOSVgoICHDt2LIJJvG7dOknDZjQa0dTUJOmJ2mw2lJaWYsmSJTj77LMxYcIEmEwmsCwbsgELDy2TkHv4s3w6Y7BE+BTDN8TR3d2NrKws7N69O+Jvwhf54MGDcDgceOqpp2L2DP73v//hpptuooQDKY+vra0Nu3fvBsdxIW1e/H5/SE82spATtLW1gef5kFAYx3FITk6mC94ll1yChx56CFdffTXmzJmDUaNGQa1WY+zYsTj77LPhdDpFxy4Mmel0Otxyyy341a9+hU2bNmHt2rW45JJLMG/ePMyYMYPmk9xuN8xmM1QqFZKSkuj1SIXBpLxgnucpuUcsFCoXViX3xOFwgOd5urjIqbZwHIeKiope78KDwSCqq6vBsqxoiJY8L6QTuzC3Sfr3ZWRkQKvVIikpSbROsqqqCllZWVQFhZRJCMta2traMH36dOTl5UUwiSdPnizKxiX3uqmpSdYTnT9/PvR6PWbMmIEHHngAhw8fxrFjx/DUU08hJSWF5miF/01NTY2Q1lMiWj8i0SN8iuEb4mhsbMTYsWNjyim99NJLMBgMkjT0b775Bv/+97+xbt06nHnmmVQNg+gaiu3yiRDxsWPHcNFFFyE/Pz+inxmRqhJ7STweD5qamtDY2EjVNITnKC4uxpgxY7B27Vo89NBDaG5uxg8//AAg9pxDbxLvXV1d2L17t6SnSzy48AW6qKgIU6ZMgV6vl2WKhgsBEMFxYX5LGCKUIxbE2u9O7BqFtZcsy8ru3knvR2G4mnhv5HkqKyuTrNMjdZJSz5Lf70cgEBDdbEh9LrzOaGG4o0eP4q9//Stmz54NnudxwQUX4OGHH4ZOp4vwwMvKymjoltyjcNUeBYkLxfANcVx66aUheRk5NDQ0RCT/p02bhuLiYni9Xqr28utf/xq7d+/GN998gw0bNiA3NxeVlZURsk1nnnkmDAYDhg0bBovFQtmTRBfR7XaD53ka6uN5PqQfIPEELRaLaM5PWDxPyAzEixRek9xi15edejTD2traioULF0Kn0yEjIwNqtRoqlQoMw8BisUQ1VMJWTsI6RrFcktxY5PrdhavtCNHbfA3JWy5fvpx6tIT0M2nSJEnPzGazwWQyRSWTiLVHAgCLxSLb8LYnYbivv/4a99xzD0pLSyUNKsdxNG0gVO1JxNCeglAohm8I49ixYzCZTBH1RmKQWzg5jsMzzzyDY8eOAYj0BHQ6HZUeI8aMGKAvv/wSo0ePjmg7Q/JgUv0ASd3WJZdcgtbWVjz33HMROTOiSCMWFhMqbJDFjtDbhUSRvibevV5vhFdXUlICr9cbMretra04fPgwVq1ahREjRmDkyJG0pCOcsq9Wq+F0OqHT6SKKqKVySceOHYPH4xHNiUrVnIV3Fxces7cMvZMnT+LIkSNobW2FXq+PKOngeV7SkKjVaowYMUKW0bl9+3bZhrjr1q2Lath6EoaT6yJCpOrCw9XCQncFiQnF8A1h/OMf/0BFRUVM3yWsPDGEhwKlPIHa2lo89dRTeOaZZ/DYY4/h9ttvD6GGCyHM60mRIQwGA+1nFl6GQRZmqbAYURchi8+9996LMWPG0HxkPBLvZEzE0yXHrK2thU6ni6g1fOaZZ5CZmQm73Q6VSgW9Xh+iqUmMvU6nQ2FhYcxlJ8ePH8esWbMwf/58UUHy0tJSZGdnRxhnwgwl92/dunX47LPP8O677+Kvf/2rbAufGTNm4KyzzkJhYSHGjh2LjIwM8DyPYcOGQavVIjU1FTabTfTeulwu0Rzf+vXr4ff7JXOjBoMBhYWFMBgMkio6VVVV+Oabb+i972uuSW4DQDx48rfwQnfF+0tcDFnDJ6cTmMhJ13hi9uzZuP/++6N+b9++fTjvvPNiKh8gC4GYp6LRaDBmzBiUlZVh9uzZmDZtGtxut2xoS450EF6oTgqk29v/jz4v91uz2UwbxBL25l/+8pe4PQPCvFp4zR1RoWloaMCnn36KefPmITs7GwsWLJDsUefz+ZCfnw+WZWXr4ZxOJ1599VUAP3YYX7RoEWbOnIlDhw5J3hu1Wh2115xarYbdbsfYsWPpOKTIIA8++CCeffZZvPLKK3j33Xfx2WefobOzEydOnADwYyE5afkUfgxCWmJZlrYpEuq7itWEkm4ifr8fr776KnQ6najqC8uymDRpEt544424aUaKGe+qqiqYzeaQsOtgofL3FENxzRxyhk9KneHYsWODQjw1HggGg3j99deh0+lw+PBhye99/vnnqK2thdlsxi233IL169dHfXGbm5tDvCbhPLrdbvz0pz9FaWkpeJ7H3LlzcfXVV8uGtrZv3y65wGdkZISQBYQsT7fbjeTkZMkwlNvtpjqhlZWV0Ov1GD16NF2Y44Fwb0Bs4SsvL4der8e1116LgwcPUsMkl8dyu92ynoZWqwXP8zjnnHNQUVGB0tJSfPfdd9i6datkCUVWVhYuvfRS6vWKIVbPXm4hF75/Vqs1YtMjJOps374dPM9HeMZEBUho2GpqamjuN7wOULggp6en49xzz4Ver4+bSLhUbnD8+PF0czCYirfDIWXYBovgdG8w5Ayf1MsqVeM0GHdjsTyopFea2IPa2dmJTZs2wWQy4ac//SltUxNL8l+usa1KpcJFF12ELVu24Ouvv8ayZcuQlJQEvV4v+pvc3FzZkJ5araatiAiVvaHh/zqmr169WlaUmOiCBgIBaDQa3HPPPfG5AQKQ503OmHEcB5/Ph5SUFKSlpckyMD0eD6099Pv9mDp1Kp577jlaO0b6FR45cgQzZsygjVvHjRsHnU4nWbeoVquxaNEi/OUvf4l5ge7q6sLq1atpLjCW0LDw/ROGp4XPJiHqrF69Gg6HQ5IJKRe1kSufWbhwYVzED6TGQwgtGo0GRqMRPp8vJsWaREM0wzZUPVhgiBk+qRciFqrzYEBfH9Rjx47h1ltvRWpqKmpqaiQ1I3uz4KjVaqoWUlFRgZEjR2LYsGEwGAy44oorMHny5Ijck9frRV5enmjHbdLXTdg0N3zR7erqQkFBgawgMll0rFYrli9fHvfdKgnNkTCiGNLS0pCVlYXx48fHpEM6evRoZGZmQqfTQa1Ww2q1Qq1Wg2VZsCyL999/Hz/72c8wceJEPPTQQ0hPT8fkyZORk5MDg8EQQbaZNm0aNmzYQM/TkwXtd7/7HUaOHIn3338/6nsi9nw0NPzYCUQY2iVsXKKn2RtPQuwaiouLkZubC6PRGHN7ot5AeG5yLURjdjCtMXLPwWD2YGPBkDJ8Ujvp1tZWyUUpUXdjYpB7UI8cOSL7oP7pT3+C0+nE7Nmz8d577/Xq/NE8FVJjVlRUBIPBgL/97W8h7MrwJq4TJkxAUVERdu/eLZmzEXpuYi9bc3MzzGazqFci9PwMBoNo/72+giwQe/bskfRcdTod0tLSsHTpUqxatUoyx1dVVYXx48dDo9GEFLcLF1diBA0GA0pLS3HmmWfij3/8I6688kqMHj0aFosFDocDGo2GMm7D83nhiitCFm44ysvLMXbs2JjmQuz5IPddaBTi4UmIRSe8Xi8KCgrw4osvSkYCOI7D//73v5jPEw65zbVOpxNt95WIHlI0r9liscgKAgyWNVMKQ8rwDWWPL1qdVnJysiyJZPLkyRFFzPEcg3Ae5eZV6E12dXXhrLPOgkqlgtVqFfU0yUsm9bKRMdXW1qKoqChCxosYmIaG/umTJlzspbwQl8uFXbt2AQj12i0WC/XoiME+cOAAZXzKGYri4mKYzWZkZGRg9OjR2LRpU0iz2W+++QYPPPAALrroIhgMBkyfPh0PPfRQiB6osJuGmNfV3d0NjuOwZs2amOais7NTlMEr7HEYb09C+Mx0d3djxYoVOOuss7Bhw4YItRqfz4ecnBxYrVbcdNNNOHLkiOSxpBBt87dkyRJR0k6iQe46srKy8Mwzzyge32DCUM3xRXtQ5XI3PM9HvOS9xbJly6L2WgNi3xV2d3dj7ty5shuTcHHicAjr+YxGIy0WN5vNEQu6y+XCk08+GbcXl7AXA4FAiBdC8lizZs2iotdCkEW2ubkZw4cPx6ZNmwAAN954I1JTU0PYonJh0cLCQjz++ON44YUX8Pbbb+OTTz7BN998g5MnT9Lvfvvtt3jssccwc+ZM6PV6XHTRRbjwwgujvg/vvvsu1Go1/vWvf0Wdh++++w5z5syhrYaExy0pKaH3t79b13R1deH888/HggULsHHjRhiNRjidThoRIEIKdrsdGo0GF198cY+Ib9E2oESKjzTMTVTEsgFRcnyDCFIEDfJwJ6J4aiw7TbEHlahhqFQqqlwSi1HqDU6cOIF//vOf0Gq1VEWFGBgxanxPdoWXX345MjMz4fP5KBU/EAhQj83n8yEvL0/y98Jwmsvloo1Md+zYETIGsjgR5Zj6+vpe3X9h6QKp2RJ2UAgEAigoKIjJUyJeVWpqKl555RV4PJ4Qjy+ahzFlyhSceeaZKC0txYQJE+BwOMBxHIYNG0ZDmURZZ968eVi0aBHNSUbbzd96662i8mLh6OzsRGlpKRYuXIgvvvgC48aNo/k7rVZLw7KFhYXYs2dPv3sS3333HUpLS7F27VocOXIETU1N2LZtG8aPHy/aOsvhcPQoRClV3lBTU5PQxi4c0QzbYBGc7g2GnOEjGAx1fGTBJurz0ZL85EEVKp3Y7XZotVp4vV7odDraz85qtYLjuD4/qO+++y6WLVsGlmWRlJQEhmGwceNGHD58GNdddx14nu+zALLBYMDHH39Md+ak1sxgMNDefCqVSjIPBSBEy9Lv94v2wquqqqIhT0Ke6U1eiYQpxRZRs9nc4wVi/PjxuOWWWyjRw+v10hyfHBGGeMM8z+N3v/sdbrrpJmzatAl1dXVYtGgRKisrMXnyZGRlZcFkMiE5ORnDhg2juUIxZGRkYNOmTWhsbERxcTFMJpPs2Nvb2zFhwgRcdtll2LhxI61hTE5OhtFohMPhwLJly6jKj1qtlmT5xtOT+OabbzBhwgSUlpbSEpHeanyGrxlDxSDEeh2JtGbGC0PW8CU6CBuQMB1j6elFHlSe5yOkrIRKHIFAALm5uZJ9xAD5h/nLL7/EzTffDKfTSaW1li9fDp/Ph/T0dHzzzTdYsGABJk6ciPfee69PiwDxaOrr60XltvLz89HV1QWPx4OioiLJuRF6xEQsWagDGr6paG//sWZOLoQajlhKFwwGQ4+7w8+ePRtPPPEEUlJSoNFowPM87e4gpMxLGYr09HRceOGFuOKKK/CrX/0Kf/jDH/Doo4/imWeewc6dO/Hee+/h888/x9GjR9Hd3S0b5uJ5HrW1tcjPzwfDMOA4DhdddBF++9vf0v54BIFAAG63GzfccAM2btwY4T2UlpYiLy8v4r6WlJQgPz+/3w0HYQQTz1ms5pMwfsWQnZ2N6upqyRDoUDEIQ+U6egLF8A0QpOrhSJ5K6iGMlRwQCASQnJyMjz76KOR7Qlklp9NJ1UWOHj2Kv//97ygtLUVycjKSk5NRVFSE+fPn09orlUqFiooKuN1uXHbZZVS7k4yrNy8PqfWSkjUj+bNo7M7u7m7YbLYQKn8gEKAF22K/yc7OhtPpjCmvJJz3eOep/H4/NmzYAJZlQ56JQCCAqVOnIjU1FS6XC2q1mpaDhDNXezrv0cJcO3bsgMViwZo1a3Dvvfeirq4OBQUFVHz8nHPOAc/zuOKKK9DR0dFjUhnHcSHdJeKN8PekN8Q3jUYTURoyVHJcpzsUwzcACNedJCCLWEZGBh544AF88803IX//4YcfMH36dNkdanNzMw3H2Ww2WgBOdqn19fURvcWysrKoGojJZMJVV12Fzz//XHRx9Pl8mDNnTlznY9myZZLXZLVaMXXq1Ii6vHD8/ve/p/nGcN1MuZ55sXp84fJkcsXq4X0Fo+E3v/kN7fMndkytVou33noLJpMJ+fn5EczVeKqRkOfk6quvhsPhwPbt20N+d/z4cdx9993Q6XSYMWMGCgoKJBnFcmVEVqsVb775Zshn8fQ8xDYnUs+zmPZneF9EMrbwzdfp6C0NBSiGbwDQ3NwcEYYkL1BWVhZUKhXMZjNGjhwJq9WK8vJylJWVQavVYurUqbJG0+/3R7zcJSUlOO+883DHHXfQ1kLhuS+tVovHHnssJIxzqujMu3btkt11L1u2TNa72bdvHyVyCOeSfM9isUSE23qa49u8eXPUWrTy8nLk5ubCZDLhmmuuCZGLk1og29raUFBQAIZhkJGRIXrurKwsWCwW3H333XHPLYmNKxgMYuzYsUhOTsZnn30W8p1//OMfSE1NxX//+188++yzmD59OtLS0kQ99mje1LJlywD0jzSWmMfX3NxMN0I2m4126iAd3YVlCGVlZUhLSxMdm8Viwd69e4esnNfpAMXwDQCECvRiUk6jRo3Cb37zGyxYsAApKSlgGAYMw1CPQGzRLSkpQW1trSz1PTk5WdKr0Gg0IdJR/U07F6KjowMsy4oaZJZl0dHRQXfh4R7ssWPHMHHiRPz+97+XzbutW7eONkglebRYWJ133303jEYjRowYgby8PFH1EYfDAb1ej7y8PKSmpiInJwe5ublISUnB//t//w8bNmyIWCA7OztxzTXXwGQywel0IikpSZYm/8ADD9DP+svLEOZG1Wo10tLSQhrJ6nQ6GI1GXH311Rg7diwmTJiABx98EMeOHcOoUaNEQ/diZUTTpk2jItX79+/vN9q8sMyFvF+k5IBlWQQCgZC5JEXoO3fuREdHB3ielxQT93g8g6ZYXUEkFMN3ikBeMJIPIS+UWOG1z+cDy7K0aWlSUhIsFgsNG4WHqdRqNcaPH08XdjF4PB5MnjwZRqNR9O9WqzXE8J1Kj6+1tZXWrgkNhHCHTQxGW1sbpk2bRheYyy+/HPPmzUN3dzeVx5JajILBHxukNjU1RR3/o48+itTUVAwfPhwXX3wxvv32W9HwYHjN1okTJ/Diiy9izZo1MJvNoiUmFRUVMJlMWLBgAf79738jKSkJKpVKlCDi8/kwb968uM21HBoaGkLYsGLz6fP5cMYZZ+C///0v2tvbsW3bNlx11VVITk7GrFmzKJOThJ2FZURkY+fxeFBfX09ZnhqNJm7PWbhAglinBxLeFHpser2elnkQhjXJafaGCaogsaEYvn5GeKhEr9dTskVBQQF90fR6PaXrEy+lsLAQdrsdJpMJ48aNiwgnBYNBPProo1Rtf+fOnbJkgpaWFqoVKfb38Bd22rRpEQ1k+2NXGwwGadjJYDDA6XRSNZbk5GSwLBtSyhAIBKDVavHwww/D6XSGiGzPnj2beii9CQU++eSTsNvtGDZsGObPnx+RZyXjjcXjOnz4MHiel8zb6fV6WgvJ8zzNq5JyFI1Gg4ULF56S8BnJO5MwerSND3l2SemJ1WrFsWPHEAwG8ac//QlmszmCqJOXlxciOm4wGKDT6ST7QFqtVowfPx7XX389du/eHUKmCocYaYvUWYpdA8uyId3a6+rqIogsU6ZMoaLiwnsdjQk62OW8Tgcohq+fIZVQz8vLC3nxyOeku3Z2djYaGxtpiDInJwdjx46NaLVSUVFBDWIwGATP87L5LKvViry8vJC/l5WVob6+PmTczz33HEaOHImCgoJ+pZ0Hg0G88847SE5ORklJCQKBAM3FCMOI2dnZWL9+Pd1EEM3KxYsXh4znkksuwa233iprmMQM13PPPQeXy4WkpCTMnDkTBw8e7PO1yYWL7XY7lZBrb2+H1+vFnDlz8J///AcMw2DEiBF4+umn+zyGnozV6XTS8TY1NcHpdIrOodgzJOw6v2bNGhQVFaG6upoKHbAsiyVLlqCzszNkA0U2fVIhXo1GQ0P9DMNAp9Nh6tSpuOaaa/DEE08gEAjgxIkToqSt7OxspKamAgi958FgMGSDGG7kSciXiISTciPy7Cse3+CHYvj6iPBFNPwF6ymFmmVZ1NbW0hfI7XbTLuRSTLz6+npqXNeuXUvVMsIFn0luw+/3UzKImDHbt28fOI7D+PHjceLEiX7JKQk94fT0dKhUKhQUFMBgMIjWrVVVVYmKAFdUVFAPtLu7G3a7XZJVKUZUWLJkCUaNGoWkpCRUVlb2SsBYan4+++wzyTINwiYN79pNujDY7fZTSpggHh/Ja5HNBSl3EdY/SkUN1Go1Pvzww5Bn3u/3Y8qUKeA4TlKEfPny5bIh6h9++AH/+te/sGTJEjgcDgwfPhwMw2DkyJFQqVQYPnw4tFqtZI6Y5L6F9zw1NZXes/ANSnjIVzgmUhs4evToISvndTpAMXy9hFgI0+v1hrxgy5Ytg91uj/itXKiE5EH8fr/kDlJOScLhcIQIPgtrpex2O5YvXy56DIIvv/wS6enp0Ol0+Oyzz+IwU+KQIjQsWbKkxwobPM9Tz5EwO2M9Z1FREVwuFz7++OMeX4MUG/G7777DbbfdhtTUVOTl5UUsoMIuEVLyV0Rh5lQsph9//DE2bNgAvV6P1NTUiJCfcDw+n08yT2yz2XDvvfeGGJFjx45RoYbCwkJaBC8kK3k8HvA8D47jYoosdHd345133sGVV16J/Px82l1eitTl9Xoj5l8oHCEMiZKoiVSIWq1WY+3atQktgaggOhTD10uEL1hiOYKioqJeJcfT09OxY8eOHi96hLxhMBiorJnQMGu1WtTW1kp6b8eOHUNhYSGMRiP+9re/9XRKejROOeUQqZowuVq/W265BTfffDNVr+nJOXsbnhIzWtOmTUNaWhpmzJiB9957T9RL53kebW1tMY2pv8Jn3d3d2LFjB+bMmYOUlBT85Cc/wTvvvCNJNNFoNFCpVHC5XJLMYLVajQ8++IBeUzAYRHV1NSoqKtDS0oKamhoqXcayLLKzs9HW1kZ/X1FRgerq6h5dazAYxJNPPin7bIR3JRF63GRzkZeXh+LiYixZsgQqlUqy/pDI4gnPr9TxDT4ohk8GUg91rKoQZDEQ07IUa75KQimEct3bHaTf7w/p50aOX1JSApZlkZWVFUKmAX5cCJcsWQKXy4VLLrmkx+fsCeRyX263u8c1YSzLIjU1FWPHjpXsJNDa2gqXyyX6t94QEuTuuV6vF31mCKO0oqICRUVFMXftjidh4vvvv8cDDzyASZMmYcyYMbjzzjtx9OhRAD/Wl0ot+DabDRqNht6f8E1eSUkJFdtetmwZPB4PDd8Kw6XBYBDPPfecpPGM1cgLvW0pHU5isKVUe8ickveUMHD37NkDnU4nGs6NZXyKMUx8KIZPBHIFtSdPnsTf/va3kGJjuYXcarVi8eLFITv+uro6aLVa5OfnR3QlnzZtWp9V3uUKwoU7XZ/Ph8LCQnR1deEXv/gFsrOz4fF4+v2FbW5ulqWwr1u3LoLEQ8JiYnmXsrIyWu945513Rpzvo48+wtlnnx1X2nxP6hzDnyeVSgWDwYDk5OSoJIn2dnHGrRSkFt2Ojg78/Oc/R1paGs4++2w8++yzIa2L3nvvPeTm5krOEal7A34MjZJegoTV6fV6cfDgQZSUlEClUslqi8ajMXS4t71q1aoI/dqysjLJkGX4PU9LS8OIESNoftPhcFCmMWFal5WVyUZg+qMQX0H/QDF8IhALYZWVlcFut2PkyJFISkqSZYURkEWroqIipN0OEaIW60oejxflv//9b0x0a7IjPuecc5Ceno6UlJQIGal44/nnn4fFYkFGRkYI+zQQCKCoqIh6oevXr6d1iTqdDhzHIS0tjeaYwvMqmzdvxsiRI6HT6XD77bejtbUVn3zyCdatWwej0Ygbb7wRGzZsiAsh4cSJE7jzzjtjZvaJPU8VFRXgeR7nnHNORP+68C4SJIcpB6lF98UXX8ScOXNgMBiwevVq7Nu3j/7m2LFjePjhh1FcXEzVV6QKtsNDyOTZXrt2LW677TY89NBDePjhh3HvvfdKEnrIvPSVFSn2vhFBAVIeJEb8Cp/f8Gsxm82i1240GqHVaqMKHvRXIb6C+EMxfGGQM2Isy+LOO+/EN998E1OOr7i4GDqdju6KMzIyRI1bPEMjjz/+uGyyP7zVisViQXJyMiZOnIhf//rXfT6/FLq7u3HzzTfDarWivr4eOTk5uPzyy2mBulqthsPhoKSgM888E7Nnz8arr76KYDCIhx56CMOHDwfLsvjss88i5svv9+OnP/0pMjMzoVKpkJKSQj0SEmrrKyHh22+/xW9/+1uYTCbKJAwvHQlf6OSeJ41Gg3vvvReZmZngeR42mw0cx4Hn+RAxao/HE+EFhT8zYosuIXGQukhyre+//z7q6+thNptx1lln4Z///CfeffddZGZmRuQkCctWbI4yMjIwZ84cLFu2DBdddBGWLFmCmTNnIj09XXT+SImOz+fDxIkTRcfL83zUeyLVaQH4MVS+bds2UeJXdnY2JdAIc4uVlZWYNGmSLEHmww8/lBwP0PscshIWHRicNoYv2gPW0dGBP//5zzjnnHNi8pbEFgiv1wuO45Cenh6iNOLz+XD++efH7QEXu5avvvoKc+bMgVqtRlFREZYtWxYRLqyqqkJ9fX2ERJrBYMDUqVNDQl/xHF8wGMSFF16IKVOmYMuWLTCbzXj//fcB/GiwxOoZZ8+eje7ubnR2dqK6uhp6vZ6WPZDwrBBnnHEGLr74YipRJbfzjmWxEX7niy++wNVXXw2O46BSqTBu3DgYDAbcf//9UQ2pmC4rgdVqRXp6OtRqNS699FIqcC0mo/Xhhx+itbUVnZ2dEZ4dYSWSqAK5rvDQNumskZKSgvnz5+PKK69ETU0NioqKoNfrI6IYra2t2LNnT4+8s2iGnmVZsCwLnU4n+b7IeUnvvvsuFi1a1GOPkVwPmb/weyZ3nxwOR9Twa1NTkySjWCx8q4RFBxZD3vBJPWDff/89XnnlFVx99dXIzc2FwWDA/Pnz8Yc//EFWBFqutIDUQoX3LYsXMy88oa/T6eD3+/HPf/4TKSkp0Gq1KC4upkaNUMSzsrKgVquRnZ2NVatWiYZzampq+jS28PGRua6trcXYsWNRXV2NAwcOwOPx4LHHHgMQfZfs9/uh1WqRnp5OvbaWlhZMnjwZ1dXVlCzy3nvvwWw2U8PRk523VGkImWOtVgutVgudTofCwkIsXboUWVlZEfJuUoZUqMsaPp5wZqFQ17KzsxOLFy+GWq2mXjmRQCMbmmAwiMbGRni9XtrDL3wRFSNx2O12TJ8+HX6/H3fddRdeeuklfPnll7BarTHrbZaXl0sap4aGhohNV2lpKQ1dk5BqSUkJNBpNTO/LK6+8glmzZiEtLQ033ngj1q1b16ewotjmrDfPDTGk4RuHaL9XwqIDiyFv+MQesNLSUuj1ekyYMAE//elP8dJLL+H48eOyv4nloXz//feRmZkp+rd4MPOkRHd5noder8e8efNEc5MVFRVgWZb2dAvf1cv1uevp+MSM6owZM3Dy5ElccMEFuOyyy+j3SXhNDA6HI6T+qq2tDdnZ2dBqtXC5XGBZNkTai2VZWK3WmEknUhui+vr6CA+0uLgY559/PsrKynD22Wfj0KFDMc0H2QgtWbJEshs8AVkg165dS0Wi7XY7bRfldrupkolUqUpLSws9llRvR4/HExIKJOjq6sKwYcNgNBojvKHwELHBYIDZbMaf//xn0esmyifE0KnVarhcrpDw4vTp07Fw4ULZ6Mr777+PLVu2oKSkBFlZWfjjH/+I7777LuT+xbOOLpb3Pvy5IfJ45J7EIlzdH6U1CnqGIW345B4wuT5sPX2pjh8/jhtuuAE6na7PNO1o1yJFPliyZIlsbtJut2PFihWwWq2ii77RaMRrr70WdQxSno3UXAcCAfA8j+uvvx55eXn4/vvvcfDgQfzyl79EWlqapDekVquxffv2kBzWtGnTUFdXJ9qB3ufz9cjjE1vkSkpKZGnxP/3pT3HixAnZOSLo6upCTU0NNBpNiPdNDJjY82Sz2XDGGWeIhqgbGhrQ2NgIu90umc8Teu3kvvv9/oh5JSLMwjG8/fbbSEpKgt/vly3jIZ83NzfDYrFg69atknMQDAaxc+dOSaaoVquV/BvP87QDxF/+8hdaCiHWQileKYRY3ns5wQHhe2W1WiXXjZaWlriW1gwEBntuckgaPnJTmpqaJJPgwgcs2otO1E/CO0YfOXIEGzZsgEajQXJyMlatWoXLL79c9MXoTYmCcFytra10wRJbKPR6vWQ/OofDAY7jsHv3bqqaL2Y8dTod9u/fHzGOaPmIw4cP45577gnJkQh/Q8g9CxYsoPJRK1asoD3NxDxyspMmOSxi9CsqKiJC0STkp1Kp4HQ6JcO58+fPp/V0YvPY2Ngo6YG4XK4eLUhiu/+KigosW7ZMkmKvUqlkN06kVU6soXibzYYdO3bQv5eUlFB2ptAbCQaDWLt2LRiG6VET3eeffx5GoxEvvPACAPH3qKmpSXROg8EgzGYzzjnnHNEOEC6XC1u3bkV3d3fU5y/ei3D4e0/SGE1NTSEKLyTUHh4x6ejogM1mw86dO0OO+9133+Huu+/GGWecEdfSmlOJoZKbHFKGL/ymkJ5fJMRCQB6w/fv3U+KE2E0U6igKmYd6vR6TJk3CiBEjoNFosHHjRtEQTGZmJg1Fhgvd9uQ69Ho9zj77bCQnJ0saco/HA47jaBiUXA+pGTSZTHA4HHA6nbJeVnJyMtUGJZCi4+fn52PSpEk0t8hxHP2OVNizpKQEBw4ciLhWOdZdUVERJUAIi77D50mtVsPpdEKn08FisUClUtGWRkuWLKFalE6nU3Ixjof4cLRIg8fjifDqysvLYbfbZUN/zc3NmDRpkqgMHvlOeD5PpVIhLS2NbniEz157+480fr1eTzcnPX0+iZLL5MmTQ567jRs34r333sN1110XMqdi2qQ6nS6ko4bZbMZ///tfmn6QCkEKiVrxXITFwpk8z9NceWFhIQ2tk/c7NTWVbuTItZHxfP755/jZz34Gi8WC8847Dy+88IJoC6rBkOMbKrnJIWX4pLyH7OzsCA9s8uTJYFkW6enp4Hke1dXVaGlpCbmJ5HhSHlJJSQk1eOEIBoNYsWJFhE5gLA+JlNFgWTZkpxieoysoKBBV1PB6vTh06BBYlsWOHTtkWYbDhg3D8OHDsW3bNrS2tuLDDz+U9DB0Oh2ef/55HDp0CK2traiurobP5+sxwYRci3BHLUQgEKBGX2hUxOapsrISPM9j+fLlIQvQlClTkJ+fT1mTUtfE83xEjq+nL3Y0JufChQtRX19PBbpZlgXP8ygoKJD15vx+P/WG5cL3ZB4WLFiA22+/HYsWLZI0luFdIkpLS5GXl4dNmzbhuuuuw2OPPYY33ngDbW1t+OGHHwDEJtfn8/lgNptxwQUXgOM4SpoR84Srqqrg8XhQXV1Na/xI81+z2SwbDu2PZrBy4Uwxz7mqqgoulwsFBQUR4yktLQXP81i1ahXNwQL9k6Psbwyl3OSQMXxyN0XYrYDjOHAcB7vdToWleZ5Heno6ZZtptVrMmjULGo1GdhHXarWYP38+VqxYgcsuuwwbN27EtddeixtvvBE333wzdDpdjx+SaHWENpsNRUVFIQQXlmVpZwOp8zU1NSEzM1P2+Gq1mnZ6V6lUSE9PR3JysqwXUl1dTedWKGkVqzLH999/jx07duDnP/85pkyZInkug8FAF0DCGpRrZyOm/CJsaFtfXy+6aAq9CLk8jRzkmJxqtRparZYSk8icE69IypgLhZTFWJNlZWXQarVwOBwRY+5Jrps8Z1qtlrYTcjgcVLzBZDJFbL6iLYb19fWwWCz0eqXGodfrUVhYSFt21dfX45VXXhEtE+ipdx5rODSW6wk/R3t7e49ZnT0dVyKgJ2pFiY4hY/jkborNZsPIkSPBcVyEBya2ezMajWAYhjahlDuuy+VCVlYWHA4HbDYbUlNTqdKD1CJus9mQm5uLOXPmYOXKlbjqqqvw29/+Fo888gjuu+8+ycacHo8HOp0OeXl5EbTzoqIi2Tqke+65R3ZxraqqwvLly+FyuUKazwYCAclFnKjShO9wSR8zqUXgP//5D2644QZMmzYNHMchPz8fP/nJT/Dkk09KejPEuPt8PrS1tWHZsmU9FiY2Go0IBAIhBk6j0YgauGAwCKfTGVK2EAuikZBqa2upGDPxbhiGodcS7gloNBrodDradYN8R8iaJOG2tWvXhnSDF0KKECPmHaWnp4d4gRUVFfB6vTjnnHOQlpYWMu9y74fL5cLevXvR1dWFlStX0rkWQ3Z2NsxmMxYvXkyZyySv2dN8bFZWFvbs2RMyn3Lh0CNHjuC///0vHnvsMaxfvz5EjjB8jFL6qVLhc7Hv9hYDbSQVjy8BEe2mkF1ltJtGduWzZ88Gz/O9CttFG4/BYMCzzz6Lf/zjH7jzzjtx/fXXY926dVi0aBHKyspkE9+kDkosHCjXfT0lJYVqKArp8Onp6dBoNFQ9ROzcdXV1ov3xpEgaGo0Gl1xyiShrUqfTYfLkyaivr8fTTz9NO5wT+TaWZUVDYTzPo7OzkxJdSL5F6vxi94UsQEIjSLyb3goSCxEMBvH3v/89hDkrDGVZLBY0NzeHHDsrKws7d+6MeFYIYcdgMKCjo0OUkBMM/ij4rNPp0NHRITs2uS4R4dct5gVyHIc///nP2LNnT0guN5YIAsuyGD58OFXekXonyHkJE9jpdGLbtm3w+/0Rz1JhYaGsrijLsigqKkJRURGmTZsW4UHPmDEDq1atoiU+NpsNLMtS/dxT5fHFgkQilCg5vgSE3E3piZtOlBrkcny9zdXF8ju/3y/ZuUGr1cLpdIr+zmg0hnhrxMsQdrFmWTbkBaqurqa5OaluAV1dXTRMaLPZ6O+kwpkZGRm44oorcMYZZ9BFheM4zJ07F1988YXkXBUVFcHtdovmPlwuF66//nqsWbMGFRUVcDqdomLIJBcabQHyeDyYOnUqzGYzLr300j69zEJ1GaPRKJuHJecn8l3kM7m8EoHYBqQ3rauI1yB2Tq/XK9oayGaz0e4FVqs1ZHNCmrOGj93lcoHneQwbNgzDhg0Dy7LgOE70ux6PB5dccglWr14NjuNgNBqhVqupeAHZ6NjtdvA8j9WrV6O6uloyx/ftt9/i6aefDjHSBO3tP6YpysvLI4xicXEx9Hq96EYvWo4vPT094tqk6vh64rkJ7xPZEAkbMJ9KDMbcpBiGlOGTuymxuunCfwt3WkJWZ6w3u7cPSVdXF23eKezcQHI9cp4k6a5O2sIsXLgQlZWVGD58OPR6PRwOR8iLJ5yXaHO0Z88eJCcnIykpCVqtVlbbcMqUKVixYgV+/vOf49prr8Wtt96KO+64A3/4wx/wxz/+EX/84x9x55134q677sJtt91G28AIx0LGSI5JGIqEEGKxWKjHRoSJ/X4/Fi1aJGsgyPGSk5Phcrlgt9uxatUqWuTf0/srVJchRAY5Q0rulXDxEj4rVqtVkt2am5tLVXj6uugQMXC1Wk09frVaHcFAFnrIxcXFqKuro2N1uVzQaDQheXQyrra2NnAchw8//BDffPMNHn74YUybNi1CfECn09FWXBqNBhaLJeL+FRcXY9asWWhubqbvQLQGttHSFMnJyZIREhJqJbwAnueRlpYGjuNgsVhoKoNce25uLrxeL7RareR4euO5kXdSTLSAREEGAgMddu0rhpThI5C6KVJFy2K1TWLHC6/j6+t45CDXuSGaJyk83yeffILk5GSkpKTg0UcfjTBs4YK/Useuq6ujjXXJoiVUrSDf9fl8yMzMxPjx47FmzRqsWbMGl112GS677DLU1dVh9erVqKmpwcKFC7FixQpceumlmD9/PvUepXJRRLx45cqVEcLQRUVFWLx4cQijMS8vjxbuh7eXqaysRE1NDebMmYORI0eCYRgMGzYMarUazz//fMzECDGGYnFxMc0fGo1G0fY2Pp9PUoy5qakJDodDVFOShErDhZj7ghdffBFJSUm4+OKLRb0nEu0QGm2tVovx48fTZ8FiscBqtYrOk1h+q7OzE3/6059ouLKkpAS1tbVRGcFqtRpmsxlTp06loen2dukGttEIb2azWfS+kjF3dHRg27Zt6OjowM6dO2kpQ21tLVQqFTIyMqhsILm3RKVmxYoVMa0/VVVVWLlyJV599VU8+eST+NOf/kQjG0Tf1mazSb4X8ZAaPB0xJA2fFMQ8MGH9USK67WKLSSyeJPkdaT47YcIEdHd3R7xA4eSVcPUJtVqN1NTUkLZAwpeP9DwTLs4vvvgi1Go1zjvvPNx///00vCm24120aBHt7t3e3k6/IyZeXFRUJOtlEq+chGX9fj/eeecd6gWLzVV3dzf+8Y9/QK/Xg2EYJCUlYdKkSdizZ4/sDj0WpmQw+KPCSW1tLfXSSEcHqZ16MBiMyKEJQ6XxkpcjuPzyy6HT6WSNjbBhMfAj+aWoqIiWfciVhxCCipgeKs/zSE5Oxp49e+jv5bw0wjKW8khjJfUUFRXBbDZHeLednZ1obGykTZqF3lVhYSHMZjN4nodOp6PEt/BzejweyhTfunUr/va3v+H3v/89rrjiCkniFmH6ms1mZGRkYPTo0Zg8eTJKSkpw1llnQaPR9Eg/WEF0nFaGjyDcmAxWt13OKBICiEqlgtlspp3JxYymmAhxcXExPB4PLXxfvHixbHlGuIKF1WrFyJEjkZKSQher/Pz8iBo54kUK9SyDwSB4no8QL5Zj8mVnZ2P27NlYuXJlVLktMTz77LPIyMigdWMMw4DjuIiyAWHOWIp9G86KJILQTU1N4HkeX331lfgN/f/vD8/zkjneeO/yx48fj/Hjx0saG7fbHeKxEW8pPI8q5v0SgevwjcP69eupwATJHRJGcrRyHqGnJ1SesVgsWL58OZ5++ml8/fXX9Hdiz7tOp4vI7ZWVlUGn0yE9PV30vpPcMalZlGtwGwgEoNFoMGnSJMybNw+XXXYZ1q1bJ6lLm5WVhUceeQTPPfccnnnmGTzxxBP461//igcffBB33303fD6fbKunwVRGkCg4LQ3fUIZUSESo2Qj8nyE4fPgwmpqaMH369JC8xeTJk/HAAw/QMGY0clBjYyNdiNrbfyxuz8/Pp+UX48aNk/TWVCoVGIahrFISOgtHMBiUZbx+8MEHlMHbU2zevBkbNmzAmjVr4PP58NRTT8mea/r06ZLyYhzHQaPRROSDH3nkEcycOVN2HI899hhlPxoMP3YCJ15svPM6P/zwA5KTk7Fu3TpZY0NKA4ixWbp0aUTpDDEwpJaTbKjCNzqk3lBoJMOjDlICDlLNcHmeD8mLchyHCRMmYM2aNXjsscfwv//9j3rfpOOFlKfe0dEhWx86depUqhcrl0e2Wq0RHTx6WwrQ2dkZtbmvgp5BMXxDCNGo5X6/Hz/88APeeecd/P73v8cFF1yAlJQUZGdno7a2Fvfeey9eeukl0Rfp888/l2RLklAMCUkKFwDS9HTUqFGwWCwhvyWLpVA4ubq6GhdccIHk4sTzfETHcuH51q9fH2HkY8HChQvx0EMP4eTJk1RxJzs7OyTPQ0Dym3PmzJHMR/r9/oiaunPPPRePPvqo6PnJRiQjIwNJSUlQq9UoLi5GY2MjOjo6+oXJ98orr0Cv1+Pxxx9HfX29qJdjMplCwsckPyl1f7RaLd58803J+jupshshY5V0Uyc5Wr1ej5SUFNEUhJh3XV9fjzfeeAO33HILZs+eDZPJBI/Hg0mTJmHq1Kmy+r3btm2TlQUkbapIZICES4UbHGKQY8nxxcrKlfKoB1sZQaJAMXxDCHJemcfjwZQpU6guaW1tLR555BFRUepwvP7668jKysKUKVNEa+zGjx8PjUaD9PT0kFY29fX1IV4ky7JYtWoVXcDEXmafz4eioiKsWbNGUp+R6HIKz0eO2dHRAaPRiM8//zzmeSP1dG+88QYA4MSJEzj//PPBcRwtx1Cr1fB6vfj444+hVqsxdepUPPLII9iwYUNIGK2mpibCIwsGg3j99dfB83zEYigWmlapVNi4cWO/08ZvvPFG8DxPNyckHCn8b319fUgud86cObDb7SgpKQkxlIFAAF6vF263GxdeeCHy8/NFvfbW1lbRz7u6umgJg8PhgMFgwOrVq9Hc3CxpRIX5VOFn4V7QyZMn8cYbb0Qwh8WOtXr1atl6PLPZTJWJyDOdl5cXEoItKytDfX296DX29p4OlTKCRIFi+BIUvck7fv3117J5OPLSyx1TeN6TJ0/ipptugsViwRNPPCH78nV2duLSSy+F3W7Hyy+/LCqrVVVVhaysLNTV1UUlRBBWo06nQ0ZGBj3Xt99+S3Nw06dPF72W+vp6rF27Nup8hRN5hISJwsJCUd1Ti8WC+vp6/OMf/0B5eTlsNhuuvPJKvPzyy7JGLTMzEyzLRixWYl5AaWmpKEM33qisrERycjK++eYbyTISYf5Wp9MhJSUFf//739HZ2Uk7RYQLuLMsS6XvYvX4iLdIOs3HyqAV83hsNhvWrFmD9957j34m3BRKsStJrpswTMM3XWKKSZWVlfB6vTAajRHsayn05Z4OVj5CLDiV16YYvlOIWG4sKWPQ6XSifdMIvvnmG7z00kv4xS9+gbKyMqSmptLFRi73IJUMF2Mwjho1Cl6vF59++mnM17FlyxbKApXaVRMquBRRxe12o7m5GQ0NDVSflIRSn332WVgsFgwbNgwXXnih6O8PHDgAk8kU1ZuVCj3J7frVanXIcd977z1cdtllMBqNmD9/Pnbs2IHu7m7Z4wuN2kBJQHV1dYFlWUyaNEk2UmC1WuF0OmEwGGAymfDPf/4TwI/9+zIzM+H3+yPyeD6fD6mpqbjwwgtFr7+wsFBUa1TMSxKOl2y6hA16xZRn9Ho91qxZg4yMDIwdOxa/+MUv8OKLL9JNofBYHo+H1jGSfo5idbS1tbWyOd/eljopGBhlGsXwnQKI3VixkFj4C0c0H6dNm4Z58+bh5z//OSoqKpCSkoLhw4cjOTkZarUaEydOxFlnnUWp1ikpKVCr1SGeUm9o39EWIyns2LFDtnWO0+nEk08+KZuPPPPMM0X72Y0fPx52ux0jRozAwoULJcdwxRVXRJAhhIhW4yWlSmOz2bBt27aIzw8fPozbb78dY8aMQU5ODn7zm99ENWoDKfq7a9cu2O122icyWiiRGLTRo0ejqqoKI0aMQHJyckxee3iEQNjRPVYviSAY/LGLh9frRV5enuzG4ocffsD/+3//D6NHj0ZSUhJ0Ol1IDSgJz+bn59OQPFl0jx07Ruto09LSQpin4VCYlX1DX3KfvYVi+E4BpJiWpPj1wIEDeOedd3D++edL0tjVajVGjBgBh8OBhQsX4s9//jMCgQC6u7tx9913w2Qyhey8A4EA8vLyUFtbG3IssYcp3p5HtIWULKZSxjZc91Csh5tGo5E1fF9++SVMJhM+++wz0b/LGR3Sv07KKMvpYnZ3d+O///0vpk+fHlW0eKA8vmAwiPXr18PtduPmm28GgBByC5HFIsX4wnGxLIvk5GRa8B9Lb0CpCEFfQlvt7e2YNGkSJk2aFGFYd+/ejfr6elitVhQWFuKPf/wjDh48iB07dqCgoAAajYaWUUyePFl20SW6oa+88gr1+MTCwYqn1zsM1DugGL5+RrgkmPCFMRgMtKmlVE6EPABOpxNvvvlmxPF//etfw+l0Su68CWNSblfdH55HfX29qN5hdna2qEyXcOG65ZZbQoyGFL3dbrfLegk//elPsWrVKtG/yb1wpGOEWI5v3LhxMV1/tBwmeaFnzpwZk75jPEC6OhDvhjSBra+vx/r166kEGOmioNFoUFhYGDLHHo8HycnJ0Gq1OPfcc2PSRY33NQijJxqNBrm5uXjxxRdx4403YuLEiXA6nbj66qvx/vvvix7j8OHDePTRR7FixQpZIgvp5LFo0SLYbDb4fD76rpHCdpfLhZKSEpw8eXJI59/6CwMV9VAMXz+DSIKJxbA9Hg8aGxvpLjotLU30GGRBEr5Q3d3d+MlPfoKxY8fi5ZdfjkrRjkZoifeu63e/+x3UanWIJiPHcbjssssijJWYoAAxGtEKmuVKFw4ePAiTyYTm5uYQbVKhSHN4fqqyspKGurxebwirc/z48UhJScHevXtjmgMxgo/P50N5eTm++OIL7Nq1CykpKaiurj4lbD0pwpHb7abF2WIsWzLH5HnIzMzEvn376DFPJc1eahPEsiwuueQSbN26FS0tLTE9s9Fym4SxabFYsHDhQtHavbKyMjidTowdOzYhuicMNige3xAFUSGRarUTDAZpM1kpNYjwBf7EiRNYuXIl8vPz8cknn+AnP/lJn1uixCvOHgwG8cYbb0CtVmPYsGFYtmwZmpqa8MYbb+Cyyy5DZmYmXn/99ZjGU1VVJdkxAgAlJhw8eFD0711dXfD5fLTxK8/zVF/UYDBg/fr1tIOE0+kUNTrhdXx//etf4XA4YioD6erqomLexKgtX74cK1asoOzHX//61+ju7u53b0FqgWlra4NWq6V1elLhaeL9EPHmcPmxU2G4o22CpkyZ0isBaKlICcuyUKlUYFkWXq9X0oMXdpgnnyk1drFDyfENQQSDQUnVBcIiI0n+6upq0d1sQUEBLYb+/vvvceGFF6KiogK33norrFYrLr74YqxatSrCe+nJw9PXBUyYh7NYLEhOTgbLshGM0CeeeAKpqan47W9/S9mPUsfz+/1QqVSyRt1qtUqGPIXeiJhnUlJSAr1ej1tvvbVHRuemm27C+PHjaS9BKezYsQMMw8Dr9YYc/4cffsDUqVMxffp0uN1u5Obm4r777sN3330X8vueGkO577e2top2Mm9oaKBdBaSiBjabDTzP0w4hUnni/g7zvfDCC5KkI4fDgby8vB4//1KLLpEmI5qechswucbHStgzOgaiRlExfP2MaEXlRUVFdBcdzoLT6/UoKCiIKDEoKChAdnY2qqqqqCxSV1cXJk+eLMrm7Al6u4CRAmhh/iMrK0uUFfrRRx8hLy8Pc+fOpcxWsfOSxXr58uWSpB9SHrF69eqI64il3ZJarY6gxEdDd3c3DQseP35c8nuFhYVISkrCQw89FPL5ZZddhlmzZuHkyZM4efIktm7dipkzZ8JsNlMd0J7Qu+Xo4EePHsW9996LgoKCiByyMKS8evVqSbo+aeHE8zzq6+tPeQjv4MGDmDNnjmweXK1W96qZcGdnJ5YtWwaVSiWqwEJSEHLPULTGxwpig1LHN4QQbdGtra2N2EWTB0Cs8zQhdWzbti3EYzp+/Di0Wi0YhsENN9xwSneasYRzw/H9999j7dq1cLlcWLZsmWwHhLa2NuTk5Ij2J2xoaIDNZoNWq40wmmTDIbf5sNvtvVqcurq6MGvWLFxyySWinmtbWxtt/is0jvfeey/OOOMMHD58OOI3H330EfWaw4k1ct6LVBE86ThB/qdWq0M2EI2NjVT8uKurC5mZmZKEJGLoxZ7Tvj5rUsc5fPgwNm7ciOHDh4NhGJx11llYv359xLVWVlZGyOERWK1W5Obm4tprr8WOHTvw/fff0+slmwW32w21Wi1aYkR6PpKogVjhu5xgteLxJSYUw3cKIJWQT01NlfTMepL07erqwty5c6FWq5GWlnbKk+tNTU2S3oJGowkR6w2HVAkHWWDJ3O3evZs2LRWyVNva2mgPNKEBi9Xj0+l0vV6cjh49iilTpuAXv/gF/YyEaFmWRUpKClQqFb0Xr732GlJTUyXZhuHjDh+r2EIabWPlcDiwdOlS3H///WhqaoLf76esTUI+Ir+97LLLwLIs9dqFcyxsiUQiE30lcwiNj8vlglarxerVq3H48GH8+te/BsuyGDZsGIxGI3bs2BHyG2FYzO/3y+bf/vGPf+DKK69Ebm4uWJZFeXk5ysvLI1ID4YadzDnZgAqbwZIuDg0NDaivrz/lOSoFfYNi+E4BpF7WcBFjId566y3JNiZWqxU5OTlYsmQJrrrqKpSWlmLy5MmiLVtOBZqamiRr1sJV6oWIZZEXzp3BYEBRURFtV0TydHl5edDr9RH1dcINh5TklVwtYCw4cOAAsrKycP/991OpM5fLFRLyzc7OxqpVq5Ceno6nnnpK9njRuoYLfx8IBHDFFVfEVEsnhMfjAcMw2L59O13UiSSZx+MR9cA8Hg9aW1vhcDhQU1MTl4W+oaEBFRUVyMvLC2HP6nQ6Wre6atUq/PDDDxG/DR+j2OayqKgITqcTubm5uPjii6HX6+HxeGjHdykiD3m2yDWFv79ETJ14h4qO5uCDYvhOIaKFhrq7u/Hqq6+GsP6kdvKVlZWorKykQsqkYwDJwZyqUAtR0pAaq5hKPUG0/Ge4B7d8+XKkpqaGLJKpqalwOp1UwkqYgwqXpiIMS9Lqx2g0Uim0vixSLS0tsFgsmDt3bkTLHeJJaLVaXH311VGPFc2DYxgGWq0WLpcLaWlpqK6u7nGozWQygWF+fPWFTWFJtEAuLB9NraUnRByDwYC8vDzRekm9Xi8bKQiHmPHJycnB8OHD4XK5UFhYSDdMzz33nORmgUi0yTV3lrpGpY5v8EAxfAmAQ4cO4bbbbkNOTg5GjRqFm266CQcOHJBknFVUVIDjOLhcrohFtqysLKouZzwQTqhgWTYiZBmu/BGOaIv8tddei6NHj4Z8t62tDatXr6aLv16vB8dxePfddyMK5IXnaW1txYcffgiVSoWpU6fG3Tvetm0bRo4cKclA1Wg02L17d0zHklK0SU1NpYLQDMNAp9PhzjvvFG0pJCU3FwwGMXLkSLAsG/I56UNXV1cny3LkOE5WY1UqjBtuFFpbW+FyuXqtkCMF4XlOnDiBtLQ0aLVa6nlrtVrqUUoZ76amJsV4DXEohm+A0N3djRdffBFLliyBXq/HkiVL8OKLL4YQJaRCKJ2dndDr9bI7b2E39P5A+OLc1tYGt9sNjuOQmZkJjUYDlUpFm5jGehyS/1SpVBg3bhzsdjvuvvtuNDc3U6MmJv9WU1NDQ1VSHSiCwWC/KY20trbCZrP1KuQbjm+//RYXXnghNTIajQYlJSXYvXs3fT5ee+01uFwuMAxDm7+KtRQihkCYl0tLS4NGo4nwaEhBf21trWSfOZVKJduoNTMzE+vWrcO2bdtw7NgxUcbp+vXr8bOf/QwjR47ssSZqTxAMBqHT6SKUcUgnBqW/3ekLxfCdYhw4cAA33XQTRo0ahZycHNx22204dOiQ7G/EdstOp1MyTJiVlYWioqJ+e4nlPDWDwYDnn38ePM8jNzcXWVlZsiEgMeO+bt06pKenIykpCdOnT0dZWRlGjx4NjuNo7WP4eYlCvtPpREZGhqin29zcjIyMDNFrisU7lrsOUhogV7NJlGOampoivIqjR4/iH//4B5YuXUq7l99www147rnnZA3yO++8I6oh2dbWBp7nqcHheR5lZWUIBAJobW2lBenCZyQ87Cl2rSTHJ2Y0Nm7ciL179+LGG29EUVEReJ7HmWeeiYqKiojcG+k72J8eX1NTk2y3db/fT1tSKXm50wuK4YsT5BbFkydPYtu2bZg3bx4MBgNWrFiBV199VbaAO9q55Dw+0m29v17iWHJzo0ePxj333AOVSgWO46Ky/8Ln76uvvkJmZiaGDRuGSZMmYc6cOTS3F36cYDBIc5yErUiuPxgM4v3338cTTzxBvaeeenyxtk3x+/3IyMgI6QBAPNI1a9ZENObV6XQoKyvDOeecA51Oh+nTp+OPf/wj/ve//4nOSU/uRUNDA/V0SLkJqXsk10D+HX58EvaUmifiPRJPU8pofPrpp5L9IXmeR0dHB3ieF83xeb1eyWsWg9g9yszMjEko3Ol09iifqGDwQzF8fYTcorh//35cf/31yMzMRG5uLu68807R+q3eoKGhgRaxh4cJifBwvA0fWYjlFkaVSkUXdrvdDo1Gg7q6Ohoq60k46cCBAzRsKlbyUF9fH9K1gRBV2traUFpaioKCAuh0OthsNtrCSSyH5fP5UFJSgj/+8Y/4y1/+gmeffRavv/46Wltb8cUXX0SlqwtVa9LS0sCyLC0tUavVsFgs0Ol0EQt8VVUVXC4XdDod1q5dS+9XZ2cnampqYioXEPO+hYXpwI/G0WKxiF6DxWIR9XSjyUi1t//Yvun555/vFXkpOzsbN9xwA0aMGAGDwRDR6f7YsWOSz4UYxMZbWFgYVcrvVJHAFCQWFMPXR4i9cNOmTcMZZ5wBo9GI1atX46233or7eQl1XqvVhqjpE0o6EcaO17mExp103g73bCorK5GWliZbl9fThaa1tVVy8eJ5XlKmra6uLqJYm0ibhYdWU1JS4PP5UFhYiJycHLhcLmqsRo4cKasWMnLkSGi12ohrLiwshN1uRyAQiDBEwmMYDAbo9XpUVFRQQ87zfEwdG6SEDoSF6cCPHlxPQ4rCeXI4HFTEmmxgqqqqMHXqVJhMJlx//fWUhCSE3AaJbAxUKhUuvfRS7N+/P0QTtSeIVqcpJszd0NCQ0Hk9hSHav1AMXx8g98LxPI8DBw7067kJo7GgoIDmgsiLTfI78XhxxIx7RUUFJk+eTFuz6HQ65Ofnx0QecTgcePzxx3Hy5Mmo5yY1guELQTAYlN3Nhy+44caHHC8QCECr1WLp0qVYv349rrnmGmzevBk/+9nPsHjxYowePVoyXJaRkYHi4mLZPBI5j5hOJvB/jXlJlw45MeTw2kahR+j1ekP+Lcw19rbOksxTc3Mz/H6/aJ3aRx99hEWLFsFut+Oee+5BV1dXyPgsFkvEpqCkpIQ2CY6H8YkWeq+pqaF6uBzHgef5COJOomAgupGfjlAMnwA92WV9+umnuOGGG/pElugLWltbkZGRIUuflyJ59ATRNAqnTp2Kc845B3V1dWhoaJAsuifz0d7+o5K+3W5HWloaVq1ahWeffZZKSYVjx44dtKURYS42NDRgx44dsrR6sW7ZwrwXuQafz0c7EyxZsgQNDQ0YO3YsLBYLVqxYgUcffVS2yH7ZsmUh3pXYNUcjAxkMBuqVRRNDXr58uWjotaqqCsuXL6ekGeFmRUjyIKSPnTt34rnnnoNarY45vyX3fuzatQtlZWUYO3Ys5s2bR8/d1dWFuro6sCxL6yeJ5xg+l73dpMUihCAceyJ7U9FCzAriA8XwIbZd1uHDh/Hvf/8ba9aswejRo5Gamgqfz9fndkC9RTAYlK2pslqtfZLjIuhJo8hoC1AgEEBJSQlKSkpoUXdaWhql15977rl45JFHQhQxJk6cCI1GA4/HQwkZFRUVVF9RypiI/a2trQ0cx8FgMIiG74qLi+Hz+fD666/j8OHDIX37xBajtWvXQqfTRWVyAtK98EiJxvbt22E2m9HR0SFrJKdMmSK72SF5zmPHjtFQpdvtpnPIsiyMRiMVPVCr1SgpKYnQqOyNceju7sbjjz8uSiAKBAJITk6G2WwW/W1fN4pDwWD0VK5OQe+hGD5IvzRLly7FddddR7ukl5eX46abbsLjjz+O8vJy6HQ6pKam9qkdUF9AdBelFkG5Jq2xoqcvo1RdnsVigdFoRFpaGh555BEAPxq2t956C7/73e8wa9Ys6PV6qFQqjBgxAqNHj0ZBQYFkVwaNRgOn0xlhTIqLi6HVamE2myPui8/nQ0pKCjQaDTQaDVpaWiKuSUhzJ5ug+vp61NfXUyPC8zzGjx9PNx5SAt1arZaWXrS1tdGO4SQfy/M81q9fj1WrVkGj0cBsNoPneXi93oixk1rFxsZGWaZiY2Mjpk2bhpqaGvzwww8IBoOorq6mGzUpcg/LsjjvvPPw3XffUe3L3oTa5DZKbrdbkuXZ14V9sMqGCTcYA9WN/HTEkDV8se5Yo4XyCgoKwPM83G439Ho98vLyYDAYkJKSgtWrV1OPYCBeOEJwCTcOPp8POTk5cRtDT3bT0XRJ33zzTaSmpkb06QN+9Bg+/vhj3HvvvZg2bVrUPnwLFy5EfX19SKG3RqPByJEjsWrVKvA8H2JoVq1aRb07YWdx4P+eF7GNTFVVFW1YS7yl0tJS7Ny5EyqVCm63W/QZIEbR6XRCrVZj1apVOHjwIBYsWACNRgO1Wk0FmsO1PQsLC0PEkHmeR2dnZ9TcJsnzqtVqDB8+nPaT0+v1CAQCIc+6MM9pMBgoCzb8eaqsrIx5ExVto7R27doIdquUykxvkMhhTCHEokxyYttarZZKrgGR1zlYrjtRMOQMX7SwZXd3Nz7//HNs2bIFmzdvxrnnniupHmG32zFlypSInbder8ff//73kO8O1INHugGQ/ndEMeWss86K6zl6atzl5uPXv/41SkpKcOLECcnft7a2IisrS/RvRGiYkEeEReEmkwmpqanUUAeDQTQ2NqK8vDxCeV+tVmPt2rWor6+HwWBAVlaWbBiR53lqnFwuFzVewvwZuWZyfJVKhRkzZkCtVuOKK64I2UR8+OGH4DhO1NASma3FixeHyNCRHGm0jt9paWm0HVFKSgrcbjf1KMTekZSUFMydO1e2J9/q1atlhdUJpLz+0tJS+P1+ZGVlhZzb7XZj/fr1ssccapDaTHq9XlG5Oo7jYLPZwHEcJk+eHEFs6q2HfrpiyBk+KQZiQUEBysvL6cJYVVWF+vp6/OlPf6K7YeFCTV52MZUQOeHlgQJZdN966y0MGzYMw4YNw+effx5XYxwv437y5ElUVFTQvoHhx2xubsasWbMkjRDLsvD7/RFC1u3t7VE7tgvP4/F4kJ+fj+zsbLS1taG6uloyjGi322mXbWKcdDod7Ha7qGZqaWkpLBYLtFot7HY7NZAqlYq2t5HTjNRoNEhJSREtI1i/fj3diFit1ojFjoRsH3zwQaxbty7E49NqtbLhTqku7IRI5XA4Ym6MK9worV69GuXl5ZKbhER8p/oL0bxiIYOWiD98/PHHaGhoAMdxyM3NDdGaLSkpwZIlS+hcDrbc5kBgSBm+aPU8//rXvyLqhLq6uuD1eulLTwgURUVFMBqNEccnxcCJrPRQXV0ds2LKQOGTTz6BwWCgXpTBYMCqVauwePFiGI1G/PKXv8S0adNEw7herxdtbW1gWTZkASBkIzFGJxCaJxESbkiYqaKiQrbubPXq1SHGJdyQCdvW6HQ6Uakuo9FIu9WTXbwYrFYr/v73v0uWEQCg+bvw84QvfETsoLi4GBzHyXb9kLr+8LmuqKhAdXW1rLEK39S89dZbcdEyHeyQy+W5XC689957NIpB7kd9fT08Hg/VUBW+0+HEJtI/8XTZSPQGQ8rw9SY5LBWWIbm99vb2kNAQefiWLFmSUIZEiPr6+qihsIFGQ0NDxILt8/nAcRwNZXIcR3un2e126PV6+P1+tLW1oaqqChaLhX6fZVlotVqMGTNGdGEn/eY6Ojoi5sPlclENUDlyDqmz6+zsRGtrKzIzM+F0Ouk5hMo2UszGkSNHQqfToaqqKiLnRhDuAcl52rGEobu6uqhcWkpKiuTGIC0tDYsXLxYlBZG6O6FKjZi3KYdoDYt37doV9RhDAXIbdNJ5w2QyIS8vD06nk0rOTZs2TfKdJsQmki8VltIoub9IDCnD11MGYiwhh6qqKsnQkNfrTTjjNxgo0dHq2kjIpri4GDqdDrt37w5pLUM8tJaWFiQnJ2PEiBGw2+148803qQdPDD+pIxN2HSd9+44dO4a6uroQDVAhg5OwMMkmp62tjXYKJxsgnU6Htra2kOsIZ14KjUV4eFOqBVFPyR6xLHAdHR3Q6XSyXm1paWlIdwZShE6e876UDZAFXIoBGw8W8mCB3Dw2NzfjmmuuwdSpU6FWq7F9+3ZZjVlhJxayaSKbRKUQXhxDyvABPXsxo3mIRLEiWq4pkTAYKNGxjpF4Ak1NTfB4PPD7/eA4Dunp6bR+zmq1guM4bNy4kRqnX//617TOz2QyibIUiawZYRgSEkxFRQUaGhrov0lndyINJpbL83g8ETllYW0feSbr6uowdepUUaNI8nXxaIwrh4aGBrjdbtFu9IWFhXQsZrMZubm52LNnD92kxGNTtXbtWrhcrogWSqtWrQLHcRH1hIMNsXpYsRLGGhoaUFBQIBki9ng8EZ1YrFYrHA6HaAcNv9+veIAYgoavJwzEWF5kubqpjIwMcByHI0eOSI7nVIcaYvGmBho9WUDT0tLw73//GzqdjjIqxWSwfD4fNBoNGIYBx3EYN24c5s+fLxtaE+bnyAKs1+vB8zxaWlpQXFxMpcCEBfNi7E2dTkdFqU0mE2VeknCm8L9iBe+BQAAsy/ZKq7In6Orqwvr168FxXIgXXFhYSIWhw2XKyHzLqcqIiRmIPffNzc1UKcfpdIZ4Iunp6aiuru6/i+9HCNMhpF4xlg4p0dYHwtqW01oVnod8Fv59EvlQq9WKB4ghaPgIYjU4Yh6iUMQ2fPdOQDw+g8GA1NRU1NfX44MPPqB/H0jNvYaGBpSXl4dcU3l5OfR6PW6++WYcP36838cQDVJKJuElB8RAabXaqHkxnU6HJ598Ev/85z9x55134tJLL6XkkfDnweVywel0Sub0WJYFx3E0r0Ik4sTuqdvtRl5eHm0VVV1djU8++QQ1NTXgeR42my3EyxXrmF5aWhq3WrZYEAwGsXPnTjz55JOyxlZoBA0Gg+xG4t///nfU555sesRY1ES+LRE2Zz1FQ8OPTXxJqyeSd45XOmTdunURovBi0YbKykq6oQkf32BXtoknhqzhixXhHqKYiO2qVaskFUTUajVlE44YMQJpaWlYvXo1Vq5cOWAP2vfffw+z2Uz7pWm1Wni9XrS0tOCss87CuHHjsHPnzpDfxMMz7ckxjhw5Ap1OR9m0hIFK8mVCQ7h9+3a6gxUakPDzpaWlITk5OUQGjWhDivWhI10mpHbSwq4MHR0dYFlWND/F8zw1yOFGhHQoEBps4TNHagf7s39ivBAMBiWb0F5wwQVwu90YM2ZMVCWjmpoa0e4Tfr9/UPbGI8ZcigsQj3SI8JnJyMgAy7JYtWoVrUElQgfZ2dlYs2ZNhCB7ouf9TzVOe8NHIFxEwxfU5uZmmM1m0fApyTFptVpkZmbS3m9SrWxOxYO2ZcsW5Ofng+d57Nq1Cx9//DEsFgveeustdHd347HHHoPdbkdtbS0OHDjQZ880Fu+WCAds374dv/nNb+hmYc+ePdi2bRv2799Pj0FyXcQYCLsLBIM/dlkQM2Y6nQ5TpkyBXq/HpEmToNFokJeXJ9rotKCgAIsXL4bdbhe9JrvdHiJALtfNm+M42v1dbNEW5viEC2MgEEBRUVHC5YnlIJdKOHjwIDiOi/rcd3Z2hnSH1+v1tNNHT5mipxJSG7utW7ciPT39lKQYwsPQDocDKpUKJpMpJD8s9PBaW1sl6zMTJe9/qqEYvhggF54hL/Q333yDd999F1u2bMEvfvGLmGrJ+gvnn38+fvvb32LEiBE0//jAAw9g8uTJlIr/+eefY926dTAajVHrwKJBLIwybdo0lJSU4Pzzz8eYMWOg0WjAsiysVisyMjIwfPhwGI3GCE3MtWvXIjk5mRIgGhoa0NnZGRJu9nq9osbMbDbjD3/4A1paWvDMM88gMzNTljY+cuTICGNG8rosy2LkyJG0dk2qtU9XVxeMRiO9NrFFmyxEHMfBZDLRWsPBoicpBjEj0BNiFQl1L1myJOb+gwMFsY3dhg0b8Oijj6KsrIze01NtXITlM+H3QrhBkRN0Vzw+BbLoSYxcLrSg1Wrxxhtv9MsYg8EgduzYQSnrwoX4u+++C1GCMBgMqK2tldyha7VazJs3DxdeeCHmzZuHuXPn4vzzz8fZZ5+N0tJS5OfnY9y4cbIvFVngi4qKsHTpUpSUlECn01HtSpfLFRLaLCsrQ3Z2tugck3CzXI5PrVYjNTUVLpcLZ5xxBlQqleRilJaWRuvafD4fJbno9XqqtJKSkkJDpfv37xetDwxveEuMvli+buXKlRg5ciRee+21Icms60lIjZSdFBYWRu0/ONCQygMTIpTT6UROTk5MvShPNYhxDG9WnGibi1MNxfDFiJ7qVUqRZioqKmC1WnHeeefhlVdeievYSKxfo9FEyFyNHz8+gvCSm5sryVhNT09HUVERCgoKaBumESNGgGVZZGRkYMKECSgvL8dZZ50VU09CORIRQXt7Oy0fEH5mMBgwdepUWCwWaDQaSbUT4flOnjyJrKwsycWIKJF0dnZiyZIlol3USXjS5/NBq9XCaDSGfIcwMaWM8Pz58/HII4/QcDLLsrReMBZPbzAWH8e6QSRGsidM0YGAnDHneR779u0D8H+GPFGFIwZr94r+gmL4eoh41Ol89913uPPOO5GVlYWSkhJs3boV3d3dvR5TtMUmEAiIemVSn5O81c0334zHH38cO3fuxEcffYRvv/1WdD5iaQIaqydgs9mg0+lCXkqr1Yrq6mp0dXVRski0Y91yyy2YOnUqbDabaOlDdXU1qqurodfro4aCAoEA9Ho97b9HWHOkuaoYHA4HVqxYgZkzZ0Kn00WEZuUWxMHchTvWBZaERROdeNGT/BgpPSD3LRGNy2DcTPUHFMPXz4gmN/WXv/wFEyZMwMSJE/HXv/61xy9JLAtHa2urpGdnsVhEWXg9KXQVM7yVlZU499xzMXfuXMqyFIOUfiYxDGILoNj5SktLUVVVhdbWVuzcuRMmkwkTJkzA0qVLsXbt2hAjYjAYwLIs0tPTaflBtMVNOE8kD1hSUiJZ6qLVajFu3DioVCrJEgCe5/HFF1/ENJ/hhjLRF7Bo4xM+t7FEAwYKgUBAVjVF7PoS/d4oUAxfQqC7uxvPPPMMiouL4Xa7cdddd9FiYgKxl6mrqwsrVqyQbUxK+q3JdSsXCiGTNifh3gYhxYS/zCdPnqTdFDiOo/kxlmUxbNgwuFwuLFq0iOqeip2fFIILF3fyt6KiIni9Xtx///20M7qYVzFu3DhanEv63F177bXUkybzt3r16ggPUK5W02g0Ys+ePZLzx/O8qLalSqWibYGkjL7VaqVzNGvWLFx33XW4++67JefKaDSis7MzJm9wMCy+xOAJRQQILT8RPKUPP/wQbrcb5eXlSn5siEExfAmGxsZGzJw5EzabDTfddBMOHTokudA1NDTQ/oByYbqqqirRPl/Cl1cuCV5WVhbSi2727Nmor69HeXk5OI6jclvJycnIyMhAbW0ttm7dGtVLKysrg06ng81mEw0LWa1WFBcX4/rrr8dFF12E3Nxc6qlNnz4dtbW1WLduHebOnRtRW0akx4SIlq8Rq0+rq6uD1+uVDGlarVYqqG2z2aDX66HVarF48WI8/PDDePXVVyUNGREkFv4vKSlJtoOB2+2OYOEKvSOhLqjT6ex3CbS+IHwDYzAYUF1dnRCyZW+//TbsdjvuuusuJT82BKEYvgTF3r17sWTJEuj1+giVj4qKCni9Xmg0GjQ2NooaK9JRgLykx44dE315hZ5cT8SjSbcEl8uFSy+9FH/729/w5ZdfSl6PlFBAZmYmkpOTJfsepqSk4Oyzz8a//vUvdHV14eTJkwgEApg/f36Ihyll+I8cOYL9+/fjn//8J/Ly8mQ1D5csWQK1Wk1lvIxGI1QqFRwOh6RHSIrWSfeH/Px8rFu3LuR7cqHL7u5utLa24o477sCcOXNgNpsla0DVajWSk5MlyToPPfQQ1qxZE9LJneQwV6xYIVqjmghItDG99NJLSE1NxeOPPx7yeaKNU0HvoRi+BAYp1pZaBFUqFV3cvF4vrdkhBa1btmyJeEnJyxseMuN5HqWlpTExNImH1Nzc3KtrCu8hJxSLJscnXce/++47PPTQQ/B6vcjIyMB1112H1atXx1ScS7wxs9kMq9UqS50nzFFS5mAymTB37lxaGhLebDYQCCA/Px9+v59uBnieh8ViiZiXnnoMdXV1kr38pO6P3W6nrFsyTqLPyLIsZcTqdDradaGmpgb79+9XFnMBnnrqKZjNZjz33HMDPRQF/QjF8CUw5AqCwzuCV1ZW0h6CxGMJl8Hq7u7G/v37sWXLFpSXl1PtP0LWKCgoiLnQtbdUczGvUkxAl+d5vPLKK3jvvffw3nvvobm5GU888QQuuOCCkDFG81K9Xi9UKhX1oqRqskjj2HHjxsFsNsNgMITUDQrzUBaLJcQz5HkeK1asQEtLiywTsS+M4Hnz5mHOnDmyudoRI0ZArVYjEAigqakJy5YtQ1lZmah+ZFtbG9UktVgsESHR08W7EV7ngw8+iLS0tH6rs1WQOFAMXwKjJ50WxAqpKysrMXPmTKxfvx4VFRUwmUxITU1FeXk5tFqtaGcCrVYboaIhVm/XW6q5nDF3u93Ytm0bgsEgbDYbPB4Pxo4dizPPPBNjxozBmDFjkJWVFVHHJxdKbG9vR3JyMg1xhhsVjUZDFWSGDx8Ok8kEq9VKDWX4eP1+fwSZRdizsa+Eh2AwiD179uCf//wnrr/+eixYsABlZWVwuVxITk6mBfdi10oYqzzP05ILjUYToYFK9CPJc1RSUoLFixejoqIC9fX1ojllKXJTokPKgIeXjHAcB7PZjHfeeWeARqrgVEIxfL3AqdoNi4UFSVgtnNoupxRz/fXX48knn8TOnTtx5MgRupiLGYyioiJK0hAqvQgXzr4s8LGUX0SjikfzGMNDiUSwWvgb4uWqVCo8/fTTAIATJ07g4MGD2Lp1a0hXh1j60cUqNN3d3Y2DBw/itddewyOPPILrrrsOy5YtQ1FRER1nTk4OZs+ejQ0bNuCOO+7Atm3bEAgEcPz48QidS3KtbW1tYFk2IiQrDBsLx0s2Tg6HAxzHwWq1gud56HS6iE7fQnITYQE3NzfH9fmP9zsVrRYykUsoFPQ/FMPXA4i9TP2xCISfR6vVhnSM4Hk+pOt3tCLb6urqiDGTFjpSBojo/5FcYLwYbVIKF8JWUNEMa0NDQ0ROsKqqCizLorGxMeRetLe3U8FrsVpDnU6Hu+66K+T4wWBQtJFsrCojJ0+eRFtbG1544QXcc889uPLKKzF//nxMnjyZGq0pU6Zg4cKFuOqqq3D//ffjpZdewueff46TJ09GncOGhgZUVFTQayWEJ5ZlY5b/ys7OppqkhFgkV7Om1+vR1NSEQCBAiVPxKKzva7G+lMGUiwIketG8gv6HYvh6AOHLJCQOOByOuKprSO1GFy9ejDvuuAM5OTk9ks4SC41NmjQJ6enpoucXy9/Fa0fu9/sxdepU1NbWwmg0wuPxQKPRgOO4kFZQcvP4xhtvhHTLMBgMsNvtMJlMET3LyGInljdbt24dkpOTsXHjxohznHvuufRY5F6TNkdi88xxHM4++2yMGTMGarUadrsdJSUluOSSS3DDDTfgr3/9K3bt2oVDhw71af4A8RxgdXU1nE6nbBhZSE4yGAwoLCxEXV0d/Y7YBoqci7S4ErZ1CgQCqKysxNq1a9HV1dWrZ6QnGrhi4xJr/BrNsL3xxhtwOp2ix00EmTQF/Q/F8MWI8Jepty9sT89DQIzYvHnzcO+992L16tUhC59YnV5lZaWkCDXP87JF2/He9XZ2dmL58uVUxV7oLRMCSVNTk+x5ycK6dOlSbN68GUeOHMG1114Lk8mEpUuXIjk5GQ6Hg4ZpxYyocHHu7u7G8OHDUVJSEnHu1tZWmEwmOsfC+rzwUKLP50N5eTn+9a9/4d1338UXX3xxykLhwlIUKS+elFp0dHSEkHksFkvE3IQ/ew0NDaL1jRaLhYojjBw5EsnJySGEH4vFgrFjx2Ly5MmorKzEhRdeiNraWmzcuBE33HAD7rjjDtx///2yxfpy8+f3++H1elFbWwu9Xg+HwwGNRoP8/Hzcc889knWXLpeL9rNTPL7TF4rhixHhDVD7K1TSk9YuwoVPzAuoqamRpfqbzWZJD0kKPd3Vk3HxPC8Z3hS7NrFjEGaiWq1GdXU1zjrrLDgcDlgsFsyaNQvjxo3Dq6++GtMYu7q6UF9fTxdrjUYDnudRX1+P77//Hlu2bIFGo4HFYkFSUhJVYBGbZ6PRiNraWnz55ZcDqrHZ0NCA7OzsiHIGQnRRq9VUSo3neVHvVUiSilZOs3r1ato/USiM3tbWRhvs2u12Oo/p6ek0RKrVajFy5EhZVRuTyYTx48fj7LPPRk1NDa677jrcfvvtmDFjBtRqNTIyMkJKNEjZiVqtlvTMNRoNnnjiCYwePTphBaUV9D8UwxcjhMauJ8apL+cRIlajGu4FyB1r586dmDFjBrRaLWw2GzQaDc4880zcfPPNeOONN3D8+HH6m97mYkg+SmoB5TgOLS0tsuQQqRIEjuNw/vnn46233gIAjBkzJubawoaGhghhgKqqKsqeJIun3W6HVquFSqUKKR4n80y81eXLl0Ov10d0wDiViykx5oTVGd5NQkhUcbvdYFk2YuNTVFSEvLw8GI1GOJ1OyYJ/u92OoqIi0WsV8xKrqqoiQsrhfRYJiGEVqtkMHz4cDMOIdtEQ5oerqqqg1WrhdrsjxlBSUoL6+nr8+c9/htVqxfjx4xU1ltMUiuHrAcgCLNcTLh6hkniGUXsieLxv3z48+uijqKurw/jx48FxHCoqKnDttddiwYIFoiGv8DF1d3fj6NGjOHDgAPbu3Que57F48WLZ9kdOp1OyHKAnJR02mw2ff/551DkJP2a4EdPpdBHXWlxcjClTpsgSJhKlr1wwGMTOnTtlx9PU1IS1a9dSD5rkV71eLyorK2ktoJzHJ6a2YzQaZX/jcrkwd+5c3HXXXaitrQ3xUMn3ysrKsHr1atx4440YN24crU/UaDRRBaOJV7d169YIOTSe5/Hqq6/CbDbDYrHQ8PZgLNNQ0Dcohq8HEIa5LBZLv4VK4qkN2Jdjff3119iyZQs2btwoueCQ/nhmsxksyyIpKQkajQapqanIyMiAyWSS9fhI2Kuzs1PUSPSkLYxWq41pAWttbYXb7Y7wYvV6PYxGo2wRv1DQWziX/RkF6A1iGQ+5fp1Oh9zcXAQCgRDSVkZGBniej/CMfT4fTCaT6LHdbrckcSQtLQ3Dhg2joWOVSkVrSYUGiuM4LFmyBBdffDEWLFhAyVjJyckxPQtWqxVPPvlkhESbzWaDSqXC6NGjUVxc3NMpVTCEoBi+XiAYDKK5uVlyEezrsYWhynjtRvtyrGgKMmeccQbsdjtGjhwJk8mEnJwcVFVVYd68edSI1NTURBTGk42CcNES/v+9e/di4cKFManJdHV1Yfjw4TH1NQwGg9DpdNTLFHZft1gsUbtdiM2lnGeq0Wjwi1/8AkeOHIk6tnihJyHzzs5O1NTUwGAw0NzcxIkTodPpkJmZCa1WS+s5jUYj6urqZEkpUn/jOA6TJ08Gy7JgGCZknoVzSmTjSJjZ6XTScHwsNaBqtRpZWVkh4XjyHVKOEd7zcahC8WjFoRi+PiJeD1YiNx+NdRE9efIkDh48iHfeeQfbtm3D5s2bYbPZ0NzcjFWrVoFlWUq04DgO9fX1aGtriyhcf/rppzFjxgxYLBa43W44HI6o3vWhQ4dgMBiiXkdraytaWlqQnJxMvVhhOLgvIUupsPKKFSuwePFimM1mXHfddfj6669jn/w+oKchczI/pOaRhDt5nkcgEAh5zsWUgkh9n5iXGH7effv2STKO1Wo1UlJSQjZKpMaQbFakcnzFxcVYtWpVyHnDw+gk0iDWwWOoIJHXk0SAYvgSBP1VHhEv9HR8pBs1obYLmYKNjY2oqKhAXV1dyKJVUVGBjIwMnHnmmVi5ciVSU1Pxy1/+Eueffz6mTp0q611//PHHyMzMlB0LKXMgAt/p6ek9ovDHWlsmNc7W1lZccsklMJlMuOqqq3Dw4MGIY8Rzh97bMHf478QMWWVlJRVGJ8cm5Smxih6IPVMlJSXQ6XSiHSrq6upQVFSEurq6iBpQ4qkSZikBMaTLli0LaXdECviHavlCoq8nAw3F8CUABoOSRE8X0ViIQIReT6S6WJbFeeedhwULFmDMmDHYvXs37rrrLkyaNAnff/+9rFF4++23MWHCBNFxiynFlJSUUKWX8DBueNE2oe73VU2E4OOPP8all14Ko9GIjRs3oqOjo1936L01puR3coZM7tjRziv1TDU3N4vm8rq6umAymWjZQnJyMmpra9HZ2Ylt27bB7XaLnsdqtcLpdNI5FUYZhmLB+mBYTwYaiuFLAPSEwDHQiGURjbX0w+FwYOrUqSHSWz6fD7m5udizZw927doFs9mMlpaWqOPasWMHSkpKIj73+/2yqjb5+fmyxJumpiZceOGFuPPOO6OOoafYv38/1q1bB6PRiPz8/F55macS/ZUvCj9uNCZvR0cHnnzySVgslpBjxNJLsqqqiurUDlVDkGhEq0SEYvgGGMFgEJdffrlk89HB+GKGF/vr9Xps3749RB1FjhKvVqvhdruh0WhwzjnnxOT1PPXUUzjvvPNCPiMkFqlNhdvtxqRJk0Q7UgiFnR9//HGcffbZPZuEHuDDDz88pSo6gwFioTqhOHu4nmowGERNTY1o54xwgW5SO5poG4t4QfH4okMxfKcIZFdLxJ8PHTqE22+/HQaDAWq1Gnl5eUMmJk9evLa2NtTX14PjOKqQotVqsWrVKpSXl9Mde/iO3+PxoLW1NaY5IL+95557sHTpUvp5d3c37rnnHqSlpcl6Ap2dnVi7di10Ol0E8YYY3CNHjkCn0+Hw4cPxnioAyg5dDGJh0HBx9vr6epSWllLtULfbDa1WC5Zl6cZJLFxss9nA8/yQJnsoOT55KIavn0FeYEKVV6vVSE9Ph1qthslkQkZGBl599dW41u4NJIghmjdvHlwul2iBss/ng9lspiUFwrwW+Tcp59i+fTvVmRRCrJ9abm4ujh49il/96lfQarW0VkyMCUh60gnH3dTUJKkXes455+Cxxx7rtzmTMs48z/ebwR0M6OjowLZt29DR0YGZM2eGeOZtbW1ITU2NyN+Wl5dj2bJlkmUVJFw6lDFU1pP+gmL4+hlk5yW2+JaVleHyyy8P+f5grbsJN0SEuCK1+LAsC47jItoLlZSUoLCwMER6i3iK1dXVlJknJWNGasSSkpJQXFwMvV5Pu5ATJiDpQt6TReBPf/oT5s2b12/3Rux6KisrkZmZifz8fLz22mtxP2ciI/x54nkeer0el1xySUixu5ywgk6ni1pWMdQxWNeT/oZi+PoRZCff3xJnAwXhSyVliIQEBCHS0tIwcuRISS9HylPU6XSYP3++ZA0Yy7LIzc3F7373O1itVrz88st05+t2u8HzfEwNY4UIb3TbHzVRUjv0H374AQ8++CDsdjuWLVuG//3vf3E7ZyJD7HmaNm0alYdrbW1FU1OTbBumN998U7ZBsYLTF4rh60eQ3M1Qy+GIUe/D8y+APIGF53nR1jHBYBAqlUqWaUlq+sQg1J788MMPQ47b253vqcyXSI3zyJEj+OlPf4qUlBRs3rwZ33//fdzPnSiIlZwR7XuBQADTp0+H3+9XvB4FIVAMXz9iqHp8Ut6dMGdGYLPZ4PV6I0K8fr9fdE4aGxthNptld/KNjY2yIslpaWlYsmRJPC414RhygUAAs2bNgsfjwVNPPRWTRNtgQ082inKRBsXDUyAFxfD1M+RyfIMx3yBnCMJbCxHquBRjUmzRqqioAMuyUWXD7HY7cnNzI+aTCCzrdLq4GKVE9da3bduGMWPGYMaMGTHVOQ4m9GSzIRYiJgoyg21DqeDUQTF8/QxhWJCwOh0Ox6DajQrDb01NTcjKyhL9nsfjgdfrpcXBwpY9YoxJqbxWfX29ZLsacmydTger1QqWZWlLnbq6OlRWVtLmq01NTXG59kTy+IQ4fvw4fvOb38BsNmPDhg0hklyDndTQW53RwXq9Ck4tFMN3ihBexzcYXlCicanT6WiNlFarle2WQMSEe0IACV+0whuqhotakwVw165dNL/odrtDNhNWqzUuhg9I/JqoL774ArW1tbBarbjrrrtQX18/6MWJFTq+gv6EYvgUiIJoXJJiYL1ej9raWkybNg0FBQWoqKiQNARutztu3tauXbuwbNkyupCH60RqtdqIzgEkxBqvzYVwEXY4HAnb0mb37t1wuVyS7Z8GIxRPTkF/QDF8CkTh9/tF287k5eUhOTkZZrM5pHN3eN+z/tZzJCDqHeEh0fr6+rien4zh1ltvxTnnnJOQi3Eih2UVKEgkKIZPQQTketKxLAuz2YzW1lbU1tYiLy+PsisHwrvo6upCUVER1ffsz5BYV1cXLr300n6t5+sLogmCv/zyy6d4RAoUJCaGMQoUhKG9vZ0xm82MzWYL+dxmszE2m405cuQIw/M8c+eddzIFBQXM+PHjmVGjRjE5OTnMxIkTmc2bN5+ysY4YMYIxm83MyZMnmWeffZZpa2tjbrnlFmbEiBFxP9emTZuYTz75hPn444+ZQCDA7Nu3j9m7dy+zadOmuJ+rN7Db7cxXX33FdHR0hHze0dHBHDx4kJk+fTozY8YM5oknnmCOHz8e8fujR48yH3zwAXP06NFTNWQFCgYGA215FSQeonUhT09PD6Hxu91ubNu2bcBCaRMnToRer+/XcwyWMKJUV4P8/HwsXLgQer0eHo8HqampaGhoQEtLSwSJqaeerJKHGxxQ7tP/QfH4FESA4zhm+fLlzIUXXki9h46ODmb58uXMokWLmO+++46x2+30887OTsbn8zEcx53ysR49epT53//+x+h0un49j5wXnJKSwrS3t/fr+WPF5s2bmYkTJzI5OTnUC58wYQJz5MgRZsqUKcxLL73EjB49mlGpVMz777/PlJWVMTabjbn33nsZi8XCfP3118zixYuZt99+O6one+LECeaKK65gHA4HM3PmTMbhcDBXXHEFc+LECcnfKF7lqUdv7tOQx0BbXgWJCSGrU1gnV1xcjLq6OgADyxgU1kdarVZoNJoea3D2BIPF4yMI393v378fWVlZuOuuuwAAW7duxYQJE2C1WiNITEQIQKfT4dFHH8WuXbvQ0dGBkydPhpyjJ2Ue/dlhXoE8Er0cZyCgGD4FkiAhMJ7naSslg8EAlUp1Sorw5UIzUlJVPe260BMM9gXkww8/hN1ux6xZs2j/OrmaTLvdjvz8fIwZMwYGgwEjRoyAzWbDxIkTqcJOrBuBwT53gxWDbcN2qpAEAAPtdSpIbBw9epRpb29neJ6nxJYjR44wdru9X8KbJ06cYDZt2sTce++9jNlsZr766iumtraW2bx5MzNixAjm6NGjjMPhYPbt2xcSeuzo6GCys7OZ2tpa5ne/+12/jeu+++5jUlJSmEOHDjE1NTV0XIMBNTU1TGtrK/P3v/+dCQaDzLnnnst8+OGHEd/Lzs5mPv/8c2bUqFHMiBEjmGHDfsyKnDhxgunq6mK+/fZbpqurSzTEa7PZGJPJxGRlZTF2u51JSUlh7rjjDubZZ59lJk2aRJ+Zjo4OJicnh2lraxuQMPnpgA8++ICZOXMmEwgEIv42atQoZsuWLczo0aMHYGQDjIG2vAoUhCOadyBH2/d4POB5vl93soOVJBC++5fzBliWFRUdlzqW8LdqtRpTp07FHXfcgT/84Q+YMmWKZAnIYOxOMpigeHziUAyfgoRCLC9qNNap2+1WFlMRiG0Y+hIyltqgnH322XC73UhKSoLZbEZ5ebnoJoa0mRrq3dAHGkqYORKK4VOQUIi1G0JBQYEkKeN03snKQWxTIWyya7VaodfrYyYJRdPTbGpqkswDGgwGFBYWwmw2K0SXfoaiexoJJcenIKEgl78j+SCGYZiMjAwmOzubaWlpYWw2G/P1118zixYtYvbu3cvk5eUxv//97wfqEhIaV1xxBbN3717mwQcfZGw2G9PR0cEsXLiQmThxIvPwww8zTzzxBFNQUNCjnBvJAYfnfOXySzabjamoqGAeeugh5uDBg8zy5cuZiRMnMj//+c9Fj6Wg75C6T6cjlDo+BQkD8mLOnz9ftIawpqaG4TiOaW9vZ1JTU5mysjLGarUyX375JWM0Gpm//vWvTHt7u1IjJgOxOr/29nZm165dzA8//MBceumlPa7z4jiOGT16dMRiKqck8/333zN33303M2LECMZmszEPPvggc9dddzEZGRlKrVk/Qeo+nZYYaJdTgYLwGi+NRoP/r717aW2iDcM4flUKUkmqGAUTTRceApksEih0obhzobgTpKRCKzZYBC2CFJlPMIuu3HkG+yXcFi3uihYxBQcLKnUWxsOiiEIw7yrzptOJPaWZJs//t0uhzZQMc+V+5n7mTqfTtb6+vtClmfo9vv3799dc110168/0m/Yb1digc/PmzR2b6NDsSTJhfzuVStVevnzZ8mMAggg+RC54cXRdtzY0NFSzLKvpDMOrV6/WEolE6IZoOgU3bqe7/sLuL/X399c+ffq05v0OHDiw6v3qx9BJMyzRGQg+RKrxwttY+Z04caK2d+/e2tDQUGi43bp1i+aWFthoM9F2NVaYzTpJw6q7+vxDnvaCViL4EKnGC2/wgnjjxo3a2bNn14Tb5OTklvefYbUo9nltpgrct29fpGOv0J0IPkSqfuF1XXfDm6vrj9sKk06na+/evdvx4+4mUe3zClaB586dW1MF1p8LW8c9XLQCXZ2IVCwWU6lU0tjYmBKJhL+FYb1pCMvLy6Hdgj9+/NDAwEDbjr8bhHV6tmOuYmOXoeM4sixLx48f18mTJ2VZllzXXbMtZbdNw0BnYh8fIletVnXnzh09fPhQS0tLSiaT/9zPl8lkdPDgQWUyGc3MzPj70erjdJaXl2nZ3oKo93m9f/9eFy5c0PPnz9Xf369sNhv6+VuWpc+fP/MZY8uo+BC53t5e3bt3T9evX1exWNTc3JwkqVgsanh4eNV+vtHRUfX09Ghubk6FQsGvUizL0uvXr6kGtiHqfV6pVErfv39XPB7XkSNH/JWAxs9/ZGREPT09evXqlf97zPjDZlHxYVeoVqu6e/euHjx4oHg8rp8/f2rPnj3KZrNaWlrypyFcunRJs7Oz/kSBxiqlUCjo69evVHwdrPHJMocPH9bk5KSePXumRCKhlZUVjY+P6/z58xodHdW1a9f069cvPX36NHSKB9AMZwd2Bdu29fbtW7mu6y9dXrlyRYODg5qdnfXDTZLS6bQ8z1MymfSrFM/z5HmeSqUSodfBHMeRbdvK5XL+l51SqaSJiQkNDAz4n+38/LzOnDmjo0eP+suh9Sf82Lat6enpiP8T7GZUfIjcRp7P2RhmYc+bvHz5sv7+/asXL17wbb8LrHe/cbPnDNCIe3yI3HodnMF7dsEuRMuyNDg4SOh1kfXuN272nAEaUfEhclv99h51FyKiQ8WH7aDiQ+Qa9/I1m8jQ7Pd42ryZtnrOABLBh10iqk3U6FycM9gqljqxq7B8ic3inMFmEXwAAKOw1AkAMArBBwAwCsEHADAKwQcAMArBBwAwCsEHADAKwQcAMArBBwAwCsEHRCQ4OZxJ4kB7EHxAm1WrVU1NTSmdTuvixYs6duyYTp8+7b9Op9OamppStVqN+lCBrkTwAW1m27YWFhZULpfluq5GRkbU29vrvy6Xy5qfn9fExATVH7ADeFYn0EbBOXLB19VqVbZt69GjR+rr69Pv379VKpXkOA5DdoEWoeID2ig4OTz4ul4NLi4uyvM8lctlLSwsyLbtKA8b6CpUfEAb/avii8fjTBUH2oCKD2ij4OTwWCymYrGo4eFhvXnzZlX1V5dMJpVIJPTly5eIjhroLlR8wDpaPei0fh/vyZMnSiQSqlQqymazWlxc1J8/f/ThwwcqPmAHEXxAE/WAevz4sQ4dOqRKpdLSRpNgoK6srOj27dv6+PGjZmZmlEwm5XmexsbGlM/nNT093YL/CgBLnUATwW0HrW40icViymQyfhUXi8V0//59FQoF5XI5nTp1SrlcTvl8Xo7jtOQ9AVDxAaGCTSh17Vp2bPXyKoD/UfEBIYLbDOra1WgSrAYBtA7BB4RIpVKqVCryPG/Vzz3P07dv35RKpSI6MgDbRfABIYLbDiT5jSbj4+NUYkAHI/iAJhzHUT6fp9EE6DI0twDroNEE6C4EHwDAKCx1AgCMQvABAIxC8AEAjELwAQCMQvABAIxC8AEAjELwAQCMQvABAIxC8AEAjELwAQCMQvABAIxC8AEAjELwAQCMQvABAIxC8AEAjELwAQCMQvABAIxC8AEAjELwAQCMQvABAIxC8AEAjPIfWez74hDROY8AAAAASUVORK5CYII=\n" }, "metadata": {} } ], "execution_count": 2, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:48:59.665Z", "iopub.execute_input": "2022-04-10T00:48:59.677Z", "iopub.status.idle": "2022-04-10T00:49:04.590Z", "shell.execute_reply": "2022-04-10T00:49:04.651Z" } } }, { "cell_type": "markdown", "source": [ "The `nodetype` argument ensures that node labels are read as integers rather than strings. Note: the package assumes nodes are labeled as integers 0, ..., n-1, where n is the network size. This way, the ith row of the data corresponds to node i.\n", "\n", "The `core` class of `networkinference` includes a function that outputs a table of network summary statistics." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "import networkinference as ni\n", "ni.core.sumstats(A)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Summary Statistics\n", "-------------------- --------\n", "# Units 1000.000\n", "# Links 2473.000\n", "Average Degree 4.946\n", "Max Degree 14.000\n", "# Isolates 11.000\n", "Giant Size 797.000\n", "Diameter 67.000\n", "Average Path Length 25.180\n", "# Components 37.000\n", "Clustering 0.578\n" ] } ], "execution_count": 3, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:04.609Z", "iopub.execute_input": "2022-04-10T00:49:04.621Z", "iopub.status.idle": "2022-04-10T00:49:06.545Z", "shell.execute_reply": "2022-04-10T00:49:06.624Z" } } }, { "cell_type": "markdown", "source": [ "The table indicates that we have a sparse network because the average degree is about 5, much smaller than the network size. The network is quite connected because there are few isolated nodes and the giant component contains a large share of nodes, a typical feature of social networks. The diameter and average path length are unusually large compared to typical social networks. Finally, the network exhibits high clustering (friends of friends tend to be friends), which is common." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "We will estimate a linear-in-means model with three regressors: average outcomes of neighbors, average covariate of neighbors, and own covariate. We instrument for the first regressor using the average covariate of 2-neighbors (friends of friends).\n", "\n", "To construct our regressor and instrument matrices, we first need to convert columns of our dataframe into arrays using the `numpy` package. `networkinference` requires node-level data to be provided as numpy arrays." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "import numpy as np\n", "Y = node_data['Y'].to_numpy() # array of outcomes\n", "X = node_data['X'].to_numpy() # array of covariates" ], "outputs": [], "execution_count": 4, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:06.558Z", "iopub.execute_input": "2022-04-10T00:49:06.568Z", "iopub.status.idle": "2022-04-10T00:49:06.582Z", "shell.execute_reply": "2022-04-10T00:49:06.631Z" } } }, { "cell_type": "markdown", "source": [ "To construct vectors of average outcomes and covariates of neighbors and 2-neighbors, we can use a function in the `utils` subpackage of `networkinference`. " ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "import networkinference.utils as nu\n", "nhbr_avg_Y = nu.nhbr_mean(Y, A)\n", "nhbr_avg_X = nu.nhbr_mean(X, A)\n", "two_nhbr_avg_X = nu.nhbr_mean(X, A, distance=2)" ], "outputs": [], "execution_count": 5, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:06.595Z", "iopub.execute_input": "2022-04-10T00:49:06.603Z", "iopub.status.idle": "2022-04-10T00:49:06.750Z", "shell.execute_reply": "2022-04-10T00:49:06.864Z" } } }, { "cell_type": "markdown", "source": [ "Finally, we construct the arrays of regressors and instruments. " ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "regressors = np.vstack([nhbr_avg_Y, nhbr_avg_X, X]).T\n", "instruments = np.vstack([two_nhbr_avg_X, nhbr_avg_X, X]).T\n", "\n", "pd.DataFrame(instruments, columns=['2 nhbr avg', 'nhbr avg', 'own']).describe() # just for printing summary statistics" ], "outputs": [ { "output_type": "execute_result", "execution_count": 6, "data": { "text/plain": " 2 nhbr avg nhbr avg own\ncount 1000.000000 1000.000000 1000.000000\nmean -0.011524 -0.004285 0.002677\nstd 0.451400 0.504014 0.983470\nmin -1.866928 -2.479348 -3.158637\n25% -0.260843 -0.309940 -0.630493\n50% -0.023416 -0.016504 -0.008020\n75% 0.220878 0.267530 0.619925\nmax 2.235487 2.860034 3.138272", "text/html": "
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
2 nhbr avgnhbr avgown
count1000.0000001000.0000001000.000000
mean-0.011524-0.0042850.002677
std0.4514000.5040140.983470
min-1.866928-2.479348-3.158637
25%-0.260843-0.309940-0.630493
50%-0.023416-0.016504-0.008020
75%0.2208780.2675300.619925
max2.2354872.8600343.138272
\n
" }, "metadata": {} } ], "execution_count": 6, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:06.764Z", "iopub.execute_input": "2022-04-10T00:49:06.773Z", "iopub.status.idle": "2022-04-10T00:49:06.795Z", "shell.execute_reply": "2022-04-10T00:49:06.867Z" } } }, { "cell_type": "markdown", "source": [ "For those unfamiliar with `numpy` arrays, note that `nhbr_avg_Y` and the other arrays should be thought of as $1 \\times n$ (i.e. horizontal) vectors. The above code stacks them vertically and then transposes the resulting matrix.\n", "\n", "Thus, `regressors` is an $n \\times 3$ matrix whose first column is the average outcomes of neighbors, and `instruments` is an $n \\times 3$ matrix whose first column is the average covariate of 2-neighbors." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "## Estimation\n", "\n", "Next we load the data into a `TSLS` object and output the two-stage least squares estimate." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls = ni.TSLS(Y, regressors, instruments, A)\n", "print(tsls.estimate)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "[0.93397924 0.52239823 2.85116605 1.00725669]\n" ] } ], "execution_count": 7, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:06.812Z", "iopub.execute_input": "2022-04-10T00:49:06.821Z", "iopub.status.idle": "2022-04-10T00:49:06.838Z", "shell.execute_reply": "2022-04-10T00:49:06.873Z" } } }, { "cell_type": "markdown", "source": [ "Given how `regressors` and `instruments` were formatted, these entries correspond to the intercept, endogenous peer effect, exogenous peer effect, and own covariate coefficient.\n", "\n", "As an aside, the `networkinference` package also includes an `OLS` class for the OLS estimator and an `IPW` class for the Horovitz-Thompson estimator. The methods discussed below are identical for all three classes." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "## Network HAC SEs\n", "\n", "We can call several different methods on the `tsls` object for constructing confidence intervals. The first uses a heteroskedasticity and network-autocorrelation consistent variance estimator due to [Kojevnikov et al. (2021)][3] with the bandwidth suggested by [Leung (2022a)][4]." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls.network_se()" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " Estimate SE 95% CI\n", "---------- ----- --------------\n", " 0.934 0.105 [0.728, 1.14]\n", " 0.522 0.046 [0.432, 0.612]\n", " 2.851 0.152 [2.553, 3.149]\n", " 1.007 0.055 [0.9, 1.115]\n" ] } ], "execution_count": 8, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:06.849Z", "iopub.execute_input": "2022-04-10T00:49:06.856Z", "iopub.status.idle": "2022-04-10T00:49:06.999Z", "shell.execute_reply": "2022-04-10T00:49:07.058Z" } } }, { "cell_type": "markdown", "source": [ "Network HAC SEs account for a general form of network dependence proposed by [Kojevnikov et al. (2021)][3] where correlation between two units decays with their network distance (measured by length of the shortest path). [Leung (2022a)][4] shows that this condition is satisfied in treatment effects models with interference and social interactions models such as linear-in-means." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "## Cluster-robust inference\n", "\n", "A common issue with HAC estimators is that, in smaller samples, they can produce CIs that are a bit too short. An alternative is to use cluster-robust inference methods, which appear to have better finite-sample performance [(Leung, 2022c)][6].\n", "\n", "However, cluster-robust methods rely on the existence of \"quality\" clusters in the network. (We will quantify what's meant by \"quality\" shortly.) As we'll see below, this is not always guaranteed, as some networks may have no quality clusters. In contrast, HAC estimators are asymptotically valid irrespective of the existence of such clusters.\n", "\n", "### Constructing clusters\n", "\n", "In order to construct the clusters, we first have to figure out how many to ask for. [Leung (2022c)][6] shows that the eigenvalues of the normalized Laplacian of the network can be used to detect the number of \"quality\" clusters in the network. The `core` class provides functions for calculating and plotting the eigenvalues." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "ni.core.plot_spectrum(A)" ], "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsQAAAEUCAYAAAAsgyAxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/IklEQVR4nO3deVhUZf8/8PcMMyAIihiQmk9pSjymhmnuYlhuwJgCmUtakvtGmGYqSm6oiSJqlvu3x699FRVRe8w9t7Dccg3UNBI1WUTAYRtmuH9/+GNyBGUUZg7MvF/X5XXNuc/M4XNgOPP25j73LRNCCBARERERWSm51AUQEREREUmJgZiIiIiIrBoDMRERERFZNQZiIiIiIrJqDMREREREZNUYiImIiIjIqjEQU5Vz7tw5DBo0CCqVCv7+/hg6dCiuXbv2XMe6cOECZsyYod8ODg5GRkZGRZWKZcuWYdasWWU+r6K/LhFRebz22mslrkmxsbEYMWIEACA6OhpxcXFPPcby5ctx4MABU5VIVKEUUhdA9Cw0Gg1GjBiBdevW4fXXXwcA7NixA8OGDcPBgwdhY2PzTMf7448/kJKSot/++eefK7ReY0n1dYmInkdISEiZz/n111/RqFEjM1RDVH4MxFSl5OXl4cGDB8jNzdW39erVC46OjtDpdNi+fTvWr18PuVyOWrVqYcGCBXB3d0dERATOnz+PnJwcCCEwZ84c1K1bF0uXLsWDBw8wZcoU/fE++ugjrFq1CnK5HLNmzcLff/+NwsJC+Pn5YeTIkbh16xYGDhyIV199Fbdv38b8+fMxYcIEtG7dGomJiRBCYMaMGWjVqpVB7deuXcOsWbOQmZkJmUyG4OBg9O7dW/+1i79unTp1zPPNJCJ6Tl988QUaN26MTz75BEuXLsX+/fuhVCpRq1YtzJs3D/v378elS5fw1VdfwcbGBm3btsXMmTORmJgImUyGTp06YcKECVAoFDhy5AgiIyMhl8vx73//G/Hx8fj+++9x8uRJbN26FXl5eXB0dMTKlSvx5ZdfIikpCVlZWahevToiIyPRsGFDDBo0CK+//jp++eUX3Lt3D4MHD8a9e/dw8uRJ5OXlYcmSJXjttdek/rZRZSaIqph169aJ5s2biy5duoiJEyeKLVu2iNzcXJGQkCDatGkj7ty5I4QQYv369WL69Oni7NmzYty4cUKn0wkhhFi5cqUYMWKEEEKIbdu2ieHDh+uP7eHhIe7duyeEEGLQoEHi4MGDQggh8vPzxaBBg8R///tfkZycLDw8PMSpU6eEEEK/vXPnTiGEEIcPHxYdOnQQGo1GLF26VMycOVMUFhaKd955R+zdu1cIIcTdu3dFp06dxNmzZ0t8XSIiqXl4eAh/f3/Rq1cv/b/OnTvrr5eTJ08Wa9asEXfu3BFvvvmmKCgoEEIIsXbtWrF//34hhBAffvih+PHHH4UQQnz++edi9uzZoqioSBQUFIjg4GCxcuVKkZGRIVq3bi0SEhKEEELExsYKDw8PkZycLLZt2ybeeust8eDBAyGEED/++KOYPXu2vsbp06eLWbNm6b/W2LFjhRBCnDt3Tnh4eOiv33PnzhVhYWGm/pZRFcceYqpyhgwZgvfffx+nTp3CqVOnsHr1aqxevRoqlQodO3bU97B+/PHH+tfUrFkTmzZtQnJyMn799VdUr179qV8jNzcXp06dQlZWFqKjo/VtiYmJaN68ORQKBby8vAyOr1KpAACdO3eGjY0Nrly5ot+flJSEgoICdOvWDQDg7u6Obt264dixY2jRokVFfFuIiCrUd999BxcXF/12bGws9u7da/Acd3d3eHp6ok+fPvD29oa3tzfatWtX4lhHjx7F//3f/0Emk8HW1hb9+vXDd999hwYNGuDVV1+Fp6cnAKBPnz6YM2eO/nWvvfYaHB0dAQA9evRA/fr1sWHDBvz11184efKkwfWza9euAID69esDADp16gQA+Ne//oWTJ09WxLeELBgDMVUpZ86cwW+//YahQ4fCx8cHPj4+mDBhgj6MymQy/XPz8/Nx+/ZtJCcnY+7cuRgyZAjeeecdNGzYEDt37nzq1ykqKoIQAps2bYK9vT0AICMjA3Z2drh//z5sbW2hUPzz6/P42OWioiKDtqKiohJfQwgBrVb77N8EIqJKQi6X43//939x8eJFnDhxAhEREWjTpg3CwsIMnvf4NbCoqAharRY2NjYQQpQ4ZjEHBwf94++//x4xMTEYOHAgVCoVnJ2dcevWLf1+W1tbg+Molcpynx9ZD84yQVWKi4sLvvnmG5w+fVrflpaWhry8PLz77rs4ceIEUlNTAQCbNm3CwoUL8fPPP8PHxwcDBgxAs2bNcODAAeh0OgAPg+yjobR429HREV5eXli/fj0AIDs7G/3798fBgwdLrSsjIwNHjx4FABw6dAhKpRIeHh76/Q0aNIBSqcS+ffsAACkpKdi7dy/at29fah1ERFVBYmIi/P398eqrr2LEiBH4+OOP9X8de/S61rFjR2zcuBFCCGg0GsTExKB9+/Z48803kZSUhMTERADA3r17kZ2dbdC5Uez48ePo06cP3n//fTRo0ACHDh3SX8uJyos9xFSlNGjQAF9//TWioqJw9+5d2NnZwcnJCbNmzYKnpycmTZqEoUOHAgBcXV0REREBtVqNiRMnQqVSwcbGBq1atcK+fftQVFSEFi1aYMmSJRgzZgy+/vprdO3aFQMGDMCKFSsQGRmJ2bNnQ6VSQaPRwN/fH7169TLokShmZ2eHHTt2IDIyEtWqVcPXX39t0EOsVCqxYsUKzJkzB8uWLYNOp8OYMWPQtm1bADD4uo8GaSKiyszT0xM9e/ZEYGAgHBwcUK1aNX3vsI+PDxYsWIDCwkKEhYVhzpw5UKlUKCwsRKdOnTBy5EjY2tpi8eLFmDx5MuRyOZo2bQqFQqH/y9yjgoODMWPGDMTGxsLGxgavv/46rl69au5TJgslE4//rYKInsmtW7egUqnw22+/SV0KEVGVolarsWLFCowbNw729va4fPkyRowYgWPHjpXaS0xkKuwhJiIiIkk4OjpCqVQiKCgICoUCCoUCS5YsYRgms2MPMRERERFZNd5UR0RERERWjYGYiIiIiKwaAzERERERWTUGYiIiIiKyalY3y8T9+zkoKuJ9hERkGnK5DLVqPX1pcEvGaywRmZKprrFWF4iLigQv1kREJsJrLBFVRSYdMqFWq+Hv71/qyl6XL19GYGAgevXqhREjRiA7OxsAkJWVhWHDhqFXr14ICgpCQkICACAnJwfjxo2DSqVC7969ER8fb8rSiYiIiMhKmCwQnz9/Hv3790dSUlKp++fOnYvx48dj586daNCgAdauXQsAWL9+PTw8PLBz506MHj0as2bN0re//PLL2LVrFxYtWoTPP//cVKUTET1VvkYrdQlERFSBTDZkIiYmBuHh4U8MrkVFRcjJyQEA5OXloWbNmqW2V6tWDQAwduxYaLUPP4Ru3bqlfz4RkTl9E3cJpxJT8ZanG0b1bip1OUREVAFMFojnzp371P1ffPEFhgwZgoiICNjb2yMmJgYAEBwcjA8++AAdO3ZETk4O1q1b90+xCgU++eQTnDhxQt9zTERkLvkaLU4lpgIATiWmYohGi2q2VncrBhGRxZFk2rX8/HxMmzYN3333HY4fP44BAwZg8uTJAIDZs2dj4MCBOH78ONatW4fQ0FB9jzEArF27Fvv370d0dDSuX78uRflEZKWmrf5V//gtTzeGYSIiCyFJIL569Srs7OzQvHlzAMAHH3yAkydPAgAOHjyIwMBAAECLFi1Qu3ZtXL9+HSdPnkRq6sOemXr16qFFixa4du2aFOUTkRXKVBfg/oMC/Xb/dxtLWA0REVUkSQLxyy+/jLt37+LGjRsAHobgZs2aAQA8PT1x4MABAEBSUhJSU1PRoEEDHD58GKtWrQIApKam4tKlS/rXEBGZmrOjHWo52QEAajnZwdnRTuKKiIioopj1733Dhg3D+PHj0axZM8ybNw+ffvophBCoXbs2IiIiAADz58/HjBkzsHr1atja2mLBggVwcnLC6NGjMW3aNKhUKtjY2GDq1KmoV6+eOcsnIivXqF5NnEpMRaN6vKmXiMiSyIQQVjWD+r17ak4aT0TPLF+jxejFR/XbKyZ4lzqGWC6XoXZtR3OWVqnwGktSqFHTHnYVMKa/QKNFdlZeBVREpmKqayzvCCEiMkI1WwXe8nTTT7nGG+qIKg87WwUmRR8p93EWhnSugGqoKuIVnYjISKN6N+VUa0REFkiSm+qIiKqiTHUBwzARkQXilZ2IyAifff0z7j8oQC0nOywa00HqcoiIqAKxh5iIqAyPzkF8/0EBMtUFZbyCiIiqEgZiIqIycA5iIiLLxkBMRGSEBnVqAADnICYiskAMxEREZVgeewFnr6YBAE4lpiJfo5W4IiIiqkgMxERET5Gv0eLs1XT99pserpxpgojIwjAQExE9g6H+/5a6BCIiqmAMxERET7Hmh9+lLoGIiEyMgZiI6Ak4XIKIyDowEBMRGYnDJYiILBMDMRHRE1SzVejnHK7lZMfeYSIiC8VATET0BMtjL+hXpbv/oIDTrRERWSgGYiKiUnD8MBGR9WAgJiIqRb5GZ7D9YTcPiSopP7VaDX9/f9y6dQsAsHnzZvj7+0OlUmHKlCnQaDQAgISEBAQGBqJ79+6YNm0atFr2iBORdWAgJiIqhbOjHZwdbQE8HD9cPJa4qjl//jz69++PpKQkAMCff/6JtWvXYtOmTdi5cyeKiorw/fffAwAmTZqE6dOnY+/evRBCICYmRsLKiYjMh4GYiKgUn339MzLVGtSsrsSiMR2kLue5xcTEIDw8HG5ubgAAW1tbfPnll3B0dIRMJoOHhwfu3LmD27dvIz8/H15eXgCAgIAA7NmzR8LKiYjMhwPiiIgek5KRi/sPHt5Ml5VTiEx1QZXtIZ47d67Bdr169VCvXj0AQEZGBjZu3Ih58+YhNTUVrq6u+ue5uroiJSXFrLUSEUmFgZiI6BHfxF3CqcRUg7ZqtjYSVWM6KSkpGDp0KAIDA9GmTRucPXu2xHNkMpkElRERmR+HTBAR/X/5Gm2JMNysoYvFzS5x/fp19O/fH3369MGYMWMAAO7u7khP/2dWjbS0NP0wCyIiS8dATET0BAobGUL7ekldRoVSq9X45JNPEBISguDgYH17vXr1YGdnhzNnzgAA4uLi4O3tLVWZRERmZVndHkRE5bDmh98NtpeGdJKoEtPZunUr0tPTsW7dOqxbtw4A0KVLF4SEhCAyMhJhYWHIyclBkyZNMHjwYImrJSIyDwZiIiJY/kIchw4dAgB8/PHH+Pjjj0t9jqenJ7Zu3WrGqoiIKgfLudoTEZXD+t2J+scKGxnGBjSTsBoikkKhtgiurk7lPk6BRovsrLwKqIjMhYGYiKze4zfTaXWiSk+1RkTPR6mQY1L0kXIfZ2FI5wqohsyJN9URkdWbsvKEwfbDVeoYhomIrAV7iInIqoVEH8WDPK1BW8TwNhJVQ0REUmAgJiKrNeyrQ9AVGbZZ2s10RERUNg6ZICKrNP9/z5QIw280eoE30xERWSF2gxCRVcnXaDF+yTFoi4RB+5seL2BsQHOJqiIiIikxEBOR1YiKOYeLNzJKtM8b3hbuLg4SVERERJUBAzERWYXhC3+CVidKtNdysmMYJiKycgzERGTRUjJyMfc/p0qEYbkMiBzTgdOrERERAzERWa4RCw+j8PE75wDUcFBgyXhvCSoiIqLKiIGYiCxGproABRodAGDaql9QMgpzvDAREZXEQExEFmHC8uPIVGue+py3PN0YhomIqAQGYiKq0vI1Wkz+Jr7EanPF5ADmDm+Lmo62XHCDiIhKxU8HIqqylsdexNmraU/cL5cBayZ3MWNFRERUFTEQE1GVk6/RYnnsRfyedL/Evn+/4oLB3TwAgMMjiIjIKAzERFSlPK1XmKvNERHR82AgJqIqI3LTb6X2CjvZ22DBqA4cI0xERM+Fnx5EVOllqgvw+Yp4aItKrjRXs7oSUeM6SVAVERFZCgZiIqrUQpcdR1ZOyenUnOxtMPOTtlxpjoiIyo2BmIgqnXyNFllqDcJW/wJdyU5h9goTEVGFYiAmokojU12A9bsTcPFGRqn75TIgckwH9goTEVGFYiAmIsmlZORi/sazpQ6NKNa0oQsm9PUyX1FERGQ1GIiJSBLFwyJmrD2JQl3RE58nlwHLQ705gwQREZkMP2GIyKzyNVp8E3fpicMiitVwUGDKh624uEYFUKvV6NevH7799lu89NJLiI+Px7x581BQUICePXsiNDQUAJCQkICwsDCo1Wq0atUKM2fOhELBjwkisnxyqQsgIuuQkpGLqJhzGL346FPDcNOGLlg8tgOWjPdmGK4A58+fR//+/ZGUlAQAyM/Px9SpU7FixQrs3r0bly5dwpEjRwAAkyZNwvTp07F3714IIRATEyNh5URE5sP/+hORSWWqCzD5mxNlDouYO6wtajracmhEBYuJiUF4eDg+//xzAMCFCxfw8ssvo379+gAAlUqFPXv2oFGjRsjPz4eXlxcAICAgAEuXLsWAAQOkKp2IyGz4yUNEJmHsjXID3/VgT7AJzZ0712A7NTUVrq6u+m03NzekpKSUaHd1dUVKSorZ6iQikhIDMRFVKGN6hGs4KDB/ZHv2BktAiJITO8tksie2ExFZA34aEVGFKKtHuHhYhJ2tDecRlpC7uzvS09P126mpqXBzcyvRnpaWBjc3NylKJCIyOwZiInomxdOlPaqsqdM4h3Dl8cYbb+DPP//EX3/9hZdeegk//PADAgMDUa9ePdjZ2eHMmTNo2bIl4uLi4O3tLXW5RERmwUBMRGUqDsHfH7ha5nRpxXijXOVkZ2eH+fPnY9y4cSgoKEDnzp3Ro0cPAEBkZCTCwsKQk5ODJk2aYPDgwRJXS0RkHvyUIqKnioo5Z3QILsYe4crn0KFD+sft2rXDzp07SzzH09MTW7duNWdZRESVAgMxEZVQ3CM8fc2v0BaVvNmqNOwRJiKiqoqfWkRkYOnW8zj3x72nPqd4FblHceo0IiKqqhiIiQiZ6gIUaHSY979nkJ1bWOpzikMwZ4kgIiJLw0BMZMWMWTxDLgMix3RgCCYiIovFQExkhYxZPMPJXoGpg1pxKAQR0TMq1BbB1dWpXMco0GiRnZVXQRVRWRiIiaxE8bCIsnqEAeCNRrUREvSGmSojIrIsSoUck6KPlOsYC0M6V1A1ZAwGYiIrMGH5cWSqnx6CmzZ0wcB3PThLBBERWR1+6hFZgOLe39JEbDiNB3naUvdxqjQiIiIGYqIqLVNdgJnrT5U5BKI0XDyDiIjoIQZioirk0Z5gY8YCP+7fr7hgcDcOiyAiInoUPxGJqojQZcefqye42JseL2BsQPMKrIiIiMgyMBATVWLFSyiHrf4FujJWUC6+Ka407BEmIiJ6Mn5CElUCxcH3Ud8fuIqLNzKe+jobOTBnKG+KIyIiKg9+ghJJKFNdgPW7E8oMvo96tCeYi2YQERGVn1GB+Pr16zh79iyCgoIwZswYXLlyBXPnzkXbtm1NXR+RRXqe2SHkMmB5qDd7gomIiCqYUZ+s4eHh6Nu3L3766Sfcv38fERERWLx4MTZv3mzq+ogsSkpG7jPPDlHDQYEpH3IJZSIiIlMxKhAXFBSgV69emD17Nnr27Ik2bdqgsLDQ1LURVVmljQmesfYkCnVFT3xNcfB9lJ2tDZwd7UxSIxERET1kVCDWaDRIT0/H4cOHsXLlSqSnp6OgoMDUtRFVSctjL+Ls1TSjn9+0oQuCff/N4EtERCQRowLxBx98AB8fH/Ts2RONGjXC22+/jdGjR5u6NqIq4dHe4A37ruD3pPtGva5pQxeM7t2UY4KJiIgkZtQn8YABA9CvXz/I5XIAwPbt21GrVi2TFkZUmaRk5JbabszUaI+Sy4C5wzhNGhERUWVi1CdyTk4OFi1ahOvXryM6OhpRUVGYPHkyqlevbur6iCSVkpFb5tjfJ3Gyt8HUQW8ZtPHGOCIiosrHqEA8Z84cuLm54d69e7Czs4NarcaMGTOwaNEiU9dHZDaP3wj3vEEYAGpWVyJqXKeKKo2IiIhMyKhAnJCQgHnz5uHIkSOwt7dHZGQk/P39TV0bkcnla7TI1+ieeXGMxz06QwRnhiAiIqpajArExWOHi+l0uhJtRFVFcU/ws4z/LR77WxoGYCIioqrNqED81ltvYeHChcjPz8exY8ewceNGtGnTxtS1EVWYTHUBCjS6574JjmN/ydLs2LEDq1atAgB4e3tj8uTJSEhIQFhYGNRqNVq1aoWZM2dCoeDNn0Rk+Yy60k2cOBGrVq2Ck5MToqKi0KlTJ067RpVecU/ws6wM9/jiGAzCZIny8vIwd+5c7NmzBzVq1ED//v0RHx+PiIgIzJkzB15eXpg6dSpiYmIwYMAAqcslIjI5owKxUqnEmDFjMGbMGFPXQ1RuKRm5z9QTXLO6El8MbMmhD2Q1dDodioqKkJeXBwcHB2i1WigUCuTn58PLywsAEBAQgKVLlzIQE5FVMCoQjxw5stT2b7/9tkKLISqvEQsPGzUzRHFPMEMwWSNHR0eEhISgZ8+eqFatGlq3bg2lUglXV1f9c1xdXZGSkiJhlURE5mNUIO7evbv+cWFhIQ4dOoTXXnvNZEURGevRqdKmrfoFT4vCTRu6YOC7HgzBZPUSExOxbds2/PTTT3BycsLEiRPx888/l3ieTCaToDoiIvMzKhD36dPHYDsgIAAffvihSQoiKsuzzBJR3BPMleGI/nH8+HG0a9cOtWvXBvDwmr527Vqkp6frn5OWlgY3NzepSiQiMqvnSghFRUVITU2t6FqInipTXWDUfMFyAHOHt2VPMNETeHp6YuHChcjNzYW9vT0OHTqE1q1bY+/evThz5gxatmyJuLg4eHt7S10qEZFZPNcY4qtXr6J169YmKYjocZnqAsxcf8qomSLkMmDN5C5mqIrIvKZOnYqIiAiDtnHjxmHZsmXPfKyOHTvi999/R0BAAJRKJZo1a4bhw4eja9euCAsLQ05ODpo0aYLBgwdXVPlERJXaM48hlslk6N+/Pzp27GiyooiAh0Mjpqz8pcwg/OhUaZwmjSxNeHg4UlJScObMGWRk/PPXEa1Wixs3bjz3cYcPH47hw4cbtHl6emLr1q3PfUwioqrqqYE4MzMTAODj41Ni34MHD+Ds7GyKmogQFXPuqUMjOEsEWYugoCBcu3YNV65cMeicsLGxQYsWLSSsjIjIcjw1ELdt21Z/l7EQAsDDHmIhBGQyGRISEkxfIVmN4pvlpq/5FdoiUepzmjZ0QbDvvxmCyWo0a9YMzZo1Q/v27fHiiy9KXQ4RkUV6aiBOTEw0Vx1k5crqEW7a0AWjezflTBFktW7evIlJkyYhKytL30EBALt27ZKwKiIiy2BUutBoNDhy5AhycnIAPFzl6ObNmwgNDTVpcWT5UjJyn9ojXMNBgfkj2zMIk9WbNWsWAgMD0aRJE84PTERUwYxKGaGhoUhOTkZaWhqaNGmC8+fPc5YJKpdMdQEmf3PiiavKyWVA5JgOHBpB9P8plUoMGTJE6jKIiCySUYE4ISEB+/btw5dffokhQ4ZACIGZM2eaujayUKHLjj9x5ojim+U4WwSRocaNG+PKlStcJZSIyASMCsRubm5QKBR45ZVXcPXqVfTs2RN5eXmmro0sTKa6AJO+/hm6UkZHsEeY6OmSk5MRGBiIunXrws7un98TjiGmqq5GTXvYcVgcScyod6CDgwN27doFT09PxMTEoGHDhvop2YiMMWH5cWSqS/YKy2XA3GFt2SNMVAbes0GWys5WgUnRR8p1jIUhnSuoGrJWcmOeNGPGDCQkJKBDhw6Qy+UYNGgQPvnkkzJft2vXLvj6+qJr167YuHFjif1HjhyBSqWCSqXCZ599pr9pr9jdu3fRunVr3Lp1y6D9ypUr8PPzM6Z0kli+RouQ6KOlhuGmDV2wZnIXhmEiI3h4eJT6j4iIys+oHuLk5GR8/vnnAIAlS5YYdeCUlBRERUUhNjYWtra26NevH9q0aYNGjRoBALKzs/HFF19gw4YNaNSoEVavXo2oqCiEhYUBAIqKijBt2jQUFhYaHDcuLg6LFi2CUqk09hxJIt/EXcKpxNQS7XIZsDzUmzNHED2D4nnhi+eBBwBXV1ccPXpU4sqIyBQKtUVwdXUq93EKNFpkZ3GYa1mMSiTLly9HeHg4AgIC8P7778Pd3b3M18THx6Nt27b61ey6d++OPXv2YOzYsQCApKQk1K1bVx+QfXx8MHToUH0gXrNmDdq3b48///xTf8wHDx7g4MGDWLx4MSZPnvxMJ0rmla/RlhqG//2KCyb18zJ/QURV3KPzwhcWFmLfvn2cK57IgikV8nIPJQE4nMRYRg2Z2Lx5M1avXo28vDz07dsXI0aMwIEDB576mtTUVLi6uuq33dzckJKSot9+5ZVXcPfuXf0F/ccff0R6ejoA4NKlS/j1119LTDHk5OSEZcuWoU6dOsadHUlmysoTJdre9HiBYZioAiiVSvj5+eHnn3+WuhQiIotg9N+sX331VUyaNAndu3fHnDlzMGHCBFy4cOGJz390JaVij04mX6NGDSxYsADTp09HUVER+vbtC6VSiby8PMyaNQtLliyBXG5UXqdK5tOlR5GdqzVoWzyWM0gQlcejNzILIXDp0iVkZ2dLVxARkQUxKhDfu3cPO3fuxPbt26HT6RAUFISVK1c+9TXu7u44ffq0fjs1NRVubm76bZ1OhxdffBFbtmwBAFy+fBn169fH6dOnkZ6ejlGjRulfN3z4cCxfvhwNGzZ85hMk8/p06bESYfhND1eGYaJyenQMMQDUrl0b06ZNk7gqIiLLYFQg7tatG7p164bw8HC0bNnSqAO3b98ey5YtQ0ZGBuzt7bFv3z7Mnj1bv18mkyE4OBhbtmyBm5sb1q1bB19fX3Tq1AmHDh3SP69Lly5YtWoVXnrppWc8NTK3xTHnkZ1reBPkG41ewNiAZhJVRGQ5OF6YiMh0jArER44cgaOj4zMd2N3dHaGhoRg8eDAKCwsRFBSE5s2bY9iwYRg/fjyaNWuGWbNmYejQodBoNGjXrp1RU7lR5bQ89gIu3bhn0PZGoxcQEtRcooqILEtRURHWrl2Lo0ePQqvVokOHDhg5ciQUCs7WQkRUXkZdSX///XcsW7YMWVlZBmODy1ohqXiO4UetXr1a//jtt9/G22+//dRjPNpbXOyll14qtZ2kka/R4uzVdIO2Zg1dGIaJKtCiRYuQmJiIjz76CEVFRdi8eTO++uorTJ06VerSiIiqPKMC8axZsxAYGIgmTZoY3BhH1i1fo0WWWoOIDacN2mtWVyK0r5c0RRFZqGPHjmHbtm36Odjffvtt9OrVi4GYiKgCGBWIlUpliSnQyLo9adGNGg4KRI3rJEFFRJZNCGGwIJGtrS0XKCIiqiBGzWvWuHFjXLlyxdS1UBXxpEU3nOwVWDLeW4KKiCyfp6cnIiIicPPmTdy8eRMRERFcupmIqIIYvXRzYGAg6tatCzu7f6bPKmsMMVmmqat+KdFmIweiQxiGiUwlPDwcc+bMQb9+/VBUVIROnTph+vTpUpdFRGQRjArEoaGhpq6DqoiQ6KN4kGc4z/CUD99E45ecpSmIyMJpNBpMnz4dXbt2xfz58wEAw4cPh42NzTPP/kNERKUzashE69atUa1aNdy4cQNeXl5QKpVo3bq1qWujSmbYgkMlwrCzox3DMJEJLV26FGq1Gi1atNC3zZ49G9nZ2Vi2bJmElRERWQ6jAnFsbCymTJmCNWvW4MGDBxg9ejRiYmJMXRtVIkPnH4LusdW4mzasjcVjO0hTEJGVOHz4MBYtWoTatWvr29zd3fHVV1/hwIEDElZGRGQ5jArEGzZswObNm+Ho6IjatWsjNjYW3333nalro0pi6IJDKHqs7U2PFzCh7xuS1ENkTZRKJapVq1ai3dHREba2ts993EOHDiEgIAA9evTAnDlzAADx8fFQqVTo1q0boqKinvvYRERVjVFjiOVyucFYtTp16sDGxsZkRVHlkZKRi6JHeoblAJZP8EY1W66ORWQOcrkcarW6xHhhtVoNrVb7hFc9XXJyMsLDw7FlyxbUrl0bH330EY4cOYLw8HBs2LABderUwYgRI3DkyBF07tz5mY5du3b5xjUXaLTIzsor1zGIiJ6VUanG2dkZCQkJ+kU5du7ciZo1a5q0MKocYo/e0D+Wy4A1k7tIWA2R9fH390dYWBgiIiLg4OAAAMjNzUVYWBi6dev2XMfcv38/fH198eKLLwIAoqKi8Ndff+Hll19G/fr1ATxcaXTPnj3PHIgj1v2C+w8KnqsuAFgY8mxfj4ioIhgViKdOnYqQkBDcvHkTHTt2hJ2dHVasWGHq2khimeoCg/mGl4dyWjUic/voo48QHh6ODh06oHHjxigqKsL169ehUqkwZsyY5zrmX3/9BaVSiU8++QRpaWnw8fFB48aN4erqqn+Om5sbUlJSKuo0iIgqNaMC8auvvoodO3YgKSkJOp0ODRo04ApJFm557EWcvZomdRlEVk8ul2P27NkYMWIEfv/9d8jlcjRr1gzu7u7PfUydTofTp09jw4YNcHBwwOjRo2Fvb1/iecV/FSQisnRGBeLly5cbbMtkMtjb26Nx48bo1InL9Fqa5bEXcPZqutRlENEjXnrpJbz00ksVcqwXXngB7dq1g4uLCwDgnXfewZ49ewzuDUlNTYWbm1uFfD0iosrOqFkmrl69is2bNyMzMxMPHjzAtm3b8NNPP2Hp0qX4+uuvTV0jmVG+RltqGH7L04030hFZCB8fHxw/fhzZ2dnQ6XQ4duwYevTogT///BN//fUXdDodfvjhB3h7c5gUEVkHoxLOvXv3EBsbqx9fNnLkSISEhGDjxo0IDAx87nFsVPk52SuwYFR7hmEiC/LGG29g6NChGDBgAAoLC9GhQwf0798fDRs2xLhx41BQUIDOnTujR48eUpdKRGQWRqWczMxMg5statWqhczMTNja2kKhYFCyJOt3J+of16yuRNQ4DokhskRBQUEICgoyaGvXrh127twpUUVERNIxashE/fr1sWjRIiQnJyM5ORlRUVH417/+hfPnz0MuN+oQVAXka7QGs0pk5RQiX/N885wSERERVRVGpdmIiAjcvn0bffr0QVBQEFJSUjBnzhxcvnwZkydPNnWNZCZTV/1qsM1xw0RERGQNjEo7Li4uWLx4cYn2AQMGVHhBJI3ITb8hU204mX7/dxtLVA0RERGR+Tw1EIeEhCA6OhoqlarU/bt27TJJUWReocuOIytHY9Dm7GgHZ0c7iSoiIiIiMp+nBuJhw4YBAKZPn26WYsj8oreeLxGGm7xSCxP7tZCoIiIiIiLzemogLp60vXXr1iX2HT161DQVkdnka7Q4/8c9g7aa1W0ZhomIiMiqPPWmukfnFx43bpzBvqioKNNURGazIu6ywXaTV2ohalxHiaohIiIiksZTA7EQQv84OTn5ifuo6lkeewGXbvzTO8yeYSIiIrJWTw3EMpms1MelbVPVUdryzOFD3pKoGiIiIiJpGd1DTJZjfPQxg23OKEFERETW7Kk31RUVFSErKwtCCOh0Ov1jANDpdGYpkCrWp0uPQqsz/I9OxPA2ElVDREREJL2nBuKrV6+ibdu2+hDcps0/wYlDJqqeTHUBsnMNl2LmanRERERk7Z6ahBITE81VB5nBrP85ZbC9eGwHDpUgIiIiq/fUMcRkOTLVBchU/7MAh1ej2gzDRERERGAgthqP9w4P7/W6RJUQERERVS4MxFagtN5hjhsmIiIieoiB2Ao4O9pBqXj4o1bYyDA+6A2JKyIiIiKqPBiIrUDosuMo1BYBALQ6gXyNtoxXEBEREVkPBmILF731PLJyOFyCiIiI6EkYiC1YvkaL83/cM2jjzXREREREhhiILdjy2IsG2+wdJiIiIiqJ6chChS47bjBUomZ1W95MR0RERFQK9hBboMUxhuOGASB8yFsSVUNERERUuTEQW5jlsRdw6YbhuGFnRzuuSkdEpVqwYAG++OILAEBCQgICAwPRvXt3TJs2DVotZ6QhIuvAQGxB8jVanL2abtDW5JVaWDy2g0QVEVFlduLECWzfvl2/PWnSJEyfPh179+6FEAIxMTESVkdEZD4MxBZk5c7LBts1qysxsV8LiaohososMzMTUVFRGDlyJADg9u3byM/Ph5eXFwAgICAAe/bskbBCIiLzYSC2EI9PsVbDQYmocZ0krIiIKrMZM2YgNDQUNWrUAACkpqbC1dVVv9/V1RUpKSlSlUdEZFYMxBZiysoTBttfBreWqBIiquy2bNmCOnXqoF27dvo2IUSJ58lkMnOWRUQkGU67ZgEy1QXIyinUbzs72vImOiJ6ot27dyMtLQ3vvfcesrKykJubC5lMhvT0f+5BSEtLg5ubm4RVUmVXo6Y97Di3PVkIvpMtQDVbGyht5CjUFcFGDiwe21HqkoioElu/fr3+cWxsLE6ePIl58+bB398fZ86cQcuWLREXFwdvb28Jq6TKzs5WgUnRR8p9nIUhnSugGqLyYSCu4r6Ju4RTian6bV3Rw/HEXJGOiJ5VZGQkwsLCkJOTgyZNmmDw4MFSl0REZBZMTVVYvkZrEIYBoFlDF4ZhIjJaQEAAAgICAACenp7YunWrxBUREZkfb6qzIAobGUL7ekldBhEREVGVwkBchWWpDZdnXhrCadaIiIiInhX/tl4F5Wu0mLLyF2Tl/BOI3/J041AJIiIioufABFXFPH4TXbEA74YSVENERERU9XHIRBVS2k10AKBUyOHu4iBBRURERERVH3uIq5DHV6MDgKYNXTCBN9IRERERPTcG4iri8dXonOxtsGBUB44bJiIioicq1BbB1dWp3Mcp0GiRnZVXARVVTkxTVcT/Hbimf2wjB6K5sg8RERGVQamQc0VBI3AMcRXw+Njh4tXoiIiIiKj8GIirgPW7Ew22OcUaERERUcVhIK7kSptZov+7jSWqhoiIiMjyMBBXcmt++N1g29nRDs6OdhJVQ0RERGR5GIgrsXyNFmevphu0RQxvI1E1RERERJaJgbgSq2arwFuebvptjh0mIiIiqnhMV5XctVtZAABnR1uM6t1U4mqIiEyLc6YSkRQYiCux6K3nkakuAABkqjXIVBdw/DARWTTOmUpEUuCQiUoqX6PF+T/u6bd5Mx0RERGRaTAQV1IhS4/rH9vIgcVjO0hYDREREZHl4pCJSigk+igKtUX6bV0ROFyCiOgZVLaxyDVq2sOunDdFV6ZaiCwNfyMqmeit5/Egz3BZ5lpOHC5BRPQsKttYZDtbRbnrqUy1ABynTZaFgbgSeXzcMADMG94W7i4OElVEREREZPkYiCuJfI0W45YcNWjzalSbYZiISEIVNfSiIlSmWogsDQNxJfBN3CWcSkw1aKvhoMT4oDckqoiIiIDKNfSiMtVCZGk4y4TE8jXaEmHYRg4sGd9JooqIiIiIrAsDscTyNTqDbSd7BVZ/3kWiaojIWixfvhx+fn7w8/PDV199BQCIj4+HSqVCt27dEBUVJXGFRETmwyETEiptqMSCUe0lqoaIrEV8fDyOHz+O7du3QyaTYejQofjhhx8QGRmJDRs2oE6dOhgxYgSOHDmCzp3553UisnzsIZZIaUMlmjV0QTXODUlEJubq6oovvvgCtra2UCqVePXVV5GUlISXX34Z9evXh0KhgEqlwp49e6QulYjILBiIJVLNVoE3PVz12wobGUL7eklXEBFZjcaNG8PLywsAkJSUhN27d0Mmk8HV9Z9rkpubG1JSUiSqkIjIvBiIJfJN3CWcvZoGAKjhoMCqST4SV0RE1ubatWsIDg7G5MmT8a9//avEfplMJkFVRETmx0AsgceHS2TnapGpLpCwIiKyNmfOnMHHH3+Mzz77DH369IG7uzvS09P1+1NTU+Hm5iZhhURE5sNALIGVOy8bbDs72nJpZiIym7///htjxoxBZGQk/Pz8AABvvPEG/vzzT/z111/Q6XT44Ycf4O3tLXGlRETmwTu4zOzx5ZlrOCixeGxHCSsiImuzdu1aFBQUYP78+fq2fv36Yf78+Rg3bhwKCgrQuXNn9OjRQ8IqiYjMh4HYzL749oTB9pfBrSWqhIisVVhYGMLCwkrdt3PnTjNXQ0QkPQ6ZMKPoreeRnVuo3+ZQCSIiIiLpMRCbyeNDJQAgYnhbiaohIiIiomIMxGYyZeUvBttverhyEQ4iIiKiSoCJzAyit55HVo7GoG2o/78lqoaIiIjo2RRqi+Dq6lTu4xRotMjOyquAiioWA7GJZaoLSgyVYO8wERERVSVKhRyToo+U+zgLQzpXQDUVj6nMhJZuPY9zj4XhZg1dMDagmUQVEREREdHjGIhNJHTZ8RLDJGo4KBHa10uagoiIiIioVLypzgRKGzNsIweWjO8kUUVERERE9CTsIa5gpU2v5mSvQHQIl0AlIiIiqowYiCtQvkZbYiW6pg1rY0LfNySqiIiIiIjKwkBcQZbHXsTZq2kl2kf3fl2CaoiIiIjIWAzEFSAq5hwu3sgo0c7p1YiIiIgqP6a1ciptNgkAeNPjBU6vRkRERFQFMBCXQ+Sm30qEYSd7GywY1YE9w0RERERVBFPbc/p06TFk5xYatNVwUHJqNSIiIqIqhoH4GeVrtBgXdRQ6Ydje5JVamNivhTRFEREREdFzYyB+BqUtxQwANavbMgwTERERlaFQWwRXVyepyyjBpIF4165d+Oabb1BYWIiPP/4YAwcONNifkJCAsLAwqNVqtGrVCjNnzoRCocCdO3cwadIk3Lt3Dw0aNEBkZCSqV6+O7OxsTJw4EcnJyXBxccGSJUvg6upqylMAAGSqC/DlupMlhkgAwBuNaiMkiPMMExEREZVFqZBjUvSR5359LSc7TA1uW4EVPWSypZtTUlIQFRWF77//Hjt27MDmzZvxxx9/GDxn0qRJmD59Ovbu3QshBGJiYgAAM2fOxIABA7Bnzx40bdoUK1asAAAsWbIErVq1wo8//oj3338fc+fONUnt+RotUjJykZKRi9BlxzFh+c8lwrBcBqyY4M0wTERERFTFmSwQx8fHo23btnB2doaDgwO6d++OPXv26Pffvn0b+fn58PLyAgAEBARgz549KCwsxKlTp9C9e3eDdgA4fPgwVCoVAMDf3x9Hjx5FYWHJXtvnlZKRi6iYcxi9+CimrPoFU1b9UuqUak72CqyZ3IUzSRARERFZAJMlutTUVIPhDG5ubrhw4cIT97u6uiIlJQX379+Ho6MjFAqFQfvjr1EoFHB0dERGRgbc3d3LXe+IyMMo1BaV+TwOkSAiIiKyLCYLxEKIEm0ymazM/WW97nFyefk7uVMycssMw00bumB076bsFSYiIiKyMCZLd+7u7jh9+rR+OzU1FW5ubgb709PT9dtpaWlwc3ODi4sL1Go1dDodbGxs9O3Aw17m9PR0vPjii9BqtVCr1XB2di5/rS4OUCrkBqG4hoMCUz5sBQCo6WjLIExERERkoUyW8tq3b49ly5YhIyMD9vb22LdvH2bPnq3fX69ePdjZ2eHMmTNo2bIl4uLi4O3tDaVSiVatWmH37t1QqVT6dgDo3Lkz4uLiMHLkSOzevRutWrWCUqmskHpXTnwbKRm5AAA7Wxs4O9pVyHGJiIiIqHIz2U117u7uCA0NxeDBg9G7d2/4+/ujefPmGDZsGC5evAgAiIyMxLx589CzZ0/k5eVh8ODBAIDw8HDExMTA19cXp0+fxqeffgoACAkJwblz5+Dn54fvv/8eM2bMqNiaXRzg7uLAMExEVm3Xrl3w9fVF165dsXHjRqnLISIyOZOOA1CpVPpZIYqtXr1a/9jT0xNbt24t8bp69ephw4YNJdqdnZ3x7bffVnyhREQE4J8pM2NjY2Fra4t+/fqhTZs2aNSokdSlERGZjNUNjJXLn3yDHhFReVX1a8yjU2YC0E+ZOXbsWKNeX7MC/sJWy6li/kpnicepTLVUtuNUploq6jiVqZbKcpyKuMaURiZKm9aBiIis0sqVK5Gbm4vQ0FAAwJYtW3DhwgWDe0CIiCyNycYQExFR1fOsU18SEVkCBmIiItJ7fErMx6fMJCKyRAzERESk1759e5w4cQIZGRnIy8vDvn379FNfEhFZKqu7qY6IiJ7s0SkzCwsLERQUhObNm0tdFhGRSfGmOiIiIiKyahwyQURERERWjYGYiIiIiKwaAzERERERWTUGYiIiIiKyalYXiHft2gVfX1907doVGzduLLE/ISEBgYGB6N69O6ZNmwatVitBlWQq0dHR8PX1hZ+fH9avXw8AiI2Nha+vL1QqFebMmaP/md+6dQsDBw7Ee++9h0GDBuH27dsljldYWIg333wT7733nv6fTqcz6zlR+WzZssXg59eyZUvMmjWr1PfFvXv3DJ7bpUsXtGjRosQxNRoNJk2ahJ49e6JPnz64fv26BGdWPs97rbxz5w4GDhyIHj16YNSoUcjJyQEAZGdnY/jw4ejZsycGDhyItLQ0s54PmUdZ75sjR45ApVJBpVLhs88+078/it29exetW7fGrVu3DNqvXLkCPz8/k9ZO0lKr1fD39y/xsweAy5cvIzAwEL169cKIESOQnZ0NAMjKysKwYcPQq1cvBAUFISEhAQCQk5ODcePGQaVSoXfv3oiPjy+7AGFF7t69K3x8fMT9+/dFTk6OUKlU4tq1awbP8fPzE7/99psQQogpU6aIjRs3SlApmcKvv/4q+vXrJwoLC0VeXp7w8fER169fF506dRIpKSlCCCHCw8PFunXrhBBCTJw4Uf/z/89//iM+++yzEse8ePGiCA4ONt9JkEldvXpVdO3aVZw7d+6J74tiOp1OfPjhh2Lnzp0ljrNmzRoxffp0IYQQJ0+eFEFBQaYvvgKV51o5fPhw8cMPPwghhFi+fLn46quvhBBCzJw5U6xcuVIIIcT27dtFSEiIeU6GzKas901WVpZo27atvm3VqlVi9uzZ+v06nU4EBwcLLy8vkZycrG/fvn276Nixo/Dx8THfyZBZnTt3Tvj7+4vXX3/d4GdfrH///uLw4cNCCCHmzZsnFi9eLIQQIioqSn+NOXjwoOjXr58QQohly5aJhQsXCiGE+OOPP0SHDh3KrMGqeojj4+PRtm1bODs7w8HBAd27d8eePXv0+2/fvo38/Hx4eXkBAAICAgz2U9XWunVr/Oc//4FCocC9e/eg0+lw4cIFeHl56Vfi8vHxwYEDBwAARUVFUKvVAIC8vDxUq1atxDEvXryIjIwM9O3bF3379sXJkyfNd0JU4b788kuEhobizp07T3xfFNu2bRvs7e2hUqlKHOfw4cPo1asXAOCtt97C/fv3cefOHdOfQAV53mtlYWEhTp06he7duxu0Aw+/J8XfK39/fxw9ehSFhYXmPTEyqbLeN0lJSahbty4aNWoEoOTv1Zo1a9C+fXvUqlVL3/bgwQMcPHgQixcvNt+JkNnFxMQgPDz8iatiFhUV6f+a8Ojn8ZPax44di08//RTAw7/21qxZs8warCoQp6amwtXVVb/t5uaGlJSUJ+53dXU12E9Vn1KpxNKlS+Hn54d27dqhefPmOH/+PP7++2/odDrs2bNHv2xtSEgI/ud//gedOnXCunXrMGzYsBLHk8lkeOedd7B582Z9mMrIyDD3aVEFiI+PR35+Pnr27AlPT88nvi8AQKfT4ZtvvsFnn31W6rFKu5bcvXvX5OdQUZ73Wnn//n04OjpCoVAYtD/+GoVCAUdHR/6uWJiy3jevvPIK7t69i8TERADAjz/+qP+9unTpEn799VcMGTLE4JhOTk5YtmwZ6tSpY4YzIKnMnTsXrVq1euL+L774AtOmTUPHjh0RHx+Pfv36AQCCg4Nx4sQJdOzYEWFhYRg/frz+NQqFAp988glGjRpV4n1VGqsKxKKUNUhkMpnR+8kyjB8/HidOnMDff/+NU6dO4bPPPsOoUaMwcOBAvPbaa1AqlQCAyZMnY9asWTh27BhmzpyJsWPHlniP9OvXD2PHjoVMJkOTJk3QvHlznD17VorTonLatGmT/qLZoEGDJ74vAODYsWNo0KABXnvtNaOPL5dXncvt814rn/UaWpW+J1S2sn7+NWrUwIIFCzB9+nQEBgbCzc0NSqUSeXl5mDVrFmbPns33BJWQn5+PadOm4bvvvsPx48cxYMAATJ48GQAwe/ZsDBw4EMePH8e6desQGhpqMC597dq12L9/P6Kjo8u8l8Oq3nnu7u4GvTypqakG3fOP709LS3ti9z1VPdevX9cPuLe3t0e3bt1w4cIFNG/eHHFxcdi0aRPq1q2L+vXrIyMjAzdu3MC7774LAOjevTvS0tJw//59g2PGxcXh5s2b+m0hhEFwoqpBo9Hg1KlT6NKlCwCgoKCg1PdFsQMHDsDX1/eJx3NzczO4aayqXUue91rp4uICtVqtv7H00fN2c3PTv0ar1UKtVsPZ2dkMZ0PmUtb7RqfT4cUXX8SWLVuwbds2NG3aFPXr18fp06eRnp6OUaNG4b333kNqaiqGDx+OGzduSHEaVMlcvXoVdnZ2+iXkP/jgA/3wxIMHDyIwMBAA0KJFC9SuXRvXr1/HyZMnkZqaCgCoV68eWrRogWvXrj3161hVIG7fvj1OnDiBjIwM5OXlYd++ffD29tbvr1evHuzs7HDmzBkAD8POo/upart16xbCwsKg0Wig0Whw8OBBtGnTBh999BHUajU0Gg02bNgAX19f1KpVC3Z2djh9+jQA4MyZM6hevTpcXFwMjnnlyhWsW7cOAHDjxg0kJCSgZcuWZj83Kp8rV67glVdegYODAwAgNze31PdFsXPnzj31z3udO3fGjh07AACnT5+GnZ0d6tata9qTqEDPe61UKpVo1aoVdu/ebdAOPPyexMXFAQB2796NVq1a8T+PFqas941MJkNwcDBSUlIghMC6devg6+uLTp064dChQ9ixYwd27NgBNzc3rFq1Cg0bNpTwbKiyePnll3H37l39f5AOHjyIZs2aAQA8PT3149CTkpKQmpqKBg0a4PDhw1i1ahWAh/8xu3Tpkv41T1SeuwKrop07dwo/Pz/RrVs3sWrVKiGEEEOHDhUXLlwQQgiRkJAgAgMDRY8ePcSECRNEQUGBlOVSBYuOjhY9e/YU/v7+YunSpUIIIWJiYoSvr6/o1q2bvk0IIc6fPy+CgoKEv7+/+OCDD8Tly5eFEEIcOHBATJ06VQghxIMHD8S4ceOEn5+f8Pf3FydOnDD/SVG5/fe//xWffvqpQduT3hdCCNG8eXORn59v0Pb999+LJUuWCCGEyM/PF59//rnw9fUVvXv3FpcuXTLtCZjA814rb926JT788EPRs2dPERwcLDIzM4UQQty/f1+MGDFC+Pr6ig8++KDUO8mp6ivrffPTTz8Jf39/0a1bNxEeHi40Gk2JY/j4+JR4fyQnJ3OWCSvw6M/+0ffN4cOHhUqlEv7+/uKjjz4SN2/eFEII8eeff4pBgwYJPz8/0adPH/Hzzz8LIR5+No8fP174+/uL9957T+zfv7/Mry0TopRBP0REREREVsKqhkwQERERET2OgZiIiIiIrBoDMRERERFZNQZiIiIiIrJqDMREREREZNUYiImIiIjIqjEQExEREZFVYyAmIiIiIqv2/wBiyifBPDmPPgAAAABJRU5ErkJggg==\n" }, "metadata": {} } ], "execution_count": 9, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:07.013Z", "iopub.execute_input": "2022-04-10T00:49:07.030Z", "iopub.status.idle": "2022-04-10T00:49:07.293Z", "shell.execute_reply": "2022-04-10T00:49:07.341Z" } } }, { "cell_type": "markdown", "source": [ "The number of eigenvalues near zero corresponds to the number of quality clusters, and we see a substantial number near zero in the histogram. (Strictly speaking, the function plots the spectrum of the giant component by default, so the plot shows the number of quality clusters in the giant.) The exercise here is similar to the use of screeplots to detect the number of principal components in PCA. \n", "\n", "The scatterplot is not as interesting, but we will see shortly below how it can be useful to see when quality clusters do *not* exist.\n", "\n", "To get a better sense of how many eigenvalues are near zero, we compute how many are below 0.05." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "eigenvals = ni.core.spectrum(A) # get vector of eigenvalues\n", "print(f'{(eigenvals < 0.05).sum()} eigenvalues are below 0.05')" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "38 eigenvalues are below 0.05\n" ] } ], "execution_count": 10, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:07.321Z", "iopub.execute_input": "2022-04-10T00:49:07.331Z", "iopub.status.idle": "2022-04-10T00:49:07.375Z", "shell.execute_reply": "2022-04-10T00:49:07.410Z" } } }, { "cell_type": "markdown", "source": [ "This suggests there could be as many as 38 \"quality\" clusters in the network. To actually construct the clusters, we can use spectral clustering by calling the following method, which stores the result in the `tsls` object for later use." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls.get_clusters(38, seed=0)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Maximal conductance: 0.12727272727272726\n" ] } ], "execution_count": 11, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:07.390Z", "iopub.execute_input": "2022-04-10T00:49:07.401Z", "iopub.status.idle": "2022-04-10T00:49:07.962Z", "shell.execute_reply": "2022-04-10T00:49:08.026Z" } } }, { "cell_type": "markdown", "source": [ "**Maximal conductance** is the key $[0,1]$-measure of cluster quality. [Leung (2022c)][6] shows that for cluster-robust inference methods to work, this value must be near zero. In practice, he recommends that the value not exceed 0.1. With the current set of 38 clusters, we're a little over that threshold, so maybe we can try asking for a few less." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls.get_clusters(30, seed=0)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Maximal conductance: 0.08771929824561403\n" ] } ], "execution_count": 12, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:07.994Z", "iopub.execute_input": "2022-04-10T00:49:08.010Z", "iopub.status.idle": "2022-04-10T00:49:08.504Z", "shell.execute_reply": "2022-04-10T00:49:08.554Z" } } }, { "cell_type": "markdown", "source": [ "We now have 30 low-conductance clusters we can use to conduct inference. \n", "\n", "Before using them, let us take a digression and look at a cautionary tale showing that not all networks have quality (low-conductance) clusters. " ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "### When not to cluster\n", "\n", "We draw a new network $G$ from an Erdos-Renyi model with the same average degree and size as $A$ and plot the spectrum." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "n, avg_degree = 1000, 5\n", "G = nx.fast_gnp_random_graph(n, avg_degree/n, seed=0)\n", "ni.core.plot_spectrum(G)\n", "\n", "eigenvals = ni.core.spectrum(G)\n", "print(f'{(eigenvals < 0.18).sum()} eigenvalues are below 0.18')" ], "outputs": [ { "output_type": "display_data", "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsQAAAEUCAYAAAAsgyAxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6z0lEQVR4nO3de1zTZf8/8NfGxikwxIBOVpoaWSqZeRbSbg8c5gFMSdPSFDwhgQdQUcuzhuKBtDT10W15G5riIU8ZCSp2a+YhDTJv85AlBxFxwBjbrt8f/tpX3MClbGPb6/l4+Hiw69o+vIfjsxfXrs91SYQQAkREREREDkpq7QKIiIiIiKyJgZiIiIiIHBoDMRERERE5NAZiIiIiInJoDMRERERE5NAYiImIiIjIoTEQk805deoUhgwZAoVCgbCwMIwYMQK//fbbAx3rzJkzmDFjhv728OHDUVRUVFulYsWKFZg1a9Z971fb35eI6GG88MILBuekrVu3Ijo6GgCwbNkypKen13iM1NRUHDhwwFwlEtUqmbULIPon1Go1oqOjsW7dOrz00ksAgO3bt2PkyJH47rvv4OTk9I+Od+HCBeTl5elvHzlypFbrNZW1vi8R0YOIjY29733++9//okmTJhaohujhMRCTTSkvL8ft27dRVlamb+vduzc8PDyg1Wqxbds2rF+/HlKpFPXr18fChQvh5+eHefPm4fTp0ygtLYUQAnPmzMGTTz6J5cuX4/bt25gyZYr+eO+88w5Wr14NqVSKWbNm4a+//kJlZSVCQ0MxatQo/PHHHxg8eDCef/55XLt2DQsWLEB8fDzatm2L3NxcCCEwY8YMtGnTpkrtv/32G2bNmoXi4mJIJBIMHz4cffv21X/vv7/vE088YZkfJhHRA0pMTETTpk3x3nvvYfny5fj2228hl8tRv359zJ8/H99++y3Onj2LRYsWwcnJCe3bt8eHH36I3NxcSCQSdOnSBfHx8ZDJZMjMzERycjKkUilefPFFZGdnY+PGjTh27Bi2bNmC8vJyeHh44NNPP8UHH3yAS5cu4datW3jkkUeQnJyMxo0bY8iQIXjppZfwww8/4MaNGxg6dChu3LiBY8eOoby8HEuXLsULL7xg7R8b1WWCyMasW7dOtGzZUnTr1k1MnDhRbN68WZSVlYmcnBzRrl078eeffwohhFi/fr2YPn26+Omnn0RMTIzQarVCCCE+/fRTER0dLYQQ4uuvvxZRUVH6Yzdr1kzcuHFDCCHEkCFDxHfffSeEEEKlUokhQ4aIb775Rly9elU0a9ZMHD9+XAgh9Ld37NghhBDi4MGDolOnTkKtVovly5eLDz/8UFRWVoo33nhD7Nu3TwghxPXr10WXLl3ETz/9ZPB9iYisrVmzZiIsLEz07t1b/y8oKEh/vkxISBCfffaZ+PPPP0Xr1q1FRUWFEEKItWvXim+//VYIIcTbb78t9uzZI4QQYvLkyWL27NlCp9OJiooKMXz4cPHpp5+KoqIi0bZtW5GTkyOEEGLr1q2iWbNm4urVq+Lrr78Wr732mrh9+7YQQog9e/aI2bNn62ucPn26mDVrlv57jRs3TgghxKlTp0SzZs305++5c+eKpKQkc//IyMZxhJhszrBhw/Dmm2/i+PHjOH78ONasWYM1a9ZAoVCgc+fO+hHWd999V/+YRx99FJs2bcLVq1fx3//+F4888kiN36OsrAzHjx/HrVu3sGzZMn1bbm4uWrZsCZlMhoCAgCrHVygUAICgoCA4OTnh119/1fdfunQJFRUV6NGjBwDAz88PPXr0wKFDh/DKK6/Uxo+FiKhWff755/D29tbf3rp1K/bt21flPn5+fvD390e/fv0QGBiIwMBAdOjQweBYWVlZ+M9//gOJRAJnZ2dERkbi888/R6NGjfD888/D398fANCvXz/MmTNH/7gXXngBHh4eAIBevXqhYcOG2LBhAy5fvoxjx45VOX92794dANCwYUMAQJcuXQAAzzzzDI4dO1YbPxKyYwzEZFNOnDiBkydPYsSIEejatSu6du2K+Ph4fRiVSCT6+6pUKly7dg1Xr17F3LlzMWzYMLzxxhto3LgxduzYUeP30el0EEJg06ZNcHNzAwAUFRXBxcUFN2/ehLOzM2Sy//v1uXfusk6nq9Km0+kMvocQAhqN5p//EIiI6gipVIovvvgCP//8M44ePYp58+ahXbt2SEpKqnK/e8+BOp0OGo0GTk5OEEIYHPNv7u7u+q83btyItLQ0DB48GAqFAl5eXvjjjz/0/c7OzlWOI5fLH/r5kePgKhNkU7y9vbFq1Sr8+OOP+raCggKUl5fjX//6F44ePYr8/HwAwKZNm/DRRx/hyJEj6Nq1KwYNGoQWLVrgwIED0Gq1AO4E2btD6d+3PTw8EBAQgPXr1wMASkpK8NZbb+G7774zWldRURGysrIAABkZGZDL5WjWrJm+v1GjRpDL5di/fz8AIC8vD/v27UPHjh2N1kFEZAtyc3MRFhaG559/HtHR0Xj33Xf1n47dfV7r3LkzvvzySwghoFarkZaWho4dO6J169a4dOkScnNzAQD79u1DSUlJlcGNvx0+fBj9+vXDm2++iUaNGiEjI0N/Lid6WBwhJpvSqFEjfPzxx0hJScH169fh4uICT09PzJo1C/7+/pg0aRJGjBgBAPDx8cG8efOgVCoxceJEKBQKODk5oU2bNti/fz90Oh1eeeUVLF26FGPHjsXHH3+M7t27Y9CgQVi5ciWSk5Mxe/ZsKBQKqNVqhIWFoXfv3lVGJP7m4uKC7du3Izk5Ga6urvj444+rjBDL5XKsXLkSc+bMwYoVK6DVajF27Fi0b98eAKp837uDNBFRXebv74/g4GBERETA3d0drq6u+tHhrl27YuHChaisrERSUhLmzJkDhUKByspKdOnSBaNGjYKzszOWLFmChIQESKVSvPzyy5DJZPpP5u42fPhwzJgxA1u3boWTkxNeeuklnD9/3tJPmeyURNz7WQUR/SN//PEHFAoFTp48ae1SiIhsilKpxMqVKxETEwM3NzecO3cO0dHROHTokNFRYiJz4QgxERERWYWHhwfkcjn69+8PmUwGmUyGpUuXMgyTxXGEmIiIiIgcGi+qIyIiIiKHxkBMRERERA6NgZiIiIiIHBoDMRERERE5NIdbZeLmzVLodLyOkIjMQyqVoH79mrcGt2c8xxKROZnrHOtwgVinEzxZExGZCc+xRGSLzDplQqlUIiwszOjOXufOnUNERAR69+6N6OholJSUVOnfsmULEhMT9bfVajUmTJgAhUKBPn36IDs725ylExEREZGDMFsgPn36NN566y1cunTJaP/cuXMxfvx47NixA40aNcLatWsBABUVFUhOTsbcuXOr3H/79u3Q6XTYuXMnFi1aVCUsExFZkkqtsXYJRERUi8wWiNPS0jBz5kz4+voa7dfpdCgtLQUAlJeXw9XVFQBw/Phx6HQ6TJo0yeD+5eXl0Gq1Ve5PRGRJq9LPYsySLKxKP2vtUkx276d12dnZUCgU6NGjB1JSUvT3y8nJQUREBHr27Ilp06ZBo2HwJyLHYLZAPHfuXLRp06ba/sTEREybNg2dO3dGdnY2IiMjAQCdO3fG5MmTDQJvv379UFxcjC5duuDtt9/GxIkTzVU6EZFRKrUGx3PzAQDHc/NtYqT43k/rVCoVpk6dipUrV2L37t04e/YsMjMzAQCTJk3C9OnTsW/fPgghkJaWZsXKiYgsxyrLrqlUKkybNg2ff/45Dh8+jEGDBiEhIaHGx6SmpiIgIABHjhzBzp07MXfuXFy7ds1CFRMRAa7OMrRu5gMAeM3fF67Odf+65Hs/rTtz5gyeffZZNGzYEDKZDAqFAnv37sW1a9egUqkQEBAAAAgPD8fevXutWDkRkeVYJRCfP38eLi4uaNmyJQBg4MCBOHbsWI2P+e677xAeHg6JRIJGjRqhVatWOHPmjCXKJSICcGe6xE/nC6xdxj9y76d1+fn58PHx0d/29fVFXl6eQbuPjw/y8vIsWisRkbVYJRA/++yzuH79Oi5evAjgTtht0aJFjY/x9/fHgQMHAABFRUU4e/YsXnzxRbPXSkQEVJ0uAdjOlIl7CWG4JJpEIqm2nYjIEVj0876RI0di/PjxaNGiBebPn4/3338fQgg0aNAA8+bNq/GxU6ZMwfTp0xEaGgqpVIr4+Hg899xzlimciBxeTEpWldsBTRrYxJSJe/n5+aGwsFB/Oz8/H76+vgbtBQUF1V4UTURkb8x+Ns/IyNB/vWbNGv3XQUFBCAoKqvZx4eHhCA8P199+7LHHsGrVKvMUSURUgxELMqC7py2q90tWqeVhtWrVCr///jsuX76Mp59+Grt27UJERASeeuopuLi44MSJE3j11VeRnp6OwMBAa5dL5LDqPeoGFzP+0V2p0UIuczLb8SvUGpTcKjfb8Wub7Q1vEBFZ0IiFhmHYVi6oM8bFxQULFixATEwMKioqEBQUhF69egEAkpOTkZSUhNLSUjRv3hxDhw61crVEjsvFWYZJyzLNdvyPYoPMfnxbYptndCIiC4hJOYh7dyFeMq4TvDxcrFPQQ7j707oOHTpgx44dBvfx9/fHli1bLFmWwzDnaJ+5R+LMPVJpayOJZJ8YiImI7qFSazBmSZZBe+tmPjYZhsn6zDnaZ+6ROEuMVBJZGwMxEdFdlm85jVMXbhi0t2ryGMaF17waDtkuc4+CmlOlRgcfH09rl1EncXSbTGWbv/1ERGYQuywLt8sNl1Jr3ewxjAtvaYWKyFJseRRULpPabO3mZsv/r2RZDMRERACGL8gw2r4yPtBmL6AjsgW2PMJty7VTVTzLE5FDK1ZWID71iEG7BMDaxG6WL4jIwZhzhNvcI7i2XDtVxUBMRA4rPvUwipVqg/aXG3sjfkCA5QsiIiKrYCAmIoc0fmkmlCqtQTunSBAROR6e9YnIoXCKBBER3YuBmIgcgkqtQeInR1FSVmnQ5+nmhGWcr0dE5LAYiInI7lW3tjAAPPqIHCkxXSxcERER1SUMxERk1+JWHMKtUsNRYQmAxTa6DTMREdUuBmIislvVXThXz12GpeMDrVARERHVRQzERGR38orKMGX1DwbtHq5SzBrRgaPCRERUBQMxEdmVqI++h0YrDNo5KkxERNVhICYiu3Al7zY+WH/caB8vnCMiopowEBORTVOpNRibkgVhOCjMC+eIiMgkDMREZLMWbzqJc5duGu3j9stERGQqBmIisjnV7TYH3BkV/pjbLxMR0T/AdwwisinxqYdRrFQb7Zsf1R5+3u4WroiIiGwdAzER2YRiZQVmfPaD0XWFufUyERE9DAZiIqrz4lYcxq1S46PCS3jRHBERPSQGYiKqs1RqDcalZEFnZAUJD1cnLH+fo8JERPTwGIiJqE5KSTuFny8WGe1r1aQBYvu3snBFRERkrxiIiahOKVZWYPLKbGiMDAt7ujlh4ehOXEGCiIhqFd9ViKjOqGkFCc4VJiIic2EgJiKrU6k1SFiVjdvlGoM+CYC1id0sXxQRETkMBmIisqrlW07j1IUbRvu42xwREVkCAzERWU3ssiyjo8KcK0xERJbEdxsisriatl5+9BE5UmK6WLgiIiJyZAzERGRR1W2yIZUAyWN54RwREVkeAzERWUReURmS1vwArZFNNuq5y7B0fKDliyIiIgIDMRFZQNRH30NjJAlLACzmcmpERGRlDMREZDYqtQZjlmQZ7fN0k2FZLEeFiYjI+hiIicgsFm86iXOXbhrt49bLRERUlzAQE1GtqmlU+OXG3hjT92Uup0ZERHUK35WIqNbUNCq8Mj6QQbgO2b59O1avXg0ACAwMREJCAnJycpCUlASlUok2bdrgww8/hEzG/zMisn9SaxdARLavWFmB4QsyjIZhTzcnrEvsxjBch5SXl2Pu3LnYsGEDtm/fjh9//BHZ2dmYNGkSpk+fjn379kEIgbS0NGuXSkRkEXyHIqKHUt26wgCwhCtI1ElarRY6nQ7l5eVwd3eHRqOBTCaDSqVCQEAAACA8PBzLly/HoEGDrFssEZEFMBAT0QNRqTUYl5IFnZF1hT3dnLAsNsjyRZFJPDw8EBsbi+DgYLi6uqJt27aQy+Xw8fHR38fHxwd5eXlWrJKIyHIYiInoH0tJO4WfLxYZ7eOocN2Xm5uLr7/+Gt9//z08PT0xceJEHDliuJW2RCKxQnVERJbHQExEJlOpNYhZmgWtzrCPo8K24/Dhw+jQoQMaNGgA4M70iLVr16KwsFB/n4KCAvj6+lqrRCIii+JFdUR0Xyq1Bos3ncSYJcbD8JJxnRiGbYi/vz+ys7NRVlYGIQQyMjLQtm1buLi44MSJEwCA9PR0BAZy4xQicgwcISaiGi3fchqnLtww2icBsDaxm2ULoofWuXNn/PLLLwgPD4dcLkeLFi0QFRWF7t27IykpCaWlpWjevDmGDh1q7VKJiCyCgZiIqhW7LAu3yzVG++ZHtYeft7uFK6LaEhUVhaioqCpt/v7+2LJli5UqIiKyHgZiIjJQrKxAfKrhRVYAR4WJiMj+MBATURXvLz+EkrJKg3YPVymmDW3LUWEiIrI7DMREBODOhXNjlmQZ7Xv0ETlSYrpYuCIiIiLLYCAmomovnJMAWMx1hYmIyM4xEBM5uOounKvnLsPS8Vx2i4iI7B8DMZGDyisqw5TVPxjt425zRETkSBiIiRxQ9EcHUWlkhw0PVymWv/+65QsiIiKyIgZiIgeiUmswdkkWhJE+XjhHRESOioGYyEGkpJ3CzxeLDNp54RwRETk6BmIiO6dSaxCzNAtGZkjgxee8MSkywOI1ERER1SUMxER2bPGmkzh36abRvtbNHsO48JYWroiIiKjuYSAmskM1bbLh4eqERWM6wdWZv/5EREQAAzGR3alpVLhVkwaI7d/KwhURERHVbQzERHaiplFhTzcnLBzNUWEiIiJj+O5IZAeq23oZ4CYbRERE98NATGTjqtt62dPNCctig6xQERERkW1hICayUTVNkeCoMBERkekYiIlsUHUXzkkArE3sZvmCiIiIbBgDMZENqWlU+OXG3ogfEGDZgoiIiOwAAzGRDVCpNViVftbo1ssAsDI+kCtIEBERPSC+gxLVcTWtIMEL54iIiB4eAzFRHaVSa5D4yVGUlFUa7eeFc0RERLWDgZioDqppVJgXzhEREdUuBmKiOqa6dYWlEmDuyPbw83a3QlVERESmq9To4OPjae0yTMZATFRH1LSCRD13GZaOD7RwRURERA9GLpNi0rLMWj9ufU8XTB3evtaPy0BMVAekpJ0yuoKEVAIkj+VcYSIiInOSmnKn//3vf9i8eTOEEBgzZgzeeOMN/PDDD+aujcjuqdQaRC363mgYrucuw2cJ3RiGiYiIzMykQDxz5ky4uLjg+++/x82bNzFv3jykpKSYuzYiu7Z8y2mMWZIFjU5UaZfgzgoSnCJBRERkGSYF4oqKCvTu3RtHjhxBcHAw2rVrh8pK40tBEdH9vb/8kNFVJKSSOytIcFSYiIjIckwKxGq1GoWFhTh48CA6duyIwsJCVFRUmLs2IruTV1SGkQszjK4t/HJjb3yWwOXUiIiILM2ki+oGDhyIrl27Ijg4GE2aNMHrr7+OMWPGmLs2IrsS/dFBVGp1Bu1SCZAax62XiYiIrMWkd+BBgwYhMjISUumdAeVt27ahfv36Zi2MyF4UKyswMfUIDKMw4Okmw7JYzhUmIiKyJpOmTJSWlmLOnDl45513UFxcjJSUFJSWlpq7NiKbF596GPHVhOFWTRowDBMREdUBJgXiOXPmwNPTEzdu3ICLiwuUSiVmzJhh7tqIbNr4pZkoVqoN2l98zhsr4wMR27+VFaoiIiKie5kUiHNychAXFweZTAY3NzckJycjJyfH3LUR2aRiZQWGL8iAUqU16Gvd7DFMigzgfGEiIqI6xKR35b/nDv9Nq9UatBERELfiMG6VGo4Ke7o5YeHoTgzCVGdkZGQgNTUVZWVl6Ny5M5KSkpCdnY358+ejoqICwcHBiIuLs3aZREQWYVKqfe211/DRRx9BpVLh0KFDiImJQbt27cxdG5HNKFZWYOTCDKNh+NFH5FgWG8QwTHXG1atXMXPmTKxcuRI7d+7EL7/8gszMTEydOhUrV67E7t27cfbsWWRmZlq7VCIiizApEE+cOBHu7u7w9PRESkoKXnjhBUyePNnctRHZhLgVdy6c01bdcE6/41xKTBer1EVUnW+//RYhISF4/PHHIZfLkZKSAjc3Nzz77LNo2LAhZDIZFAoF9u7da+1SiYgswqQhK7lcjrFjx2Ls2LHmrofIpoxYmIF7dl4GwOXUqG67fPky5HI53nvvPRQUFKBr165o2rQpfHx89Pfx9fVFXl6eFaskIrIckwLxqFGjjLZ/8skntVoMka3IKyrDlNU/GO1r1aQBV5CgOk2r1eLHH3/Ehg0b4O7ujjFjxsDNzc3gfhKJxArVERFZnkmBuGfPnvqvKysrkZGRgRdeeMFsRRHVVSq1BuOXHYLm3vkRuDMqvHB0R84VpjrvscceQ4cOHeDt7Q0AeOONN7B37144OTnp75Ofnw9fX19rlUhEZFEmvXP369evyu3w8HC8/fbbZimIqK5avOkkzl26abSPo8JkS7p27YqEhASUlJTgkUcewaFDh9CrVy+sXr0aly9fxtNPP41du3YhIiLC2qUSEVnEAw1l6XQ65Ofn13YtRHWSSq3BmCVZRvskAD6OD+SoMNmUVq1aYcSIERg0aBAqKyvRqVMnvPXWW2jcuDFiYmJQUVGBoKAg9OrVy9qlEhFZxAPNIT5//jzatm1rloKI6pLlW07j1IUbRvtebuyN+AEBli2IHNrUqVMxb968Km0xMTFYsWLFPz5W//790b9//yptHTp0wI4dOx6qRiIiW/SP5xBLJBK89dZb6Ny5s9mKIqoLYpdl4Xa5xqCdm2yQpc2cORN5eXk4ceIEioqK9O0ajQYXL160YmVERPahxnf04uJiAHfmm93r9u3b8PLyMkdNRFZVrKxAfOoRo31LxnWCl4eLhSsiR9e/f3/89ttv+PXXX6sMUDg5OeGVV16xYmVERPahxkDcvn17/bI7Qty5ql4ikUAIAYlEgpycHPNXSGQhKrUGiZ8cRUlZpUGfBMDaxG6WL4oIQIsWLdCiRQt07NgRjz/+uLXLISKyOzUG4tzcXEvVQWRVnCtMtuDKlSuYNGkSbt26pR+kAICdO3dasSoiIttn0iRItVqNzMxMlJaWArizqPuVK1cQFxdn1uKIzE2l1iBhVbbRucIAsJIrSFAdMmvWLERERKB58+bcNIOIqBaZ9E4fFxeHq1evoqCgAM2bN8fp06e5ygTZvJpGhV98rj4mRXJuJtUtcrkcw4YNs3YZRER2x6RAnJOTg/379+ODDz7AsGHDIITAhx9+aO7aiMzm/eWHqp0rzHWFqa5q2rQpfv31V+4USkRUy0x61/f19YVMJsNzzz2H8+fPIzg4GOXl5eaujajWqdQajEvJgs5w52XOFaY67+rVq4iIiMCTTz4JF5f/W+2Ec4iJiB6OSYHY3d0dO3fuhL+/P9LS0tC4cWP9kmxEtiIl7RR+vlhk0M5RYbIVvG6DiMg8pKbcacaMGcjJyUGnTp0glUoxZMgQvPfee/d93M6dOxESEoLu3bvjyy+/NOjPzMyEQqGAQqHAhAkT9BftlZSUICoqCsHBwRg8eDAKCgoAAIWFhRg1ahTCwsIwcOBAnDx58p88V3JQKrUGIxdlGA3D9dxlWJvYjWGYbEKzZs2M/iMioodjUgq4evUqJk+eDABYunSpSQfOy8tDSkoKtm7dCmdnZ0RGRqJdu3Zo0qQJgDuhNzExERs2bECTJk2wZs0apKSkICkpCUuXLkWbNm2wevVqpKenY+7cuVi6dCkWLFiA5s2b45NPPsHVq1cxbNgw7Nq1C66urg/27MnuVTcqDHCTDbI9f68N//da8ADg4+ODrKwsK1dGRGTbTBohTk1NRbdu3ZCamoq8vDyTDpydnY327dvDy8sL7u7u6NmzJ/bu3avvv3TpEp588kl9QO7atSsOHDgAADh48CAUCgUAICwsDFlZWaisrEROTg6Cg4MBAA0bNoSXlxdHickolVqDqEXfVztFYl1iN4Zhsjm5ubnIyclBbm4uzpw5g+TkZPTp08faZRER2TyTAvFXX32FNWvWoLy8HAMGDEB0dLQ+vFYnPz8fPj4++tu+vr5VwvRzzz2H69ev6zf/2LNnDwoLCw0eK5PJ4OHhgaKiIjRv3hzffPMNAOD8+fO4cOGC/jFEf0tJO4UxS7KgMXLl3Pyo9txxjuyCXC5HaGgojhwxvs04ERGZzqRADADPP/88Jk2ahBUrVuDmzZuIj4+v8f5376L0t7sXkq9Xrx4WLlyI6dOnIyIiAr6+vpDL5dUXKpViypQpuHz5MhQKBf7973+jXbt2NT6GHIspo8J+3u6WL4yolhQXF+v/3bx5E4cOHUJJSYm1yyIisnkmzSG+ceMGduzYgW3btkGr1aJ///749NNPa3yMn58ffvzxR/3t/Px8+Pr66m9rtVo8/vjj2Lx5MwDg3LlzaNiwIYA7o8mFhYV4/PHHodFooFQq4eXlhby8PMyePRseHh4AAIVCgWeeeeafPWOySzXNFeZyamQv7p5DDAANGjTAtGnTrFwVEZHtMykQ9+jRAz169MDMmTPx6quvmnTgjh07YsWKFSgqKoKbmxv279+P2bNn6/slEgmGDx+OzZs3w9fXF+vWrUNISAgAICgoCOnp6Rg1ahR2796NNm3aQC6X44svvsBjjz2GESNG4PDhw6isrIS/v/8DPG2yJyMXZUCrM2yXSoDUOC6nRvbj7ylmRERUu0xKCpmZmfpRWVP5+fkhLi4OQ4cORWVlJfr374+WLVti5MiRGD9+PFq0aIFZs2ZhxIgRUKvV6NChg34pt9jYWCQmJiI0NBSenp5ITk4GAERFRWHChAnYvn07HnnkEaSmpkIqNXnWB9mZYmUF4lONz5/kqDDZI51Oh7Vr1yIrKwsajQadOnXCqFGjIJPxjz4ioodh0ln0l19+wYoVK3Dr1q0qc4PvtzvS32sM323NmjX6r19//XW8/vrrBo/z8vLCJ598YtDu7e2N9evXm1Iy2bnqtl7mqDDZs8WLFyM3NxfvvPMOdDodvvrqKyxatAhTp061dmlERDbNpNQwa9YsREREoHnz5lUujCOyNI4KkyM7dOgQvv76a/3FxK+//jp69+7NQExE9JBMCsRyuRzDhg0zdy1ENYpPPYJiZYVBO7deJkchhKiyso6zszNX2iEiqgUmTcBt2rQpfv31V3PXQlStRRt/MhqGufUyORJ/f3/MmzcPV65cwZUrVzBv3jxu3UxEVAtM3ro5IiICTz75JFxc/m93r/vNISZ6WMXKCkxemW2wyYYEwGJuvUwOZubMmZgzZw4iIyOh0+nQpUsXTJ8+3dplWUS9R93gwj98ichMTDq7xMXFmbsOIgPVXTjn4eqE5e8HWaEiIutQq9WYPn06unfvjgULFgC4s+qOk5PTP14ByFa5OMswaVmm2Y7/USzPKUSOzKQpE23btoWrqysuXryIgIAAyOVytG3b1ty1kYMqVlZg+IIMo2HYSQqGYXI4y5cvh1KpxCuvvKJvmz17NkpKSrBixQorVkZEZB9MCsRbt27FlClT8Nlnn+H27dsYM2YM0tLSzF0bOaD3lx+qcRWJNZO7WbgiIus7ePAgFi9ejAYNGujb/Pz8sGjRIhw4cMCKlRER2QeTAvGGDRvw1VdfwcPDAw0aNMDWrVvx+eefm7s2ciAqtabaUeF67jKsjA/kkmrksORyOVxdXQ3aPTw84OzsbIWKiIjsi0lziKVSaZV5ak888QScnJzMVhQ5lsWbTuLcpZsG7bxwjugOqVQKpVJpMF9YqVRCo9FYqSoiIvth0gixl5cXcnJy9Jty7NixA48++qhZCyP79/eosLEw/PdyagzDREBYWBiSkpJQVlambysrK0NSUhJ69OhhxcqIiOyDSSPEU6dORWxsLK5cuYLOnTvDxcUFK1euNHdtZMdS0k7h54tFRvuWcFSYqIp33nkHM2fORKdOndC0aVPodDr873//g0KhwNixY61dHhGRzTMpED///PPYvn07Ll26BK1Wi0aNGnF3JHogKrUGMUuzoNUZ9nm6OWEZlz4iMiCVSjF79mxER0fjl19+gVQqRYsWLeDn52ft0oiI7IJJgTg1NbXKbYlEAjc3NzRt2hRdunQxS2Fkf1K3/oyfzhcY7eOoMNH9Pf3003j66aetXQYRkd0xKRCfP38eJ0+eRM+ePeHk5IRvv/0WTz31FPbs2YMzZ87wIzu6r+RNJ/GLkbnC3GSDiIiIrM2kQHzjxg1s3boVPj4+AIBRo0YhNjYWX375JSIiIhiIqVo1TZFo1aQBYvu3snxRRERERHcxKRAXFxfrwzAA1K9fH8XFxXB2doZMxr3lybiallP7OD4Qrs587RAREZH1mZRIGjZsiMWLF2PAgAEAgC1btuCZZ57B6dOnIZWatHIbOZBiZUW1u815usmwLDbQwhURERERVc+kNDtv3jxcu3YN/fr1Q//+/ZGXl4c5c+bg3LlzSEhIMHeNZCNUag3iVhyuNgy3atKAYZiojlm4cCESExMBADk5OYiIiEDPnj0xbdo0bvpBRA7DpBFib29vLFmyxKB90KBBtV4Q2ablW07j1IUbRvs83ZywcHQnTpEgqmOOHj2Kbdu24fXXXwcATJo0CXPmzEFAQACmTp2KtLQ0nueJyCHUmFBiY2OxbNkyKBQKo/07d+40S1FkW+JWHMatUrXRPi6nRlQ3FRcXIyUlBaNGjUJubi6uXbsGlUqFgIAAAEB4eDiWL1/OQExEDqHGQDxy5EgAwPTp0y1SDNkWlVqDpWmnjYZhbrJBVLfNmDEDcXFx+OuvvwAA+fn5VS6e9vHxQV5enrXKIyKyqBoDsbe3NwCgbdu2Bn1ZWVnmqYhsQnVbL0slQPJYjgoT1WWbN2/GE088gQ4dOmDr1q0AACGEwf0kEomlSyMisooaA/HYsWOxbds2AEBMTAxWrFih70tJSUFgIC+QckQjF2Vw62UiG7Z7924UFBSgT58+uHXrFsrKyiCRSFBYWKi/T0FBAXx9fa1YJRGR5dQYiO8eMbh69Wq1feQYVGoNxiwx/snAo4/IkRLDbbyJbMH69ev1X2/duhXHjh3D/PnzERYWhhMnTuDVV19Feno6Bz2IyGHUGIjv/rjs3o/O+FGa41CpNViVftboFAkJgMW8cI7ILiQnJyMpKQmlpaVo3rw5hg4dau2SiIgswuQRYnJMNS2nNj+qPfy83S1cERHVpvDwcISHhwMA/P39sWXLFitXRERkeTUGYp1Oh1u3bkEIAa1Wq/8aALRarUUKJOtQqTVI/OQoSsoqjfavS+xm4YqIiIiIzKPGQHz+/Hm0b99eH4LbtWun7+OUCftV06iwBMBahmEiIiKyIzUG4tzcXEvVQXWASq1Bwqps3C433K5VKgHmjuQUCSIiIrI/3EuXAFS/rjAA1HOXYel4Xm1ORERE9omB2MGp1BrELM0yuq4wN9kgIiIiR8BA7MA4KkxERETEQOyQahoV9nB1wqwR7TkqTERERA6DgdjB1LSCRKsmDRDbv5WFKyIiIiKyLgZiBxK7LMvoChIerk5YNKYTXJ35ciAiIiLHwwTkAIqVFYhPPWK0j6PCRERE5OgYiO1YTbvNSQB8HB/IUWEiIiJyeExDdqqmucIvN/ZG/IAAyxZEREREVEcxENuhuBWHcKvUcFQYAFZyVJiIiIioCiYjO6JSazB55REoVVqDvhefq49Jka9YoSoiIiKiuo2B2E5UN0WCc4WJiIiIasaUZONUag0SVmUbXU7N002GZbHcbY6IiIioJgzENqymrZe5nBoRERGRaRiIbZBKrcH4pYeg0QmDPk6RICIiIvpnmJpsTE2jwlxOjYiIiOifYyC2ERwVJiIiIjIPJigbwE02iIiIiMyHgbiOi12WZXQFCakESI3jqDARERHRw2KaqqOKlRWITz1itI+jwkRERES1h4G4Dnp/+SGUlBluvcxRYSIiIqLax2RVh9Q0KlzPXYal47nJBhEREVFtYyCuI+JWHMatUrVBuwTA4nGd4OXhYvmiiIiIiBwAA7GVqdQajEvJgpHV1DgqTERERGQBDMRWtHjTSZy7dNOgnaPCRERERJbDQGwFNc0V9nSTYVksR4WJiIiILIWB2MKqmysMAK2aNEBs/1YWroiIiIjIsTEQW0hNc4U93ZywcHQnLqdGREREZAVMYBZQ09bLSzhXmIiIiMiqGIjNrLqtl194pj4SBr1ihYqIiIiI6G4MxGZS04VzrZs9hnHhLS1cEREREREZw0BsBtVtvSwB8HE8t14mIiIiqkuk1i7AnqjUGgxfkGE0DL/c2BtrE7sxDBNRnZCamorQ0FCEhoZi0aJFAIDs7GwoFAr06NEDKSkpVq6QiMhymM5qSXWbbADASo4KE1Edkp2djcOHD2Pbtm2QSCQYMWIEdu3aheTkZGzYsAFPPPEEoqOjkZmZiaCgIGuXS0RkdhwhrgXDF2QYDcMvPlcf6zgqTER1jI+PDxITE+Hs7Ay5XI7nn38ely5dwrPPPouGDRtCJpNBoVBg79691i6ViMgimNQeQk0XznFUmIjqqqZNm+q/vnTpEnbv3o0hQ4bAx8dH3+7r64u8vDxrlEdEZHFMbA+ouh3nPN2csCyWHzESUd3322+/ITo6GgkJCZDJZPj999+r9EskEitVRkRkWQzED2Dkwgxojew4x002iMhWnDhxAuPHj8fUqVMRGhqKY8eOobCwUN+fn58PX19fK1ZIRGQ5nEP8DxQrKzBigfEwvC6xG8MwEdmEv/76C2PHjkVycjJCQ0MBAK1atcLvv/+Oy5cvQ6vVYteuXQgMDLRypURElsERYhNVN0WCO84Rka1Zu3YtKioqsGDBAn1bZGQkFixYgJiYGFRUVCAoKAi9evWyYpVERJbDQGyCkYsyoNUZtrdq8hhi+3PHOSKyLUlJSUhKSjLat2PHDgtXQ0RkfQzE9zF+aabRMMztl4mIiIjsAwNxDUYszIDunvnCnm4yLBzdkUuqEREREdkJprpqDF+QYdDW/Ln6mBjJ+cJERERE9oSrTBjxnpEw7OXhzDBMREREZIc4QnyPEQsycO+qalxfmIiIiMh+MRDfZcTCDNx7/Ry3YCYiIiKyb5wy8f8tSTtlcAHda/6+DMNEREREdo5pD4BKrcHZi0VV2jhNgoiIiMgxcITYiNbNfBiGiYiIiBwEAzEAV2cZ6nveCcBeHs4YF97CyhURERERkaUwEOPOlImbtysAAMVKNVRqjZUrIiIiIiJLYSDGnRHi1/x9AfBCOiIiIiJHw+T3/43u+zKGqTUMw0REREQOhiPEd2EYJiIiInI8DMRERERE5NAYiImIiIjIoTEQExEREZFDYyAmIiIiIofGQExEREREDo2BmIiIiIgcGgMxERERETk0BmIiIiIicmgMxERERETk0BiIiYiIiMihMRATERERkUNjICYiIiIih8ZATEREREQOjYGYiIiIiByaWQPxzp07ERISgu7du+PLL7806M/JyUFERAR69uyJadOmQaPRAAD+/PNPDB48GL169cLo0aNRWloKACgpKUFUVBSCg4MxePBgFBQUmLN8IiIiInIAZgvEeXl5SElJwcaNG7F9+3Z89dVXuHDhQpX7TJo0CdOnT8e+ffsghEBaWhoA4MMPP8SgQYOwd+9evPzyy1i5ciUAYOnSpWjTpg327NmDN998E3PnzjVX+URERETkIMwWiLOzs9G+fXt4eXnB3d0dPXv2xN69e/X9165dg0qlQkBAAAAgPDwce/fuRWVlJY4fP46ePXtWaQeAgwcPQqFQAADCwsKQlZWFysrKWqtZpdbU2rGIiIiIyDaYLRDn5+fDx8dHf9vX1xd5eXnV9vv4+CAvLw83b96Eh4cHZDJZlfZ7HyOTyeDh4YGioqJaqXdV+lmMWZKFVelna+V4RERERGQbzBaIhRAGbRKJ5L7993vcvaTSh38KKrUGx3PzAQDHc/M5UkxERETkQMwWiP38/FBYWKi/nZ+fD19f32r7CwoK4OvrC29vbyiVSmi12irtwJ1R5r8fo9FooFQq4eXl9dC1ujrL8Jr/ne/xmr8vXJ1lD31MIiIiIrINZgvEHTt2xNGjR1FUVITy8nLs378fgYGB+v6nnnoKLi4uOHHiBAAgPT0dgYGBkMvlaNOmDXbv3l2lHQCCgoKQnp4OANi9ezfatGkDuVxeK/WO7vsyVsYHYnTfl2vleERERERkG8w6QhwXF4ehQ4eib9++CAsLQ8uWLTFy5Ej8/PPPAIDk5GTMnz8fwcHBKC8vx9ChQwEAM2fORFpaGkJCQvDjjz/i/fffBwDExsbi1KlTCA0NxcaNGzFjxoxarZkjw0RE918yk4jI3pg1ASoUCv2qEH9bs2aN/mt/f39s2bLF4HFPPfUUNmzYYNDu5eWFTz75pPYLJSIiAP+3ZObWrVvh7OyMyMhItGvXDk2aNDHp8Q0aeJi5QiKi2udwQ6JSafUX6BERPSxbP8fcvWQmAP2SmePGjTPp8R+nncQtZUWt1zV1eHvU93Sp9ePezZaPb8u1m/v4rN2+jv+oh3lqlghjyzoQEZFD+vTTT1FWVoa4uDgAwObNm3HmzBnMnj3bypUREZmPWbduJiIi2/JPl74kIrIHDMRERKR3vyUziYjsEQMxERHp3W/JTCIie+RwF9UREVH17l4ys7KyEv3790fLli2tXRYRkVnxojoiIiIicmicMkFEREREDo2BmIiIiIgcGgMxERERETk0BmIiIiIicmgOF4h37tyJkJAQdO/eHV9++aVBf05ODiIiItCzZ09MmzYNGo3GClWSuS1cuBCJiYkAgMzMTCgUCigUCkyYMAGlpaUAAKVSiQkTJqBv377o27cvzp07Z3CcyspKtG7dGn369NH/02q1Fn0u9PBWr16Nnj17QqFQYNWqVQCqf11cuHABkZGR6N27N4YMGYJr164ZHE8IgYULF6JXr14ICQnBiRMnLPp8asODniv//PNPDB48GL169cLo0aP1P7eSkhJERUUhODgYgwcPRkFBgUWfD1nG/V431f1eVff6KCwsxKhRoxAWFoaBAwfi5MmTFn0+ZDlKpRJhYWH4448/DPrOnTuHiIgI9O7dG9HR0SgpKanSv2XLFv17OgCo1WpMmDABCoUCffr0QXZ29v0LEA7k+vXromvXruLmzZuitLRUKBQK8dtvv1W5T2hoqDh58qQQQogpU6aIL7/80gqVkjllZ2eLdu3aiYSEBHHr1i3Rvn17/etg9erVYvbs2UIIIaZOnSo++ugjIYQQmZmZon///gbH+vnnn8Xw4cMtVzzVuiNHjoiwsDBx+/ZtodFoRHR0tNi3b1+1r4u3335bZGZmCiGE2Lhxo4iPjzc45p49e8TIkSOFVqsVFy9eFP/6179EZWWl5Z7UQ3qYc2VUVJTYtWuXEEKI1NRUsWjRIiGEEB9++KH49NNPhRBCbNu2TcTGxlrmyZDF3O91U9P5trrXx4QJE8SyZcuEEEJcuXJFvPHGG6K8vNyCz4os4dSpUyIsLEy89NJL4urVqwb9b731ljh48KAQQoj58+eLJUuWCCGEUKlU4qOPPhIBAQEiISFBf/+0tDTx/vvvCyGEyM3NFV26dLlvDQ41QpydnY327dvDy8sL7u7u6NmzJ/bu3avvv3btGlQqFQICAgAA4eHhVfrJ9hUXFyMlJQWjRo0CAFy6dAlPPvkkmjRpAgDo2rUrDhw4ACEE9u/fj6ioKABAYGAg5s2bZ3C8n3/+GUVFRRgwYAAGDBiAY8eOWe7JUK345Zdf0LlzZ3h4eMDJyQldunRBWlqa0dcFAKxfvx6BgYHQ6XT4888/Ua9ePYNjZmZmIiQkBFKpFI0aNcKTTz5pUyNbD3qurKysxPHjx9GzZ88q7QBw8OBBKBQKAEBYWBiysrJQWVlp2SdGZnW/101151ug+tdHTk4OgoODAQANGzaEl5eXTf0ukWnS0tIwc+bManfF1Ol0+k8TysvL4erqCgA4fvw4dDodJk2aZHD/8vJyaLXaKveviUMF4vz8fPj4+Ohv+/r6Ii8vr9p+Hx+fKv1k+2bMmIG4uDh9iHnuuedw/fp15ObmAgD27NmDwsJC3LhxA87Ozvjiiy/Qt29fDB061OhUCIlEgjfeeANfffUVPvjgA8TFxaGoqMiiz4kezksvvYTDhw+juLgYFRUVyMjIgE6nM/q6AACZTIaSkhIEBgbiP//5DwYMGGBwzHu3O/bx8cH169ct84RqwYOeK2/evAkPDw/IZLIq7fc+RiaTwcPDg78rduZ+r5vqzrf3Pvbu10fz5s3xzTffAADOnz+PCxcuVNlanOzD3Llz0aZNm2r7ExMTMW3aNHTu3BnZ2dmIjIwEAHTu3BmTJ082CLz9+vVDcXExunTpgrfffhsTJ068bw0OFYiFkT1IJBKJyf1k2zZv3ownnngCHTp00LfVq1cPCxcuxPTp0xEREQFfX1/I5XJotVoUFhbi0UcfRXp6OqKjozF27FiDY0ZGRmLcuHGQSCRo3rw5WrZsiZ9++smST4seUocOHRAeHo4hQ4ZgxIgRePXVVyGXy42+Lv5Wr149HD58GEuWLMHo0aMN/lgydi6RSm3ndPug58p/eg61pZ8J3d/9/v+rO99WRyqVYsqUKbh8+TIUCgX+/e9/o127djU+huyPSqXCtGnT8Pnnn+Pw4cMYNGgQEhISanxMamoqAgICcOTIEezcuRNz5841er3H3Rxq62Y/Pz/8+OOP+tv3juL4+flV+cuzoKCg2uF7sj27d+9GQUEB+vTpg1u3bqGsrAzz5s3DgAEDsHnzZgB3Ju43bNgQ9evXh0wmQ1hYGACgU6dOKCsrw40bN9CgQQP9MdPT09G6dWs888wzAO68IfBkbVuUSiW6d++OYcOGAbgzJeLpp5/G448/bvC6AO68joKDgyGRSBAYGAiVSoVbt27B29tbf0w/P78qF43Z2rnkQc+V3t7eUCqV0Gq1cHJyqvK8fX19UVhYiMcffxwajQZKpRJeXl4We05kfvd73Wi12mp/r6p7feTl5WH27Nnw8PAAACgUCv35lhzD+fPn4eLiot9CfuDAgVi2bFmNj/nuu++QkpICiUSCRo0aoVWrVjhz5gyeeuqpah/jUH+ed+zYEUePHkVRURHKy8uxf/9+BAYG6vufeuopuLi46K8IT09Pr9JPtm39+vXYtWsXtm/fjvHjx6Nbt25ITEzE8OHDkZeXByEE1q1bh5CQEDg7O6Njx476j+pOnToFNzc31K9fv8oxf/31V6xbtw4AcPHiReTk5ODVV1+1+HOjB/fHH39g7Nix0Gg0uH37NjZv3oyQkBCjrwsAWLduHb799lsAwA8//ID69etXCcPAnTnnO3fuhFarxeXLl3Hp0iW0aNHC4s/tQT3ouVIul6NNmzbYvXt3lXYACAoKQnp6OoA7f1S0adOGfzzamfu9biQSSbW/V9W9Pr744gts2rQJAHD48GFUVlbC39/f4s+NrOfZZ5/F9evXcfHiRQB3wu79zqf+/v76+elFRUU4e/YsXnzxxZq/0YNeEWirduzYIUJDQ0WPHj3E6tWrhRBCjBgxQpw5c0YIIUROTo6IiIgQvXr1EvHx8aKiosKa5ZKZfP311/orUr///nsRFhYmevToIWbOnCnUarUQQoi8vDwRHR0tQkNDRZ8+fcSpU6eEEEIcOHBATJ06VQghxO3bt0VMTIwIDQ0VYWFh4ujRo9Z5QvRQUlNTRXBwsOjRo4fYuHGjEKL618Vvv/0mIiMjRe/evcXgwYPF+fPnhRBVXxc6nU4sWLBAhISEiJCQEHHo0CHrPLGH8KDnyj/++EO8/fbbIjg4WAwfPlwUFxcLIYS4efOmiI6OFiEhIWLgwIFGryQn23e/1011v1fVvT5u3Lgh3n33XREWFiYGDhxosNoJ2ZeuXbvq/+/vft0cPHhQKBQKERYWJt555x1x5cqVKo+7+z1dCCEKCgrEqFGjREhIiAgLCxM7d+687/eWCGFk0g8RERERkYNwqCkTRERERET3YiAmIiIiIofGQExEREREDo2BmIiIiIgcGgMxERERETk0BmIiIiIicmgMxERERETk0BiIiYiIiMih/T9bI+X+vkz4mAAAAABJRU5ErkJggg==\n" }, "metadata": {} }, { "output_type": "stream", "name": "stdout", "text": [ "1 eigenvalues are below 0.18\n" ] } ], "execution_count": 13, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:08.518Z", "iopub.execute_input": "2022-04-10T00:49:08.527Z", "iopub.status.idle": "2022-04-10T00:49:08.920Z", "shell.execute_reply": "2022-04-10T00:49:08.967Z" } } }, { "cell_type": "markdown", "source": [ "The histogram shows only one eigenvalue near zero, while the scatterplot reveals a \"spectral gap\" between the first and second eigenvalues. Note that the first eigenvalue is always zero and uninteresting. The next largest eigenvalue is around 0.19, which is rather far from zero. This suggests there are no good clusters in this network. Indeed, if we implement spectral clustering and ask for a just 2 clusters, we already obtain a large maximal conductance. Asking for more will only result in a worse value." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "for num_clusters in [2,6]:\n", " clusters = ni.core.spectral_clustering(num_clusters, G, seed=0)\n", " # alternate way of implementing spectral clustering without \n", " # having to instantiate TSLS object\n", " print(f'Maximal conductance with {num_clusters} clusters:\\\n", " {ni.core.conductance(clusters, G)}')" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Maximal conductance with 2 clusters: 0.34775888717156106\n", "Maximal conductance with 6 clusters: 0.5071283095723014\n" ] } ], "execution_count": 14, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:08.945Z", "iopub.execute_input": "2022-04-10T00:49:08.956Z", "iopub.status.idle": "2022-04-10T00:49:09.497Z", "shell.execute_reply": "2022-04-10T00:49:09.612Z" } } }, { "cell_type": "markdown", "source": [ "Due to the lack of quality clusters, using cluster-robust inference methods with the network $G$ can result in severe size distortion [(Leung 2022c)][6], so an alternative method should be used instead.\n", "\n", "**Summary:** We can use `ni.core.plot_spectrum()` to assess the number of quality clusters in the network and then construct the clusters using the method `get_clusters()` of the `tsls` object." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "### Clustered SEs\n", "\n", "Returning to our original network $A$ with 30 quality clusters, we have two options for constructing confidence intervals. The first is conventional clustered standard errors. " ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls.cluster_se()" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " Estimate SE 95% CI\n", "---------- ----- --------------\n", " 0.934 0.078 [0.781, 1.087]\n", " 0.522 0.029 [0.466, 0.579]\n", " 2.851 0.136 [2.585, 3.118]\n", " 1.007 0.043 [0.923, 1.091]\n" ] } ], "execution_count": 15, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:09.519Z", "iopub.execute_input": "2022-04-10T00:49:09.533Z", "iopub.status.idle": "2022-04-10T00:49:09.563Z", "shell.execute_reply": "2022-04-10T00:49:09.618Z" } } }, { "cell_type": "markdown", "source": [ "Note that, after calling `get_clusters()`, the result is stored in the `tsls` object, and that result is subsequently used by the `cluster_se()` method." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "### t-Statistic based test\n", "\n", "What if we only have a very small number of quality clusters, say 8? Then we have two alternatives. The simplest is the t-statistic based procedure due to [Ibragimov and Mueller (2010)][2]." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls.get_clusters(8)\n", "tsls.trobust_ci()" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Maximal conductance: 0.0374331550802139\n", " Estimate CI\n", "---------- --------------\n", " 0.934 [0.721, 1.179]\n", " 0.522 [0.355, 0.628]\n", " 2.851 [2.386, 3.22]\n", " 1.007 [0.85, 1.14]\n" ] } ], "execution_count": 16, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:09.581Z", "iopub.execute_input": "2022-04-10T00:49:09.591Z", "iopub.status.idle": "2022-04-10T00:49:09.775Z", "shell.execute_reply": "2022-04-10T00:49:09.809Z" } } }, { "cell_type": "markdown", "source": [ "### Approximate randomization test\n", "\n", "A potentially more powerful but slightly more computationally demanding alternative is to use approximate randomization tests [(Canay et al., 2017)][1]." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls.arand_ci(-4, 4)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " Estimate CI\n", "---------- --------------\n", " 0.934 [0.747, 1.173]\n", " 0.522 [0.373, 0.64]\n", " 2.851 [2.347, 3.2]\n", " 1.007 [0.853, 1.12]\n" ] } ], "execution_count": 17, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:09.787Z", "iopub.execute_input": "2022-04-10T00:49:09.793Z", "iopub.status.idle": "2022-04-10T00:49:11.985Z", "shell.execute_reply": "2022-04-10T00:49:12.019Z" } } }, { "cell_type": "markdown", "source": [ "These CIs are a bit narrower, indicating higher power. They are more computationally demanding because they are constructed by inverting a hypothesis test. The parameters -4 and 4 are the left and right endpoints of a grid, and a parameter `grid_size` controls how many grid points to potentially evaluate. \n", "\n", "Both the t-statistic and approximate randomization tests work by computing $L$ TSLS estimates separately across the $L$ clusters, so when using more clusters, the results can sometimes be more unstable since there's less data within each cluster. Here, though, we're only asking for 8 clusters, and the results are similar to clustered SEs." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "## Dependence-robust test\n", "\n", "The previous CIs use network data to account for network dependence, but in some situations, we may have insufficient knowledge of the dependence structure. The network data may be low quality with many missing or mismeasured links, or the analyst may be unsure about the type of dependence (network dependence? spatial dependence? cluster dependence? some combination of these?).\n", "\n", "For these scenarios, we can resort to a resampling method due to [Song (2016)][7]. [Leung (2022b)][5] shows that this method works for regular estimators (like OLS, TSLS, IPW) whenever they are $\\sqrt{n}$-consistent, a requirement satisfied by most forms of weakly dependent data. Since the method does not require knowledge of the dependence structure, it can be applied to time series data when time periods are unobserved, spatial data when spatial locations are unknown, or clustered data when cluster memberships are unknown. \n", "\n", "We compute dependence-robust CIs using the following method." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls.drobust_ci(-4, 4, seed=0)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " Estimate CI\n", "---------- -------------\n", " 0.934 [0.16, 1.6]\n", " 0.522 [0.16, 0.907]\n", " 2.851 [1.333, 3.84]\n", " 1.007 [0.64, 1.44]\n" ] } ], "execution_count": 18, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:12.000Z", "iopub.execute_input": "2022-04-10T00:49:12.304Z", "iopub.status.idle": "2022-04-10T00:49:20.240Z", "shell.execute_reply": "2022-04-10T00:49:20.283Z" } } }, { "cell_type": "markdown", "source": [ "This method does not use the network data supplied to the `tsls` object. For use with non-network data, we can invoke this method by instantiating the `tsls` object without supplying the network argument at all." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "code", "source": [ "tsls_no_network = ni.TSLS(Y, regressors, instruments) # no network argument\n", "tsls_no_network.drobust_ci(-4, 4, seed=0) # same results" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " Estimate CI\n", "---------- -------------\n", " 0.934 [0.16, 1.6]\n", " 0.522 [0.16, 0.907]\n", " 2.851 [1.333, 3.84]\n", " 1.007 [0.64, 1.44]\n" ] } ], "execution_count": 19, "metadata": { "collapsed": true, "jupyter": { "source_hidden": false, "outputs_hidden": false }, "nteract": { "transient": { "deleting": false } }, "execution": { "iopub.status.busy": "2022-04-10T00:49:20.254Z", "iopub.execute_input": "2022-04-10T00:49:20.264Z", "iopub.status.idle": "2022-04-10T00:49:28.086Z", "shell.execute_reply": "2022-04-10T00:49:28.103Z" } } }, { "cell_type": "markdown", "source": [ "This method is the most computationally demanding because it resamples the data to construct test statistics and then constructs the CI via test inversion. It has two other important disadvantages.\n", "\n", "1. It requires a larger amount of data to ensure adequate coverage and power. In particular, it is much less powerful than the other procedures, yielding wider CIs. Of course, since this is a method that relies on little prior knowledge about the dependence structure, it is unsurprising that it can deliver weaker conclusions.\n", "\n", "2. The output is random, which is why we set the seed in the `drobust_ci()` method. However, the method contains a parameter `L`, the number of resampling draws, which can be made arbitrarily large to make the randomness of the output arbitrarily small. This only requires more computation time." ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "## Other methods\n", "\n", "Invoke\n", "\n", " help(tsls)\n", " \n", "for information on the other available methods, which are scalar equality dependence-robust and approximate randomization tests.\n", "\n", "Finally, all of the methods discussed above are accessible as independent functions through the 'core' class for use with custom estimators. See\n", "\n", " help(ni.core)" ], "metadata": { "nteract": { "transient": { "deleting": false } } } }, { "cell_type": "markdown", "source": [ "## References\n", "\n", "* [1]: Canay, I., J. Romano, and A. Shaikh, \"Randomization Tests Under an Approximate Symmetry Assumption,\" *Econometrica*, 2017, 85 (3), 1013-1030.\n", "* [2]: Ibragimov, R. and U. Mueller, \"t-Statistic Based Correlation and Heterogeneity Robust Inference,\" *Journal of Business and Economic Statistics*, 2010, 28 (4), 453-468.\n", "* [3]: Kojevnikov, D., V. Marmer, and K. Song, \"Limit Theorems for Network Dependent Random Variables,\" *Journal of Econometrics*, 2021, 222 (2), 882-908.\n", "* [4]: Leung, M. \"Causal Inference Under Approximate Neighborhood Interference,\" *Econometrica*, 2022a, 90(1), 267-293.\n", "* [5]: Leung, M. \"Dependence-Robust Inference Using Resampled Statistics,\" *Journal of Applied Econometrics*, 2022b, 37(2), 270-285.\n", "* [6]: Leung, M., \"Network Cluster-Robust Inference,\" *arXiv preprint arXiv:2103.01470*, 2022c.\n", "* [7]: Song, K. \"Ordering-Free Inference from Locally Dependent Data,\" *UBC working paper*, 2016.\n" ], "metadata": { "nteract": { "transient": { "deleting": false } } } } ], "metadata": { "kernel_info": { "name": "python3" }, "language_info": { "name": "python", "version": "3.9.9", "mimetype": "text/x-python", "codemirror_mode": { "name": "ipython", "version": 3 }, "pygments_lexer": "ipython3", "nbconvert_exporter": "python", "file_extension": ".py" }, "kernelspec": { "argv": [ "/usr/local/opt/python@3.9/bin/python3.9", "-m", "ipykernel_launcher", "-f", "{connection_file}" ], "display_name": "Python 3 (ipykernel)", "language": "python", "metadata": { "debugger": true }, "name": "python3" }, "nteract": { "version": "0.28.0" } }, "nbformat": 4, "nbformat_minor": 0 }