diff --git a/scripts/finding_gravity_model_exponent.ipynb b/scripts/finding_gravity_model_exponent.ipynb new file mode 100644 index 0000000..7079707 --- /dev/null +++ b/scripts/finding_gravity_model_exponent.ipynb @@ -0,0 +1,221 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from src.model import PyTradeShifts\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from src.utils import get_distance_matrix\n", + "from scipy.stats import linregress\n", + "import os\n", + "\n", + "data_dir = os.path.join(os.path.dirname(os.getcwd()))\n", + "os.chdir(data_dir)\n", + "\n", + "\n", + "def _get_matrices():\n", + " \"\"\"Get trade and distance matrices for the global wheat trade in 2018.\"\"\"\n", + " Wheat2018 = PyTradeShifts(\n", + " \"Wheat\",\n", + " 2018,\n", + " region=\"Global\",\n", + " testing=True,\n", + " make_plot=False,\n", + " )\n", + " Wheat2018.load_data()\n", + " Wheat2018.remove_net_zero_countries()\n", + " Wheat2018.prebalance()\n", + " Wheat2018.correct_reexports()\n", + " np.fill_diagonal(Wheat2018.trade_matrix.values, 0)\n", + "\n", + " dm = get_distance_matrix(\n", + " Wheat2018.trade_matrix.index, Wheat2018.trade_matrix.columns\n", + " )\n", + " return Wheat2018.trade_matrix, dm" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Loaded data for Wheat in Y2018.\n", + "Removed 0 countries with no trade or production.\n", + "Prebalanced trade matrix.\n", + "Corrected re-exports.\n" + ] + } + ], + "source": [ + "tm, dm = _get_matrices()\n", + "# a little melting trick, this will produce:\n", + "# a df with three columns: country A, country B, volume of trade from A to B\n", + "tm = tm.stack().reset_index()\n", + "# a df with three columns: country A, country B, distance from A to B\n", + "dm = dm.stack().reset_index()\n", + "# index order should match because tm and dm should have the same index/columns\n", + "# so we can just take the values and roll with it\n", + "# otherwise we should do a formal join here\n", + "y = tm[0].values\n", + "x = dm[0].values\n", + "# filter out zeroes\n", + "x = x[y != 0]\n", + "y = y[y != 0]\n", + "# sort values by distance\n", + "y = y[x.argsort()]\n", + "x = x[x.argsort()]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use the approach from [Karpiarz et al.](https://arxiv.org/pdf/1409.5963.pdf) (Fig. 1), since it is beyond the scope of our work to assess and validate the (well-established) gravity model.\n", + "\n", + "We do not trouble ourselves with questioning whether it is or isn't a \"proper\" power-law etc., we only focus on obtaining comparible quantities to the work present in the above paper.\n", + "\n", + "For a more advanced approach to power-law fitting see [Clauset et al.](https://arxiv.org/pdf/0706.1062.pdf)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/luka/.conda/envs/pytradeshifts/lib/python3.12/site-packages/numpy/lib/function_base.py:1242: RuntimeWarning: divide by zero encountered in divide\n", + " a = -(dx2)/(dx1 * (dx1 + dx2))\n", + "/home/luka/.conda/envs/pytradeshifts/lib/python3.12/site-packages/numpy/lib/function_base.py:1243: RuntimeWarning: divide by zero encountered in divide\n", + " b = (dx2 - dx1) / (dx1 * dx2)\n", + "/home/luka/.conda/envs/pytradeshifts/lib/python3.12/site-packages/numpy/lib/function_base.py:1244: RuntimeWarning: divide by zero encountered in divide\n", + " c = dx1 / (dx2 * (dx1 + dx2))\n", + "/home/luka/.conda/envs/pytradeshifts/lib/python3.12/site-packages/numpy/lib/function_base.py:1250: RuntimeWarning: invalid value encountered in add\n", + " out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]\n", + "/home/luka/.conda/envs/pytradeshifts/lib/python3.12/site-packages/numpy/lib/function_base.py:1259: RuntimeWarning: divide by zero encountered in scalar divide\n", + " out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABM8AAAIKCAYAAAA04mjzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAABcSAAAXEgFnn9JSAABbI0lEQVR4nO3dd3yV5f3/8ffJOCcDAmRBQiBgwlKQPQUX1gJFq4xa62qt1taqddcO66zt136riKu2X62j1oVbQX4qgspeKmgYIYSQwcqCkORk3b8/MJFwzp2c5Iz7nJzXs488LNe9Puc+uUjOm+u+LpthGIYAAAAAAAAAuIiwugAAAAAAAAAgWBGeAQAAAAAAACYIzwAAAAAAAAAThGcAAAAAAACACcIzAAAAAAAAwAThGQAAAAAAAGCC8AwAAAAAAAAwQXgGAAAAAAAAmCA8AwAAAAAAAEwQngEAAAAAAAAmCM8AAAAAAAAAE4RnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABgIsofJ62rq9M333yjLVu2aM+ePSorK1NTU5NSUlI0evRonX322YqJiXF77Jo1a7RixQrt27dPkZGRGjhwoGbMmKGsrCx/lAoAAAAAAACYshmGYfj6pCtXrtSLL74oSUpLS1NaWppqa2uVl5en2tpa9e7dWzfffLO6d+/e6rhFixZp2bJlio6O1rBhw9TQ0KBt27ZJkq666iqNGjXK16UCAAAAAAAApvwy8iwyMlLTpk3T9OnTlZqa2tJeWVmpJ554Qnv37tVrr72mK6+8smXb9u3btWzZMsXHx+u2225rOS4vL08LFizQCy+8oMGDBysuLs4fJcOEYRjas2ePMjMzZbPZrC4HCEn0I8A36EuA9+hHgPfoR4D3Qq0f+WXOs0mTJuniiy9uFZxJUo8ePXTRRRdJkr744gs1NDS0bPvoo48kSTNnzmx13EknnaSpU6eqpqZGq1at8ke5AAAAAAAAgFsBXzAgIyNDktTQ0KCjR49Kkurr67V9+3ZJ0ujRo12OGTNmjCRpy5YtAaoSAAAAAAAAsCA8O3TokKRjj3Y2P4K5b98+NTQ0qFu3burVq5fLMf369ZMkFRUVBa5QAAAAAAAAhD2/zHnWlk8++USSdPLJJys6OlqSVF5eLklugzNJcjgcio2NVXV1tWpra01X6jyeH9ZBCEuGYbR8Aegc+hHgG/QlwHv0I8B79CPAe1b1o87OrxbQ8Gzr1q1atWqVIiMjdd5557W0O51OSZLdbjc91uFwqKamRk6n06PwbM+ePd4XDBmGobKyMkmd/yYDwh39CPAN+hLgPfoR4D36EeA9q/rRgAEDOnVcwMKzkpISPfvsszIMQxdeeGHL3GeSZ6PEOppGZmZmdrhGuGq+76GyAgYQjOhHgG/QlwDv0Y8A79GPAO+FWj8KSHhWXl6uxx9/XNXV1Zo+fbrOPvvsVtubR5LV1dWZnqN5m8Ph8OiaoXDzQ4XNZmv5AtA59CPAN+hLgPfoR4D36EeA90KpH/l9wYCqqiotXLhQZWVlmjx5subMmeOyT/NcZ81zn53I6XSqpqZGsbGxHj2yCQAAAAAAAPiCX8Oz2tpaPfbYY9q/f79GjRqlSy65xG2i2Lt3b0VFRamqqsptgLZ3715JUt++ff1ZLgAAAAAAANCK38Kz+vp6/eMf/1BBQYFOPvlkXXnllYqIcH85u92uIUOGSJI2b97ssn3Tpk2SpBEjRvirXAAAAAAAAMCFX8KzpqYmPfPMM9qxY4eys7P1i1/8QlFRbU+vNn36dEnSkiVLdODAgZb2vLw8ff7554qJidGUKVP8US4AAAAAAADgll8WDFi+fLm+/PJLSVJ8fLxefvllt/vNmTNH3bp1kyQNHTpUZ511lj755BM98MADGjp0qBobG5WTkyPDMPSzn/1M8fHx/igXAAAAAAAAcMsv4Vl1dXXL/28O0dz5wQ9+0BKeSdL8+fOVkZGhFStWaNu2bYqMjNTQoUM1Y8YMZWdn+6NUAAAAAAAAwJRfwrPZs2dr9uzZnTp28uTJmjx5so8rAgAAAAAAADrOr6ttAgAAAAAAAKGM8AwAAAAAAAAwQXgGAAAAAAAAmPDLnGfwrdKjpXp4xULlHtql7OQs3XTGDUqKT/Lb8dV11Vqc84F2HszVoORsDe92si9eRqfrMTtGkmnbjoM7FRcdpwGJmRqedopmDZuhOHucquuq9dbWd7Q891PZJJ2RfYYuGH6eJH33mlOyW/YP9Os8XvP7kLN/m5wNTjmiHBrWe6hmDZuhmvoa03N7et3CikL9cfHdKqwoVGREpLrHdNfAxIGSpL0Ve12OLawo1B3v/VEF5QWKsEVqTL/R+sM5v233ur64Dy9uelmLvnhDNfU1ckQ65Ih2yNngVO/uqUrtnqq80jxV1hxWY1Ojesb00EnJA7XvyP52781Px1+mp1Y/rc1FX6ipqUmxUTGqqj8qwzDUv1c//XX2/dp3eL9ueus2NamppSZHpEPORqdssh17X1KH6lD1IRVXFrvcm+q6ai3+5gNtzv9CfQv7aum2/6eK2gr1cByrs7CySA2NDYqOtCs5Lklf7/tGjWpsuVZURJTG9hvT5r2+YPj5unfpn1VWU6YIW4TGZIzWdVN/pQUrHtPm4i8kSZG2SM0fOUdnZZ+hPy65W+U15eoV20t//N7v9NbWd7TtwHY1NDYoMiJKNpu07/B+GTKO1WCLUqPRqMiISPWM7anBKYOUENNdSfFJWrx1icqdFS31RtoiJUkxUQ7FOeJU56xTZf3hlu3REdGSpF5xvXTukOnaU75XeaV5iomKUYQtQlXOKvWM66mJmRO169Au5R7KVUNjg2y2CBlGkxoaG1VeWy5JirBF6MIRP9SPRs3VI58+7vI+NjQ1tPpeirBFKCk+SdMGnqZPd32msuoyRUZEKVKRqm2qbfP7MFKRGtZnqAorixQVEammpiaV11S03CNJskfaldajjyqrKlVRV9lynKSW99QR4VBUZKScDXVyRNkVa49VtbNa1Q01bV7/RNG2aDki7KpqPNqqxuO/d45nk00RipBhM9RkNLndx5dssskeYVd9U32rvtPW/vfO+JPOHnyW32sDgGBUdfSo3n1viQ4fOaKhQwZr0oRxcjgcVpcFAPiWzTAMo/3dYJXSo6Wa9+zFqmusa2mzR9q16KcveRRAdPT46rpq/fK167SrNK+lrX9CP/3fxU8p3uFdmNSZesyOibJFyRZhU31jfZttzbKTs/TQDx/Ub968RbvL8lttG5CYKZtsrdqzk7P05LxHOx2gefu+uXsfjq+3qLK41etsPrckj65bWFGoHz9/WasP/u40H1tTX6OLnr/UZXt0ZLRe/+nLptf914+e0NWvXuvVfbjy5V9ob0Vhu/u29xrc1ehP0ZHR+s8l/9bv3/+T2/exM+czu9cAfMsmmyIVqQY1uGyLUpTbdqtER0Trf8//q8b1H2t1KQgRhmFoz549yszMlM1ms7ocSKqtrdXv/3Sf6uu/+90uPT1NN99wLQFakKIfAd4LtX7EyLMg9/CKhS4fkusa67RgxaO6b9bdPj9+cc4HLh/0Cw7v1ZJtH2jeyDkdrt/besyOaTAadOIAC3dtzXIP7dLDKx51Cc4kKb9sj9v9l+Qs1dyRF5q+lrZ4+765ex/aqrf53IYMj677x8V3txucHX+sWXhV31jf5nVvfecOr++DN8HZ8ddzV6M/1TfW687Fd/skOGs+nxWvAwhHhgzTgCyYgjNJqm+q12/eusXqMrz2g2EzdeMZ13s96hsIRW+/u6RVcCZJxcUlWrt+o06fOsWiqgAAxyM8C3K5h3Z1qN3b43cezHW//0HPrufretrb1hG7OnienYfc3wtPePu+mb0P7V3TLBA78brFh4s7dN7S6tJOXbe8psKjesx05j6YXc+TsNDXig6X+PR8Vr0OAPC393OW6P2cJVaXERAj+gzXn39wT4emMEDXlpef77a9sMjz39cAAP5FeBbkspOz3I68yU7O8svxg1Ky3Z8nxbPr+bqeto7pqKzkLBVU7PV4/0HJ7u+FJ7x938zeh/auacjw6LrpCekeh4PZyVlyVDhM92/rur1ie+rg0UPt1mOmM/fBnbZq9Ke+CWna4UUIeyKrXgcAwHe27Nuq85+ea3UZAZHWPU0LLvybMnpmWF1KUOverZvb9oy+6QGuBABghvAsyN10xg1auXu1y5xRN55xvV+OnzVsht79+v1WI4P6J/TTzKEzOvkKvKvH7JjOzHl20xnXa0/5HuWV7m61bUBipiJsEa3as5OzNHPY9zv1Gs1q7sj75u59OL7e4soS03N7ct37Z93t8ZxnN55xfZtznrV13f89/69u5zzryH14/as3VVDueehp9hrc1ehP0ZHRum/W3frD4rt8MnqyrXsNAEAwKjlS4vb3h1ATExWjp+Y/pmwf/aPeiVJTU7R9x85WbVFRURp16nC/XA8A0HEsGBACSo+WasGKR1tW1rvxjOs7vNpmR46vrqvWkpyl2nkoV9nJWRoef7KGZA3x2SR+nXk97o6RZNp24mqbM4d9v2W1zbe3vqvluZ9Kks7MPl0//Ha1zebXPCg5u2X/QL/O4zW/Dyeutjlz2PdVU19jem5Pr1tYUag7F9+jvcettnlS4gDJZtPe8r0uxxZWFOp3792pPeUFirBFaGy/Mfr9Obe3e11f3If/bnpFi754XdX1NYqJdMhhj5Gzvla9u/duWW2zoqby29U2eyrruNU227o3V4y/VP9c/bQ2Na/SGB2rqroqGYahzMT++ssP7nO72mZMpEO1x622eXLvYTpUfUhFFcUu96Z51dTNu79Q35S+WrrtQ1XUHFttM6vVapvRx1bb3J9zbP6+b0VHRGlsv7Ft3usfDj9P9y59oGW1zbEZo/Xrqb/UIyse06ZvV9uMiojSvFMv1FnZZ+jOJfeorKZcibG99Ifv3aG3t7573GqbkbLZbNp3eJ+aTlhtM6p5tc3UwUpwHFtt8/2vl7Ssfikdt9pmdIzi7M2rbVYe93qOrbaZGNdL5w45x2W1zSPOI+oV10sTMydo16G841bbtMkwDDU0NqjsuNU254z4oeaPmquFnz7u8j66W20zOT5J0wZO1Yq8z1R2tEyREZEer7Z5cp9hKqwsVGRElJqMJlVUl7fcI+lYUJveI00VRytV8e0KpCeuthkTGaPIyEg5651yRDsUG9282mZ1m9c/UbQtWo5Ih6oaqlraohSphnZW25RNajRMJof0IZtsckTYVefhapsAEE5unHq95o/5bvThY0/+yyU8k6T5cy9gzrMgFWoTnQPBKNT6EeEZ2hRq39BAMKIfAe4VVhTqt+/+QfnlrguhWCni2/+5WxwgWtGql+sIZwDorKGFpyiuPt6lvTbxqP56873MjxeE+N0O8F6o9SMe2wQAAJbI6JmhFy97zuP9Q+2XLH/ZULBRN791e8uIRgChrSGyXu4y+QO2g7rgmfl668rXCNAAwGKEZwAAACFkXP+x+vSGj60uw2ulR0t121u/0/bSHVaXAliqxl6jhNqerdqq7dUq63ZITUaT3xeYuG/GXTp78Fl+vQYAhDrCMwAAAARcUnySnrnkn1aXERCFFYX6zZu3aN+R/VaXgiBk2Fxn0cntvU1NEYGZM/LOD+7RnR/cE5Br3X7mLfrhqecF5FoA4EuEZwAAAIAfZfTM0Os/e8XqMiT5//Hn3IO5uual61SrthdiwXdshuv7kL1/qHakfROwAC1QHlz+dz24/O8u7TbZ9D+z/6zTTmKBBADBifAMAAAAgE9kp2Tr4xs+sLoMr7226XUt+PzRgFwrti7WpS2uLk6JVck6lHAgIDVYzZCh29/7fYePu3byNbpk/MV+qAgAWiM8AwAAAIDjzB8zV/PH+GausdyDubrqpV+ZrtQb1Rjttj3OGeeT63dlT6x+Sk+sfqpTx0ZFROnxOQs0PH24j6sC0BURngEAAACAn2SnZGv5DR+abn/syX9p+46dLu3Vjmp/lhX2GpoadM2i61zabbJpwQX/q3H9x1pQFYBgRXgGAAAAABZJT+vjEp41r7aJwDNk6Ddv3eJ228BeA/TPi55QbLTro7YAujbCMwAAAACwSFSU60eyJ+79uxwOh1+ut+Sbpbr/o7/45dxd3e7yfH3vH7M6dMyDsx9gIQSgCyA8AwAAAIAg4q/gTJJmnvx9zTz5+347//FKj5bqp/+5SmXO8oBcLxh5shDChSf/ULeec1MAqgHQWYRnAAAAAACfS4pP0rvXvOl2W3Vdtf7y0d+0LPeTAFcVfN785m1JIkADghjhGQAAAAAgoOLscbpv1l26T3d5fMzW4q1uJ/nvCt785u2WEE2Ssntl67nL/s/CigAcj/AMAAAAABD0hqcP18oblnf6+Cc+/Yde/OJl3xXkR7nluTpt4Zku7RP7jdf9P7hHcfa4wBcFhDHCMwAAAACwiGEYVpcQNq49/Ze69vRfurSHUqi2du96fe8fsxQVEaXH5yzQ8PThVpcEhAXCMwAAAABA2DIL1STp7sX36cPcjwNcUfsamhrcPsLK4gOAfxCeAQAAAADgxt2z7tTdurNVm2EY2rNnjzIzM2Wz2Vpte3bN8/rXumcCWWIrx8+dRpAG+A7hGQAAAAAAPvDTSZfrp5MuN90eyEUPWMUT8J0IqwsAAAAAACAcDE8frqfmPRaw6x2/gieAzmPkGQAAAAAAAeJu1dDSo6W69Nmf6nDjEZ9f78RVO+0Rdr1w6TPK6Jnh82sBXRXhGQAAAAAAFkqKT9KSX7/rdtvNb9yutYXrfHatuqY6XfT8pS7tQ5IH67F5CxRnj/PZtYCugvAMAAAAAIAg9dCcByVJ/1r5tJ7d+ILfrrP90A597x+zWrVdO/kaXTL+Yr9dEwgVhGcAAAAAAAS5q0/7ua4+7eet2lbmrdLt7/3eb9d8YvVTemL1U23u88NTztN1037FiDV0aYRnAAAAAACEoNNOmqKVNyzXT/97lXYeyrWkhre/fldvf+3+kdPOGJo6RA+e94CS4pN8dk7AW4RnAAAAAACEsPtn3e12HrNQtO3Adp3/9Fy32+yRdj164UManj48wFUh3BGeAQAAAAAQwjJ6ZuiVy/+jn73wC1Ub1VaX4zd1jXW6ZtF1rdp4bBSBQHgGAAAAAECIy+iZoQ+vX+x221//39/07rb3A1xRYLT32Cij1eALhGcAAAAAAHRhd5x7m+449zaX9q4cqjVrHq321LzHCNDQaYRnAAAAAACEoRNDNX+v3mmlEx/3PN7I9FN138y7WKQApgjPAAAAAABAy+qd7Xlx/Ut6YvVT/i8oQL4s/krnPz1XMVExeuSC/2WEGlwQngEAAAAAAI9dMv5iXTL+Yp+d77VNr2vB54/67HydVdtQ63aEWu9uqVo45yFl9MywoCoEA8IzAAAAAABgmflj5mr+mLlutxVWFOo3b96ifUf2B7iq7+yvOqCLnr9UkhRhi9Bff3C/TjtpimX1IPAIzwAAAAAAQFDK6Jmh13/2Squ2ZTs+0Z0f3GNJPU1Gk25/7/d6cPYDBGhhhPAMAAAAAACEjLMHn6WzB5/ldlvp0VJd8PR8NanJrzXc/t7v9c7PX2eRgTBBeAYAAAAAALqEpPgkvfXz13Tf0ge0vnCjX691/tPfPWo6NHWIHjzvAcK0LorwDAAAAAAAdBlJ8UlaMOfvbe5TWFGoW97+rQori3xyzW0Htuv8p+fquYv/T9kp2T45J4JHhNUFAAAAAAAABFJGzwy9csWLWnnDcj017zFFK9on5/3pS1er9GipT86F4MHIMwAAAAAAELaGpw/X8hs+bNWWezBX1716k440HunQuQwZWrDiUd03624fVgirEZ4BAAAAAAAcJzslWx/8+l1J0oaCjfrNW7d4fOymvV/4qSpYhcc2AQAAAAAATIzrP1Y3Tr3e4/0rnBU6beGZmv7EDG0t3urHyhAojDwDAAAAAABow/wxczWszxD95o1bVdtU69ExtQ21umbRdZKk7o7uemzOwywmEKIIzwAAAAAAANoxPH24Pr7ug5Y/n/PkDNXUexakHXEe0RUvXeXSHhURpb+f/z8a13+sz+qE7/HYJgAAAABYxDAMq0sA0EkZPTK8PkdDU4N+89Yt2lCw0QcVwV8IzwAAAAAAADrofh+uqHnbu7/z2bnge4RnAAAAAAAAHZTRM0OvXP4fJcYmen2uusY6H1QEfyE8AwAAAAAA6ISMnhl69+o39M7PX9fw1JM7fR57pN2HVcHXCM8AAAAAAAC8kBSfpKd+/IQ+/OVinTv4ex0+vrGpUTe/fbtKj5b6oTp4i/AMAAAAAADAB+Lscbprxh+08obleufnr2tknxEeHddoNGrtnnU6/+m5unfpn1VdV+3nStERUVYXAAAAAAAA0NUkxSfpiR896tL+0/9epZ2Hck2PW7r9Qy3d/qEkqVdsLy244G/KTsn2W51oH+EZAAAAAABAgBQfLvZ43/Kacl3x0lWKtEUqKT5JF4w4T/NHzlWcPc6PFeJEPLYJAAAAAAAQIOkJ6R0+ptFo1IGqA/rn6qf181eu4bHOACM8AwAAAAAACJD7Z93t1fEF5Xv19tZ3fVMMPEJ4BgAAAAAAECAZPTP0yuX/Uf+e/Tp9jne3vu/DitAev815VlBQoJycHOXn5ys/P1+VlZWKiorSwoUL3e7/3nvvafHixabnO/fcc3XBBRf4qVoAAAAAAIDAyOiZoZcuf0HVddV69NMn9c43HRtJtqeiQKVHS5UUn+SnCnE8v4Vnixcv1ldffdXh47KyspSSkuLS3r9/f1+UBQAAAAAAEBTi7HH67Tm36Lfn3KKtxVt145u3qaaxxqNj71v6gBbM+bufK4Tkx/DspJNOUkZGhjIzM5WZmak77rjDo+OmTJmiyZMn+6ssAAAAAACAoDM8fbg++vUSSVLuwVzd/NZvVVpTarr/xqLNgSot7PktPDv33HP9dWoAAAAAAIAuKzslW+9c/bpyD+bqipeucruPYRgBrip8sWAAAAAAAABAEMpOyVaPmB5utyXGJwa4mvDlt5FnnbV9+3YVFhaqvr5evXr10imnnMJ8ZwAAAAAAICwtvPDvbkef/el7v7egmvAUdOHZunXrWv353Xff1ejRo3XZZZcpJibG4/MwfNE3DMNo+QLQOfQjwDfoS4D36EehgfcnuNGPEGg9Y3spUhFqVFOr9lvfvUOLrnhZSSE4As2qfmSz2Tp1XNCEZykpKZozZ45OOeUUJSYmqrq6Wrm5uXrzzTe1efNmNTU16ZprrvH4fHv27PFjteHDMAyVlZVJ6vw3GRDu6EeAb9CXAO/Rj4JP5eFKlzY+ywQ3+hEC7ZH1j7sEZ5JU31ivBz74q24Y92sLqvKOVf1owIABnTouaMKziRMntvqzw+HQ+PHjNXjwYN1///368ssvtWvXLmVlZXl0vszMTH+UGXaaU+DMzEx+MACdRD8CfIO+BHiPfhR8Er7c6tLGZ5ngRj9CoJV8WmK6bW3xev3s/V/oRyPn6fLxlyrOHhvAyjov1PpR0IRnZnr06KHJkyfro48+Uk5OjsfhWSjc/FBhs9lavgB0Dv0I8A36EuA9+lFwcfc28N4EP/oRAik7OUt7KwpNt9c11uk/m/6rFXmf6pkf/1Nx9rgAVtd5odSPQmK1zdTUVElSZaXrkGYAAAAAAICu6qYzbvBov70VhVqSs9TP1YSnkAjPqqurJR17lBMAAAAAACBcJMUnacqAyR7tm7N/m5+rCU9BH54ZhqEvvvhCktS/f39riwEAAAAAAAiwO6bf6tF++w/v93Ml4SkowrOqqiqtWbNG9fX1rdpra2v10ksvKT8/XwkJCRo1apQ1BQIAAAAAAFgkKT5JT817rN39vtrnuggJvOe3BQO2bNmiJUuWtGprbGzUgw8+2PLnmTNnasSIEXI6nXr++ef16quvqk+fPkpMTFR1dbX27t2ro0ePKjY2VldffbXsdru/ygUAAAAAAAhaw9OH65XL/6Pb3vmdCir2ut2noakhwFWFB7+FZ1VVVcrPz2/VZhhGq7aqqipJUnx8vM4991zt3r1bBw8eVGFhoSIiIpSUlKRJkyZp+vTp6tmzp79KBQAAAAAACHoZPTP00uUv6NwnZ+lofbXL9nh7vAVVdX1+C88mT56syZM9m9AuJiZGF1xwgb9KAQAAAAAA6DIe+uGDumbRda7t5/+PBdV0fUEx5xkAAAAAhCPDsLoCAKFoePpwPTXvMXWzd5NNNnWzd9NT8x7T8PThVpfWJRGeAQAAAAAAhJi0Hmk6NX2EYqJj1Gg06rkNL6r0aKnVZXVJhGcAAAAAAAAhpPRoqeY8c5FW5a9WTX2NauprtCp/tS54ej4Bmh8QngEAAAAAAISQh1csVIPhurJmk5r0wIfMe+ZrhGcAAAAAYBGbzeoKAISi3EO7TLdtLNwcwErCA+EZAAAAAABACMlOzjLdVt9Ur8KKwgBW0/URngEAAAAAAISQm864oc3tFz1/KQGaDxGeAQAAAAAAhJCk+CQ9Ne+xNvf58fOXaWvx1gBV1LURngEAAAAAAISY4enDlZbQx3S7IUPXLLpOS75ZGsCquibCMwAAAACwiGFYXQGAUDY0dUi7+9z/0V+0bMcnAaim6yI8AwAAAAAACEHtzX3W7M4P7lHuwVw/V9N1EZ4BAAAAAACEoKT4JP1swuUe7fvzl3+p/256RdV11X6uqushPAMAAAAAAAhRPxnzY6X3SGt3vwajQY9//qSuee3XBGgdRHgGAAAAAAAQouLscXru4qd10aj5Hu2fV7pbS3JYRKAjCM8AAAAAAABCWJw9Tjec/mu98/PX1S26W7v7bynZGoCqug7CMwAAAAAAgC4gKT5Jb/78VY3vN7bN/faUFQSooq6B8AwAAAAAAKCLiLPH6YEf3KeTkgaa7nPYeTiAFYU+wjMAAAAAAIAuJM4ep6fmP67MXv3dbt93ZL/OfOwcXf7ilSqsKAxwdaGH8AwAAAAAAKCLibPHqW+Pvqbb65satKs0Txc9f6lyD+YGsLLQQ3gGAAAAABYxZFhdAoAubG/FXo/2u/rVa1VdV+3nakIX4RkAAAAAAEAXlJ2c5dF+dY11WpKz1M/VhC7CMwAAAAAAgC7opjNuUJQtyqN9dx7i0U0zhGcAAAAAAABdUFJ8kt648hWNSR/V7r6DkrP9X1CIIjwDAAAAAADoopLik/TovAV65fL/KDspSxFuoqDs5CzNHPZ9C6oLDZ6N3QMAAAAAAEDIyuiZoecueVqSVF1XrSU5S7XzUK4GJWdr5rDvK84eZ3GFwYvwDAAAAAAAIIzE2eM0d+SFVpcRMgjPAAAAAAAAwkzp0VI9vGKhcg/tUnZylm464wYlxSdZXVZQIjwDAAAAAIvYZLO6BABhqPRoqeY9e7HqGuskSXsrCrVy92ot+ulLBGhusGAAAAAAAFjEkGF1CQDC0MMrFrYEZ83qGuu0YMWjFlUU3AjPAAAAAAAAwkjuoV0dag93hGcAAAAAAABhJDs5q0Pt4Y7wDAAAAAAAIIzcdMYNskfaW7XZI+268YzrLaoouBGeAQAAAAAAhJGk+CQt+ulLOjv7TPXv2U9nZ5/JYgFtYLVNAAAAAACAMJMUn6T7Zt1tdRkhgZFnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABggvAMAAAAAAAAMMGCAQAAAABgFcPqAgCEq+q6ai3O+UA7D+ZqUEq2Zg2boTh7nNVlBSXCMwAAAAAAgDBSXVetX752nXaV5rW0vfv1+3py3qMEaG7w2CYAAAAAAEAYWZzzQavgTJJyD+3SkpylFlUU3AjPAAAAAAAAwsjOg7nu2w+5bw93hGcAAAAAAABhZFBKtvv2ZPft4Y7wDAAAAAAAIIzMGjZD2clZrdqyk7M0c9j3LaoouLFgAAAAAAAAQBiJs8fpyXmPaknOUu08lKtBydmaOez7LBZggvAMAAAAAAAgzMTZ4zR35IVWlxESeGwTAAAAAAAAMEF4BgAAAAAWMQzD6hIAAO0gPAMAAAAAAABMEJ4BAAAAAAAAJgjPAAAAAMAiNpvN6hIAAO0gPAMAAAAAAABMRFldAAAAAAAAAAKvuq5ai3M+0M6DuRqUkq1Zw2Yozh5ndVlBh/AMAAAAAAAgzFTXVeuXr12nXaV5LW1PrHxKV4y/VPNHziVEOw6PbQIAAAAAAISZxTkftArOJMnZ4NQ/Vz+tq1/9larrqi2qLPgQngEAAAAAAISZnQdzTbfll+3Rff/vAQK0bxGeAQAAAIBFDMOwugQAYWpQSnab2z/N+1y/WnQ9AZoIzwAAAAAAAMLOrGEzlByf3OY+uYd26e2t7waoouBFeAYAAAAAABBm4uxxeubHTykuuu2FAZ5f/5+wH31GeAYAAAAAABCGkuKT9PLlLygpLsl0n8POI1qSszSAVQUfwjMAAAAAAIAw1RygXTf1V4qOiHa7z85D5osLhAPCMwAAAAAAgDAWZ4/TxWMu0i+mXOV2+6DkthcX6OoIzwAAAAAAAKALhp+n7OSsVm3ZyVmaOez7FlUUHKL8deKCggLl5OQoPz9f+fn5qqysVFRUlBYuXNjmcWvWrNGKFSu0b98+RUZGauDAgZoxY4aysrLaPA4AAAAAAACdF2eP05PzHtWSnKXaeShXg5KzNXPY9xVnb3tRga7Ob+HZ4sWL9dVXX3XomEWLFmnZsmWKjo7WsGHD1NDQoJycHOXk5Oiqq67SqFGj/FMsAAAAAAAAFGeP09yRF1pdRlDxW3h20kknKSMjQ5mZmcrMzNQdd9zR5v7bt2/XsmXLFB8fr9tuu02pqamSpLy8PC1YsEAvvPCCBg8erLi48E47AQAAAAAAEDh+C8/OPffcDu3/0UcfSZJmzpzZEpxJx0K4qVOnavny5Vq1apXOOeccn9YJAAAAAAAAmAmKBQPq6+u1fft2SdLo0aNdto8ZM0aStGXLloDWBQAAAAAAgPAWFOHZvn371NDQoG7duqlXr14u2/v16ydJKioqCnRpAAAAAAAACGN+e2yzI8rLyyXJbXAmSQ6HQ7GxsaqurlZtba1iYmLaPadhGD6tMVwZhtHyBaBz6EeAb9CXAO/Rj4KPu/eC9ye40Y8A71nVj2w2W6eOC4rwzOl0SpLsdrvpPg6HQzU1NXI6nR6FZ3v27PFZfeHMMAyVlZVJ6vw3GRDu6EeAb9CXAO/Rj4LPkSNHXNr4LBPc6EeA96zqRwMGDOjUcUERnnmSNHY0jczMzOxsOThO833PzMzkBwPQSfQjwDfoS4D36EfBp3v3L13a+CwT3OhHgPdCrR8FRXjWPJKsrq7OdJ/mbQ6Hw6NzhsLNDxU2m63lC0Dn0I8A36AvAd6jHwUXd+8D703wox8B3gulfhQUCwY0z3XWPPfZiZxOp2pqahQbG+vRI5sAAAAAAACALwRFeNa7d29FRUWpqqrKbYC2d+9eSVLfvn0DXRoAAAAAAADCWFCEZ3a7XUOGDJEkbd682WX7pk2bJEkjRowIaF0AAAAAAAAIb0ERnknS9OnTJUlLlizRgQMHWtrz8vL0+eefKyYmRlOmTLGqPAAAAAAAAIQhvy0YsGXLFi1ZsqRVW2Njox588MGWP8+cObNlNNnQoUN11lln6ZNPPtEDDzygoUOHqrGxUTk5OTIMQz/72c8UHx/vr3IBAAAAAAAAF34Lz6qqqpSfn9+qzTCMVm1VVVWtts+fP18ZGRlasWKFtm3bpsjISA0dOlQzZsxQdna2v0oFAAAAAAAA3PJbeDZ58mRNnjw5YMcBAAAAAAAAvhY0c54BAAAAAAAAwYbwDAAAAAAAADBBeAYAAAAAAACYIDwDAAAAAAAATBCeAQAAAIBFDMOwugQAQDsIzwAAAAAAAAAThGcAAAAAAACACcIzAAAAAAAAwAThGQAAAAAAAGCC8AwAAAAAAAAwQXgGAAAAAAAAmCA8AwAAAAAAAEwQngEAAAAAAAAmCM8AAAAAAAAAE4RnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAAAAJqKsLgAAAAAAgEBxOp1as26DioqK1bdvuiZNGCeHw2F1WQCCGOEZAAAAACAsOJ1OPfTI4you2dfStmrNOt18w7UEaABM8dgmAAAAACAsrFm3oVVwJknFxSVau36jRRUBCAWEZwAAAACAsFBUVOy2vdCkHQAkwjMAAAAAQJjo2zfdbXuGSTsASIRnAAAAAIAwMWnCOKWnp7VqS09P08TxYy2qCEAoYMEAAAAAAEBYcDgcuvmGa7V2/UYVFhUro2+6Jo4fy2IBANpEeAYAAAAACBsOh0OnT51idRkAQgiPbQIAAAAAAAAmCM8AAAAAAAAAE4RnAAAAAAAAgAnCMwAAAACwiGEYVpcAAGgH4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAgEVsNpvVJQAA2kF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAAAAJgjPAAAAAAAAABOEZwAAAABgEcMwrC4BANAOwjMAAAAAAADABOEZAAAAAAAAYILwDAAAAAAAADBBeAYAAAAAAACYIDwDAAAAAAAATBCeAQAAAAAAACYIzwAAAAAAAAAThGcAAAAAAACACcIzAAAAALCIYXUBAIB2EZ4BAAAAAAAAJgjPAAAAAAAAABOEZwAAAAAAAIAJwjMAAAAAAADABOEZAAAAAAAAYILwDAAAAAAAADBBeAYAAAAAAACYIDwDAAAAAAAATBCeAQAAAAAAACYIzwAAAAAAAAAThGcAAAAAAACAiSirCzjeww8/rJ07d5pu//Wvf61TTjklgBUBAAAAAAAgnAVVeNZs9OjRcjgcLu09e/YMfDEAAAAAAAAIW0EZns2ZM0dJSUlWlwEAAAAAAIAwx5xnAAAAAAAAgAnCMwAAAAAAAMBEUD62uWrVKlVVVSkiIkKpqakaOXKkEhMTrS4LAAAAAAAAYSYow7MlS5a0+vMbb7yhmTNnatasWR6fwzAMX5cVlgzDaPkC0Dn0I8A36EuA9+hHQcjNe8H7E9zoR4D3rOpHNputU8cFVXiWnZ2tKVOmKCsrSwkJCSovL9fmzZu1ZMkSvffee4qJidHZZ5/t0bn27Nnj52rDg2EYKisrk9T5bzIg3NGPAN+gLwHeox8Fn8NHjri08VkmuNGPAO9Z1Y8GDBjQqeOCKjw777zzWv25d+/emjFjhvr376/HHntM77//vqZOnSq73d7uuTIzM/1VZlhpToEzMzP5wQB0Ev0I8A36EuA9+lHw6d69u0sbn2WCG/0I8F6o9aOgCs/MnHzyyerfv78KCgq0e/duDRkypN1jQuHmhwqbzdbyBaBz6EeAb9CXAO/Rj4KLu3eB9yb40Y8A74VSPwqZ1TZTU1MlSYcPH7a4EgAAAAAAAISLkAnPqqurJUkOh8PiSgAAAAAAABAuQiI8O3LkiHbt2iVJ6tevn8XVAAAAAAAAIFwEzZxnu3fvVl1dnQYPHtzqedfS0lI9++yzcjqdOvXUU9WrVy8LqwQAAAAAAEA4CZrwbN++fXrhhRfUo0cPpaamKiEhQRUVFSooKFB9fb3S0tJ0ySWXWF0mAAAAAAAAwkjQhGcDBgzQtGnTlJ+fr5KSEu3atUsOh0MZGRkaM2aMpk2bJrvdbnWZAAAAAAAACCNBE56lpaXp4osvtroMAAAAAAgYw+oCAADtCokFAwAAAAAAAAArEJ4BAAAAAAAAJgjPAAAAAAAAABOEZwAAAABgEZvVBQAA2kV4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAAAAJgjPAAAAAAAAABOEZwAAAAAAAIAJwjMAAAAAsIhhdQEAgHYRngEAAAAAAAAmCM8AAAAAAAAAE4RnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAgFUMw+oKAADtIDwDAAAAAAAATBCeAQAAAAAAACYIzwAAAAAAAAAThGcAAAAAAACACcIzAAAAAAAAwAThGQAAAAAAAGCC8AwAAAAArGKzWV0BAKAdhGcAAAAAAACACcIzAAAAALCKYVhdAQCgHYRnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMRFldAAAAAAAAXYHT6dSnn69Syb79yuzfT5MmjJPD4bC6LABeIjwDAAAAAMBLTqdT9/75QR0+ckSStH7DJq1as04333AtARoQ4nhsEwAAAAAAL32+ck1LcNasuLhEa9dvtKgiAL5CeAYAAAAAFjEMw+oS4CM7cnPdthcWFQe4EgC+RngGAAAAAICXUpKT3bZn9E0PcCUAfI3wDAAAAAAAL508bKhLW5/eqZo4fqwF1QDwJcIzAAAAAAC85HDYXdquvvIKFgsAugDCMwAAAAAAvGSTzaUt2h5tQSUAfI3wDAAAAAAAL9kiXMMzsR4E0CUQngEAAAAA4CWbzTU8YzVVoGsgPAMAAAAAwA8Mhp4BXQLhGQAAAAAAXopwN/KsifAM6AoIzwAAAAAA8Ja78MyCMgD4XpTVBQAAAAAAEOrcLBcg+WDOM6fTqfc/+H/6+pttampqUmREpPr0SdVlP7lIDofD6/MDaB/hGQAAAAAA3nI78sy78Ozw4SO67y9/U21tbav2/QcO6OChUt18w7UEaEAA8NgmAAAAAABe8sdqm2+9855LcNasuLhEa9dv9Or8ADxDeAYAAAAAgJf8EZ6t37i5ze2FRcVenR+AZwjPAAAAAADwkvvwzL/XzOib7t8LAJBEeAYAAAAAgNf8MfKsLenpaZo4fqzfzg/gO4RnAAAAAAB4yU125rehZ1FRUSwWAAQQq20CAAAAAOAldyPPmvwUniV070ZwZoGDB0v15jvvqrbWqZEjhmvSxHG8D2GC8AwAAAAAAC/Z5G7omX907949YNfCMVVVVbr/r39TU1OTJGln7i6tXL1Gt9x4HQFaGOCxTQAAAAAAvOUmO/PXnGcREa4f5Z1Op5Z/+rn+899XtOLTlXI6nX65drh6+bU3W4KzZiX79mvl6rUWVYRAYuQZAAAAAABecjvyzN/LbX7L6XTqbw8t1P4DByVJa9dv1MrVa3XLjb9mVJSP7Mrb7bb9iy+36OwzT/f79bfv2KmNm79UaWmpevbsqTOmnabeqSlas26DioqK1bdvuiZN4DFSfyE8AwAAAADAWwFebfN4a9ZtaAnOmpXs26e16zfq9KlTAlJDV2cYTW7bTxyN5g9Op1P/fPo51dXVtbTlbNuu2JhYHTj43fu+chWPkfoLj20CAAAAAOAH/orOTlycoKBgr9v99pi0o+Ps0Xa37T179PD7tdes29AqOJOkI0eqWgVnEo+R+hPhGQAAAADA55xOp1Z8tlL/ffk1rfgsPOfgCtTIs7r6erft9Sbt6LjIyEi37TEx/h/lVVRU7PG+X3y5xY+VhC8e2wQAAAAA+NSJc3BJ4flIWaDCM7NRUfbo6IBcPxzEx8fpUGmpS3v//v38fu2+fdP9fg20jZFnAAAAAACf+nzVGjdzcIXfI2WBCs/6988wae94sMOIQfeSEhNd2tLT0zRx/Fi/X3vShHGKjYlp1RZtEoyOGjnC7/WEo6AbeVZfX6+lS5dqw4YNKisrU3x8vE4++WTNnj1bvXr1sro8AAAAAEA7zB4dC9TKhEEjMNmZJk0Yp0VvvO3S3tFgx+l06i9/e1ilpWUtbZ99vkq33XxDWI0YdCciwnVBiJtvuDYg98XhcGjo0MHa/MVXLW3jxozS6rXrW+3XvXs3nTZ5ot/rCUdBNfKsvr5ejzzyiBYvXiyn06lTTz1VvXr10urVq/WXv/xFB0+YDA8AAAAAgGBl+Dk9O3z4sP7v3y/of/7+iNvtHQ12Pln+aavgTJL2Hziovy94TIcPH+50nV1Bk5tRhIEMFO321o/mupuDbdKEcWEfcvpLUI08++CDD5SXl6eBAwfq+uuvV8y3wxI//vhjvf7663rhhRd08803W1wlAAAAALR28OAh/fPpZ1VaVq7kpCRdfeXlSklJbve4QD3WF2ijRo5Q/p4Ct+3hxN3b63Q69fEnn2r37nwNHpyt06dO6VTgUVFRoT/d+xeffg99vPxTt+0l+/brrvv+qnvuvEMJCQk+u14gHT58WK+98baKi0uUnp6m+XN+2KHXYnVfjYxoPfapobHRZZ/GxqZAlRN2giY8a2xs1PLlyyVJP/7xj1uCM0maPn261qxZo9zcXBUUFKh///4WVQkAAAB0bYcOlWn9xo0qKytXRkZfRjJ44ODBQ7r3gQdb/lyyb5/ufeBB/en3t3sUoPmC0+nUmnUbVFRUrL590y1/36ZOmaS33nm/VVtERERQP1LWfA8LC4uUnJSoM8+YprLyCj35z6dVVXVUqSnJ+vlPL1NycpLb4+vczA1mGK3DDKfTqb/+7WEd+nZ017YdO7Vs+af63W03dSjIsdls+vfz/2030Hn62RdUXFyixMREVVRWqqKiQklJifr5FZcpJSXZ5fumttZ8frOGhgYtevMdXXnFpR7XGSwqKyv1p3v/oqamY+/HgYOHtPXrHLdhoGEYeuud91Wyb7+OVB1Rj4Qemj/nAtN7feDAQT3x1P+p8vARpSQn6eorr+hwv/ek/0acMNKsoaHB7XngHzbD6vj0W9u3b9cjjzyilJQU3XPPPS7bFy9erPfee0+zZs3S7NmzLagwPBmGoT179igzM1M2m+sz3qHG6XTq81VrWuZgGDVyhKZOmcQvhPCrUOtHzT+8Cwr2qq6+XvZou/r3z/DrL+HB9gv/icz+7nA6nXr5tTeUtztfhmEoKTFRY8eMcvv3SjC/Rm//JdbT4725B06nU2vWbtC2HTs0dPBgTZrY+tjOnPvEY0aPHKHVa9dr1ep1qquv00kDB+iieRd6dC+OP1fv3qlqaGjQ1q9zJLn/WeN0OrX808+1ctVaHa2uVnx8nE6bPFFnnj7VtG53r1E6Nin3ps1fqqrqqLp376bRo05tuZ6798bprNPTz72g0tKyVh+gpGPv5YsvL1LurjxJhhITE9U3PU0DB2S2e0+dTqc+WrZcn362StU1NS3tkZGRGjd2tHbl7dahQ9+tUta9ezddccnF+uiTFdqZu0uGYSg5OUm9U1O0Zes3rc5ts9kUGxururo6NTY0yBETo2t/caW6deumf/37Oe3ff7DlA1F7bDabIiMjFB1tV81xdTaz26OVkJCg8vJyNTU1yTCOfeD39PySFBUVpZ49e+jw4cOqq6uXzWZTn94pKtl3oN1jT7xWXFysoiIjdfhIVbuvy92v9Xa7XY2NDYxGAIAwFB0draysgRrQv5++ydmuqqqjio+PU2xcrAoKClVbW9uyb2RkpOzR0YqIiFD37t0UHR2tyMhIv31mDrXPSEETni1btkyLFi3SmDFjdNVVV7ls37Jli5588kmNHDlS11xzjQUVhqdQ+4Zui9Pp1F33/kVHq6tbtXfv3k13/eG3QfMhFl1PKPUjp9OpO+95wO0Hyj69e+vWm3y/vLzT6dRDjzyu4pJ9LW3p6WkBm4C1PU6nU3ff91dVHT3aqr179246YvJhNq1Pb91y43f3Kphf4+HDh/WHu+5v1RYREaH77vq9R6HRpi++1L+fe7FVW1RUlMu/5HpzD9o7tjPndneMO5EREbq3nXvhybmO/55wOp36w133u/3X4T69U3XrTde7DV9PvEafPr1lNDW5rGbXvO2C82bpH//6d6t2m839I0R/+v3tWrNuvf7fR5+Yvoa27qnT6dQDDz6ksrJy0+MBAEBoOvF3W18Ipc9IUhAtGFBWdmzYas+ePd1ub15ps3m/9hiGwRdfrb7WrN3gEpxJ0pEjVVqzboPl9fHFVzB8rVm7wW1wJkn79u/3S19Zs3aDS+hQXFwSNP1yzdoNLsFZ898dZkr2tb5XwfwaX3vddWWupqYmvfbG2x4d/8nyz1yOb2hocDnem3vQ3rGdObe7Y9xp9OBeeHKu478n1qzdYPpYxb79B9zW7e4a+/btdxucNW97/c13XdoNk38yffq5/2jj5i/bfA1t3dM1azcQnAEA0EWd+LttKH91VtDMedb8S+SJK0g0a2739BnePXv2+KawMGcYRktgGQppcFu27dhhum379h3K7Nc3gNUgnIRSP2qrn0j+6Stm1wyWftnePTFzfP3B/Br3FLhO5ixJe/bs9ehn6cGDhzw63pt70N6xnTl3R97X9u6Fp+dqr94T9+vMNY5XUVnh8b4HDx6SJ387md3TzvYTAAAQGnz9e6tVn5EGDBjQqeOCJjxr5qublpmZ6ZPzhLvmZDZUhlK2ZejgwS3zz5xoyJDBfM/Ab0KpH7XVTyT/9BWzawZLv2zvnpg5vv5gfo2Z/furvGKLa3tmP49qS0lJ1lE3q6mdeLw396C9Yztz7o68r+3dC0/P1V69J+7X2Xqb9ezRUwcPuQ83T5SSkqza2lo5S9se4W92TzvbTwAAQGjw9e+tofQZSQqixzaPnxfGnbq6ulb7tcdms/HFV6uvSRPHKSmxl8v3SlJSoiZNGGd5fXzxFQxfkyaOU6KbfiJJaWl9/NJXJk0cp/T0tFbXSk9PC5p+eezvjkSX+5GY2Es2m/sf9Gl9ereqP5hf4/y5P3RZ+jwyMlLz5/zQo+PnzfmhIk44PioqyuV4b+5Be8d25tzujnHHk3vhybmO/56YNHGcevXs4Xa/Pr1T3dbt7hp9+vRW79QU0+tdfumPXd4bs+/Zn19xqeZeeH6br6Gtezpp4jglJblffQ4AAIS2E3+3DeWvzmLBALTJMEJrEr/2OJ1OrVy9ttWKeadNnmj5hN3o2kKtHzmdTq1dv1F7Cvaqvr5e9uho9e/fTxPHj/Xraptr129UYVGxMvqm+/VanWH2d4fT6dQri97Urrx8yTCUlHRstU13f68E82s8fPiwFr35joqKS9Q3PU3zLjy/w6ttenK8N/egeaXJ7dt3aMiQwS4rP3bm3CceM+rU4VqzboNWrV4nZ12dsgYO0I/mXeDxapvN5+rjZrXNE78nnE6nVny2UitXrVXV0aPqFh+v0yZP1Bmnn9bmIgcnvkZJWrl6bctqm926x2vMqJEt13P33jiddXrmuf/oUGmZkpMSdeUVl7ZabfO/ryzSztw8SVJSYq9jq20OHNDuPXU6nfr4kxX69LNVreYYjYyM1Pixo7UrL7/VSLiW1TaXr9DOncdW20z5drXNr05YbTPCdmy1Tee3q23GxMToV9+utvl//35e+/YfkNHUJE9+qbXZbC0rilW7mePRYbere0J3lZdXfLvapqHIiAg1dmq1zSOqq6tThM2mPr1TVbxvf7vHnrjaZnxcnCIjI3X4yBGPr388e3SUGpsMNTY2dup4AEDoio6OVnbWScdW29x2bLXNuPhYxcXGqaBgr2q+XW3TJimiebXNyAgldO+u6OgoRUZG+e0zc6h9Rgqa8Gz79u165JFHlJKSonvuucdl++LFi/Xee+9p1qxZmj17tgUVhqdQ+4YGghH9CPAN+hLQtoMHD5kGo80C3Y86G5p7G+p7UteDDy3UgeMW3UhNTdHtN98gh8Ph0b10d87m19qzR4K+2vK1SsvK3R7vdDr1t4cWmi76IUnz516g06dO6fBrO772nj0SdPDQITU2fhfIRkZG6t4//c6n99Mdp9OphxY+oeLikpa241ft7Wi43577HnhQB9zMw5mamqKRI4brw4+/W0143NjRuuLSi+V0OvU///tIq2A/JTlZv731N6Y1StJrb7ytPXv2KjOzn84+Y5r+89Krbq8tSXMvOE+nnDxMDzz4kBoaGjx6LdKxoP/O392mhITuWrl6rT79bJVKj1s8Lzk5SXfceqMkacVnK/XxJ5+q2s3iaDO/f46ioqL07vsfeHxtd+649Ub17Zvu1Tk6y5u/D47vD4m9emrkiOEq2b9f+fkFqqmpUWxsrAYOzFR21kktfz+9896SVt8vkvTgX+7VuvUbtegN14WW5s+9QMOGDO70964nf08ePHhI//jXMzp4qFSGYbj8Q4t07B9feqem6Oorr/D42lYKtd/rgiY8a2ho0G9/+1vV1NTod7/7nfr169dq+5///GcVFRXp9ttv7/QEb+i4UPuGBoIR/QjwDfoS4D360XesHhHcfP1lyz9VqZv5BidPmqCfXDTP6+v4O4hsSyDv8dPPvtAyQvx4o0edquSkpNbh2ZhRuuKyn3S6RrN+dP1Nt7vs+/DfHlBUVFTL+5CzbYdqvx3tc7zY2FgNHTLI9H1atvxTvfn2e62OSevTW7fceF1L0PeHu+53Oe+f7/mjEhISWkYWf/3NtjZfmztWBmdWePvdxfpo2fJWbX//n/tlGEabgXAgWf33ly+E2s+joFkwICoqSmeccYY++OADvfLKK7r++utb3vyPP/5YRUVFysrKIjgDAAAA4DWHw9GpkV2+vr5hGG5Hs2T4KKxISEjQlVdc6pNzdVQg7/H8OT/Ulq3ftHpEOSoqSvMuPF/LP1vZat/jR48Eqsbm92HFZyvdvt+xMY4236ddebtd2kr27dfK1Wt19pmna7Ob4FCSvvhqq06fOkUJCQn65dVXtoR4O3PzVFVV5VHt4RScSZJhuD6mb4uIkD0qSjffcG1QhFZW//0VjoImPJOkmTNnatu2bcrLy9Pdd9+trKwslZWVKT8/X/Hx8brsssusLhEAAAAAfGbShHFatWady2iW5nkN4ZmEhATd+6ffuR1lZ1PrUS1WPnw1acI4Lf/0cx06VNqqPT4+vs3jioqK3bZ/8eUWnX3m6abbC09obw7x/vvya1q9dn0HKg8fTW6+Pxrq6xUdFUVoFcaCKjyLjo7WjTfeqKVLl2rDhg366quvFBcXp0mTJmn27NlKdLPaGQAAAACEKofDETSjWUKd2Si7YHoizOFw6I5bb9Std9zZqv3E1ZE7ymx0mNkIxnAbTdYR9fWuc9MtePQfuvk3gX88E8EjqMIzSbLb7TrvvPN03nnnWV0KAAAAAPgdo1n8LYjSM6lTAUx6eppKy8pd2keNHCGp4yMYJ00Yp89XrdE+D1YBDjfu7klxSYnWrt9IPw1jQReeAQAAAADgKyeOPPPHY5tOp9O0PSqq/Y/d7U2YnnXSQG3Z+k2rtj59euu0yRMldXwEo8Ph0K03XqeVq9e2LLQwauQIl0UJwpHZXHAnPgKL8EJ4BgAAAAAIH36Y8uzzVWvctq9es07nTD/L6/OfGMClJCfr1m9X2mzW0RGMDodDZ595us4+8/SWNsIzKS4+zm27rxbxQGjy7sFqAAAAAACC2gkLBvghPfvCZLXLL7d87fNrSVLPnj2Yf8tPUpOTXdpYxAOEZwAAAACALiuYFgzoNOsWCA077hZvuPkGFgsId4RnAAAAAIAu7MRJz3x/heGnDHPbfsrJQ31/MbU/Rxo6z92ceARnIDwDAAAAAHRZLgsG+CE9M1sUwJPFAhBc/LGgBEIf4RkAAAAAoMs6cZSWP7KR/fsPuG8/cND3F4NfNbn5Bjl8+LAFlSCYEJ4BAAAAAOCFviYrMfZNT/PJ+f0xWg7uOZ1Ol7a77vsrAVqYIzwDAAAAAHRhJ84P5vsgatKEcerZs4dL+4RxY3x+LfhXQUGhS1tDQ4MWvfmOBdUgWBCeAQAAAAC6rBPnPDt4qNTt6CJvOBwOnTv9bLftnmABgODR0Njgtr2ouCTAlSCYEJ4BAAAAALqshobGVn8uKdmnhxY+4fMALTraj4sD8NRmwJw0cIDbdl89govQRHgGAAAAAOiyCouKXNqKi0u0dv1GC6rxDUaq+c+P5l7gskpqVFSU5l14vkUVIRiwbi4AAAAAoMs6fPiI2/bCouIAV4JQkJCQoHvuvEOL3nxHRcUl6puepnkXnq+EhASrS4OFCM8AAAAAAF2WWeiRYbJCJpCQkKArr7jU6jIQRHhsEwAAAADQZX3/nLNcVsJMT0/TxPFjLaqo4wwmPQMsxcgzAAAAAECXNWBApv54x61au36jCouKldE3XRPHj/V4JcxgxIxnQGARngEAAAAAujSHw6HTp06xugwAIYrHNgEAAAAACGKGwWObgJUIzwAAAAAAsJDN1sEHMTu6PwCvEJ4BAAAAAGChpqYmq0sA0AbCMwAAAAAAAsTpdLq0FZeUuG0HEBwIzwAAAAAACJA16za4tNXV1Wvt+o0WVAPAE4RnAAAAAAD4gbu5zIqKit3uW2jS7v68nS4JQCcQngEAAAAAECB9+6a7bc8waQdgPcIzAAAAAAC85eFwsEkTxrm02e3Rmjh+rMeXKisvZ440IIAIzwAAAAAA8JZheLSbw+FwaUtPS3Pb3qy+vr7Vn/fvP6iHFj5BgAYECOEZAAAAAAAWioyMbHN7wd5Cl7bi4hIWGQAChPAMAAAAAIAgVll52G17RxYZANB5hGcAAAAAAASx1JQUt+0sMgAEBuEZAAAAAABB7KL5Fyo9Pa1VW3p6WocWGQDQeVFWFwAAAAAAAMw5HA7dfMO1Wrt+owqLipXRN10Tx49tc5EBAL5DeAYAAAAAgIVstvb3cTgcOn3qFP8XA8AFj20CAAAAAOAtTxIwACGJ8AwAAAAAAAAwQXgGAAAAAAAAmCA8AwAAAADAW4bhxcE88gkEM8IzAAAAAADgltPptLoEwHKEZwAAAAAAwK2HFj5BgIawR3gGAAAAAICFmpqarC7BVHFxidau32h1GYClCM8AAAAAAAgQd6O4CouKgnp0V2FRsdUlAJYiPAMAAAAAwFs210n/bW7aPl+1xqWtrq5eK1ev9UtZvpDRN93qEgBLEZ4BAAAAABAgX3y5pUPtwWDYkMFWlwBYivAMAAAAAACYenfxB1aXAFiK8AwAAAAAgAAZNXJEh9oDKT4+zm17UXFJgCsBggvhGQAAAAAAATJ1yiT16dO7VVufPr112uSJFlX0nZMGDHDb3jc9LbCFAEEmyuoCAAAAAAAIeYbh0W4Oh0O33nid1q7fqMKiYmX0TdfE8WPlcDj8XGD7fvyjOcq5b4caGhpa2qKiojTvwvMtrAqwHuEZAAAAAAAB5HA4dPrUKVaX4SIhIUH33HmHFr35joqKS9Q3PU3zLjxfCQkJVpcGWIrwDAAAAAAASDoWoF15xaVWlwEEFeY8AwAAAADAWzab1RUA8BPCMwAAAAAA/MBGoAZ0CYRnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAwFuGYXUFAPyE8AwAAAAAAAAwQXgGAAAAAIC3bDarKwDgJ4RnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJiIsroASdqxY4cWLFhgun3AgAG6/fbbA1cQAAAAAAAAoCAJz5qlpKQoKyvLpT05OdmCagAAAAAAABDugio8y8rK0uWXX251GQAAAAAAAIAk5jwDAAAAAAAATBGeAQAAAAAAACaC6rHNAwcO6K233tLRo0fVrVs3ZWVl6eSTT1ZEBBkfAAAAAAAAAi+owrO8vDzl5eW1aktPT9cvfvELpaamduhchmH4srSwZRhGyxeAzqEfAb5BXwK8Rz8CvGfaj9z0K/oa4J5VP49sNlunjguK8Cw2NlbnnHOORo8e3RKSFRYW6p133tHu3bu1cOFC/eEPf1BsbKzH59yzZ4+/yg0rhmGorKxMUue/yYBwRz8CfIO+BHiPfgR4z6wfHSotddmXz6WAe1b9PBowYECnjvNJePavf/1LxcXFHTrmiiuuaCm6X79+6tevX6vtQ4YM0S233KIFCxYoNzdXK1as0IwZMzw+f2ZmZofqgXvNKXBmZia/YAGdRD8CfIO+BHiPfgR4z6wf7dt/0GVfPpcC7oXazyOfhGelpaXav39/h46pq6trd5+IiAh973vfU25urr755psOhWehcPNDhc1ma/kC0Dn0I8A36EuA9+hHgPfc9iM3fYp+BpgLpZ9HPgnP7rjjDl+cxq3mxzgPHz7st2sAAAAAAOCN4P/4D6Czgn4Zy+rqakmSw+GwuBIAAAAAANxjaQCg6wr68Gzz5s2SpP79+1tcCQAAAAAAAMJNUIRnn332maqqqlq1GYahzz77TMuWLZPNZtO0adMsqg4AAAAAAADhyidznnlr6dKlevXVV5WWlqbExERJUlFRkUpLS2Wz2TRv3jxGngEAAAAAACDggiI8mz59unJyclRSUqLt27ersbFRCQkJmjBhgs4880wNGDDA6hIBAAAAAAAQhoIiPDvrrLN01llnWV0GAAAAAAAA0EpQzHkGAAAAAAAABCPCMwAAAAAAvGSzugAAfkN4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAOAlw+oCAPgN4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAAAAJgjPAAAAAADwks3qAgD4DeEZAAAAAAAAYILwDAAAAAAAADBBeAYAAAAAAACYIDwDAAAAAMBLhtUFAPAbwjMAAAAAAADABOEZAAAAAAAAYILwDAAAAAAAADBBeAYAAAAAgJdsbtqcTmfA6wDge4RnAAAAAAB4qb6+3qXtoYVPEKABXQDhGQAAAAAAXsrbne/SVlxcorXrNwa+GAA+RXgGAAAAAICXysor3LYXFhUHthAAPkd4BgAAAACAl3okdHfb3qd3aoArAeBrhGcAAAAAAACACcIzAAAAAAC8VF5R6ba9uLgkwJUA8DXCMwAAAAAAvGQzaa9zswongNBCeAYAAAAAgJeSkhLdtkdHRwe4EgC+RngGAAAAAICXMjP7u2/v3y/AlQDwNcIzAAAAAAC8NGnCOKWnp7VqS09P08TxYy2qCICvRFldAAAAAAAAoc7hcOjmG67V2vUbVVhUrIy+6Zo4fqwcDofVpQHwEuEZAAAAAAA+4HA4dPrUKVaXAcDHeGwTAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAAAAJgjPAAAAAAAAABOEZwAAAAAAAIAJwjMAAAAAAADABOEZAAAAAAAAYILwDAAAAAAAADBBeAYAAAAAAACYIDwDAAAAAAAATBCeAQAAAAAAACYIzwAAAAAAAAAThGcAAAAAAACACcIzAAAAAAAAwAThGQAAAAAAAGDCZhiGYXURAAAAAAAAQDBi5BkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAAAAJgjPAAAAAAAAABOEZwAAAAAAAIAJwjMAAAAAAADARJTVBSD4PPzww9q5c6fp9l//+tc65ZRTAlgRELoqKyv14YcfauvWrSovL1d0dLSSkpI0ZMgQzZkzx+rygKC2Y8cOLViwoN39Zs+erVmzZvm/ICBE7d69Wx9++KHy8vJUVVWlmJgY9evXT9OmTdOYMWOsLg8ICbt379bSpUu1a9cuOZ1OJSYmasyYMZoxY4bsdrvV5QFBo6CgQDk5OcrPz1d+fr4qKysVFRWlhQsXtnncmjVrtGLFCu3bt0+RkZEaOHCgZsyYoaysrABV3jbCM5gaPXq0HA6HS3vPnj0DXwwQgvLy8vT444+rpqZGaWlpGjFihJxOp0pKSrRs2TLCM6AdCQkJmjRpktttTU1NWrdunSQFzS9VQDDatGmTnn76aRmGof79+2vQoEGqrKzUjh07tH37dp177rm64IILrC4TCGrr1q3T888/r6amJvXv31+9evVSQUGBPvjgA23dulU333yzYmJirC4TCAqLFy/WV1991aFjFi1apGXLlik6OlrDhg1TQ0ODcnJylJOTo6uuukqjRo3yT7EdQHgGU3PmzFFSUpLVZQAhqaKiQo8//rgaGhr0i1/8wuUv/Pz8fEvqAkJJnz59dPnll7vd9vXXX2vdunXq1auXBg0aFODKgNDQ2Niol19+WYZh6Morr9S4ceNatuXl5emRRx7Rhx9+qNNOO00pKSkWVgoEr/Lycr344otqamrSZZddpsmTJ0uS6uvr9dxzz2nTpk168803dfHFF1tcKRAcTjrpJGVkZCgzM1OZmZm644472tx/+/btWrZsmeLj43XbbbcpNTVV0rGfUwsWLNALL7ygwYMHKy4uLhDlm2LOMwDwg7feeks1NTW68MIL3f5LyYABAwJeE9CVNI86Gz9+vCIi+HUGcGf//v2qqqpS7969WwVn0rEPN8OGDZNhGCooKLCoQiD4rVmzRvX19Ro6dGhLcCZJ0dHRuuiii2S327Vq1SpVVVVZWCUQPM4991zNnj1bI0aMUEJCQrv7f/TRR5KkmTNntgRn0rGfU1OnTlVNTY1WrVrlt3o9xW+bAOBj1dXV2rRpk2JjY3XaaadZXQ7Q5TidzpbHASZMmGBxNUDwiory7CGT+Ph4P1cChK7mcHnw4MEu27p3764+ffqosbFRX3/9daBLA0JefX29tm/fLunYtFEnap6Xc8uWLQGtyx0e24Sp5n9BiYiIUGpqqkaOHKnExESrywKC3q5du9TQ0KChQ4cqMjJSmzZt0q5du9TY2KjevXtr7NixHv0rDAD3vvjiCzmdTvXr10/p6elWlwMEreTkZCUnJ2v//v3auHGjxo4d27ItLy9POTk5SkpKUnZ2toVVAsHN6XRKkukjY83hc2FhoSZOnBiwuoCuYN++fWpoaFC3bt3Uq1cvl+39+vWTJBUVFQW6NBeEZzC1ZMmSVn9+4403NHPmTFY0A9pRUlIi6di/Rv7973/X7t27W21/5513dNlll7HCGdBJzY9sMuoMaFtERIQuv/xyPfnkk3r66af14YcfKiUlRZWVldq1a5cGDBigK664wuMRakA46t69uySprKzM7fbmdrPtAMyVl5dLktvgTJIcDodiY2NVXV2t2tpaSxfm4CclXGRnZ2vKlCnKyspSQkKCysvLtXnzZi1ZskTvvfeeYmJidPbZZ1tdJhC0qqurJUlr165VdHS0Lr30Up166qlyOp1avny5Pv74Yz377LNKTU1VRkaGxdUCoaWyslLbt29XRESEyxxOAFxlZ2frpptu0lNPPaWCgoKWR9BiYmI0dOhQVlEH2jFo0CCtX79eGzZs0OzZs1uFzbt379b+/fslSbW1tVaVCISs5pGddrvddB+Hw6Gamho5nU5LwzPmPIOL8847TxMnTlRycrLsdrt69+6tGTNm6JprrpEkvf/++6qrq7O4SiB4NTU1tfx37ty5mjJlirp166akpCTNnTtXo0ePVkNDgz788EOLKwVCz4YNG9TU1KShQ4eqR48eVpcDBL3169frwQcfVGJiom6//XY9/PDDuvvuuzVu3DgtWbJEjzzyiBobG60uEwha48ePV2JiosrKyvSPf/xDxcXFqq2t1ddff61//etfLYvW2Gw2iysFQo9hGD7ZJxAYeQaPnXzyyerfv78KCgq0e/duDRkyxOqSgKDU/C8iNpvN7dwXU6ZM0ebNm7Vz585AlwaEvOZHNplXBmjfgQMH9Pzzz6t79+669tpr5XA4JEmpqan6yU9+osrKSm3ZskWrV6/W1KlTLa4WCE4Oh0O/+tWv9OSTT+qbb77RN99807ItMTFR06dP14cffmg6JxoAc82fm9oanNO8rflnmFUIz9AhqampKigo0OHDh60uBQhazQtrJCQkKDo62nT7kSNHAloXEOpKSkq0d+9eORwOjRw50upygKC3YcMGNTY26pRTTnH7oWPMmDHasmWLduzYQXgGtKFv377605/+pE2bNqmgoECNjY3KyMjQ+PHjW+aJTktLs7hKIPQ0z3XWPPfZiZxOp2pqahQbG2vpI5sS4Rk6qHkuJ6tTXyCYNa8KU11dLcMwXIbxHz16VBL9COio5lFno0aNanNuDADHVFRUSJLpB47Y2FhJ3/1+B8Cc3W7XpEmTNGnSpFbt27ZtkyQNHjzYirKAkNa7d29FRUWpqqpK5eXlLgsH7N27V9KxANtqzHkGjx05ckS7du2S9F04AMBV3759lZSUpPr6euXn57tsb35ck34EeM4wDK1fv14Sj2wCnkpISJAk7dmzx+325vbmEdEAOmbHjh3au3ev0tLSlJWVZXU5QMix2+0t00Ft3rzZZfumTZskSSNGjAhoXe4QnqGV3bt3a/v27S6T8pWWluqf//ynnE6nTj31VNOlZAEcc+6550qSXn31VVVVVbW0FxQU6KOPPpIkTZs2zZLagFCUm5ursrIy9ejRg3/dBzx06qmnSjrWfz799NNW23bv3q1ly5ZJOvb4JgBze/fudVlYo6CgQP/+979ls9n0ox/9yKLKgNA3ffp0SdKSJUt04MCBlva8vDx9/vnniomJ0ZQpU6wqr4XNCJalCxAUVq9erRdeeEE9evRQamqqEhISVFFRoYKCAtXX1ystLU033nijunfvbnWpQFBramrSM888o02bNik+Pl4DBw6U0+nU7t271dDQoNNOO02XXHKJ1WUCIePFF1/UypUrdc4552jOnDlWlwOEjDfeeKPlH23S0tKUlpamyspK5eXlyTAMTZ06VT/5yU8srhIIbg8//LBKSkqUkZGhbt26qbS0VPn5+bLZbPrxj3/MnIHAcbZs2dIyF6Cklr6SmZnZ0jZz5sxWo8lee+01ffLJJ7Lb7Ro6dKgaGxuVk5MjwzB01VVXafTo0QF9De4w5xlaGTBggKZNm6b8/HyVlJRo165dcjgcysjI0JgxYzRt2jTmmQE8EBERoSuvvFKDBg3SqlWrtGPHDtlsNvXv31/Tpk3jsTOgA+rr61uG8k+YMMHiaoDQMmfOHJ100kn67LPPVFBQoP379ysmJkaDBg3SaaedpvHjx1tdIhD0JkyYoHXr1qmwsFA1NTXq1q2bxo0bp3POOYdpOIATVFVVuUxdYxhGq7bjn8yRpPnz5ysjI0MrVqzQtm3bFBkZqaFDh2rGjBnKzs4OQNXtY+QZAAAAAAAAYII5zwAAAAAAAAAThGcAAAAAAACACcIzAAAAAAAAwAThGQAAAAAAAGCC8AwAAAAAAAAwQXgGAAAAAAAAmCA8AwAAAAAAAEwQngEAAAAAAAAmCM8AAAAAAAAAE4RnAAAAAAAAgAnCMwAAAAAAAMAE4RkAAAAAAABggvAMAAAAAAAAMEF4BgAAAAAAAJggPAMAAAAAAABMEJ4BAAAAAAAAJgjPAAAAAAAAABP/HzqqxfHhQnvqAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# integrate over distance\n", + "vr = np.cumsum(y)[::-1]\n", + "# convert to log-log\n", + "x = np.log(x)\n", + "vr = np.log(vr)\n", + "plt.plot(x, vr, \".\")\n", + "# this part kind of needs to be done by hand (eye)\n", + "# we use the second derivative as the guide\n", + "# the second derivative close to zero indicates (local) linearity\n", + "d1 = np.gradient(vr, x)\n", + "d2 = np.gradient(d1, x)\n", + "# scale the derivative to make comparison easier\n", + "d2 /= np.max(d2[~np.isnan(d2)])\n", + "d2 *= np.max(vr)\n", + "plt.plot(x, d2, \".-\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2804 0.8947335558777068 0.8304317218117174\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABMsAAAIICAYAAABwwWmCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAABcSAAAXEgFnn9JSAABUj0lEQVR4nO3deXzdVZ0//vfNctOmULqkINgF2hQoFpFNKVVZRGgri1IYQBEUcVB+WGXccB0cncFBESgoIoICw4CCCgJt+YJYHBAERIVChaZbWovQnbZJs937+6O3IWmTNMvNvTfJ8/l45NHccz7L+yb9NO2rZ0mk0+l0AAAAAABRlO8CAAAAAKBQCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAECGsAwAAAAAMoRlAAAAAJAhLAMAAACADGEZAAAAAGQIywAAAAAgQ1gGAAAAABnCMgAAAADIEJYBAAAAQIawDAAAAAAyhGW0kk6nY9myZZFOp/NdCvRZniPoOc8R9JznCLLDswQ919eeI2EZAAAAAGQIywAAAAAgQ1gGAAAAABnCMgAAAADIEJYBAAAAQIawDAAAAAAyhGUAAAAAkCEsAwAAAIAMYRkAAAAAZAjLAAAAACBDWAYAAAAAGcIyAAAAAMgQlgEAAABAhrAMAAAAADKEZQAAAACQISwDAAAAgAxhGQAAAABklOS7AFp79dV/xoKXFkayNBnJZGmUJkubP0+WJre9btFWWpqM0tKSSCQS+S4dAAAAoM8TlhWY6pX/iN8+MLfL57UK05qDtNJIJrcHbaVRmkxGskVbc3/m19LSbW3r1q2N0tJkJMvKMn3b+ouLi3vhHQMAAAAUDmFZgWmor+/WefX1DVFf3xCxJcsFtVBUVNQiXGsZvm0b4dZqxFsnRsQlW7Vte11SUhJFRWYHAwAAAPkhLCsw9Q0N+S6hXalUKrZu3Rpbt27t1fuUlpbuMFLuzRFwpS1Gym0fDbd95FvrkXI7jKTboa24uNjUVQAAAGAnwrICM2L48HjbQQdGfX1DNDQ0RH19fdRnfm2ob4i6+vpIpVL5LrNXNTRse+9boqbX7lFUVJQJ31qOhnszVNspqGs1ei4TzrUaPdcyyHtzmqtRcgAAANC3CMsKzDsOOTjeccjBHR7T1NS0LUSrb4j6hoZoaBGotRmyNTRkpmnW73D8Dv0N286vq6uLxsbGHL3j/EilUrG1ri621tX16n1KSkraDN9ajohrvZ5cy2mrbUxzbT7+zVCupMQGDwAAAJAtwrI+qLi4OAYPHhyDBw/O+rXT6XQsX748xo4dG42NTdHQUN8qSGto8Xl9w7bRbq36G94M8baPhqtvaDuoa6ivj8ampqy/h0LS2NgYjY2NURO1vXaPRCLR5gYPba4l19EGDx1MeS3NTF0FAACA/k5YRpuaA5hkaQwZ0nv3aWpq2jayrUWA1jpca3+U3I4j4pqDupZBXqYtnU733pvIs3Q6HXV19VFX173NITqrpLi4dYC2w5pw7Y6G68KmD6WlRskBAACQX8Iy8qq4uDiKi4tj0KBBvXaPdDodjU1N0VBfH3U7jnZrdzRc56a5tjy2v09dbWxqisba2qit7b1RchHR5ii5liPidgrq2hol19E012TSKDkAAADaJSyj30skElFaUhKlJSVRXl7ea/dJpVI7j3ZruT5c87TVFlNRm0fD7Ry+NQd1Lc6tq6/v16PkIqJ5dGBs6b17FBUVvRm+tRWydbhuXBtTXnfY4KG0tKTff58AAAD6K2EZZElRUVGUlZVFWVlZr90jnU5v2+ChVZC2w5pwLUfJtbvxQztTXhvePLY/S6VSsXXr1ti6dWuv3men8K2dNeF2HA3XeqTczlNeW7YVFxebugoAAJBFwjLoQxKJRJSUlERJSUlEL2zwsF0qlYrGxsb2R8nV13dr04cdw7xUKtVr76EQNDRsCx63RE2v3aOoqGinHVLb3F211bTVHYK6DkbJbZ/mWlRU1GvvAQAAoJAIy4CdvDlNMRkRvbfDQ1NTU/PIt12tD7fjBg877q6646YPLUO8/iyVSsXWurrYWlfXq/cpKSlpJ3zr7LpxLXdqbT3ldXvAV2KUHAAAUACEZQVo7Za1cfVjs6NqzeKorJgQlx4zK0YOGdlr59bU18SchfNi0eqqmFhRGZN3Oygbb6Pb9bR3XkS02/bK6kVRXloe+44YF5P3flvMmDQtypPlUVNfE/cu+G3Mr/pDJCLimMpj4qQDTojfVz227f2Oqmw+Nl/vdbvt34eFr/096hrroqykLCbtdWDMmDQtahtq2712Z++7csPK+Pqcy2PlhpVRXFQcuw/aPfYbsV9ERKzYsGKnc1duWBmXPfD1qF5fHUWJ4jhszKHxtRO+vMv7duXrUFxcHIMHD47BLUbJ1dTXxB3P3RX3/PXXUdtQG2XFZVFWWha19bVRUlwcg5ODI9WUii0NNZFIJOJtex0UyZJkm+9h7Za18YP518bi1Uti9O6jI5oiVqxbGanGVBSliyPVkIqtDVsjkUpExeCKOLHyffHHxX+KVRtejaJ0URSliqIoXRQlURrRlI6idHEUR3EMLhocTU2pSDemoihdFMXpkoh+vkRZY2NjNDY2Rk307gYPqaJUpBOpSBelI1WUisZ0YzQVNUUqkYpEUcTew/eOulRdbKjfEKlEU6SLIprSjbGlaUs0JZoiVZSKVCIV6UQ6EsWJqBg6MvYYMjQWrVscTdEYg5JlsaluczRGxxtyDB80LPYbuW8sWbcsSoqKI5VKxfraDZFu8Y0uThTHW4a+JWrqamLD1tZ9idgW+hUlEpFKp6O0uDRGlA+PccPHRfWG6nj1jX/u8mtRnCiKw0YfGt848atd+rMEAADomUTaKtQFZe2WtXHGz8+J+qb65rZkcTLu+didu/zHUnfOramviU/dfUksXrukuW3s0DHx03NujCFlPQuQuvte2jqvJFESiaJENDQ1dNi2XWXFhPjBaVfGZ3/z+Vi6blmrvmRxadS3OKeyYkLccMZ1PQrMevJ9i2j7+7DdviPGxT82rmr1PrdfOyI6dd+VG1bG2bd9tNU/5tuy/dzahto467Zzd+ovLS6NX33srnbve9O//Cg++cuLe/R1uOCuf40VG1bu8thdvYe2auxNpUWlcfOZN8Z/zP3PqF67IorSxa0Ct9afF7dub9GXyHy+37D9oiSKo7Zua6zdtDYS2/vDyKueSkWqOVRLJVp+3hTpRHt92183ddDX8jrpyOa3qry4PIoSidjc2LmdLxKRiKIoinQiHal07093LikqidIoja2pra3+nCkpKokR5SPigwefEmceMrNLf86m0+lYvnx5jBs3zohD6CbPEWSHZwl6rq89R0aWFZirH5u90z/u65vq45rHrotvz7g86+fOWThvp4Cm+o0VMffv8+KMQ07vcv09rae98xrTjRFNscu27arWLI6rH7tup6BsWw0NOx07d+FDMfOQD7Vb06705PsW0fb3Ybtl65bv1Lb92ulId+q+X59z+S6DspbnthdWNTQ1dHjfL/z2sh5/HXoSlLW8X1s19qaGVEN853f/FVVvLI4o7fn1Xk2ujIc+9UB8fc6/x/NVf9nWmI5IpBMtQrbidsK4lq+3hXaJdsK5bZ8Xt+5L9+/1yYpi2/vvTelItwrS0jsFcE0dBG5t9BWlojGRimSirFVo114gl450NEVTzkY8NqYa2xyt15hqjNc3vx4/efLm+MmTNze3lxaVxvdP/W4cMfbw3BQIAABdICwrMFVrFnepvafnLlpd1fY5q3d9v96opzP9nbW4C9dZtKbtr0Nn9eT7FtH+92FX92wvANvxvqveWNWl666tWdut+66v3dCpetrTna9De/frTDiYbf9449WsXWtL/bYRRK2+domIdGJbCNIUTRHRS+uxpaNTo+ESOwV1xe2EcTsGdcX9fpRcIhJRnC6O4nTxtoZ2gv2e2nGU3M4j45raGf3WeiRc+33brpntb1VDqiE+e+/ns3vRHiiKRKR2+DOjpKgkDh9zWKvp5wAADAzCsgJTWTGhzZE1lRUTeuXciaMq277WqF3frzfq6ei8rppQMSGqN6zo1LETK9r+OnRWT75vEe1/H3Z1z3SkO3XffYbu0+lAsLJiQpRtKGv3+I7uO3zwsFi9Zc0u62lPd74Obemoxt701qF7xys9DF63G5LctrFCtp6HLknEtqAkUhHFvXSPdGaqYLsBXPtTVbsyos4ouexoDtR2Mf204ymsbwZz6Tb60pHdqatden9thOuNqcb40/Kn49SbZ+a0lsl7vS2u/tD3srKWJgAA3SMsKzCXHjMrnlj65E5rPn3umM/0yrkzJk2L+198sNXolbFDx8T0A6d18x30rJ72zuvOmmWXHvOZWL5+eSxZu7RVX7I42eralRUTYvqkk7r8/nZVc2e/bxFtfx+223fEuFi18dV2r92Z+35nxuWdXrPsc8d8psM1yzq67/dP/W6ba5Z15evwq+d/E9XrOxdydvQe2qqxN5UWl8a3Z1weX5vz71kZHXnVqd+NiLZ/b/ULicxUweLto+R6SRuj5BJtTT/tbDDX5qi54n49Si4iMuFj746Sa5662qnRcLs6runNkXY79BX6t2rBay/G+388I6vXPPgtk+M/P/AtI+QAADrJAv8FaO2WtXHNY9c17yb4uWM+06XdMLt6bk19Tcxd+FAsWlMVlRUTYvKQg+KACQdkZdG97r6Xts6LiHbbdtwNc/qkk5p3w7xvwf0xv+oPERFxbOV748QDToj5VX+IRWu27f65/dh8vdfttn8fdtwNc/qkk6K2obbda3f2vis3rIxvzPlWrGixG+b4EftGJBKxYv2Knc5duWFlfOWBb8Ty9dVRlCiKw8ccFl894Uu7vG82vg7/+9wv4p6//ipqGmpjUHFZlCUHRW1dTRQXl7TYDXNLJBKJmPyWt23bDbON99CyljHDRkckErF07dJobGqMRCIRDY2N8Ub9G5FOp2PciLHx5eM+Hz996mfxzMo/t6ppUHFZbG2qi0QkoqykLA7aa1KsqVkT/9iwaqevzfZdTf+y9K/x1lFvjYf+/nBsqN0Qe5TtERMq9ouVG//RfP/6hoZYX7e+1b2GJIfED07975i8z+Tmth2/psdVHhPffviKqG+qj2RxMr7x/q/Egy/NjWeq/xxN0RRFiaKoGDIyPnTwaXHYW98R35j7rVhXuz5GDB4eX3v/ZXHfgvvj76+/HI1NjdGUaop1teubg9ShyaGxtXFrNKQaoqSoOIYNHhb777l/DC3bPYYOGhq/XfBA1DZu2xEzEYkoThRFOiIGlQ6KQSWDYn3N+m0j0jL9iUhEebI8Dt77bbF2y7qoaaiJ8SP3i8qKCfHksj/FhtoNMbx8eLxr3Dtj8ZolUbWmqvnr09jUFBtrNzavhVVWUhYnT5oRq7esaXVcQ2NjbNi6ofm+2+9dUlQSb91jnxhRPiIWvv73SKVSUZ4sj021m3a9G+bg4TF+xL6xZN3SKC4qiVQ6FRtq1rcagVScKI6993hL1GytifW1G7etJ5cJ0oozIVtxFEcilYjSRGnsXjI0KgaPjPVbNsTmms3tBHWtw7mWU12bp3TSI28GbE0djJLbYf24Fn07TXVtY0RdPkfJ9cTXT/hKTD+oZ/95RGHpa4spQ6HyLEHP9bXnSFhGK33tNzAUIs8RPbVg1YL47K+/EFtTWyMiYreS3SKRjtjcsKXtnVbb2LShePvmDulE+xs6pFqMtNt+jejfU1dzodUGD7vc0GF7ANfJ9eX60Ci59uxetntcf/rVUZmlqfe0z88jyA7PEvRcX3uOTMMEgAIzeZ/J8btL5uXl3k1NTVHf0BAN9Q1RX18f9Q3bfm3I/Fpf35Dp397XEPUN9Tsc3xAbt2yIl/+5KJoaU20Gev156mqrDR5Suz6+u3bc4KFlqLbzZg87BnXtrEO3w4i6dCL7o+Q21W2K8++8sMNjEpGI/z75P2Pq+KOze3MAgE4QlgEAzYqLi2NwcXEMHjQo69euWl0Vn/v1F2P91vXbRry1M/20o00bErsYUdeqzwYPPbbzKLnWo9zaWheuwxF17YR2OwZy6UjHlx74aq++tx2VJIrjqtOujCPGHp7T+wIAhUdYBgDkROWoynjgot906thsDNVPpVKZEXHbRr/V13dmRNzGeKLqyVi9aU0kUhFFqeI3g712Nn4wSq7ndhwlt/PIuE5OU+0gnEvvYupqY7opPnvv53v8Xt6y+1vi2g99P0YPG93jawEA+SEsAwD6paKioigrK4uysrIunfcv8aFOH5tOp6OxqWlbCNcylGs1bTUTyrVoa+5v2BbQ/XXF81FXV9fuiDqj5LKjrU0bdrU2XFdHzf1z0z+bd5S2PhsA9E3CMgCAbkokElFaUhKlJSVR3vONldu1bZRcY4vgLRPE1TdEQ8MOQV19Q6xYuyLuf2FORCranuq60xTWbcf051FyEZEJHzM7yzb1zj3S0To8++73rul84NbOiLq2pru29a0qLSqN60+/utWOygBA1wnLAAAK3LZRcskoK0t2+pyPnn7OLo9ZuWFlXHzPZ2NtzdptDTtNOW1rDbniDgK3dqaqtvi8eHtY1U8loiiK00VR3Eth3HY7rh+3/fX3rruuw00bWu/CuvM6dC2PS8e2DR6Glg2N607/gRFyAAwYwjIAgAFq9LDR8dsLf5WTez1b/edta4KlIxIdjnBrMf20C6HdTteMfj51tXl6bkkvjpJ7c4OHK6+8bpej4dKJVFQMHRkfmDw9hu02LJKlpZFMJqM0Wdr8ebK0NEqTyShLJiOZLI3S0tIoKurf3ysA+h5hGQAAve6IsYfHE7Pmt9u/dsva+OK9X4mX176SnRumE1GUTuy8S2p74Vw701QTbawfZ4OHDmyMeGjFo126T6vgbccNHXYI5tLtrS2XOXb4kGHxtWmXxZ5DR2VCum2hXElJSbc3CwFg4BGWAQCQdyOHjIxbPvKTnN/3iSV/jC898NXunZyON3dLbSt0azOAK95paurOId4OAV9/3+Ahy+/xxh/+fKe2lqPkUkWpGDKoPPbaY88YVDYoM8JtW6jWejTcm22lmVFxyWSmv8VIuZaj54qL+/c0Y4CBQlgGAMCANXX80R2OeOuMNwO3hqzUtJN0tBm47WqaaqLD0K6NteUG0Ci5xs1N8Y/Nr2b9PsXFxZlwbYcwbXvQVprcoa9FUNeyf8egrkVbaWmpUXIAvUxYBgAAPbBj4Fa1uiouvuuzsSW9JTs3SGSmKkYqorcGLqW3BUpt7ZLa5tpwHWzo0NGIuv4+Sq6pqSlqm5qiduvWXr1PW8Fa82i3nUbD7dC/Y4DXYi25lqPnjJIDBjJhGQAAZFHlqMr4f595MCf3WrlhZXz2N5+Pf256rWcXSmybqthU3BRNvbVjQEQbo+RaBnDt7LTayfXlWgZz/XmUXEREfX1D1Nc3RGQpj21LUVFR6xFwbY2G2x7UtRoNt330XMsgr+3Rc6WlJTZ4AAqSsAwAAPqo0cNGx68+/os2+9LpdCxfvjzGjRvX5Wl7a7esjW88eHn87Z8vZKPMN+VilFxEZoOHrm/o0OFU1zZG1PVnqVQqtm7dGlt7eZRcaWlJlJYm46BJB8T5557Tq/cC6CxhGQAA0MrIISPjR/9yXaeOve/5++PK+Vf1ckVdlEhndszs3VFyiQ4Dt86tLzcsOSwqR0yIpsZU1NfXR31DQ9TX10dDfUPUNzREU1MvvocC0NDQGA0NjdtGygEUCGEZAADQbae9/ZQ47e2n5O3+z1b/OT577+dzf+NERDqRiqZI9Xji6l9Sf44jxhwWV5zynShPlrfqa2pqivqGhm3hWcswLfNrfSZUa2jua4j6hvodjm+IhobMsc3nbjtuW19DpNPpHr6LnkkmS/N6f4CWhGUAAECfdcTYwzu1o2lBjoBr4dl/PBf/cutH4pfn39EqMCsuLo7BxcUxeNCgXrt3Op2OxsbG1kFbizCuZdDWOmR7M6xrORpue0jXMsBrqK+Pxg5GySVLhWVA4RCWAQAA/d6uRsBdPufb8XDV73JY0c7W166PuQsfipmHfCin900kElFaum3HzCFRvusTuqmpqWnbqLYWAdr2cG233XbrtfsCdJWwDAAAGPAun/GNuDy+0fx6waoF8el7Zm3bjCCHFq2pyun9cqm4uDiKi4tjUC+OkgPIBmEZAADADibvMzn+b9aj7fY/seSP8aUHvpr1+97/4oNx/4sPxqCSQXHjmddH5ajKrN8DgI4JywAAALpo6vij210r7e7nfhXXPN653UTbs7Vxa5x/54VdOue48cfEd07+Vo/uC4CwDAAAIKvOPGxmnHnYzHb7e2uzgd8veSymzj52p/aLp1wUHznynKzfD6C/Ksp3AQAAAANJRxsN9IYfPXljTJ19bPPH3Jceyun9AfoaI8sAAAAGkO88ckV855Er2uz7yDvOjovf+6kcVwRQWIRlAAAAOfbJd14QNz19S77L2Mkdf70r7vjrXa3aJowYH1d/6HsxcsjIPFUFkFvCMgAAgBz72FHnRUQUZGC2o8XrlsSpN7e/Blt7vj3t3+P4/Y/rhYoAepewDAAAIA8+dtR5zaFZRMSjr/w+vjGv/+xmuf29CMyAvkZYBgAAUACO3/+4LgVLNfU18eXffi2eW/WXXqyqZ7457z+EZUCfIywDAADog8qT5XHdGVfv1H7HM3fGj568MQ8V7Swd6Zg6+9hWbYNKBsW1H/x+TN5ncn6KAtgFYRkAAEA/8pEjz4mPHHlO8+snlvwxvvTAV/NYUWtbG7fGRfdcEjeecb3ADChIwjIAAIB+bOr4o+OJWfPb7b/v+fvjyvlX5a6gjIvuuaT588Elg+PHZ14XlaMqc14HwI6EZQAAAAPYaW8/JU57+ymt2q79/XXxyxd+lbMaahtr4/w7L+zwmA8ddFp84YRLc1QRMJAJywAAAGjls8d9Jj573Gc6ffyO65L1ht+8dF/85qX7On38J995QavdRgE6qyjfBQAAANC3fenYz+e7hJ3c9PQtMXX2sXHqTafH2i1r810O0IcYWQYAAECPbJ/GmY+1z3Zlbe26OPXmmc2vrzz5v2Lq+KPzWBFQ6IRlAAAA9Fhba59FRCxYtaDVYv759qUHvhpfOvbzbdYKEGEaJgAAAL1o8j6T48Yzro/SKM13Kc2unH9VTJ19bNz3/P35LgUoQEaWAQAA0Ksm7zM55s96uFVbrnfcbMuV86+KeS/9v7jh7OvyWgdQWLISllVXV8fChQtj2bJlsWzZsti4cWOUlJTE7Nmz2z1nw4YNMXfu3HjppZdi48aNUVxcHHvuuWccccQRceyxx0ZpaeH8rwMAAADZ1ZkdN7//yNVd2gGzO55//YWYOvvY+PoJX4npB53Uq/cC+oashGVz5syJ559/vtPHv/baa3HVVVfF5s2bo6KiIg4++OCoq6uLxYsXx29+85t4/vnn43Of+1wUFxdnozwAAAD6oC+ccGl84YRLO3Xso6/8Pr4x71vdvtd3HrkivvPIFRERMahkUNx45vVROaqy29cD+q6shGXjx4+P0aNHx7hx42LcuHFx2WWXdXj8vffeG5s3b45jjjkmzjzzzCgq2rZ02qZNm+Kqq66KxYsXx9NPPx1TpkzJRnkAAAD0c8fvf1wcv/9xERFxxzN3xo+evLHb19rauDXOv/PCuPWcn8aEignZKhHoI7ISlp144oldOr6qqioiIqZPn94clEVE7L777vHe97437rnnnli2bJmwDAAAgC77yJHnxEeOPCciIr5471fij9VPdus65995YfPnwwcPjx+feV2MHjY6KzUChSsvu2GWlOw6oxsyZEgOKgEAAKA/+94Hr4iLp1zU4+usr10fZ912bqzcsDILVQGFLC+7YU6aNCmeeuqpmDdvXpxxxhmtpmH+4Q9/iKKiojjyyCPzURoAAAD9zPaRZpf99uvxf8se79G1zrrt3Dbbh5SWxy3n/MTIM+gH8hKWnXbaabF8+fKYP39+LFiwIMaMGRP19fVRVVUVQ4YMiYsuuij23nvvLl0znU73UrUDSzqdbv4AusdzBD3nOYKe8xzBzq445dvxbPVz8bn7Pp/1a29pqNkpSPvsuz8TZx56etbvBX1NPn4mJRKJ7p+b7oVKL7744igpKYnZs2e3e8zmzZvjlltuib///e9vFpNIxLve9a6YMWNGVFRUdOmey5Yt6265tJBOp2PdunUxYsSIHv3GgoHMcwQ95zmCnvMcQcdeWbsovvXEf+b0nhdMPj/eN/64nN4TCkE+fibtu+++3T43LyPLVq5cGT/60Y+iqKgoPvWpT0VlZWXU19fHc889F/fdd1+89NJL8fnPfz5GjRrV6WuOGzeuFyseOLZnp+PGjfOXKugmzxH0nOcIes5zBB0bN25cvP+wEyIiYsGqF+NTv7qk1+95y4Jb4+9vvBJXnvZfvX4vKCR97WdSzsOypqam+OlPfxobN26Myy67LMaMGRMREeXl5XH88cdHOp2OX/3qV3H//ffHBRdc0Onr9oUvdl+RSCSaP4Du8RxBz3mOoOc8R9A5B791ctx4xvVx0T29H5j9sfrJuPVPt8fHjjqv1+8FhaQv/UzK+W6YS5cujddffz0qKiqag7KWDjvssIiIWLRoUa5LAwAAYICavM/keGLW/OaPT76z84M3uuqmp2/ptWsDPZfzsGz9+vURETFo0KA2+7e3b9myJWc1AQAAQEsfO+q8Xg3MgMKV82mYQ4cOjYiI119/PbZu3bpTaLZ8+fKIiBg5cmSuSwMAAIBmHzvqvDj/XR+N5cuXt7nW0vcfuTp+89J9eaoO6C05H1k2fvz42H333aOuri5+8YtfRENDQ3Pfhg0b4p577omIiEMPPTTXpQEAAECnfeGES1tN3fzSsZ/v9Lmn//SsqKmv6cXqgO7KysiyF154IebOnduqrampKa688srm19OnT4+DDz44SktL45xzzomf/vSn8ac//SlefvnlGDt2bDQ0NMTSpUtj69atMWbMmDjxxBOzURoAAADkxGlvPyVOe/spza87Gnn2Ws1r8f4fz4iIiCP2OTyuPeOqnNQI7FpWwrLNmzfHsmXLWrWl0+lWbZs3b27+/B3veEd8+ctfjkceeSQWLVoUL774YpSUlMSoUaPisMMOi+OPPz6SyWQ2SgMAAIC8+O3CBzp13LOr/hxn3vLhuPuC/+3lioDOyEpYNmXKlJgyZUqXzhkzZkx8/OMfz8btAQAAoOA0pZs6feyqzati6uxjm19/6KDT4gsnXNoLVQG7kvM1ywAAAGAgKE4Ud/vc37x0X5x/+4VZrAboLGEZAAAA9IIrPvDtHp1ftb4qbnri5ixVA3SWsAwAAAB6wdTxR8eVJ/9Xj67x8z/fnqVqgM4SlgEAAEAvmTr+6Hhi1vx4Ytb8mL7/Sd26xqOv/D7LVQEdEZYBAABADnx92lfiiVnz412j39ml874x71sxdfaxzR8fvu28WLtlbS9VCQjLAAAAIId+cPqVcfGUi7p9/vIN1XHqzTNjwaoFWawK2E5YBgAAADn2kSPPaZ6e+cSs+VGU6Po/zy+655JYuWFlL1QHA5uwDAAAAPJsZPmIbp131m3nmpIJWSYsAwAAgDz7/qnf7fa5p918hsAMskhYBgAAAHlWOaoybj3np906Nx3p+Mw9l2a5Ihi4hGUAAABQACpHVcYTs+bHxw7/aJfPXb6xOn7+1G29UBUMPCX5LgAAAAB40yenfiI+OfUTrdqmzj52l+fd9PQt8bGjzuulqmDgMLIMAAAAClx3RpsB3SMsAwAAgAL3yamfEJhBjgjLAAAAoA/45NRPxPT9T2q//50X5LAa6L+EZQAAANBHfH3aV9oMzD75zgusVwZZYoF/AAAA6EO+Pu0rcfT4o+Ib877V3HbT07fETU/f0vx6+v4nxdenfSUf5UGfZ2QZAAAA9CGPvvL7VkFZW+a+8lB85NaP5aYg6GeEZQAAANCHfHPef3TquGUbl8XMm8/u5Wqg/xGWAQAAQB+SjnSnj/3nln/G1NnHxtTZx8Ydz9zZi1VB/yEsAwAAgD4kEYlunfejJ29sDs5+/tRtWa4K+g9hGQAAAPQh/zHtmz2+xk1P3yIwg3YIywAAAKAPOX7/4+Lb0/69x9dpuXsm8KaSfBcAAAAAdM3x+x8Xx+9/XPPr7z9ydfzmpfvyWBH0H8IyAAAA6OO+cMKlUbHbyC6PFps6+9jmzz/37s/EmYfNzHJl0PeYhgkAAAD9wMeOOi9uPOP6bp9/zePXxXHXnxgrN6zMYlXQ9xhZBgAAAP3E5H0mxxOz5kdExKOv/D6+Me9bXTq/PlUfZ912bkREDCkdEj8649qoHFWZ7TKhoBlZBgAAAP3Q8fsfF0/Mmh9DkuXdOn9Lw5Y4/84LY+rsY2Pq7GPj0Vd+n+UKoTAZWQYAAAD92D5D94lFa6p6fJ1vzPtWmyPV3jXmyPjOB74V5d0M5aDQGFkGAAAA/dh3Zlzeq9f/04pn4v0/nhF3P/erXr0P5IqRZQAAANCPjR42On5x3v/E+bddGFtja6/d55rHr4trHr+u3X67bdJXCMsAAACgnxs9bHT8bta8Vm03PXFz/PzPt+eshl2FaV8/4Ssx/aCTclYPtMc0TAAAABiAPjn1E/HErPkxff/CCKi+88gV8aM//DjfZYCwDAAAAAayr0/7Sjwxa348MWt+fOQdZ+e1ljv+epe1z8g7YRkAAAAQEREXv/dTzcHZtR+8Ki81XPP4dTF19rFx+i3/Eis3rMxLDQxs1iwDAAAAdnLE2MPjiVnz2+yb+9JD8Z1HrujV+7+2+fU467Zz4xfn/U+MHja6V+8FLRlZBgAAAHTJ9INOioc/NSdmHDi91+911m3nxtTZx8ZJPz45qlZX9fr9wMgyAAAAoMvKk+XxtRO/HF878cvtHvPzp26Lm56+JSv321y/Oc6/88KoKB8ZPzzjWqPN6DXCMgAAAKBXfOyo8+JjR53XZt8dz9wZP3ryxi5fc03NWtMz6VXCMgAAACDnPnLkOVFTXxM///Pt3Tr/rNvOjYiIokRRXH3a9+KIsYdnszwGMGuWAQAAAHnxyamfiF+c9z8xsnxkt6+RSqfis/d+PqbOPjaO/9FJsWDVgixWyEAkLAMAAADyZvSw0fHbC38V137wqh5fq66xLi665xKBGT0iLAMAAADy7oixh2clMIuI+PxvL8vKdRiYrFkGAAAAFIQjxh4eT8ya36rt3bOPi3Sku3SdzfWbY8GqBTF5n8lZrI6BwsgyAAAAoGD998n/2a3zLrrnkjj2+vebkkmXGVkGAAAAFKyp44+OK0/+r/jSA1/t8rkNqYa46J5Lml8nIhFHjDksvnHiV2PkkO5vKkD/JiwDAAAACtrU8Ue3mp65dsva+Or934gFr7/UpeukIx3PrPhzzPzZ2fGrj98lMKNNpmECAAAAfcrIISPjxrN/FMdVHtOt8xtSDXHNY9dluSr6C2EZAAAA0CddesysSESiW+dWrVmc5WroL4RlAAAAQJ80csjIuO8T98SYoaO7fG71hhXxgZ98MKpWV/VCZfRl1iwDAAAA+qyRQ0bGXR/7n5j70kPxnUeu6NK5G7ZuiPPvvLD59UkHvD++cNylUZ4sz3aZ9CFGlgEAAAB93vSDToonZs2PX5z3P7HnbqO6dY2HXn44zrn9vKipr8lydfQlRpYBAAAA/cboYaPjNxfc3fz6mOvfF42ppk6fv2bLmnj/j2dERMRuyd3ihzOvicpRlVmvk8JlZBkAAADQbw0fPLzb526u3xzn33mhdc0GGGEZAAAA0G99/9Tv9vgaX/ztV7JQCX2FsAwAAADotypHVcat5/w0hiWHdfsa62rXZ68gCp41ywAAAIB+rXJUZTz4qXubX3d158wRPZjKSd9jZBkAAAAwoGzfOfPKk/+rU8d/79TOB2v0fcIyAAAAYECaOv7oDgOz3ZK7xa3n/NRumAOMaZgAAADAgDV1/NHxxKz5za9r6mvipqduiV8/f29srt8cF98zK35w2pUxeZ/J+SuSnDKyDAAAACC2BWXn/e8F8cu/3hONqcaIiNjSUBMX3XNJPFv95zxXR64IywAAAAAiYs7CefHqG/9ss+/Se78YNfU1Oa6IfBCWAQAAAETEotVV7falIhVzFz6Uw2rIF2EZAAAAQERM3MVC/ovWtB+m0X8IywAAAAAiYsakabHHoKHt9t//4oMxdfax8R8P/acpmf2YsAwAAAAgIsqT5XH7R34Wg4sGd3jcQy8/HKf+dGas3bI2R5WRS8IyAAAAgIyRQ0ZGUXFil8fVNtbGR++4wAizfkhYBgAAANBCTUNtp47buHWjRf/7IWEZAAAAQAtDkuWdPtai//2PsAwAAACghatO/e9OHzuxouMdNOl7hGUAAAAALUzeZ3LceMb1URqlHR43eo+3xvRJJ+WoKnKlJN8FAAAAABSayftMjvmzHm5+vWDVgvi3+74cWxq2RElRSZzx9g/FJ476eJR3YcomfYOwDAAAAGAXJu8zOf7fpx/MdxnkgGmYAAAAAJAhLAMAAACADGEZAAAAAGRYswwAAACgEx595ffxjXnfan69e9nucf3pV0flqMo8VkW2GVkGAAAAsAs7BmUREZvqNsX5d14YVaur8lQVvUFYBgAAALAL35z3H+32nX/nhfGBn3xQaNZPCMsAAAAAdiEd6Q77N2zdYJRZP5GVNcuqq6tj4cKFsWzZsli2bFls3LgxSkpKYvbs2R2e19DQEPPnz48///nP8frrr0c6nY499tgjJkyYEKecckoMGzYsG+UBAAAA9EgiErsMzCIivvjbr8RvPnF3Diqit2QlLJszZ048//zzXTpn06ZNce2118aqVati6NChccABB0RExOrVq+PJJ5+MKVOmCMsAAACAgvAf076505plbVmzZW0OqqE3ZSUsGz9+fIwePTrGjRsX48aNi8suu6zD41OpVPz4xz+OVatWxbRp0+IDH/hAFBcXN/evWbMmBg0alI3SAAAAAHrs+P2Pi4jYZWCWilRMnX1s8+sTKo+PL5/whShPlvdmeWRRVsKyE088sUvHP/XUU7F06dI49NBD49RTT92pv6KiIhtlAQAAAGTN8fsfF2OHj4nz77yw0+c8UvVo/O3VF+J/P3qrwKyPyMsC/48//nhERLzvfe/Lx+0BAAAAuqVyVGXces5Pu3TO6i2r43+f+0UvVUS25Tws27p1ayxfvjwGDRoU++67byxZsiTuvffe+N///d+YO3durFq1KtclAQAAAHRa5ajKmFhR2aVzfvb0rXHc9SfGglULeqkqsiUr0zC74tVXX410Oh2jRo2Ku+++Ox577LFW/Q888EC8733vi9NPP71L102nd70jBbuWTqebP4Du8RxBz3mOoOc8R5AdniXa8+3pl8fZt5/bpXPqU/Vx0T2XxI9nXh+T93lbL1VWePLxHCUSiW6fm/OwrKamJiIi/vGPf8SKFSvihBNOiGOOOSbKysrib3/7W9x9993xyCOPREVFRbz3ve/t9HWXL1/eWyUPKOl0OtatWxcRPfuNBQOZ5wh6znMEPec5guzwLNGRCyafH7csuLXL5/3bfV+Mm2bc0AsVFaZ8PEf77rtvt8/NeViWSqWafz3iiCNajSCbOnVqNDY2xi9+8YuYN29el8KycePGZb3WgWh7yjtu3Dg/CKCbPEfQc54j6DnPEWSHZ4mOXDDuY3HAuAPiyw9+tUvn1TZuHVA5Rl97jnIelg0aNKj586OPPnqn/ilTpsQvf/nL2LBhQ7z++uux5557duq6feGL3VckEonmD6B7PEfQc54j6DnPEWSHZ4mOvHvC0fHErPnNr5+t/nN89t7Pd3jOkOSQAff7qS89Rzlf4H/kyJHNn48YMWKn/mQyGbvttltERGzatClndQEAAAD01BFjD49rP3hVh8dcdep3c1QN3ZHzsGzEiBExZMiQiIjYsmXLTv2pVCpqa2sjIqKsrCyntQEAAAD01BFjD48nZs2Pb0/791btg0rK4sYzro/J+0zOU2V0Rs6nYUZEHHzwwfHUU0/FK6+8Evvtt1+rviVLlkRjY2OUlpbGW97ylnyUBwAAANBjx+9/XBy//3H5LoMuyvnIsoiI97///VFUVBSPPPJIrFixorl906ZNcffdd0fEtvXMSkrykuUBAAAAMEBlJY164YUXYu7cua3ampqa4sorr2x+PX369Dj44IMjImLvvfeOmTNnxt133x3f+973Yr/99ouysrJYsmRJ1NTUxJgxY+K0007LRmkAAAAA0GlZCcs2b94cy5Yta9WWTqdbtW3evLlV/3HHHRd77bVXPPLII7F8+fJoaGiIUaNGxfHHHx8nnHBCJJPJbJQGAAAAAJ2WlbBsypQpMWXKlC6fd9BBB8VBBx2UjRIAAAAAoMfysmYZAAAAABQiYRkAAAAAZNhuEgAAAKAXVa2uii/89rJYX7s+hg8eHt8/9btROaoy32XRDiPLAAAAAHpJ1eqqOP/OC2P1ljXRmGqK1VvWxPl3Xhg/mH9t1NTX5Ls82iAsAwAAAOgln7/vy222/+r538TH7/qkwKwACcsAAAAAesm62nXt9q3c8I+4+2+/zmE1dIawDAAAAKCXJIuTHfbf/swdRpcVGGEZAAAAQC8545CZHfbXNtbGSTeeHAtWLchRReyKsAwAAACgl5x/5Lmx99C3dHhMKp2Ki+65JGb+7OxYuWFljiqjPcIyAAAAgF5SniyP2z58S4wdNmaXx/5z0z/jrNvOFZjlmbAMAAAAoBeVJ8vj+pnXdPr4L9//td4rhl0SlgEAAAD0spFDRsa1H7yqU8cuW7/c6LI8EpYBAAAA5MARYw+PG8+4vlPHfmPOt3q5GtojLAMAAADIkcn7TI6HPzUnZhw4vcPjXlmzKKpWV+WoKloSlgEAAADkUHmyPL524pfj3MM/3OFx5995ocAsD4RlAAAAAHlw/pHnxtjhHe+S+cXffiVH1bCdsAwAAAAgD8qT5XHzWTfG0EFD2z1mXe36HFZEhLAMAAAAIG/Kk+Vx+OhD2+0fMXh4DqshQlgGAAAAkFeXHjMritqJaL5w3KU5rgZhGQAAAEAejRwyMu4877ZIFiV36vvyA1+LlRtW5qGqgUtYBgAAAJBnTy1/OupT9Tu1pyMd35jzrTxUNHAJywAAAADybNHqqnb7/vHGqhxWgrAMAAAAIM8mjqpst2+v3fbMYSUIywAAAADybMakaTGoZFCbfSOHjMxxNQObsAwAAAAgz8qT5TGopKzNvpdfeznH1QxswjIAAACAAjA4ObjN9jfqN9kRM4eEZQAAAAAF4JS3faDdPjti5o6wDAAAAKAAnHnIzEhEos0+O2LmjrAMAAAAoACUJ8tj/Mj92ux769B9clzNwCUsAwAAACgQ//WB/9hpdFkiEvHtGf+ep4oGHmEZAAAAQIEYPWx03HXe7bF/xcQYkhwS+1dMjLvOuz1GDxud79IGjJJ8FwAAAADAm0YPGx0/+/BN+S5jwDKyDAAAAAAyhGUAAAAAkCEsAwAAAIAMYRkAAAAAZAjLAAAAACBDWAYAAAAAGcIyAAAAAMgoyXcBAAAAALxp5YaV8dUHvxkrNqyIQSWDY+YhH4wPH3Z2lCfL813agGBkGQAAAECBWLlhZZx127mxeO2SqG9qiDfq3oifPX1bfOIXn4qa+pp8lzcgCMsAAAAACsTX51zeZnv1+uqYu/Ch3BYzQAnLAAAAAArEqjdWtdu3aE1VDisZuIRlAAAAAAVin6H7tNs3saIyh5UMXMIyAAAAgALxnRmXt9m+zx77xPRJJ+W2mAFKWAYAAABQIEYPGx03nnF9JCLRqn3N5jVR21Cbp6oGFmEZAAAAQAG56693RzrSrdrqm+rjmseuy1NFA4uwDAAAAKCAVK1Z3KV2sktYBgAAAFBAKismtNn+j42r4qdP3RI19TU5rmhgEZYBAAAAFJBPHf3JNtub0k3xs6dvi4/c8TGBWS8SlgEAAAAUkKeWP91h/+ubXo9bn/mfHFUz8AjLAAAAAArIotVVuzzmnr/9KgeVDEzCMgAAAIACMnFU5S6PqW9qyEElA5OwDAAAAKCAzJg0LcYOH9PhMcMHD89RNQOPsAwAAACggJQny+Pms26Mj7/z/CiJ4jaPmTHppBxXNXAIywAAAAAKTHmyPC486uMx7aC2Q7GVG/+R44oGDmEZAAAAQIEaN2Jsm+2/r3osVm5YmeNqBgZhGQAAAEDBSrTb87UH/z2HdQwcwjIAAACAArV83fJ2+6o3rMhhJQOHsAwAAACgQE0cVdlu36CSQTmsZOAQlgEAAAAUqBmTpsXwwcPa7DvjkA/ltpgBQlgGAAAAUKDKk+Vx64dvjuGDh7dqHzt8bJxz2Fl5qqp/K8l3AQAAAAC0b+SQkfHL8++IuQsfikVrqmJiRWVMn3RSlCfL811avyQsAwAAAChw5cnymGnaZU6YhgkAAAAAGcIyAAAAAMgQlgEAAABAhrAMAAAAADKEZQAAAACQISwDAAAAgAxhGQAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAECGsAwAAAAAMoRlAAAAAJAhLAMAAACAjJJ8FwAAAABAx6pWV8UXfntZrK9dH8MHD4/vn/rdqBxVme+y+qWshGXV1dWxcOHCWLZsWSxbtiw2btwYJSUlMXv27E5f49prr42XX345IiKuuOKK2GOPPbJRGgAAAECfVrW6Ks6/88Lm16u3rInz77wwbjzj+pi8z+Q8VtY/ZSUsmzNnTjz//PPdPv/JJ5+Ml19+ORKJRKTT6WyUBAAAANAvfP6+L7fZfvGvPhvzLro/ypPlOa6of8vKmmXjx4+PGTNmxKc//en47ne/26VzN23aFL/+9a9j0qRJMXz48GyUAwAAANBvrKtd12Z7U7op5i58KMfV9H9ZGVl24okndvvce+65J+rr6+Pss8+Oa6+9NhvlAAAAAPQbyeJkbG2sa7Nv4Wt/z3E1/V9ed8N86aWX4plnnolp06bFqFGj8lkKAAAAQEE645CZ7fbVtROi0X15C8vq6+vjzjvvjL322ive//7356sMAAAAgIJ2/pHnRnlp2+uSlZWU5bia/i9vYdn9998fa9eujXPOOSdKSrIyGxQAAACg3ylPlsfH33V+m32T9jowx9X0f3lJqaqrq+P3v/99HHXUUbH//vtn5Zp20cyOdDrd/AF0j+cIes5zBD3nOYLs8CxRKE572ykxb+FDsXjtkua2CSPHx7QDTyr435/5eI4SiUS3z815WJZKpeKOO+6IwYMHx+mnn5616y5fvjxr1xrI0ul0rFu3bZeNnvzGgoHMcwQ95zmCnvMcQXZ4ligkX33Xl+IPKx6P5RurY9weY+O9Y94dq199Pd9l7VI+nqN999232+fmPCx79NFHY8WKFXHuuefGbrvtlrXrjhs3LmvXGsi2p7zjxo3zgwC6yXMEPec5gp7zHEF2eJYoNAdMOCDfJXRZX3uOch6WvfDCC5FIJOKpp56KP/3pT6363njjjYiI+MlPfhIlJSVxyimnRGVlZaeu2xe+2H1FIpFo/gC6x3MEPec5gp7zHEF2eJYoFGu3rI3v/u778Zd//DUiIg596zvisvd9IUYOGZnfwjqhLz1HeVmzLJ1OR1VVVbv9S5cujYiIzZs356okAAAAgIK1dsvaOP2Ws6Ix3djc9sdlT8bMn58dv/rYXX0iMOsrch6WXXrppe32ff3rX49169bFFVdcEXvssUcOqwIAAAAoXFc/NrtVULZdQ1ND/Pfvvh9XnnpFHqrqn4ryXQAAAAAAHatas7jdvieWPRlrt6zNYTX9W1ZGlr3wwgsxd+7cVm1NTU1x5ZVXNr+ePn16HHzwwdm4HQAAAMCAUlkxIVZsWNlu/zWPXRffnnF57grqx7ISlm3evDmWLVvWqi2dTrdqs/4YAAAAQPdcesys+L/FT7Q5FTOi45FndE1WwrIpU6bElClTenyd73znO1moBgAAAKB/GTlkZPz6gl/Ex+/8ZKytWbdT/9Dk7nmoqn+yZhkAAABAHzByyMj40Rmz2+xb8PpLMfelh3JcUf8kLAMAAADoI55a/nS7fd955AoL/WeBsAwAAACgj1i0uqrD/mseuy5HlfRfwjIAAACAPmLiqMoO+y3033PCMgAAAIA+YsakabHn7qPa7d9n6N45rKZ/EpYBAAAA9BHlyfK44yO3xnvHv6fN/qeqn44nlvwxx1X1L8IyAAAAgD6kPFkexUXtRzpfeuCrUbWLtc1on7AMAAAAoI/Z1dpkX/ztV3JUSf8jLAMAAADoYyorJnTYv652fY4q6X+EZQAAAAB9zKXHzIpEJNrtHzF4eA6r6V+EZQAAAAB9zMghI+Ou826PZKK0zf7vnXpFjivqP4RlAAAAAH3Q6GGj48GL7otzD/9wDCopi6JEUVQMqYhbz/lpVI6qzHd5fVZJvgsAAAAAoHvKk+Xx6an/Gp+e+q/5LqXfEJYBAAAA9HFrt6yN7/3+6njh1QUxuHRQnPK2D8SZh8yM8mR5vkvrc4RlAAAAAH3Y2i1rY+bPzo6GVENERGyojfjJkzfH/3v5kbjpX24QmHWRNcsAAAAA+rCrH5vdHJS1tGzd8pi78KE8VNS3CcsAAAAA+rCqNYvb7Vu0piqHlfQPwjIAAACAPqyyYkK7fRMr7IrZVcIyAAAAgD7s0mNmRWlR6U7twwcPi2Mr35uHivo2YRkAAABAHzZyyMj4n3N/FuWlg1u1r6/dEJ+79wtRU1+Tp8r6JmEZAAAAQB/31PKno6ahdqf2JWuXWuS/i4RlAAAAAH3cotXtL+Rvkf+uEZYBAAAA9HETR7W/kL9F/rtGWAYAAADQx82YNC3Gj9xvp/bxI/eL6ZNOykNFfVdJvgsAAAAAoGfKk+Vx45k/jPsW3B/zq/4QERHHVr43Tpt8SpQny/NcXd8iLAMAAADoB8qT5XHOYWfFOYedle9S+jTTMAEAAAAgQ1gGAAAAABmmYQIAAAD0Eys3rIyvPvjNWLFhRQwqGRwzD/lgfPiws61b1gVGlgEAAAD0Ays3rIyzbjs3Fq9dEvVNDfFG3Rvxs6dvixk/OS2qVlflu7w+Q1gGAAAA0A98fc7lbbY3pBri/DsvjA/ffn6s3LAyt0X1QcIyAAAAgH5g1RurOuxfvn55nH3bRwVmuyAsAwAAAOgH9hm6zy6PSUc6vvbgv+egmr5LWAYAAADQD3xnxuWdOm7JuqVRU1/Tu8X0YcIyAAAAgH5g9LDR8Yvz/ieGDx7W4XGpdCrmLnwoN0X1QcIyAAAAgH5i9LDR8cvz/zfGDh/b4XGL1tgdsz3CMgAAAIB+pDxZHjef9eP4t2M+GxMrKts8pr12hGUAAAAA/U55sjxmHvKh+NEZs6OyYkKrvsqKCTF90kl5qqzwleS7AAAAAAB6R3myPG4447qYu/ChWLSmKiZWVMb0SSdFebI836UVLGEZAAAAQD+2fZQZnWMaJgAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAECGsAwAAAAAMoRlAAAAAJAhLAMAAACADGEZAAAAAGQIywAAAAAgQ1gGAAAAABnCMgAAAADIEJYBAAAAQIawDAAAAAAyhGUAAAAAkCEsAwAAAICMknwXAAAAAEDvqamviTkL58Wi1VUxcVRlzJg0LcqT5fkuq2AJywAAAAD6qZr6mvjU3ZfE4rVLmtvuf/HBuOGM6wRm7TANEwAAAKCfmrNwXqugLCKias3iuG/B/XmqqPAJywAAAAD6qUWrq9psv+2Z/4ma+pocV9M3CMsAAAAA+qmJoyrbbH+jblPMXfhQjqvpG4RlAAAAAP3UjEnTYmjZ7m32LVrT9qizgU5YBgAAANBPlSfL46NHnttm38SKtkedDXTCMgAAAIB+7IOTT4nKigmt2iorJsT0SSflqaLCVpLvAgAAAADoPeXJ8rjhjOti7sKHYtGaqphYURnTJ50U5cnyfJdWkIRlAAAAAP1cebI8Zh7yoXyX0SeYhgkAAAAAGcIyAAAAAMgQlgEAAABAhrAMAAAAADKEZQAAAACQISwDAAAAgAxhGQAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAECGsAwAAAAAMoRlAAAAAJAhLAMAAACAjJJsXKS6ujoWLlwYy5Yti2XLlsXGjRujpKQkZs+evdOxqVQqlixZEs8//3xUVVXFmjVrYuvWrTFs2LA48MAD48QTT4yKiopslAUAAAAAXZKVsGzOnDnx/PPPd+rYNWvWxA9+8IOIiBg2bFiMHz8+EolELF++PB5//PF49tln4+KLL47KyspslAYAAAAAnZaVsGz8+PExevToGDduXIwbNy4uu+yydo9NJBJx0EEHxUknnRQTJ05sbm9oaIg777wznnrqqfj5z38e3/rWt6K4uDgb5QEAAABAp2QlLDvxxBM7feyoUaPikksu2am9tLQ0zjnnnPjb3/4W69ati8WLF8f++++fjfIAAAAAoFMKaoH/0tLS2HPPPSMiYuPGjXmuBgAAAICBpqDCslQqFevWrYuIiKFDh+a5GgAAAAAGmoIKy5599tnYtGlT7LbbbjF+/Ph8lwMAAADAAJOVNcuyYd26dXHPPfdERMTJJ58cpaWlXTo/nU73RlkDTjqdbv4AusdzBD3nOYKe8xxBdniWoOfy8RwlEolun1sQYVldXV385Cc/ic2bN8chhxwS733ve7t8jeXLl/dCZQNPOp1ungrbk99YMJB5jqDnPEfQc54jyA7PEvRcPp6jfffdt9vn5j0sa2pqiptuuimqq6tjwoQJ8fGPf7xb1xk3blyWKxuYtqe848aN84MAuslzBD3nOYKe8xxBdniWoOf62nOU17AslUrFz372s3jppZdi9OjR8elPfzqSyWS3rtUXvth9RSKRaP4AusdzBD3nOYKe8xxBdniWoOf60nOU1wX+77rrrnjuuedizz33jM985jNRXl6ez3IAAAAAGODyFpbdd9998fjjj8eIESNi1qxZsfvuu+erFAAAAACIiDyFZb/73e/ioYceiqFDh8asWbNixIgR+SgDAAAAAFrJypplL7zwQsydO7dVW1NTU1x55ZXNr6dPnx4HH3xwrFixIn79619HRMTIkSNj3rx5bV7z6KOPjsrKymyUBwAAADCg1dTXxJyF82LR6qqYOKoyZkyaFuVJy2G1JSth2ebNm2PZsmWt2tLpdKu2zZs3R0REbW1t8y4IS5cujaVLl7Z5zYkTJwrLAAAAAHqopr4mPnX3JbF47ZLmtvtffDBuOOM6gVkbshKWTZkyJaZMmdKpY/fff//40Y9+lI3bAgAAALALcxbOaxWURURUrVkccxc+FDMP+VCeqipced0NEwAAAIDetWh1Vdvta9puH+iEZQAAAAD92MRRbS9zNbHC8ldtEZYBAAAA9GMzJk2LyooJrdoqKybE9Ekn5amiwpaVNcsAAAAAKEzlyfK44YzrYu7Ch2LRmqqYWFEZ0yedZHH/dgjLAAAAAPq58mS5xfw7yTRMAAAAAMgQlgEAAABAhrAMAAAAADKEZQAAAACQISwDAAAAgAxhGQAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAECGsAwAAAAAMoRlAAAAAJAhLAMAAACADGEZAAAAAGQIywAAAAAgQ1gGAAAAABnCMgAAAADIEJYBAAAAQIawDAAAAAAyhGUAAAAAkCEsAwAAAIAMYRkAAAAAZAjLAAAAACBDWAYAAAAAGcIyAAAAAMgQlgEAAABAhrAMAAAAADKEZQAAAACQISwDAAAAgAxhGQAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyCjJdwEAAAAA9L6a+pqYs3BeLFpdFRNHVcaMSdOiPFme77IKjrAMAAAAoJ+rqa+JT919SSxeu6S57f4XH4wbzrhOYLYD0zABAAAA+rk5C+e1CsoiIqrWLI65Cx/KU0WFS1gGAAAA0M8tWl3VdvuattsHMmEZAAAAQD83bsTYttuHt90+kAnLAAAAAPq9RJutz696IWrqa3JcS2ETlgEAAAD0c8vXLW+z/Q9LHo9P3/MZgVkLwjIAAACAfm7iqMp2+yz035qwDAAAAKCfmzFpWlRWTGi3f+Frf89hNYVNWAYAAADQz5Uny+OGM66LiRVtjzDbUr8lxxUVLmEZAAAAwABQniyPCRXj2+z7yz/+at2yDGEZAAAAwAAxaa8D22zfVLfZumUZwjIAAACAAWLGpGkxtGz3NvsWranKcTWFSVgGAAAAMECUJ8vjo0ee22Zfe+uZDTTCMgAAAIAB5IOTT9lpZ8zKigkxfdJJeaqosJTkuwAAAAAAcmf7zphzFz4Ui9ZUxcSKypg+6aQoT5bnu7SCICwDAAAAGGDKk+Ux85AP5buMgiQsAwAAABiAauprYs7CebFodVVMHFUZMyZNM7oshGUAAAAAA05NfU186u5LYvHaJc1t97/4YNxwxnUDPjCzwD8AAADAADNn4bxWQVlERNWaxTF34UN5qqhwCMsAAAAABphFq6vabl/TdvtAIiwDAAAAGGAmjqpsu72i7faBRFgGAAAAMMDMmDQtKismtGqrrJgQ0yedlKeKCocF/gEAAAAGmPJkedxwxnUxd+FDsWhNVUysqIzpk04a8Iv7RwjLAAAAAAak8mR5zDzkQ/kuo+CYhgkAAAAAGcIyAAAAAMgQlgEAAABAhrAMAAAAADKEZQAAAACQISwDAAAAgAxhGQAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAECGsAwAAAAAMoRlAAAAAJAhLAMAAACADGEZAAAAAGQIywAAAAAgI5FOp9P5LgIAAAAACoGRZQAAAACQISwDAAAAgAxhGQAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAEBGSb4LIP+uvvrqWLRoUbv9/9//9//F2972thxWBH3Xxo0b4+GHH44FCxbE+vXro7S0NEaOHBkHHHBAnH766fkuDwrWK6+8Etdcc80ujzv55JNjxowZvV8Q9GFLly6Nhx9+OJYsWRKbN2+OQYMGxZgxY+I973lPHHbYYfkuD/qEpUuXxkMPPRSLFy+Ourq6GDFiRBx22GExbdq0SCaT+S4PCkZ1dXUsXLgwli1bFsuWLYuNGzdGSUlJzJ49u8PznnrqqXjsscfin//8ZxQXF8d+++0X06ZNiwkTJuSo8o4Jy2h26KGHRllZ2U7tw4YNy30x0ActWbIkfvjDH0ZtbW3svffecfDBB0ddXV28+uqr8eijjwrLoANDhw6No446qs2+VCoVTz/9dEREwfwFCgrVc889FzfffHOk0+kYO3ZsTJw4MTZu3BivvPJKvPzyy3HiiSfGBz/4wXyXCQXt6aefjttuuy1SqVSMHTs2hg8fHtXV1TFv3rxYsGBB/Nu//VsMGjQo32VCQZgzZ048//zzXTrnnnvuiUcffTRKS0tj0qRJ0djYGAsXLoyFCxfGhRdeGO94xzt6p9guEJbR7PTTT4+RI0fmuwzokzZs2BA//OEPo7GxMf71X/91pz/gly1blpe6oK94y1veEuedd16bfS+++GI8/fTTMXz48Jg4cWKOK4O+o6mpKe66665Ip9NxwQUXxBFHHNHct2TJkrj22mvj4YcfjqlTp8aoUaPyWCkUrvXr18cdd9wRqVQqPvrRj8aUKVMiIqKhoSFuvfXWeO655+I3v/lNnHPOOXmuFArD+PHjY/To0TFu3LgYN25cXHbZZR0e//LLL8ejjz4aQ4YMiS9+8Yux5557RsS2n1PXXHNN3H777bH//vtHeXl5LspvlzXLALLg3nvvjdra2vjQhz7U5v+E7LvvvjmvCfqL7aPKjjzyyCgq8lcXaM9rr70Wmzdvjr322qtVUBax7R8zkyZNinQ6HdXV1XmqEArfU089FQ0NDXHggQc2B2UREaWlpXHWWWdFMpmMP/7xj7F58+Y8VgmF48QTT4yTTz45Dj744Bg6dOguj3/kkUciImL69OnNQVnEtp9T7373u6O2tjb++Mc/9lq9neVvnAA9VFNTE88991wMHjw4pk6dmu9yoF+pq6trHtr/zne+M8/VQGErKencpJEhQ4b0ciXQd20Pk/fff/+d+nbfffd4y1veEk1NTfHiiy/mujTo8xoaGuLll1+OiG3LQO1o+7qaL7zwQk7raotpmDTb/j8kRUVFseeee8YhhxwSI0aMyHdZUPAWL14cjY2NceCBB0ZxcXE899xzsXjx4mhqaoq99torDj/88E79Lwuws7/+9a9RV1cXY8aMiX322Sff5UBBq6ioiIqKinjttdfiz3/+cxx++OHNfUuWLImFCxfGyJEjo7KyMo9VQmGrq6uLiGh3Ctj2sHnlypXxrne9K2d1QX/wz3/+MxobG2O33XaL4cOH79Q/ZsyYiIj4xz/+kevSdiIso9ncuXNbvf71r38d06dPt+sY7MKrr74aEdv+t/Gqq66KpUuXtur/7W9/Gx/96EftQAbdsH0KplFlsGtFRUVx3nnnxQ033BA333xzPPzwwzFq1KjYuHFjLF68OPbdd984//zzOz0CDQai3XffPSIi1q1b12b/9vb2+oH2rV+/PiKizaAsIqKsrCwGDx4cNTU1sXXr1rxupOEnJVFZWRlHH310TJgwIYYOHRrr16+Pv/zlLzF37tx44IEHYtCgQXH88cfnu0woWDU1NRER8ac//SlKS0vj3HPPjbe//e1RV1cX8+fPj9/97nfx85//PPbcc88YPXp0nquFvmPjxo3x8ssvR1FR0U7rLwFtq6ysjEsvvTRuvPHGqK6ubp5SNmjQoDjwwAPtcg67MHHixHjmmWfi2WefjZNPPrlVuLx06dJ47bXXIiJi69at+SoR+qztIzeTyWS7x5SVlUVtbW3U1dXlNSyzZhlxyimnxLve9a6oqKiIZDIZe+21V0ybNi0uuuiiiIh48MEHo76+Ps9VQuFKpVLNv86cOTOOPvro2G233WLkyJExc+bMOPTQQ6OxsTEefvjhPFcKfcuzzz4bqVQqDjzwwNhjjz3yXQ70Cc8880xceeWVMWLEiPjSl74UV199dVx++eVxxBFHxNy5c+Paa6+NpqamfJcJBevII4+MESNGxLp16+LHP/5xrFq1KrZu3Rovvvhi3HTTTc0bzSQSiTxXCn1POp3OyjG5YGQZ7TrooINi7NixUV1dHUuXLo0DDjgg3yVBQdr+Px6JRKLNtSuOPvro+Mtf/hKLFi3KdWnQp22fgmlNGOic119/PW677bbYfffd4+KLL46ysrKIiNhzzz3jwx/+cGzcuDFeeOGFePLJJ+Pd7353nquFwlRWVhaf/vSn44YbboiXXnopXnrppea+ESNGxPve9754+OGH213TDGjf9n83dTQYZ3vf9p9h+SIso0N77rlnVFdXxxtvvJHvUqBgbd8IY+jQoVFaWtpu/6ZNm3JaF/Rlr776aqxYsSLKysrikEMOyXc50Cc8++yz0dTUFG9729va/EfGYYcdFi+88EK88sorwjLowFvf+tb45je/Gc8991xUV1dHU1NTjB49Oo488sjmdZ733nvvPFcJfc/2tcq2r122o7q6uqitrY3BgwfndQpmhLCMXdi+FlO+U10oZNt3bampqYl0Or3TsPwtW7ZEhOcIumL7qLJ3vOMdHa5rAbxpw4YNERHt/gNj8ODBEfHm3++A9iWTyTjqqKPiqKOOatX+97//PSIi9t9//3yUBX3aXnvtFSUlJbF58+ZYv379Tgv9r1ixIiK2Bdb5Zs0y2rVp06ZYvHhxRLwZBgA7e+tb3xojR46MhoaGWLZs2U7926dfeo6gc9LpdDzzzDMRYQomdMXQoUMjImL58uVt9m9v3z7iGeiaV155JVasWBF77713TJgwId/lQJ+TTCabl3f6y1/+slP/c889FxERBx98cE7raouwbIBbunRpvPzyyzstord27dr4yU9+EnV1dfH2t7+93a1dgW1OPPHEiIj45S9/GZs3b25ur66ujkceeSQiIt7znvfkpTboa6qqqmLdunWxxx57+J976IK3v/3tEbHtGfrDH/7Qqm/p0qXx6KOPRsS26ZhA+1asWLHTRhjV1dXxs5/9LBKJRPzLv/xLniqDvu9973tfRETMnTs3Xn/99eb2JUuWxOOPPx6DBg2Ko48+Ol/lNUukC2WrAfLiySefjNtvvz322GOP2HPPPWPo0KGxYcOGqK6ujoaGhth7773jc5/7XOy+++75LhUKWiqViltuuSWee+65GDJkSOy3335RV1cXS5cujcbGxpg6dWp85CMfyXeZ0Cfccccd8cQTT8QJJ5wQp59+er7LgT7l17/+dfN/0uy9996x9957x8aNG2PJkiWRTqfj3e9+d3z4wx/Oc5VQ2K6++up49dVXY/To0bHbbrvF2rVrY9myZZFIJOLss8+25h+08MILLzSv5RcRzc/KuHHjmtumT5/earTY3XffHb///e8jmUzGgQceGE1NTbFw4cJIp9Nx4YUXxqGHHprT99AWa5YNcPvuu2+85z3viWXLlsWrr74aixcvjrKyshg9enQcdthh8Z73vMdaMdAJRUVFccEFF8TEiRPjj3/8Y7zyyiuRSCRi7Nix8Z73vMdUMuikhoaG5mH573znO/NcDfQ9p59+eowfPz7+7//+L6qrq+O1116LQYMGxcSJE2Pq1Klx5JFH5rtEKHjvfOc74+mnn46VK1dGbW1t7LbbbnHEEUfECSecYFkN2MHmzZt3WoomnU63ams58yYi4swzz4zRo0fHY489Fn//+9+juLg4DjzwwJg2bVpUVlbmoOpdM7IMAAAAADKsWQYAAAAAGcIyAAAAAMgQlgEAAABAhrAMAAAAADKEZQAAAACQISwDAAAAgAxhGQAAAABkCMsAAAAAIENYBgAAAAAZwjIAAAAAyBCWAQAAAECGsAwAAAAAMoRlAAAAAJAhLAMAAACADGEZAAAAAGQIywAAAAAgQ1gGAAAAABnCMgAAAADI+P8B0s92xe3+Yj8AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# find the best fit in a reasonable range assessed from the above plot\n", + "# in theory we could run this for each possible segment (computationally insane)\n", + "# but a) we want to start from index 0 (we're interested in highest trade volume the most)\n", + "# and b) we need to make sure we do not take too few points, becasue we risk overfitting\n", + "best_r = 0\n", + "best_a = 0\n", + "best_b = 0\n", + "best_end_of_lin = 0\n", + "for end_of_lin in np.arange(100, 4001, 1):\n", + " slope, intercept, r, p, stderr = linregress(x[0:end_of_lin], vr[0:end_of_lin])\n", + " if r**2 > best_r:\n", + " best_r = r**2\n", + " best_a = slope\n", + " best_b = intercept\n", + " best_end_of_lin = end_of_lin\n", + "\n", + "# print the result, we print a+1 because the slope fitted to the integral\n", + "# is a - 1 where a is the coefficient in the trade data\n", + "print(best_end_of_lin, best_r, best_a+1)\n", + "# plot the fit alongside data\n", + "plt.plot(x, vr, \".\")\n", + "plt.plot(x[0:end_of_lin], best_a * x[0:end_of_lin] + best_b, \"-\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The fit isn't perfect but it almost never is (and honestly shouldn't be, this distribution is not exactly a power-law).\n", + "\n", + "What's interesting is that the ceofficient is much lower than for all trade (which is what usually is considered).\n", + "\n", + "We're finding a = 0.83 whilst in Karpiarz et al., it's in range [1.2, 1.6]" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}