{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# LASSO when conditioning on signs and active set\n", "\n", "One of the first works in this line of conditional inference\n", "is [Lee et al.](projecteuclid.org/euclid.aos/1460381681) which\n", "considers the LASSO (squared-error loss) and conditions\n", "on the active set and their signs.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np, pandas as pd\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "%matplotlib inline\n", "\n", "from selectinf.tests.instance import gaussian_instance # to generate the data\n", "from selectinf.algorithms.api import lasso\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will know generate some data from an OLS regression model and fit the LASSO\n", "with a fixed value of $\\lambda$. In the simulation world, we know the\n", "true parameters, hence we can then return\n", "pivots for each variable selected by the LASSO. These pivots should look\n", "(marginally) like a draw from `np.random.sample`. This is the plot below." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "np.random.seed(0) # for replicability\n", "\n", "def simulate(n=500, \n", " p=100, \n", " s=5, \n", " signal=(5, 10), \n", " sigma=1): \n", "\n", " # description of statistical problem\n", "\n", " X, y, truth = gaussian_instance(n=n,\n", " p=p, \n", " s=s,\n", " equicorrelated=False,\n", " rho=0., \n", " sigma=sigma,\n", " signal=signal,\n", " random_signs=True,\n", " scale=False)[:3]\n", "\n", " sigma_hat = np.linalg.norm(y - X.dot(np.linalg.pinv(X).dot(y))) / np.sqrt(n - p)\n", " L = lasso.gaussian(X, y, 2 * np.sqrt(n), sigma=sigma_hat)\n", " soln = L.fit()\n", " active_vars = soln != 0\n", " \n", " if active_vars[truth != 0].sum() == s: # ensure we have screened for ease of interpretation\n", " projected_truth = np.linalg.pinv(X[:, active_vars]).dot(X.dot(truth))\n", " S = L.summary(truth=projected_truth)\n", " S0 = L.summary()\n", "\n", " pivot = S['pval'] # these should be pivotal\n", " pvalue = S0['pval']\n", " return pd.DataFrame({'pivot':pivot,\n", " 'pvalue':pvalue})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a look at what we get as a return value:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Index(['pivot', 'pvalue'], dtype='object')" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "while True:\n", " df = simulate()\n", " if df is not None:\n", " break\n", "df.columns" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "dfs = []\n", "for i in range(200):\n", " df = simulate()\n", " if df is not None:\n", " dfs.append(df)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAHSCAYAAADfUaMwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1iV5RvA8e8roDgQFy4UR6GZKxO3\npgnuPSFzZcNylWlqpblN0zRnmubItMy9N6KVObCs3CaOHD8V3IoDeH5/PB4PR0hR4bznHO7PdXHF\n+7xvnBtF7vOs+zGUUgghhBDCPGnMDkAIIYRI7SQZCyGEECaTZCyEEEKYTJKxEEIIYTJJxkIIIYTJ\nJBkLIYQQJnM364Vz5MihChYsaNbLCyGEEHa1Z8+eSKWUT2L3TEvGBQsWJDw83KyXF0IIIezKMIyT\n/3VPhqmFEEIIk0kyFkIIIUwmyVgIIYQwmSRjIYQQwmSSjIUQQgiTmbaaOimuXbvGhQsXuHfvntmh\niEfw8PAgZ86cZM6c2exQhBDCKTlsMr527Rrnz5/H19eX9OnTYxiG2SGJRCiliI6O5syZMwCSkIUQ\n4ik47DD1hQsX8PX1JUOGDJKIHZhhGGTIkAFfX18uXLhgdjhCCOGUHDYZ37t3j/Tp05sdhkii9OnT\ny3SCEEI8JYdNxoD0iJ2I/F0JIcTTc+hkLIQQQqQGkoyFEEIIk0kyTkGDBg3CMIwHH3nz5qVFixYc\nO3YMgI4dOxIQEJDsr3vkyBEGDRrElStXkv1rCyGESH4Ou7XJVXh7e7Nu3ToAIiIiGDBgAIGBgezf\nv58BAwYQHR2d7K955MgRBg8eTMeOHcmSJUuyf30hhBDJ67HJ2DCMmUBD4IJSqkQi9w1gPFAfuAV0\nVEr9ntyBOit3d3cqVqwIQMWKFfHz86NatWqsWbOGVq1amRydEEIIR5CUYerZQN1H3K8H+N//eAf4\n+tnDcl1ly5YF4MSJEzbD1MePH8cwDFavXm3zfGxsLLlz56Z///4P2kJDQ6lQoQKenp7kypWLLl26\ncOPGDQDCwsJo1KgRAIUKFcIwDAoWLGiH70wIIVzH1atwcuf/ICbGLq/32GSslNoGXHrEI02A75S2\nA8hiGEae5ArQ1Zw4cQKA3Llz27QXKlSI8uXL89NPP9m0b926lfPnzxMSEgLA/v37qVu3Ljly5GDx\n4sUMHjyY+fPn07JlSwBefvllxowZA8CSJUv47bffWLp0aQp/V0II4SKU4u7mn9njH0Leivk59tVK\nu7xscswZ+wL/xrs+fb/tXDJ8bRuOsJVVqSf/f2Luv7OKiIigS5cueHl5ERQUxObNm22eCwkJYfDg\nwdy5c4d06dIBsGDBAooXL06JEnqGYOjQoRQoUIAVK1bg5uYGQLZs2QgODua3336jUqVKFC1aFIAy\nZcpIr1gIIZLi+nX4/nvUlCmwbx+LgcLAyX5TyNi2GQ/1n5KdXVdTG4bxjmEY4YZhhF+8eNGeL22a\nqKgoPDw88PDwoGjRokRERLBgwQLy5Ek4eNC6dWuuXbv2YMFXTEwMS5YsITg4+MEzu3btolmzZg8S\nMUCLFi1wd3fnl19+SflvSAghXM2KFVC4MHTpgrFvHz2BKcAWwD//HXJnuZ3iISRHMj4D5I93ne9+\nWwJKqW+UUgFKqQAfH59keGnH5+3tze7duwkPD+f06dOcOHGCevXqJfqsr68vVatWZcGCBQBs3ryZ\nyMjIB0PUAOfOnSNXrlw2/5+bmxvZs2fn0qVHzSYIIYSwcfs2dO8OTZpAZOSD5q6kpwOvcqHNn+SL\n2AaenikeSnIMU68AuhmG8SNQAbiqlEr2IWp4uiFis7m7uz/RXuLg4GD69etHdHQ0CxYsoEyZMvj7\n+z+4nydPngQHMsTGxhIVFUW2bNmSLW4hhHBpBw5ASAj8/TcAN4HPycwVhjCXjtRt7c1Hc+03PfrY\nnrFhGD8AvwFFDcM4bRjGm4ZhvGsYxrv3H1kDRAD/ANOBLikWbSrQqlUroqOjWbp0KUuXLrXpFQNU\nqFCBpUuXEhsb+6BtyZIlxMTEULVqVQDSpk0LwO3bKT+0IoQQTuXCBejZE15++UEivg6UJzvDucFk\nyhNQ05vvvoM0dpzIfWzPWCn12mPuK6BrskWUyuXMmZMaNWrQu3dvrly5QuvWrW3u9+/fnzJlytC0\naVPee+89Tp8+Td++falTpw6VKlUCeLCAa9q0aYSEhJAhQwZKlixp9+9FCCEcxuXLMGYMjB8PN28+\naL6UNh0lYvJyLu4UMJ+XXqrE0qVwfw2t3Ug5TAcUEhLCuXPnqFixYoLV0MWLF2ft2rVcuHCB5s2b\n079/f1577TUWLVr04JkCBQowZswYlixZQpUqVR7sOxZCiFTnwAE9L1ygAIwYYZOIL5Ypy3PG85yL\nOw0sJE+eYFavhsyZ7R+moUyaiA0ICFDh4eH/ef/gwYMUK1bMjhGJZyV/Z0IIhxAXB0uXwqRJEBaW\n8H6pUqghQwmakI3Q0EbAXDw8GrJ1K9wfYEwRhmHsUUoluohIalMLIYRwDUrpbUr9+8O+fQnvv/AC\nDBzInSZNmDApPaGhAMeBLEyenLKJ+HEkGQshhHBucXGwcSN89hns2mV7z81Nb13q0gVq1uTsuXOU\n9y/L2bMfAW8AWXj3XXj7bTMCt5JkLIQQwjkdOQJz58L338P9UsMPZMyo54q7doV8+QD4999/KVeu\nJufP/w94HoAqVfSaLrNJMhZCCOE8bt2CH3+EadMS9oJBL4Pu0gX69YOcOR80Hz9+nHLlahIVdQnY\nAFSiWDFYvBju7wY1lSRjIYQQju/oUZg6FWbN0tuUHpYlC7RtC337PugJW1y9epWyZatz+fINYDMQ\nQKlSemQ7Xr42lSRjIYQQjmvHDhg1CpYtS3jPwwMaNIB27fR/E9kcrBRMnuzN5csfAjWAl3j5Zdiw\nAbJnT+ngk06SsRBCCMeiFKxbp5Pw1q0J7xcuDO+9Bx07Qo4c//ll/vprP0OH3mLRonLABwBUqKC/\ndJYsKRP605JkLIQQwnEcOaLnfB86YhaA+vWhWzeoU+extSp37fqTatWCuHs3J/AX4EbNmnr7sRlF\nPR5HkrEQQgjz3bkDI0fqKll371rb3d3h9dfho4+gePEkfamwsHBq1apNTExGYBngRnAwzJlj/zKX\nSSXlMFPQoEGDMAzjwUfevHlp0aIFx44ds8vrG4bBpEmT7PJaQgjx1DZuhNKlYdAgayJOk0ZvSzp2\nDGbPTnIiXrlyB4GBgcTEeAPbAH/efx/mz3fcRAzSM05x3t7erFu3DoCIiAgGDBhAYGAg+/fvJ2PG\njCZHJ4QQJjpyBHr3hpUrbdvLldNbl8qUeaIvd/gwvPbaZOLifIBQwI9Ro3Sn2l5HIT4tScYpzN3d\nnYoVKwJQsWJF/Pz8qFatGmvWrKFVq1YmRyeEECa4eBE+/xwmToSYGGu7l5cepn7vPV056wn89pui\nUSODmzdnAFdwd8/Ft99C+/bJG3pKkWFqOytbtiwAJx6uFnNfoUKF+OijjxK0t2rV6sF5xTdv3qRb\nt24ULVqUDBkyUKhQIbp27cq1a9ce+doFCxakd+/eNm2zZ8/GMAxu3LjxoO3SpUu888475MqVC09P\nTypXrszOnTuf5NsUQghbSsFvv+ltSPnywbhxtom4Y0fdte3W7YkS8Z078M47G6hatSJRUZFAOjJk\nyMXKlc6TiEGSsd1ZknDu3LkTvd+6dWsWLlxo03bjxg1Wr15NSEgIALdu3SI2Npbhw4ezdu1ahg4d\nSmhoaLL0tO/cuUNQUBCbNm1i9OjRLFu2DB8fH4KCgvjf//73zF9fCJHKHDsGY8dC2bJQubIuXRl/\ngVbVqhAerot55MmT5C+rFCxYAH5+q5k+vRFxcXcARY4csGUL1K2b/N9KSpJhajuIuf/uLyIigi5d\nuuDl5UVQUFCiz4aEhPDFF1+wY8eOB8PbK1eu5O7duw+SrY+PD19//bXN1y9UqBBVq1bl1KlT+Pn5\nPXWs33//Pfv27WP//v34+/sDEBQURNGiRfnyyy8ZPXr0U39tIUQq8fff8NNPulBHYqcngd7w27s3\ntGjxxBO6+/fDm2/Czp3LgNZAKWADRYtmY8UKKFLkWb8B+3OuZOwIM/BPeP5zVFQUHh4eD679/PxY\nsGABuXPnfpCkQa98dnNzo0yZMhQpUoQFCxY8SMYLFiygevXq5MqV68Hzc+fOZezYsRw9epSb8Q7L\nPnLkyDMl402bNlG2bFkKFSpkE1/16tV51PnTQohU7tw5vWR57lz488/En/H0hDZt9D7i+1N2T+q7\n7+DddyE6ehXQCggga9a1DB6chc6dHaPO9NNwrmTshLy9vdm0aROGYZA7d27y5s2LYRiEhYXx6quv\nPniuevXqhN0/BDs4OJiZM2cyduxYrl+/zrp165g4ceKDZ5cuXUr79u157733GDFiBNmyZePcuXM0\na9aM27dvP1O8kZGR7Nixw+YNhMVzzz33TF9bCOGCLOUqV6zQRxk+LF06qF1bH2PYrBlky/ZULxMd\nDT16wIwZlpaXSZOmDd27T2TQoMwOV1HrSUkyTmHu7u4EBAQkaC9btiy7d+9+cO3l5fXg8+DgYIYO\nHcovv/zC8ePHiYuLo3nz5g/uL1y4kAoVKjBlypQHbVsTKxn3EE9PT+7Gn6sBLj9UcD1btmwEBATY\nDINbpHPkTXpCCPtRCtav10U6Evvd4+mpk2/r1joRZ8r0TC8XEaFHs/fuBdgCVOOFF/KyaNGcpG4/\ndnjOlYyfcIjYkXl5eSWapAGKFy9OiRIlWLBgAcePHycoKIjs8SqaR0dHJ0iM8+bNe+xr5suXj4MH\nD9q0bdiwweY6MDCQDRs24OfnR05HOc5ECGGumzchLAx279bHFu7eDZGRCZ+rXl2vlm7ZEry9k+Wl\nt27ViTgqCuAboDNlyoxm27bez5rjHYpzJeNUJDg4mPHjx3P16lWmT59uc69WrVp07dqV4cOHU6FC\nBdasWcPmxOq4PqRZs2Z0796dESNGUK5cORYvXsz+/fttnmnfvj1Tp06lRo0a9O7dm8KFCxMVFcWu\nXbvInTs3PXv2TNbvUwjh4JYtg86d4cKFxO+7u+t54I8+ghIlkvWlp0/X08t6+cokoDslStTn11+7\nkT59sr6U+ZRSpnyULVtWPcqBAwceed8ZDBw4UGXPnv2p/t+jR48qQKVLl05duXLF5l5MTIzq1auX\n8vHxUV5eXqp58+Zqx44dClArV6588BygJk6c+OD67t27qmfPnipXrlwqS5YsqkePHmratGkKUNev\nX3/w3JUrV1SPHj1Uvnz5lIeHh/L19VXNmjVTv/zyyyNjdoW/MyHEfZcuKdW2rVJ6TDLhR5YsSvXo\nodTJk8n+0vfuKfX++/Ff7ksFqFdeaaJu376d7K9nL0C4+o+caCiThn4DAgLUo1bnHjx4kGLFitkx\nIvGs5O9MCBcQG6t7wz16wNmz1vY8efQccLly+uP55x97ctLTUErX//juO0vLGdKkKUq9evVZunRe\nootLnYVhGHuUUonOT8owtRBCCF2icuZMmDoVHq4Q2LYtTJgAWbOmeBhffRU/EUPz5r706fMbZcsW\nw93ddVOW635nQgghHu/MGfjss4SVsQBy5tQHNjRtapdQQkP11DMooD8VKuRh4cJupElT0i6vbyYp\nhymEEKnRrVswdKguVzVzpm0izp4d+vTRpa7slIhPnNCj4LGxCugNjKBUqf0YhuvsonkU6RkLIURq\ncvo0rFkDw4bBv//a3qtQQS9fbt1a7xW2k1u3dD2QqCgF9AAm0alTd6ZNG4/hCJUX7UCSsRBCuLK4\nOFi7Flatgs2b4ejRhM+ULq1PUYpXFdBerl2D4GDYu1cB7wHTeO21XsyYMTrVJGJw8GSslEpVfxnO\nzKxV+UKI/3D3rp4HHj0aDh1K/JmcOWH4cHjjjSc+Pzg5RERAo0Zw4ACAARSlbt2PmTdveKr73e+w\nydjDw4Po6GgyZMhgdigiCaKjo516y4EQLuPWLb0ieuxYvTjrYZ6eUKUK1Kmji3lkzmz/GIlfWSsG\nOAK8SP/+PRk61JRwTOewyThnzpycOXMGX19f0qdPn+reJTkLpRTR0dGcOXPG5lQpIYSd3bmjT1EY\nNgwePns8c2Z4+21o2BAqVrTrfPDDYmL0Lqm+fSEm5h7wOrCOiRMP061b0s8zdjUOm4wz33+3dvbs\nWe7du2dyNOJRPDw8yJUr14O/MyGEHZ09CytXwuefw8mTtvdy54aePXUPOJlqRT+LPXvgnXfg998B\n7gDBwHJ69PgyVSdicOBkDDohyy94IYSI5+ZN2LABNm3SG3MTmw/Omxf699dzwSb2gi2uX4cBA2Di\nRMspi7eBFsAaBg+eyGefdTM3QAfg0MlYCCEE+pCGlSth+XLYuBH+69zyHDng44/hvfdwhJMUlNKV\nNbt3t52+dnefQGzsWiZNmkaXLu+YF6ADkWQshBCOSCnd850yRSfh2NjEn0ubVi/IatBAjwHHOxvd\nTKdOQbdu+j1EfLVrw/jxPTlzpiyBgYHmBOeAJBkLIYQjuX4dZs3SSfjw4cSfKVFCL8YKCoLKlR2i\nF2yhlA69b189om7h43ONYsXeZ86ckeTOnYsXXpBEHJ8kYyGEcAQXLuhlxpMnw5UrCe9Xrqz3AjVp\nAs89Z//4kuDaNXjrLVi40La9Y8cr/P13XbZv38OePS1p0KCBOQE6MEnGQghhpn//hZEjdX3oh+eC\nvbygQwc9B/zii+bEl0R//QUtW9oW+CpeHMaMuUT//rX566+/WLRokSTi/yDJWAghzHD+PIwYoQt0\nPHxa0vPPw4cfQrt2kCmTOfE9gXnzdI84/nuJLl3g448v0rBhLQ4dOsSyZcuoX7++eUE6OEnGQghh\nTxcu6EN7x4/X1bLiK1tWT7Y2b25KecqnsXKlfs9gqYibMSNMnw6vvQYXLijc3d1ZuXIltWrVMjdQ\nByfJWAghUppSsH27Xtm0cCE8XMioQgUYMgRq1QInqja4d69OupZE/OKLsHgxZM16nnv3spEzZ052\n7dpFmjRyWu/jyJ+QEEKklNhYmD8fXnoJqlbVn8dPxKVL667lb7/pPT9OlIjPntWHPFhWTBcqBGFh\nkCHDKapUqULnzp0BJBEnkfSMhRAiuVmqXQwYAPv3J7xfuTJ88IFeHe2EyermTWjcWB+NDLr09apV\ncOPGcWrWrMnly5cfJGORNJKMhRAiOYWG6nnf8HDb9gwZ4PXX9cqml14yJ7ZkEBcH7dvrOtOgp7YX\nLQIPj6O88kpNbt26xebNmylbtqy5gToZScZCCJEc/vkHevfW1bLiy5RJ94I//BCyZjUntmQ0bBgs\nWWK9njQJataMpWTJJty+fZvQ0FBKly5tXoBOSpKxEEI8i6tXdYYaP952PjhdOujaFfr1Ax8f8+JL\nRitWwMCB1uv334d33wVwY+bMmXh5eVG8eHGzwnNqzjdZIYQQjiA2FqZNA39/GDPGNhF36KB7yl9+\n6TKJ+NAhaNvWel2zJrRtu5dJkyYBULFiRUnEz0B6xkII8aRCQ/XQ899/27ZXrqz3EJcrZ05cKeTq\nVWjaVJfNBihQAD7+OJzatWuTKVMm2rdvL8fdPiPpGQshRFIopffu1K0LgYG2idjPD378EX75xeUS\nsWXBluXMivTpYdCg32jRIhBvb2+2bdsmiTgZSM9YCCEeJS5OL8oaNQp27rS9lzGjnhPu1cuhTk5K\nTkOG6Llii969t9G9ewNy585NaGgo+fPnNy84FyLJWAgh/su+fdCxo3Ufj0WaNLoG5IgRkDevKaHZ\nw7JlMHiw9bpXL8if/zD58uVj8+bN5HXh793eDGWpY2ZnAQEBKvzhfXhCCOEIYmJg9GgYNMj2EId0\n6XRy7t1bH+bgwg4ehPLl4cYNfV2jxjU2bsyMuzvcvn0bT09PcwN0QoZh7FFKBSR2T+aMhRAivkOH\ndOnKTz6xJuJ06XQhjxMn9ClLLp6ILQu2LIk4V65V7N1bkPDwHQCSiFOADFMLIQTo8/9GjoTPP7ft\nDZcrB3PmQLFi5sVmRzExegvTkSP6Om3apURFBfPSS6UpUqSIucG5MEnGQggRGgrvvWfNQAAeHnrC\n9KOPwD11/KqMi4M339R1prWfiIlpQ/ny5Vi3bh3e3t5mhufSUsdPmBBCPOzyZVi9Gn76SZ+cFF+F\nCvpQ3pIlzYnNBErpstnffWdp+RXDeI0qVaqwevVqvLy8zAzP5UkyFkKkHrdvw/ffww8/wNatuopW\nfJkz62Hqzp31CQiphFK6dPa0ada2Tp0q4u8/gu7du5ExY0bzgkslJBkLIVzf1avw9de6Otb584k/\n06qVri+dJ499YzOZUvDpp/qPRptLs2av8s03+XBz62tmaKmKJGMhhOu6eVP3dCdOhGvXEt4vX14v\nG27aNNUs0IpPKb1ofORIS8sE4H3y5OmGm9tEEyNLfSQZCyFc06pV+tSkU6ds2319dV3p117Tn6dS\nSukiHuPGWVpGA31o0qQZ48Z9aWJkqZMkYyGEazlzBnr0sD10F6BoUb1X+PXXIW1ac2JzEHFx0K2b\nHrnXhgP9adUqmHnz5uLh4WFidKmTJGMhhGu4ckUfZfjVV3p42iJ7dn2UYbt2uoxlKhcdrQcMZs16\n0IK39080bNiO2bNn4p5KtnE5GvlTF0I4txs3YMIEXb7yyhXbe506wRdf6IQs2LhRb6c+dgxAATG0\naZOe8eO3kjWrF26paAW5o5FkLIRwTseP6704M2ZAVJTtvZIl9aKt6tXNic3BnD+vty7Nn29pUUAv\n8uf/h2+/XYynZxYToxMgyVgI4UzOn4ft2+Hbb2HNGr0KKT5/f101KzhYhqTv++03aNgQLl2ytMSR\nNm0P7t6dTNOmPUiXTtKAI5C/BSGEY4qNhb17YcsW2LEDdu9OuDLaokABGDAAOnRINaUrk2LrVmjQ\nIP4UehzPPdeZY8dm0Lt3b7744gsMwzAzRHGf/NQKIRzHtWt6LHX9eggLSzgH/LC6dXUNx/r1U1XF\nrKTYtAkaN9YLtgB8fKBKlZ4sWzaDTz/9lKFDh0oidiCSjIUQ5jt/Xi/CmjxZV8v6L56eUKaMngt+\n802XP8rwaa1ZA82bw507+jpPHti8GW7dak+lSr706dPH3ABFApKMhRDmiYjQ25FmzrRmjvhy54bA\nQHjlFV0tq3hxfZqSSJRSek1bjx5w755u8/W9R9++yylWrCVQlrJly5oao0hckpKxYRh1gfGAGzBD\nKTXyoft+wBwgy/1n+iml1iRzrEIIV7F3L4wapU9Miouzvefvr4ee69SBF14AGUpNkhs39PkW1hXT\nUKDAHfz9W9OjxwoqVdpNQECAeQGKR3psMjYMww2YDNQCTgO7DcNYoZQ6EO+x/sBPSqmvDcN4EVgD\nFEyBeIUQzuzQIV2Kcv36hPfKloV+/aBZM5n/fUIHDkDLlnDwoLWtZMlocuRowaZNa5k8ebIkYgeX\nlLX/5YF/lFIRSqm7wI9Ak4eeUUDm+597A2eTL0QhhEtYuBDKlUuYiGvV0quNdu/WGUUS8RP5/nv9\nxxo/EXfseAsfn8aEha1j+vTpdOnSxbwARZIkZZjaF/g33vVpoMJDzwwCNhiG0R3ICAQlS3RCCOd3\n757u8Y4da21Lk0Yn3j59dI9YPLHbt/UgQ/wziNOn1/WmfXzCaNIkjFmzZtGhQwfzghRJllwLuF4D\nZiulvjQMoxIw1zCMEkopm8kgwzDeAd4B8PPzS6aXFkI4rDNn9OlIP/9sbXv+eVi0CEqXNi8uJxcR\noY9f/v13a1vRovDTT4pSpQygPocPH6Zw4cKmxSieTFKGqc8A+eNd57vfFt+bwE8ASqnfAE8gx8Nf\nSCn1jVIqQCkV4OPj83QRCyEcX2wsTJqkzwiOn4gbN9bD0ZKIn9ry5fDyy7aJODgYNmy4TPfuNdi4\ncSOAJGInk5RkvBvwNwyjkGEYaYEQYMVDz5wCAgEMwyiGTsYXkzNQIYST+OMPqFQJuneH69d1W5o0\nMGIELF0KWaQO8tO4dw9694amTa1bsT089HueSZOiaNo0kN9++41oS5UP4VQeO0ytlIoxDKMbsB69\nbWmmUmq/YRhDgHCl1AqgFzDdMIye6MVcHZV6uGisEMKl3bgBn30G48fbblcqWhS++UbvFRZP5cwZ\n3fv99VdrW4ECek1cgQIXqFkziCNHjrB8+XLq1atnXqDiqSVpzvj+nuE1D7V9Fu/zA0CV5A1NCOE0\nli/Xp9WfPm1tS5cOPvkE+vbVn4un8vPPuppWZKS1rWFDmDMHDOMyVarU4MSJE6xatYqgIFk766yk\nApcQ4ulFROiz+ZYvt22vWVMv6y1SxJy4XMTSpXr9m6U4mWW0/6OP9Odxcd4EBQXRokULqstxkU5N\nkrEQ4sn9VwUtHx+9hen116Vy1jOaOhW6drX+8ebMqf+4q1eHkydPopSiYMGCTJgwwdxARbKQZCyE\nSLpff4Vhw2DduoT33npLJ+hs2ewflwtRCgYNgiFDrG3PP69rpRQuDBEREbz66qtkz56dPXv2yMlL\nLkKSsRDi8fbsgf79E0/CtWrBwIFQRZaNPKuYGF2We/p0a1u5crBqle4ZHzlyhJo1axIdHc3SpUsl\nEbsQScZCiMQppfcEjxoFS5bY3jMMXXVCKmglm+hoCAmBFfE2jtapo+ujZMoEBw4cIDAwkNjYWLZs\n2UKpUqXMC1YkO0nGQghbJ07ogsdz58KRI7b3DAPatoUBA/TpSiJZXLoEjRrB9u3Wtnbt4NtvrSdG\n9unTB6UUYWFhvPjii+YEKlKMJGMhhHbypC52vGxZ4vdbtNATmZIIktXp07oHfCDeOXgffQQjR+oV\n0xZz584lMjISf3kT5JKSUoFLCOHKYmLgyy91kn04EWfKBB066DnjRYskESezq1cTJuKxY+GLL3Qi\n3rVrF8HBwdy+fZusWbNKInZh0jMWIjXbvl2vGPrzT9v2evX0OGmTJpAhgzmxubiYGGjd2pqIPTx0\nIY/XXtPX27dvp27duuTIkYOoqCh8fX3NC1akOEnGQqRGf/6p531XrrRtL1FCn8lXubI5caUSSkGP\nHrBhg7Vt1ixrIt62bRv162Owzb8AACAASURBVNcnb968hIaGSiJOBWSYWojUQildrCMkBF56yTYR\np0+vJyl//10SsR1MmKALlFl89pmukwIQGhpK3bp18fPzY+vWreTLl8+cIIVdSc9YCFcWE6MLdSxf\nrueDjx+3vW8YOjkPHw6FCpkTYyqzcCH07Gm9DgnRRT4scuTIQYUKFViwYAE5c+a0e3zCHJKMhXBF\nN27AjBl6YVb8wxvia9IEhg6FkiXtG1sqde8efPopjB5tbatYUQ9PGwYcPHiQF154gVKlShEaGioF\nPVIZGaYWwpVERupqWH5+uvv1cCL28tJdsZ07dU9ZErFdnD2rz86In4ife07/FXh6wuLFiylVqhQz\nZswAkEScCknPWAhXcPs2fPWVHm6+ccP2no8PtGype8I1ashxhna2eTO0aQMXLljb6teH776D7Nnh\nxx9/pG3btlSoUIHWrVubF6gwlSRjIZyZUrB4sa4SceKE7b3ChXV7hw56gZawq7g4fabGoEH6rwn0\n3uFhw/QRz2nS6EIeHTt2pGrVqqxatQovLy9TYxbmkWQshDO6d0+PcX71lW0NRYBixfTy3JYtwV3+\niZvh4kVdNTT+1qVcueDHH/XgBOhjEDt16kSNGjVYsWIFGTNmNCVW4RjkX6oQzuTMGX2kzzffwLlz\ntveyZ9flKt95R5KwiXbu1JVDz5yxttWoAT/8ALlzW9sKFCjAmjVrqFq1Kull5CLVk3+xQjiDa9f0\nfPBXX8Hdu7b33N2hWzfdG86a1Zz4BAChofrAh1u3rG2ffAKDB1vfH02cOJFChQrRsGFDatWqZU6g\nwuHIamohHFlsrO4J+/vrgsXxE3GePHrl9IkTMG6cJGKTrV8PDRpYE3G2bLB6tX4PZUnEo0ePpkeP\nHvzwww/mBSockvSMhXAkMTG6WPHu3fpjy5aExxiWLw+9e0PTptbz9YSpVq7UU/SW90q+vnoVddGi\n1meGDRvGgAEDCAkJYc6cOeYEKhyWJGMhHMGNG7pG4pdf6sNtE5M/P4wapfcJyz5Uh7FkCQQH6/dR\nAAUK6OHqwoX1tVKKgQMHMnToUNq1a8esWbNwc3MzL2DhkCQZC2Gm6GiYOhU+/1wvwU1MhgzQrx/0\n6iUnKDmYvXv1HmJLIn7uOd0jLlDA9rlLly7x5ptvMm3aNEnEIlGSjIUww40bei74yy9tl92CXnJb\nqRKUK6c/ypeHzJnNiVP8p6tXoVUruHNHXxcponvElgOWlFJcuHCBXLlyMWHCBADSpJFlOiJxkoyF\nsKfISJg4UX9cvmx7z89Pr4ju0EG2Jjk4peCtt+Cff/R1pkywYoU1EcfFxdGtWzdWrFjBH3/8gY+P\nj3nBCqcg/+KFsIebN/Vq6DFjbPe9gO4Jf/opvP22lKp0EhMnwqJF1usZM6yLtWJjY+ncuTPffvst\nffr0IUeOHOYEKZyKJGMhUlJcHMybp+d8z561vSflKp3Szp16MbtFly56ARdATEwMnTp1Yu7cuQwY\nMIDBgwfLoQ8iSSQZC5FSdu+Grl31f+MrVUpXgmjRQoajncxff0GzZroaKUDZsjB2rPX+yJEjmTt3\nLkOHDqV///7mBCmckvwmECK5Xb0K/fvD5MnWEwJAD0d//jm0b69PCRBOJSxMH3x17Zq+9vaGhQtt\nZxZ69OhBgQIFaNeunSkxCuclvxGESC5K6YnEYsVg0iRrIk6XTveEjxyBjh0lETuhRYugTh1rIs6c\nGZYvh0KF4M6dO3z22WfcvHmTzJkzSyIWT0V6xkI8q5s3Yf58mDJFbzyNr25d3UO2VIAQTmfqVD0v\nbHlvlTs3rFsHpUtDdHQ0zZo1Y/369ZQtW5YmTZqYG6xwWpKMhXhaZ8/qilhz5uih6fhy54bx4/VG\nVFnA47RWrrRNxEWK6BrUBQvCzZs3ady4MVu2bGHGjBmSiMUzkWQsxNNYuRLeeAOiomzb06fXG1CH\nDIEsWcyJTSSLw4f1mcSWRFyunD74wccHrl+/ToMGDfj111+ZM2eODE2LZybJWIgncfu23o40aZJt\nu7+/7kJ16CCnJ7mAa9f0ORyWOWI/P2siBjh//jwRERHMnz+fYMu+JiGegSRjIZLqjz/0Aqy//rK2\n5cunJxXr1ZOFWS4iLg7atYNDh/R1+vSwbJlOxDdv3iRDhgw8//zzHDlyhAxSK1wkE/ntIcTjHDqk\nqzq8/LJtIm7WDP78Ux9iK4nYZQwdqktbWsyYAWXKQGRkJFWrVn2wf1gSsUhO8htEiP9y7JjuCRcv\nDj/9ZG339ISvv4bFi/UJ8sJlzJ0LgwZZr3v10qcynT9/nldffZVDhw5RrVo10+ITrkuGqYV42B9/\n6FXSCxfqMcv4mjaFESP0XmLhUtatg06drNeBgTByJJw9e5bAwEBOnjzJqlWrCAwMNC9I4bIkGQth\nsX27XgW9fn3Ce7Vrw7BhekmtcDm7dunqpJZziUuW1IU+lLpHUFAQp0+fZt26dbzyyivmBipcliRj\nIU6e1Ac5/Phjwnu1a+sTleSXsMs6ckRP+1sO0/Lz071kvTPNgyFDhpA3b14qV65sZpjCxUkyFqnX\njRt6OHrMGL1lySJNGl2so08fvWhLuKw9e/Q6vMhIfZ09ux4YiY4+xooV+2ncuDEtW7Y0N0iRKkgy\nFqnP4cN6Adbs2QkrZ7VuDcOHw/PPmxKasA+l9I60Dz6Au3d1W/r0sGoVGMZhXnmlJnFxcQQGBpIx\nY0ZzgxWpgiRjkTpcuaKrNsyaBZs3J7xftix89RVUrWr/2IRd3bgBnTvrcuIWmTPrBfOZMx+gevWa\nKKXYtGmTJGJhN5KMhes6exaWLtUVG8LCrKtz4nv+eT0nLMcapgoHDkDLlnDwoLWtdGm9WOvWrb+o\nUSMINzc3tmzZQjFZMS/sSJKxcC03bsCSJXrD6ObNtucJW6RJA40aQdeuev+KJOFU4fvvdY/YslAL\ndBnxCRP0EPWgQUtImzYtoaGhFClSxLxARapkqMR+WdlBQECACg8PN+W1hQs6dw4GDoR582x/28ZX\nrpw+Hb5dO71kVqQKt2/D++/DN99Y29Kn18sGOnSAmJgY3N3dUUoRGRmJj6UAtRDJzDCMPUqpgMTu\nSc9YOLfoaBg3ThfiuHnT9p5hQI0aelyycWNdR1qkKpcu6d1pe/ZY24oW1fVcSpaEX3/9lTfeeIOV\nK1dStGhRScTCNJKMhXOKjtYTfQMG6H3C8RUvrnu/bdpA/vzmxCdMd/cuNG9um4iDg2H6dPDygrCw\nMBo2bIivry+ZMmUyL1AhkGQsnElcHGzdqueDFy2C69dt75coAWPHQlCQ7hWLVEspePdd/eNiMWEC\ndOumfzQ2bdpE48aNKVSoEJs2bSJPnjzmBSsEkoyFs1i3Dnr0gKNHE97LkUMftfPWW+AuP9JC13KZ\nNct6/fnn0L27/nz79u00bNiQIkWKsGnTJnLmzGlOkELEI8tIhWM7dw5CQvR5wQ8nYn9/nYSPHtXd\nIEnEAn2Y1scfW6/feAP69rVely5dmk6dOrFlyxZJxMJhyG8v4Zhu3YJvv4X+/eHaNWu7tze0bavn\nhMuXl+FoYeP33/WPhkX16rrSlmHAxo0bqVixIl5eXkyZMsW8IIVIhPSMhWM5ehQ+/BB8ffWwdPxE\n3K6druo/aRJUqCCJWNi4cEGfcBkdra/9/fWW87Rp4YcffqBevXoMGDDA3CCF+A/SMxaOYe9e+OQT\nWLs24T1/f70pVM6RFf/h3j1dVvzff/W1t7euM50tG8yZM4dOnTpRrVo1hg0bZm6gQvwH6RkLc50/\nD2+/rU9HejgRFy6sT1T66y9JxOKReve2rpw2DF13ukgRmDFjBm+88QY1a9ZkzZo1soVJOCzpGQtz\nXLwIM2boZa7xtygZhj5ctksXqFNHSlWKx5o9W29bshg6FOrXhxs3bjB48GDq1q3LkiVL8PT0NC1G\nIR5HkrGwn9u3YcUKvU943bqEBzfUr697wlKgXyRBbCxMm6aXGFg0b65nO5RSZMqUiZ9//pk8efKQ\nLl068wIVIgkkGYuUd/48jB+v532vXEl4v1gxXayjbl37xyac0p9/6kMfdu60tr34ou4ljxo1krNn\nzzJ+/HgKFixoVohCPBEZAxQp59gxeO89KFBAD0c/nIgrV9ZD1X/+KYlYJMmlS/DRR/r46fiJuEgR\nWL5cMW7cED7++GMiIyOJjY01L1AhnpD0jEXyu3gRPvtMH5MTF2d7r1AhfVRO27bw3HPmxCeczvnz\nevBkyhR9SqZF2rS6wEffvophw/ozYsQIOnTowLfffoubm5t5AQvxhCQZi+Rz967eAzxkCFy9anuv\nXDno108fYSi/JEUSXb+u39dNm2bdP2xRo4Yu6FG0KHzyyad8/vnnvP3220ydOpU0svBPOBlJxuLZ\n3bih95KMGZOwZGXNmvDpp/Dqq1KkQzyREyegUSPYt8+2vXhxfVhX69bWH6lKlSrx/vvvM3bsWEnE\nwikZSilTXjggIECFh4eb8toimRw6pBdlzZ5tWykL9CTel1/qbUqShMUT+uUXaNYMIiOtbWXL6uqo\njRvrHW9xcXHs3r2bChUqmBeoEE/AMIw9SqmAxO7JW0jxZKKidAKuXFmvgp4wwTYRZ8kC48bB339D\nw4aSiMUTmzVLD6hYEnHatDBzJuzerctdpkkDsbGxvPXWW1SuXJm//vrL3ICFSAYyTC0eTykIDdXz\nwatX69qDD/P314U6OnbUCVmIJ6SUXm4waJC1LWdOWLpUv/eziImJoWPHjsybN4+BAwdSsmRJu8cq\nRHKTZCz+W2ys/k04ahQkNqXg7q57v1266HKVMlcnnpJSelX0qFHWtlKldI2YAgWsbffu3eP1119n\n4cKFDB8+nE8++cT+wQqRAiQZi4Tu3IHvvoPRoxMuyAJ9dGG7dhAcDD4+9o9PuBSloGdPXRfGok4d\nWLQIHi4lvWjRIhYuXMiYMWPo1auXfQMVIgVJMhZW167pPSTjxsG5c7b3PD31Ke3vv6/3kgiRDOLi\noGtXvUXJolEjWLgQEqtgGRISgp+fH1WqVLFfkELYgYwrCr0neOBA8PODPn1sE7G3ty72e+KErrgg\niVgkkytXoEUL20TcsqXuEcdPxNHR0bRp04Z9+/ZhGIYkYuGSpGecmt28CRMnwhdfwOXLtvfy5NEV\n+N95BzJnNic+4bJ+/x1atYKICGtbmzYwZ45eimBx8+ZNGjVqRFhYGA0aNKBEiRL2D1YIO0hSz9gw\njLqGYRw2DOMfwzD6/cczrQ3DOGAYxn7DMOYnb5giWd25o7ckPfecXjUTPxHrQ2Dh+HF9SKwkYpGM\nlNJVUitXtk3EH3yglynET8TXr1+nXr16bN26le+++47XX3/d/gELYSeP7RkbhuEGTAZqAaeB3YZh\nrFBKHYj3jD/wMVBFKXXZMIycKRWweAb37umux5Ah8O+/tvcKF4bBg+G116RcpUgxQ4fqGRELLy+9\nh7hlS9vnrl69Sr169di1axfz588nODjYvoEKYWdJGaYuD/yjlIoAMAzjR6AJcCDeM28Dk5VSlwGU\nUheSO1DxDJSCxYt1L/iff2zv+frq4r9vvAEeHubEJ1KF8HD9PtCiVCm9UKtIkYTPpk2blqxZs7Jw\n4UKaNWtmvyCFMElSkrEvEL8bdRp4uP5cEQDDMH4F3IBBSql1yRKheDbHj+vlqmvX2rb7+OiFWe++\nq1dKC5GC7t6FTp301nWAatVg/XpIn972ucjISNzc3MiaNSurVq3CkApuIpVIrgVc7oA/UAPIB2wz\nDKOkUsrmAFvDMN4B3gHw8/NLppcWibp3T9eGHjLE9ribLFn0iunu3RNu4hQihXz+ua6QCjoBz5qV\nMBGfP3+ewMBAcuTIwZYtWyQRi1QlKcn4DJA/3nW++23xnQZ2KqXuAccNwziCTs674z+klPoG+Ab0\nQRFPG7R4jIgIPQn3xx/WNsPQlbKGDZNylcKu/v4bhg+3Xo8YkfAo67NnzxIYGMipU6eYMGGCJGKR\n6iRlNfVuwN8wjEKGYaQFQoAVDz2zDN0rxjCMHOhh6wiE/a1erY+3iZ+IS5eGHTt0bWlJxMKOYmL0\n8LSlnHmlSnpQJr5///2X6tWrc/r0adatW0fNmjXtH6gQJntsMlZKxQDdgPXAQeAnpdR+wzCGGIbR\n+P5j64EowzAOAFuAj5RSUSkVtEhEbKxeiNWwoa6mAPq4m9Gj9cqZ8uXNjU+kSiNHWsuap00L336b\ncLF+hw4duHDhAhs2bKBatWr2D1IIByDnGTu7K1d0yaLp02HXLmt7/vx6BXW5cubFJlK1efOgbVvr\n9YgRekH/wyIiIrh06RIBAYke8yqEy5DzjF1RWJguYZQ7N7z9tm0iDgrSJY4kEQuTbNigT9O0qFZN\n15CxOHToEP369SMuLo7ChQtLIhapniRjZ3P7tj6s4dVXdY/4zh3rPTc36N8f1q2DHDnMi1GkauHh\n0Ly5ni8GKFECli+3bmPft28fNWrUYNasWZw58/BaUCFSJ6lN7UwOHYKQEPjzT9v2MmX0kYavvaZ7\nykKY5J9/oH59XfYc9GzJ2rWQNau+/vPPPwkKCsLDw4PQ0FDy58//319MiFREkrEzuH1br3zp0wdu\n3bK2N2qkN3AWL25ebELc9++/eobk4kV9nTWrLuyRL5++3rNnD7Vq1SJjxoyEhobi7+9vXrBCOBhJ\nxo4sIkKfLzdzJkTFW5yeLp0u6NGli94/LITJzp2DmjXh5El9nT49rFoFxYpZn4mKiiJXrlysWbOG\nQoUKmROoEA5KkrEjioyEzp1h6VJdVzq+YsXgxx91YV8hHMDFi7pHbCl77uEBS5bok5lAl7jMkSMH\ntWvX5u+//8bdXX7tCPEwWcDlaPbv13uClyyxTcQFClg3bUoiFg7i8mWoXRsO3D82xs0NfvoJ6tbV\n12FhYRQuXJjFixcDSCIW4j/IvwxHsnq1XoR1/bq1rW5dfdBDvXpytKFwKDdv6sVae/fqa8OA77+H\npk319caNG2nSpAmFChWiSpUq5gUqhBOQZOwI4uL0HHDfvtbecMaMMH8+NG786P9XCBPcuwfBwbrK\nqsXMmXqxP8CaNWto3rw5L7zwAhs3bsTHx8ecQIVwEpKMzfbzz9CzJ+zZY20rUABWrJDhaOGQlNJL\nGlavtrZNmGAt8nH06FGaNm1KqVKl2LBhA9myZTMlTiGcicwZm+XECWjdGl55xTYRV6miq2lJIhYO\n6tNP9RGIFh9/bHv4g7+/P19//TWbNm2SRCxEEkkyNsOGDfDii7BwobXN01NXz9q8GXLmNC82IR5h\n/Hi9td3ijTesxyP++OOP/P777wC8+eabZJETwoRIMknG9rZzJzRrBtHR1raQEDh8GIYO1XuIhXBA\n48bBBx9Yrxs0gG++0Qu3Zs+eTZs2bRgxYoR5AQrhxGTO2J4OHNDLTy1VtPz84IcfrBsyhXBQI0fa\nnrhUqZLewuTuDt988w2dO3emVq1afPfdd+YFKYQTk56xvZw6BXXqwKVL+jpHDj1cLYlYODClYNAg\n20Rctao+iyRDBpg0aRKdO3emQYMGrFixggwZMpgWqxDOTJKxPURE6MoIp0/r60yZdPX8okXNjUuI\nR4iJgY8+gsGDrW01a+pEnDkzxMXFsX79epo2bcqSJUvw9PQ0L1ghnJwMU6eku3f1/uEhQ/RhDwBp\n08KyZSDntwoH9r//6fozYWHWtrp1dWG49OkhOjqa9OnTs3DhQtzc3PCwnI8ohHgq0jNOKb/8Ai+/\nDJ98Yk3Ebm4wbx4EBpobmxCPEBamT+WMn4ibNNHvIT09FYMGDaJy5cpcvXoVT09PScRCJANJxsnt\n5k296bJaNV1n2uKll2D7dmjZ0rzYhHiMceP0e8X//U9fG4Yepl68GNKmVXz66acMHjyYl156iUyZ\nMpkbrBAuRIapk9Mvv+gyRMeOWdsyZtRblrp310tPhXBQS5bAhx9ar318dEXWoCBQStG7d2/Gjh1L\n586dmTJlCmnSyHt5IZKL/GtKDrdvQ69euppW/ETcsCEcPKjLXUoiFg7s2DFdwMOiUiX44w+diAFG\njBjB2LFj6d69O19//bUkYiGSmWSIZ/W//+ljanbutLZ5e+tSRe3b63E+IRzY7dvQqhVcu6avCxbU\ndaezZrU+0759e9zc3Ojbty+G/EwLkezk7e2z+P13KFfONhHXqQP79kGHDpKIhVP44APdCwa92H/h\nQp2IY2NjmTFjBrGxseTPn59+/fpJIhYihUgyflqLF+vqB5a9w2nS6NUva9dCvnzmxiZEEs2bB9Om\nWa/HjdO77mJiYmjfvj1vv/02q1atMi9AIVIJGaZ+UrGxennp0KHWNm9vXRuwdm3z4hLiCf39tz4K\n0SI4GN57D+7du8frr7/OwoULGTFiBE2aNDEvSCFSCUnGTyIqCtq00WUsLZ5/HlauhBdeMC8uIZ7Q\npUt6qcPNm/q6SBGYPh3u3r1DSEgIy5Yt48svv+TD+MurhRApRoapkyo8XBfxiJ+IAwP1fLEkYuFE\nYmN1da2ICH2dMaPe1uTlBQcOHGDDhg1MnDhRErEQdiQ946SYMQO6dtXlLS0+/lgPVbu5mReXEE/h\n009t31N+9x288EIs4EaZMmX4559/yJMnj2nxCZEaSc/4UaKj4c034e23rYk4c2ZdF3DECEnEwuks\nWACjRlmvP/0Uate+Qa1atZh2fyWXJGIh7E+S8X85fhyqVIGZM61tJUvCnj26UK8QTuavv6BTJ+t1\n/frQq9c16taty9atW/Hy8jIvOCFSOUnGiVm+HMqWtW6+BGjbFnbs0Au2hHAyUVF6wdatW/ra3x+m\nTLlCvXq12blzJz/++CNt2rQxN0ghUjFJxvFdvqyrZjVtqj8H8PCAyZP1xJocnC6cUEyMXrB1/Li+\nzpQJfvrpDi1aBPH777+zaNEiWrVqZW6QQqRykowt1q6FEiVg7lxrm68vbNsGXbpINS3htD7+GDZu\ntF7PnQsvvZSO1157jWXLlsk+YiEcgKymBvjsM9siHgCvvw4TJkC2bObEJEQy+OEHGDPGev3hh//D\n1/dfoBy9evUyLS4hhC3pGX//vW0i9vHRmy6//14SsXBqq1fbLtgKCjrDqlXVadasGbdv3zYvMCFE\nAqm7Z/z773rbkkWtWrpYr4+PeTEJkQxmzdI/2rGx+rpw4VMcO1aTyMgLrF27Fk9PT3MDFELYSL09\n44sXoVkzfX4cQLFisGiRJGLh1JSCzz/XPWJLIvb1Pc6dO69w6VIkGzdupEqVKuYGKYRIIHUm43v3\noHVrOHVKX3t760IemTObG5cQzyAmBnr0gE8+sbaVLg1BQWOIjr7O5s2bqVChgnkBCiH+U+obplYK\nevaEsDB9bRh6aLpIEVPDEuJZnDsHISF68b/Fq6/C0qWQPv04Tp78AH9/f/MCFEI8UurrGY8YofcN\nWwwbBg0amBePEM8oNBReesk2Edepsw+lArl79yJp06aVRCyEg0tdyXjaNOjf33rdurXehCmEE4qN\nheHD9brDCxd0W5o00KXLXsLDa3DkyCEuW4rXCCEcWupJxosW6ZPTLYKCdFUtKeYhnNDx43oYun9/\niIvTbTlzwqRJ4fzwQ00yZMjA1q1bKSLTL0I4hdSRjDdv1kU8lNLXAQF6L3G6dObGJcQTUkqfXVKq\nFPz8s7W9WjWYPTucfv0C8fb2Ztu2bTwvddSFcBqun4z379dbmCxHIBYtCmvW6JPUhXAiV69C8+b6\nVM8bN3SbmxsMHKjnjcuUyUe1atXYtm0bBQsWNDVWIcSTce3V1BcvQqNGcP26vvb11aeqy15i4WTO\nnYN69eDPP61tRYroOtNp0+4FSpA7d25WrVplWoxCiKfnuj3jO3d0N8JyVE3GjLo+oJ+fuXEJ8YQO\nH4ZKlWwTcdeu+oTPy5fXU6lSJQYOHGhegEKIZ+aayVgpePdd+OUXfW0YMH++roAghBPZuROqVIGT\nJ/W1m5sudTlpEoSGrqJx48a88MIL9OzZ09xAhRDPxDWT8dixMHu29XrkSGjc2LRwhHhSt27pLfE1\na0JUlG7LkAFWrICOHWHp0qU0b96cUqVKsXnzZnLkyGFqvEKIZ2MoywpjOwsICFDh4eHJ/4Xv3NFz\nwpZ54o4d9fJT2cIknEBcnD4w7NNP4fRpa3v27HqWpUIFuHr1KoUKFaJo0aKsW7cOb29v8wIWQiSZ\nYRh7lFIBid1zvQVcly5ZE3GWLDB1qiRi4RT++QeCg/VhYvEVLw6LF+uNAADe3t5s2LCBokWL4iW7\nAoRwCa43TH3njvXzzJllL7FwCgcPwiuv2CbiXLl00bi9e3UinjlzJuPHjwcgICBAErEQLsS1k7Ek\nYuEE/voLqlfX25cAPD1hwAA4ehTeeQfc3WHq1Km8+eabrF27lljL2YhCCJfhesPUkoyFE9mzB2rX\n1rMrAJky6bnhV16xPjNhwgTef/99GjRowKJFi3BzczMnWCFEipGesRAmCQ2FwEBrIs6cWdekiZ+I\nx4wZw/vvv0+zZs1YsmQJnp6e5gQrhEhRkoyFsLO4OH1yZ61ausQlQNasuoR6pUq2z6ZLl47g4GAW\nLFhA2rRp7R+sEMIuJBkLYUeRkfr47AEDrKct5cqle8kB9zc8KKU4ceIEAN27d+eHH37Aw8PDnICF\nEHYhyVgIO/n7b3j5ZVi3ztpWrZpeQf3SS/paKcXHH39MiRIlOHz4MACGbM0TwuVJMhbCDs6c0Qc9\n/Puvta1vX90jzptXXyul6NWrF6NGjaJdu3b4+/ubE6wQwu5kNbUQKezmTV2N9cwZfZ05s66y1aiR\n9Zm4uDh69OjB5MmT6dGjB1999ZX0iIVIRVy7ZywLXoTJ4uKgbVtrMQ83N11NK34iBpg9ezaTJ0+m\nd+/ekoiFSIWkZyxECurXD5Yts15PmQJBQQmfa9++PRkyZCA4OFgSsRCpkGv3jCUZCxNNmQKjR1uv\ne/XSFbUsYmJi6NOndvq6CAAAIABJREFUD+fOncPd3Z2QkBBJxEKkUq6XjO/etX4uyViYZPx46NrV\net24MYwaZb2+e/cuISEhjB49mtWrV9s/QCGEQ3G9ZCw9Y2GyL76ADz6wXpcvD/Pm6fligDt37tCq\nVSsWL17M2LFjeeutt8wJVAjhMGTOWIhkohQMHQoDB1rbqlSBNWt0zWmA6OhoWrRowdq1a5k8eTJd\nunQxJ1ghhEORZCxEMnk4EdeoAStXWhMxwK1btzh9+jTTp0+XHrEQ4gFJxkIkg6lTbRNx7dqwdClk\nyKCvb9y4Qdq0acmePTvh4eFSZ1oIYUPmjIV4RkuWQPzR5jp1YPlyayK+evUqderUoW3btiilJBEL\nIRKQZCzEM9i6Fdq00fPFAOXKwaJFYDnp8PLly9SuXZtdu3bRunVr2bokhEiUDFML8ZT+/huaNLH+\nyPn7w+rV1jniqKgoatWqxb59+1i0aBFNmjQxL1ghhENLUs/YMIy6hmEcNgzjH8Mw+j3iuRaGYSjD\nMAKSL8QnJMlY2MG5c1C/vvU84ty5Yf168PHR10opmjdvzoEDB1i+fLkkYiHEIz22Z2wYhhswGagF\nnAZ2G4axQil14KHnvID3gZ0pEWiSSTIWKezWLV3E4/Rpfe3lBWvXQqFC1mcMw+CLL77g+vXrBCVW\n/1IIIeJJSs+4PPCPUipCKXUX+BFI7G3+UGAUcDsZ43tykoxFCoqLg3btIDxcX7u56Tliy3nEZ86c\n4ZtvvgGgQoUKkoiFEEmSlGTsC8Q7hZXT99seMAzjZSC/Usr8un6SjEUK+uQTvXraYuJEvY0J4NSp\nU1SvXp3evXtz7tw5cwIUQjilZ15NbRhGGmAs0CsJz75jGEa4YRjhFy9efNaXTpwkY5FCZs60rS/9\nwQfw3nv684iICF555RUiIyPZuHEjefLkMSdIIYRTSkoyPgPkj3ed736bhRdQAggzDOMEUBFYkdgi\nLqXUN0qpAKVUgI9lpUtyk/OMRQpYs8b2xKWGDWHMGP350aNHqV69OtevXyc0NJQKFSqYE6QQwmkl\nZWvTbsDfMIxC6CQcArSx3FRKXQVyWK4NwwgDeiulwpM31CSSnrFIZjt3QqtWEBurr0uXhvnzrQc/\n7Nixg7t377JlyxZKlSplXqBCCKf12J6xUioG6AasBw4CPyml9huGMcQwjMYpHeATk2QsktGhQ9Cg\ngV5BDVCwoF457eWlT18CaNeuHUeOHJFELIR4akmaM1ZKrVFKFVFKPaeUGn6/7TOl1IpEnq1hWq8Y\nJBmLZHPmjC5tGRWlr3Pk0HuJ8+SBP/74A39/f7Zu3QqAt7e3iZEKIZyd65XDvHvX+rkkY/GUTp6E\nwEA4dUpfZ8igq2sVKQK7du2iZs2aGIZBvnz5zA1UCOESXC8ZS89YPKO//4bKleHwYX3t7g6LF0P5\n8rB9+3aCgoLImjUr27Zt47nnnjM3WCGES5BkLEQ827ZBtWpw9qy+TptWL9aqWxcOHjxI7dq1yZ07\nN9u2baNAgQLmBiuEcBmulYyVkmQsntrSpbqAh6XetKXMZatW+rpo0aL07NmTrVu3yvC0ECJZuVYy\njomxnmXn5mbdeyLEYyxfrpOu5b1c7ty6l1yzJmzatIlTp06RJk0ahg4dKgU9hBDJzrWSsfSKxVNY\ntw5at7buI/b3h+3bdb3plStX0qBBAz788ENzgxRCuDRJxiJVCwuDZs2si/Cffx62btUnMC1evJjm\nzZtTunRppk+fbmqcQgjXJslYpFrbt+uylrfvnzNWoABs3qz3Ef+/vTuPr6q4/z/+GpIgjShQpFSj\ngCKKiFb5IlpFCGGLrGUVEQShBQSqgKJWK+6CVbC2BVwgZRWQRYhRQHYtLaA1lB+bCIgsKiBQZEnM\nNr8/JulNIMtF7n7fz8eDR8+cnHPz6THkzcyZM2fWrFncfffdNGrUiGXLllGlSpXgFisiEU1hLFFp\n7Vo3Q/rUKde+7DIXxDVqQG5uLq+99hp33HEHS5Ys0YIeIuJ33qxNHT4UxuKF1atdj7ggiH/xCxfE\ntWtDXl4eMTExLF68mAoVKnDhhRcGtVYRiQ7qGUtU+egjuOsuTxBXr+6CuG5dmDBhAh06dODHH3+k\natWqCmIRCRiFsUSN1FRo395zj/iyy9xkrfr14c9//jNDhgyhXLnI+ishIuEhsn7zKIylGEeOQL9+\n0LGjZ9Z0jRruOeJrr4U//elPDB8+nC5dujBv3jwu0M+OiARY5IZx+fLBq0NCgrUwdaobgv773z37\nr7rKBXHt2jB27Fgee+wxevTowezZsymvnxsRCYLIDWP1bqLav/8NiYnQty98/71nf5cubiZ1wbLS\nzZo1Y8iQIcyYMYPY2Miazygi4UNhLBFl717o3RsaNnS93wI1asD778O8eVC9umXZsmUANGjQgL/9\n7W/EaOlUEQkihbFEhNxcGDXKvW94xgzP/thYePhh2LrVPc5krWX48OG0atWKJUuWBK9gEZFCImtc\nrmB2DiiMo0hmJtx7LyxYUHT/b34DY8a4SVrgniEeOnQoEydOZNiwYbRu3TrwxYqIFCOywlg946jz\n3/+6WdKFh6RvuQXGjnXvJS6Ql5fHgAEDmDx5Mo8++ihjxozBGBP4gkVEiqFhaglbBw64wC0cxMOG\nwbp1RYMY4J///CcpKSk89dRTCmIRCTnqGUtY+vRTNzN63z7PvldecfeHi8vZxo0b8/nnn3PTTTcF\nrkgRES+pZyxhxVoYPx7uuMMTxLGxMG0aPPJI0SDOysqiZ8+e/5uopSAWkVClMJawceIE9OwJQ4dC\ndrbbV6mSe2Spd++ix/7444907dqVWbNm8eWXXwa+WBGRc6BhagkL+/ZBq1awfbtn3803w9y5biWt\nwjIyMujcuTNLlixhwoQJPPDAA4EtVkTkHKlnLCHvu++gefOiQTxgAPzzn2cHcWZmJu3bt2fp0qVM\nmjRJQSwiYUE9Ywlp338PLVpAwUhzXBxMnnz2sHSB8uXLc80119CnTx96l3SQiEiIURhLyPrvf93Q\n9JYtrh0TA3PmQKdOZx97/Phxjh07Rq1atZgwYUJgCxUROU8KYwlJJ07AXXdBerprGwPTpxcfxMeO\nHaN169YcO3aMLVu26M1LIhJ2FMYSck6fdutIr1vn2Td5Mtxzz9nHfv/997Rq1YotW7Ywb948BbGI\nhKXIDWP9Ug5LmZluTenCq2qNHw/333/2sYcOHaJFixZ8+eWXLFq0iOTk5MAVKiLiQ5EbxuoZh52s\nLOjWDfLfbgjAq6/C4MHFHz9y5Eh27txJWloazZs3D0yRIiJ+oEebJCTk5ECvXpCW5tn3/PNuecuS\nvP7666xatUpBLCJhT2EsIWHkSLeAR4E//AGefPLs4/bs2UO/fv3IyMigcuXK3HrrrYErUkTETxTG\nEnTTp8Of/+xpDxsGL7549gsfdu3aRdOmTXnvvffYvXt3YIsUEfGjyArjrCzPtsI4LHz+uVtNq0Cn\nTu5dxGcG8RdffEGTJk04deoUK1eu5Prrrw9soSIifqQJXBI0hw65mdOZma5drx5MnQrlzvgn4tat\nW0lKSsJay6pVq7jhhhsCX6yIiB9FVs9YYRw2fvgB7r7b8xrESpVg4UK46KKzj83Ly6N69eqsXr1a\nQSwiEUk9YwmozZthwgR3n/jkSbfPGHjnHahTp+ix+/fvJyEhgfr165Oenk65M7vMIiIRIrJ+uymM\nQ9aXX0JSEtxwA0yc6AlicI8wtWlT9Pj169dTv359xo0bB6AgFpGIpp6x+N2GDdC2rXsDU2H16sEj\nj0DfvkX3r127lrvuuotq1arRtWvXgNUpIhIskRPG1iqMQ9DixdC1q1tvGiA21s2YHjIEmjQ5e9b0\n6tWradeuHQkJCaxcuZKEhITAFy0iEmCRE8bZ2Z7t2Nizp+RKwE2bBv37u9W1AKpWhQ8+gJLW6Thy\n5AgdOnSgZs2arFixgl/+8peBK1ZEJIgiJ7HUKw4Z1sILL0CfPp4grlkT1q4tOYgBqlatysyZM1m9\nerWCWESiSuT0jBXGISEjw/WGZ83y7LvxRjdcfdllxZ+zaNEiADp27Ej79u0DUKWISGhRz1h85ptv\noGnTokGclOReh1hSEM+bN4+uXbsyduxYrLWBKVREJMQojMUnNm2CW26BTz/17HvgAViyxC3oUZxZ\ns2bRo0cPbr31VtLS0jBnzuYSEYkSkRnG5csHr44otHs3tGrlesYAMTEwfrxb3CMurvhzpk6dSq9e\nvWjcuDFLlizh4osvDlzBIiIhRveM5bwcOgStW8PBg65dqRLMmwctWpR+3saNG0lKSmLRokXEx8f7\nv1ARkRCmMJaf7MQJt3LWzp2uXaECpKVB48Yln3P8+HEqVarEuHHjyMrK4gL9txIRidBhav2C97us\nLOjSBf79b9cuVw5mzy49iF977TXq1avH3r17McYoiEVE8imM5Zzk5MDMmXDTTbBsmWf/G29Ax44l\nn/fyyy8zYsQIbr/9di699FL/FyoiEkYiZ5g6K8uzrTD2uexst6LW6NGwa1fRrz37LPzudyWf+/zz\nzzNq1Cjuuecepk2bRmxs5PzYiYj4QuT8VlTP2G/27oXu3WH9+qL7L7oInnwSHn205HMnT57MqFGj\n6NOnD5MnTyYmJsa/xYqIhCGFsZRq8WLo1QuOHvXsq1IFhg2D3//ebZeme/fuHD16lIcfflivQRQR\nKUHk/HZUGPtUTo7r9bZp4wni2Fg3JP311zBqVMlBbK3lr3/9KydPnuSiiy5i5MiRCmIRkVKoZyxn\nOXYM7r676ASthASYMwfuuKP0c/Py8hgyZAhvvPEGcXFxDBo0yL/FiohEAIWxFLFjB7Rv7/63QMuW\nbgZ1tWqln5ubm8uAAQNISUnh8ccfZ+DAgf4tVkQkQkTO2KHC+LwtX+5ecVg4iEeNcveNywrinJwc\n+vbtS0pKCk8//TQvvfSS1poWEfGSesYCuLWkH3oIcnNd+2c/g6lToVs3787/9ttvWb58OS+++CJP\nPPGE/woVEYlACuMol53tQnjiRM++hARYtAj+7/+8OT+b2NhYrrjiCrZs2cLPf/5z/xUrIhKhNEwd\nxY4eheTkokF8yy2wYYN3QZyZmUnnzp0ZOXIkgIJYROQnUhhHqe3b3f3hlSs9+3r0gDVr4LLLyj4/\nIyODjh07kpaWRp06dfxXqIhIFIjMMNb7jEu1ZIkL4oK3LQE8/zy88467V1yWU6dO0bZtW5YtW0ZK\nSopmTYuInCfdM44i1sJf/gIjRkBentsXH+/WnO7SxdvPsHTs2JE1a9Ywbdo0evXq5b+CRUSihMI4\nSmRlwdCh8Pbbnn2XXw6pqXDzzd5/jjGGwYMHM2DAALp37+77QkVEopDCOAocPw6dOsGqVZ59t94K\nCxfCL3/p3WccPXqUDRs2kJycTOfOnf1TqIhIlFIYR7hvvoG77oJNmzz77r0XJk2CChW8+4zvv/+e\nli1b8uWXX/LVV19RrawVQERE5JxE5gQuhTEAX3wBt99eNIhfeAGmT/c+iA8ePEhiYiLbt29nwYIF\nCmIRET+InJ5xVpZnW2HMhg3ujUtHjrh2TAxMngx9+nj/Gd988w3Nmzdn7969fPDBByQlJfmnWBGR\nKBc5Yaye8f/s2OFe7vDDD64dHw/z5rnh6nMxa9Ys9u/fz5IlS7jzzjt9X6iIiABeDlMbY5KNMV8Y\nY3YaYx4v5usjjDFbjTGbjDErjDE1fV9qGRTGAGRkuPWkC4K4alW3sMe5BLG1FoARI0awadMmBbGI\niJ+VGcbGmBhgPHAXUA+4xxhT74zD0oGG1tobgXnAn3xdaJkUxgA8+KDnHvEFF8BHH7mZ097auXMn\njRo1Ytu2bRhjuPLKK/1TqIiI/I83PeNGwE5r7W5rbRYwG+hY+ABr7Spr7en85jrgct+W6QWFMdOm\nuVnSBV5/HRo08P787du306RJE/bs2cOPha+niIj4lTdhnADsK9Ten7+vJP2BxedT1E8S5WG8eTMM\nGuRp33svDBhwLudvJjExkdzcXFatWsVNN93k+yJFRKRYPp3AZYzpBTQEmpbw9QHAAIAaNWr48ltH\ndRh//bVbzjIjw7Wvuw7eeAOM8e787du306xZM+Li4li5ciV169b1X7EiInIWb3rGB4ArCrUvz99X\nhDGmBfAk0MFaW+wYp7X2LWttQ2ttQ58/rxqlYbx2rXvt4Y4drl0wc7piRe8/o0aNGrRu3Zo1a9Yo\niEVEgsCbMP4UqGOMudIYUx7oAaQWPsAYczPwJi6ID/m+TC9EYRhPnQpJSXD4sGvHxbn7xvXOnF5X\ngvT0dI4fP058fDwzZszQqxBFRIKkzDC21uYAQ4GlwDbgXWvtFmPMc8aYDvmHvQJUBOYaYzYaY1JL\n+Dj/sLbooh8R/gpFa+Hxx6FvX8//7WrV3NrT3r596R//+AdNmjRhyJAhfqtTRES849U9Y2vth8CH\nZ+wbVWi7hY/rOjeFgzguDspFziqfxXnpJXj5ZU/7hhvc25dq1fLu/NWrV9O2bVuuuOIKXi78QSIi\nEhSRkVqFh6gjvFf87rvwxz962u3bu/vG3gbxsmXLaNOmDbVq1WL16tUkJJQ2MV5ERAIhMpbDjJL7\nxevXF11bulkzN1nL239/ZGdnM3jwYOrUqcPy5cv10gcRkRChMA4TX38NHTtCZqZrX3MNzJ9/bgMB\ncXFxLFmyhMqVK1O1alX/FCoiIucs8oapIzCMv/sO2rWDgwdd++c/hw8+gCpVvDt/7ty5PPjgg1hr\nqV27toJYRCTEKIxDXHo6NGrkVtgCNz/tvffg6qu9O3/mzJn06NGD9PR0MgpWBRERkZCiMA5hCxZA\n48awL38x0nLlICUFmjTx7vwpU6bQu3dvmjZtyuLFi4mPj/dfsSIi8pMpjEOQtfDii+6Z4dP5r9+4\n+GI3NN2rl3efMWnSJO6//35atGhBWloaFc9lSS4REQmoyAjjws8ZR0AYT5hQ9PGl2rVh3TpITvb+\nMy699FI6depEamqqesQiIiEuMsI4gnrGn34Kw4d72omJ7pGm667z7vxt27YB0LZtWxYsWECFChV8\nX6SIiPiUwjiEHDsG3btDdrZr33wzLF4M3k5+Hj16NPXr1+fjjz/2X5EiIuJzCuMQYa1b0GPPHteu\nVAnmzgVvOrbWWp599lmeeOIJ7rnnHm6//Xa/1ioiIr6lRT9CxKuvwvvve9p//7u7V1wWay1PPvkk\no0ePpm/fvkyaNImYmBj/FSoiIj6nnnEImD8f/vAHT3vECOjUybtzly1bxujRoxk4cCCTJ09WEIuI\nhCH1jINswgQYOtQNUwP8+tcwZoz357ds2ZL333+ftm3bYozxT5EiIuJX6hkHibXu8aUhQzxBfO21\n7j5xXFzp5+bl5TFy5Ej+85//YIyhXbt2CmIRkTCmnnEQ5OTAwIFuNa0Ct94KaWlwySWln5ubm8tv\nf/tbpkyZQpUqVfjVr37l32JFRMTv1DMOMGvhoYeKBnGbNrBiRdlBnJOTw3333ceUKVN45pln+EPh\nG80iIhK2Iq9nfC7vFAyCv/7V3Scu0KcPvP122UPT2dnZ3HvvvcydO5eXXnpJQSwiEkEiL4xDuGf8\nwQdFV9e6+27XQy7nxfhEXl4eJ06cYOzYsYwYMcJ/RYqISMApjANk0ybo0QPy8lz7ttvcs8RlBXFm\nZiYZGRlUqVKFtLQ0PbokIhKBdM84AL79Ftq1g5MnXbtmTVi4EH72s9LPO336NB06dCA5OZmcnBwF\nsYhIhFLP2M8OH4YWLTzvJL7oIjdrunr10s87efIk7du3Z82aNaSkpBAbGxn/qURE5GyR8Rs+RMP4\n2DFo1Qq2bnXtmBh4912oX7/083744QfatGnDv/71L2bMmEHPnj39X6yIiASNwthPfvjBvX9440bX\nLlcOZs707p3EAwcOZP369cyePZtu3br5t1AREQm6yLhnnJXl2Q6BMD51yt0j3rDBsy8lxc2e9saY\nMWNYuHChglhEJEpERhiHUM84Nxd69oRPPvHsmzjRPU9cmsOHD/Pcc8+Rl5dHzZo1adu2rX8LFRGR\nkKEw9rHHHoPUVE977FgYNKj0c7777jsSExMZPXo0WwtuMIuISNTQPWMfeustF74FHnnEvQ6xNAcO\nHCApKYn9+/fz4YcfUr+s2V0iIhJxFMY+snw5DB7saf/mN2W/CnHv3r0kJSVx6NAhli5dSuPGjf1b\npIiIhCSFsQ9s2wZdu7r7xQANGsCMGe5RptLs3r2bU6dO8dFHH3Hbbbf5v1AREQlJCuPzdPgwtG0L\nx4+79mWXuXvGF15Y8jknT56kYsWKJCYmsmvXLuLj4wNTrIiIhCRN4DoPmZnQqRN89ZVrx8e71bUS\nEko+Z9u2bdStW5fp06fnn6MgFhGJdgrjn8ha6N8f1q51bWNg1iy4+eaSz9m8eTOJiYnk5OTQoEGD\ngNQpIiKhL/LCOEDvM37+eXjnHU/71VehQ4eSj9+4cSOJiYnExsayZs0arr/+ev8XKSIiYSHywjgA\nPeN334Wnn/a0Bwwo+p7iMx08eJCkpCTi4+NZs2YN1157rd9rFBGR8BH+YZyXB9nZnrafe8a7d8Nv\nf+tpt2gBf/ubG6YuSfXq1Xnuuef4+OOPufrqq/1an4iIhJ/wn01deF3quDj3RgY/yc52S12eOOHa\ntWvD3Lnu2xbnk08+oUKFCtxyyy0MHTrUb3WJiEh4C/+ecQCHqJ99Ftavd9uxsW7CVuXKxR+7YsUK\nkpOTGTZsGNZav9YlIiLhTWHspdWr4aWXPO0XXoBbbin+2KVLl9KuXTuuuuoqFixYgCltDFtERKKe\nwtgLR45Ar17ucSaA5s1h5Mjij01LS6NDhw7UrVuXVatWUb16db/UJCIikUNhXAZrYeBAOHDAtatW\nhWnTSr41PXXqVG688UZWrFjBJZdc4vN6REQk8kTWBC4/hPHMmTB/vqedkuKWvDxTTk4OsbGxzJgx\ng8zMTCpVquTzWkREJDKFf8/4uutcIJ84AevW+fSj9+2DwpOgBw4sfmGPGTNm0KhRI44cOcIFF1yg\nIBYRkXMS/mFsjHu2qGJF8GEI5uVBv36eF0BcdZVbZetMKSkp3HfffVSuXJkKFSr47PuLiEj0CP8w\n9pOJE907isHl/bRpLu8Le/PNN+nfvz8tW7YkLS2NC0t7VZOIiEgJFMbF2LGj6GzpRx+FO+4oesy0\nadMYNGgQbdu2ZdGiRXr7koiI/GQK4zNs3gytW0NGhmvfcINb7ONMzZs358EHH2TBggUanhYRkfOi\nMC4kLQ1+/WvYs8e14+Jg+vSik7QXLlxIbm4uCQkJvP7665QP0FuiREQkcimMcc8Sv/KKmyl98qTb\nV7EivPce/OpXBcdYnnnmGTp16kRKSkrwihURkYgT/s8Z+8Ajj8C4cZ52rVqQmuqGqMEF8RNPPMGY\nMWO4//776devX1DqFBGRyBT1PePFi4sG8Z13woYNRYP44YcfZsyYMQwaNIhJkyYRExMTnGJFRCQi\nRXUYHzkC/ft72u3auceZqlXz7Nu5cydvvvkmDz74IBMmTKCcH1/RKCIi0Smqh6mHDIFvv3Xb1au7\npS4L5mNZazHGUKdOHTZu3MjVV1+tty+JiIhfRG03b/ZsmDPH0377bU+PODc3l379+jFhwgQA6tSp\noyAWERG/icowPnAABg/2tPv3h/bt3XZOTg69e/dmypQpHDlyJDgFiohIVIm6YeqsLOjdG44dc+1a\ntTwTuLKysujZsyfz589nzJgxPPbYY0GrU0REokdUhXFeHtx/P6xa5drGwJQpcPHFkJeXR7du3UhN\nTWXcuHEMHz48qLWKiEj0iKowHjkS3nnH0372WWja1G2XK1eOpk2b0qpVK4YMGRKcAkVEJCpFTRi/\n+mrR54kHDoQ//hFOnz7Njh07uOmmmxgxYkTwChQRkagVFRO4Zswo+hamTp1g/Hg4deokbdq0oVmz\nZhwruIksIiISYBHfM969GwYM8LTvvNMNVZ869QNt2rRh3bp1TJ8+nSpVqgSvSBERiWoRHcbWwgMP\neF6HeN11bs3pjIxjJCcn8/nnnzNnzhy6dOkS3EJFRCSqRfQw9ezZ8NFHbrtg5nTlyjB27FjS09OZ\nP3++glhERILOWGuD8o0bNmxoP/vsM799/tGjrid86JBr//738Je/uO3s7GzS09Np1KiR376/iIhI\nYcaYf1trGxb3tYjtGT/2mCeIExJg6NBv6dy5MwcPHiQuLk5BLCIiISMi7xl/8glMmuRpP/30ftq1\nS+Kbb75h165dVK9ePXjFiYiInCHiwjgryz1DXKBly68ZMyaJw4cPs3TpUm6//fbgFSciIlKMiAvj\n11+Hbdvcdnz8V2zdmsipUz+wfPlyDU2LiEhIiqh7xvv3uyUuCzz6aDy1al3BihUrFMQiIhKyIiqM\nH34YTp0C+Irrr8/miSeq88knn9CgQYNglyYiIlIir8LYGJNsjPnCGLPTGPN4MV+/wBgzJ//r640x\ntXxdaFmWL4d33wX4f8Ct1KkzjLg4MMYEuhQREZFzUmYYG2NigPHAXUA94B5jTL0zDusPHLPWXg28\nBrzs60JLk5UFQ4cCpAPNiI8vz8svPxTIEkRERH4yb3rGjYCd1trd1tosYDbQ8YxjOgJT87fnAc1N\nALukr70GX3zxKZCEMReyYsUarrnmmkB9exERkfPiTRgnAPsKtffn7yv2GGttDnAcqOqLAsuybx88\n++yPQGegCk899TG33VY7EN9aRETEJwL6aJMxZgAwAKBGjRo++cwFCyAj4wJgLnXrXs5TT13uk88V\nEREJFG96xgeAKwq1L8/fV+wxxphYoBJw5MwPsta+Za1taK1tWK1atZ9W8RkeegiWLYO6dW/jrbcu\nJzbinpwWEZFI5010fQrUMcZciQvdHkDPM45JBfoA/wK6AittAN9A0aIFbN4MMTGB+o4iIiK+U2YY\nW2tzjDFDgaWEOmlCAAAEy0lEQVRADJBird1ijHkO+MxamwpMBqYbY3YCR3GBHVAKYhERCVdeDepa\naz8EPjxj36hC25lAN9+WJiIiEh0iagUuERGRcKQwFhERCTKFsYiISJApjEVERIJMYSwiIhJkCmMR\nEZEgUxiLiIgEmcJYREQkyBTGIiIiQaYwFhERCTKFsYiISJApjEVERIJMYSwiIhJkCmMREZEgUxiL\niIgEmbHWBucbG3MY+NqHH3kJ8L0PPy9a6TqeP13D86dreP50Dc+fr69hTWttteK+ELQw9jVjzGfW\n2obBriPc6TqeP13D86dreP50Dc9fIK+hhqlFRESCTGEsIiISZJEUxm8Fu4AIoet4/nQNz5+u4fnT\nNTx/AbuGEXPPWEREJFxFUs9YREQkLIVdGBtjko0xXxhjdhpjHi/m6xcYY+bkf329MaZW4KsMbV5c\nwxHGmK3GmE3GmBXGmJrBqDOUlXUNCx3XxRhjjTGa1VoMb66jMaZ7/s/jFmPMO4GuMdR58fe5hjFm\nlTEmPf/vdJtg1BmqjDEpxphDxpjNJXzdGGP+kn99NxljGvilEGtt2PwBYoBdwFVAeeA/QL0zjhkM\nvJG/3QOYE+y6Q+mPl9ewGRCfv/2AruG5X8P84y4CPgbWAQ2DXXeo/fHyZ7EOkA5UyW//Ith1h9If\nL6/hW8AD+dv1gD3BrjuU/gBNgAbA5hK+3gZYDBjgNmC9P+oIt55xI2CntXa3tTYLmA10POOYjsDU\n/O15QHNjjAlgjaGuzGtorV1lrT2d31wHXB7gGkOdNz+HAM8DLwOZgSwujHhzHX8HjLfWHgOw1h4K\ncI2hzptraIGL87crAd8EsL6QZ639GDhayiEdgWnWWQdUNsZc6us6wi2ME4B9hdr78/cVe4y1Ngc4\nDlQNSHXhwZtrWFh/3L8KxaPMa5g/lHWFtfaDQBYWZrz5WbwGuMYYs9YYs84Ykxyw6sKDN9fwGaCX\nMWY/8CHw+8CUFjHO9XfmTxLr6w+UyGGM6QU0BJoGu5ZwYowpB4wD+ga5lEgQixuqTsSN0HxsjLnB\nWvvfoFYVXu4Bplhrxxpjfg1MN8bUt9bmBbsw8Qi3nvEB4IpC7cvz9xV7jDEmFjcscyQg1YUHb64h\nxpgWwJNAB2vtjwGqLVyUdQ0vAuoDq40xe3D3mVI1iess3vws7gdSrbXZ1tqvgB24cBbHm2vYH3gX\nwFr7L6ACbs1l8Y5XvzPPV7iF8adAHWPMlcaY8rgJWqlnHJMK9Mnf7gqstPl34QXw4hoaY24G3sQF\nse7Rna3Ua2itPW6tvcRaW8taWwt3372Dtfaz4JQbsrz5+7wQ1yvGGHMJbth6dyCLDHHeXMO9QHMA\nY8x1uDA+HNAqw1sqcF/+rOrbgOPW2m99/U3CapjaWptjjBkKLMXNIkyx1m4xxjwHfGatTQUm44Zh\nduJuyvcIXsWhx8tr+ApQEZibP/dtr7W2Q9CKDjFeXkMpg5fXcSnQyhizFcgFRlprNdKVz8tr+DDw\ntjFmOG4yV191UDyMMbNw/+C7JP+++tNAHIC19g3cffY2wE7gNHC/X+rQfxMREZHgCrdhahERkYij\nMBYREQkyhbGIiEiQKYxFRESCTGEsIiISZApjERGRIFMYi4iIBJnCWEREJMj+P50h0JPWRc/hAAAA\nAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results = pd.concat(dfs)\n", "import statsmodels.api as sm\n", "thresh = 0.001 # POSSIBLE BUG? several very small pivots -- fine for pvalues\n", "grid = np.linspace(0, 1, 101)\n", "fig = plt.figure(figsize=(8, 8))\n", "plt.plot(grid, sm.distributions.ECDF(results['pivot'][results['pivot'] > thresh])(grid), 'b-', linewidth=3, label='Pivot')\n", "plt.plot(grid, sm.distributions.ECDF(results['pvalue'])(grid), 'r-', linewidth=3, label='P-value')\n", "plt.plot([0, 1], [0, 1], 'k--')\n", "plt.legend(fontsize=15);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "cell_metadata_filter": "all,-slideshow", "formats": "ipynb,Rmd" }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }