Test Statistic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Statistic for Random Forest Independence Test\n",
    "In this notebook, we will work on generating a test statistic using our random forest estimate of conditional entropy. Our test is setup as so:  \n",
    "\n",
    "Given $(Y, X)$, where $Y$ is binary.  \n",
    "$H_0$: $X$ is independent from $Y$  \n",
    "$H_a$: $X$ is dependent on $Y$  \n",
    "\n",
    "This can be translated to:  \n",
    "$H_0$: $MI(X, Y) = 0$  \n",
    "$H_a$: $MI(X, Y) \\neq 0$  \n",
    "\n",
    "There are two methods we explore in this notebook:\n",
    "1. Permutation test\n",
    "2. Compute the cdf of H(Y) evaluated at H(Y|X; under alternative)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Permutation test\n",
    "Under the null, $X$ and $Y$ are independent. This means you should be able to permute around the Y labels and still get the same mutual information / coditional entropy. So, for the permtuation test, we will do exactly this:  \n",
    "\n",
    "1. Permute the Y labels\n",
    "2. Calculate Mutual information\n",
    "3. Repeat 1000 times\n",
    "4. How many times out of total trials did we get mutual information = 0?\n",
    "\n",
    "We will focus on the same example as the previous notebooks.\n",
    "Consider $(X, Y)$. Let $Y \\sim Bernoulli(1/2)$ ($Y$ = -1 or 1) and $X \\sim \\mathcal N(y\\mu, 1)$  \n",
    "\n",
    "Note, we are focusing on one dimensional case for now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampling Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAERCAYAAABy/XBZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd8VfX9+PHXuTt7D8gghCxGCCMYZAgiLuzS4qhabQv6dbVW/epP0WoLttoK+BVXFayDKqUOUFtR3CCgEPYII4QQEkhyswc3d53z+yMmJmQn9+bcm/t5Ph48yD33jPe5633OZ0qKoigIgiAIPkejdgCCIAiCOkQCEARB8FEiAQiCIPgokQAEQRB8lEgAgiAIPkokgCFKrcZdQ7FRmTvPaSi+XoL3EAlAJb/85S9JT09v/Td27FimT5/O73//e44fP95u3QcffJBrrrmm1/tesWIF7733XrfrtN1ncXEx6enpbNq0qe8n0sbRo0e59dZbWx+/9957pKenY7VaB7Tfnpz7WmZkZDBp0iSuu+66dufU19fRZrPxpz/9ia1bt7okznNfj9zcXO6///7W55999lmmT5/er31mZWVhsVg6XefSSy91yfvrKe644w6WLl3aYfnBgwe5/vrrmTBhAhdffDFvv/22CtF5F53aAfiyadOmcffddwNgt9spLS3l1VdfZf78+axZs4aMjAyg+QPflx/R559/nj/+8Y/drtPXffbGJ598wqFDh1ofz549m7Vr12IwGFx6nM60fS0VRaGhoYHVq1dz22238fbbbzN27Ng+77O8vJy33nqLCy+80CUxnvt6vPPOOxQVFblk31arlS1btjB37tx2y48ePUphYaFLjuEJli1bxueff05ycnK75Wazmd/85jdMnDiRZ555hu+++45HHnmEsLCwDq+J8AORAFQUGhrKhAkT2i2bO3cuV111FY899hhr164FIDEx0eXHdsc+zxUeHk54eLjbjwOdv5ZTpkzhggsuYO3atSxevHhQ4uiOO1+PCRMm8Nlnn3X4sdu4cSNpaWkcPXrULccdLKWlpSxZsoTNmzdjMpk6PP/WW29hMBhYsWIFBoOBWbNmUVVVxYsvvigSQDdEEZCH8fPzY8GCBezZs6e1KOjcoosvv/ySK6+8kvHjxzNjxgwef/xxmpqaAEhPTwfgj3/8I7/85S8BmDNnDk8//TRXXnklWVlZvP/++50WhxQUFHDdddeRmZnJlVdeybZt21qf66x4YtOmTaSnp1NcXMyzzz7Lc889R0VFBenp6Xz33XcdijxkWebNN9/kiiuuYPz48Vx++eW8++67rftrKYr6+uuvufnmmxk/fjwXXngh//rXv/r1WppMJpKSkjh9+nSnzzc0NPDnP/+ZCy+8kPHjx3P99deze/fu1lguuugiAG655RYefPDBDtu/9tprTJw4EYfD0brsZz/7GZdeemnrY1mWycnJ4e233273ejz44IOsW7eOvXv3tr6GLd577z3mzp3bGlN+fn6P53rxxRfz5Zdf4nQ62y3fuHFju3haHD9+nIULFzJhwgRycnJYsmRJ62cImu+iVq1axbx58xg3bhzZ2dncddddlJWVta6Tnp7O+++/z1133cWECROYPn06zz33XLdxti2qO/dfd8WW//d//0dJSQlvvfUWERERHZ7ftm0bM2fObHe3OWfOHA4cOEBNTU23MfkykQA80Pnnnw/Anj17Ojx38uRJfve733HeeeexcuVK7r77bt59912effZZgNa7hoULF/LYY4+1bvfKK6/w05/+lKVLlzJ16tROj7t06VJycnJ47rnniI2N5ZZbbuHYsWO9ivnqq69m/vz5hIaGsnbt2k6LXJ566imeeOIJ5s2bxwsvvMCMGTNYtGgRb775Zrv1HnroIc4//3xeeuklJk2axGOPPdarH8FzORwOSkpKiIuL6/Cc0+nk17/+NRs2bODOO+9kxYoV+Pv7c9NNN3Hw4EGio6Nbf8weeugh7rjjjg77mDFjBmfPnm0t9qqrq+PIkSMUFhZSVVUFwKFDh6ipqWHmzJnttr3jjjuYNWsWqamprF27lujoaACqqqpYuXIl//u//8vSpUs5ffo09957b4/nOnfuXGpra8nNzW1dVlRUREFBAXPmzGm3bllZGTfccAONjY0sX76cBx54gP/+97888MADreusXLmS5557jhtuuIF//OMf3HvvvXz77bf87W9/a7evxx9/nMTERF588UXmzZvHs88+y9dff91lnGvXru3y3+zZs7vc7pZbbmHdunWMGzeu0+cLCwsZMWJEu2Xx8fGtr4PQOVEE5IFarnAqKys7PHfgwAFsNhsLFiwgOjqanJwcDAZD61VoSzFIfHw8KSkprduNGzeOX/3qV90e96qrruKee+4BYPr06cydO5fVq1f3qvgkNjaW2NhYdDpdh6IYaP5hW716NXfeeSe333470PwD2tjYyIoVK7j22mvbxXHbbbcBkJWVxccff8ymTZvanc+5FEVpfQ1kWaa0tJS///3vVFZWMn/+/A7rf/nll+zbt4+33nqLyZMnAzBz5kzmzZvH888/zwsvvMDo0aMBSE5O7rTILCUlhdjYWLZv38748ePZuXMnw4YNo66ujl27djF37ly2bt1KamoqsbGx7bZNTEwkPDycmpqadq+XLMs888wzpKWlAc1l24sXL6a+vp6goKAuzz82NpbMzEw+++wzcnJygOY6mZycHIKDg9ut+/rrryNJEitXriQwMBCA4cOH86tf/Yq8vDxGjx5NeXk5v/vd77jhhhsAOO+88ygoKODzzz9vt68ZM2a0Jo6pU6e2vlezZs3qNM7OPhu9MWrUqG6fb2hoICAgoN2ylseNjY39OqYvEAnAy4wfPx6DwcA111zDj370I+bMmcOPf/xjNJrub+ZGjhzZ474vvvji1r91Oh0zZsxg165dA44ZYN++fdjtdi677LJ2y+fNm8e7775LQUEB/v7+QPOPfgt/f3+Cg4M5e/Zst/vfsGEDGzZsaLcsPDycxYsXk5mZ2WH9nTt3EhUV1frjD6DVarn00ktb76J6Y/r06ezYsYOFCxeyY8cOpkyZQmlpabsEcO7Vf3cCAgJaf/zhh6vYnhIANL9///rXv3j44YcB+PTTT7nqqqs6rLdjxw4mT56MyWRqTZotj7/77jtGjx7NI488AkBFRQXHjx8nPz+fXbt2Ybfb2+2r7XslSRIxMTHdvldti8vOpdVqkSSp23PsiqIoXW7b3336ApEAPJDZbAZoLRZoKyEhgVdffZW///3vvPbaa6xcuZK4uDgee+yxLq+6gE7LTXtaJywsjPr6+j5G37na2loAIiMjOz1mQ0NDawIwGo3t1tFoND22l58xYwa///3vW9cPDg4mPj6+yy9/XV1dh1igOWn05YpxxowZPProozidTnbs2ME111zDmTNn2Lp1K01NTezcubNd09ie+Pn5tXvcEr8syz1uO3fuXJYtW0ZeXh5hYWEcOHCAF154oV3ZPkBNTQ379u3rtJiuvLwcgPz8fB5++GH27NlDQEAAY8eOxWg0dngf+vpeddca64knnug0YfVGYGBgh/et5XFPidOXiQTggbZv3w7AxIkTO30+OzubVatW0djYyDfffMPf//537rnnHr799tsBNbmsq6tr97iysrK11YokSR0qGHu6Km8rJCQEaL6ibPuFrKioaPd8fwUHB3d6pd/d+i3HbquysrJPsUybNo3GxkZ27drFoUOHmDRpEmVlZaxcuZKtW7ei0+nIzs7u9f4GIjk5mZSUFD7//HNCQ0OZNGkSkZGR7SqYofkH8dJLL+WWW27psI/IyEhkWeb2228nNjaWDRs2MHLkSCRJ4qmnnhpwefo777zT5XMtdzv9kZSU1OE8i4uLkSSpQ92A8ANRCexhbDYbr776KlOmTOn0g7t+/Xouuugi7HY7AQEBrV/kxsZGGhoaAHosDupK2w5PNpuNTZs2tRaR+Pv7U1dX1+5Hf+fOne227+6448ePR6/X8/HHH7dbvmHDBsLCwkhKSupXzP01efJkzGZzu3NwOp1s3LixtZxaq9X2uJ/Q0FDGjRvHqlWrCA4OZtSoUUyYMAFZlnnppZc477zzukzK/X2fujN37ly++OILPvvsMy655JJO15k0aRIFBQWMGTOGzMxMMjMziYmJYfny5a0V2EVFRfziF78gOTkZSZKQZZlt27b16k6kOy3H6+xfWFhYv/c7depUNm/ejM1ma132xRdfMHbs2NZ6DqEjcQegopqamtaWPg6Hg+LiYv75z39SWlrK008/3ek2kydPpqKignvvvZfrrrsOi8XCiy++yKRJk1qv1oODg9m+fTsTJ05s7UzWG2+99Rbh4eEkJyfz+uuvY7FY+PWvfw00F3X89a9/5Q9/+ANXX301u3fv5qOPPmq3fXBwMLW1tXz11Vcd7l7Cw8O5/vrref7555FlmQkTJrBp0ybee+89Fi1a1KsfW1eaPXs248aN4+677+aee+4hMjKSN998k5KSEp566ingh6KDzZs3Ex8f36HzUYsZM2bwwgsvtNah+Pv7M2bMGPbs2cMf/vCHLmMIDg6mqKiIbdu2dXm311eXXHIJL730EjqdjieffLLTdW6++WbWrVvHnXfeyXXXXYfD4WhtwpuRkUFoaCjDhg3jlVdeISAgAFmWWbNmDYcOHeq0Db4nuP7661m9ejW33347N910E9u3b2f9+vWtreOEzok7ABVt3bqVa6+9lmuvvZZf//rXPPfccyQlJfH222+Tmpra6TYJCQm8+OKLlJWVcdddd/HQQw+RmprKihUrWte5/fbb+eqrr9oNM9AbDz/8MO+//z533HEHDQ0NvP76661NKNPS0li8eDG7d+/m1ltvJTc3t/WHssW8efNISUnhrrvu4ptvvumw/wcffJA777yTd955h9tuu42tW7fyl7/8hZtuuqlPcbqCTqfjlVde4YILLuBvf/sbd999NxaLhTfeeIPx48cDzeXKv/nNb1i7di1//etfu9zXjBkzANpVKLcU+3RXAXzNNdcQGBjIrbfe2q4H9UCMHTuW4cOHM3r06A4tj1okJCTwz3/+E5vNxu9//3sWLVrEsGHDeOONNwgLC0OSJFasWIFGo+F3v/sdjz76KIGBgSxfvhyLxcKRI0dcEqsrxcTE8Morr1BbW8tdd93Fp59+yp///Od2DRuEjiQxJaQgCIJvEncAgiAIPkokAEEQBB8lEoAgCIKPEglAEATBR4kEIAiC4KNEAhAEQfBRHtsRrLq6EVkeWAvViIhAKisbXBSR5xnK5yfOzTuJc1OHRiMRFhbQ84rn8NgEIMvKgBNAy36GsqF8fuLcvJM4N+8hioAEQRB8lEgAgiAIPkokAEEQBB8lEoAgCIKPEglAEATBR4kEIAiC4KM8thmo4BscRXuRa8vQxqaiiUhE0gzuxDCC4MtEAhBUocgy1u1vY9+3oXWZJmw4/j99BMngr2JkguA7RBGQMOgURaHp8xew79uAfswcAn7xFMYLfo1cU0rTV68g5igShMEhEoAw6JxFe3CcyMWQfRWmGTehCYrCkDELY87VOAp3Yt+/Ue0QBcEn9CoB2O12Fi9eTE5ODjk5OSxduhRZlnvcbuHChdxzzz0DDlIYOhTZifW7fyOFxGKYMK/dc/rMy9AlTcK6/d846ipVilAQfEevEsDy5cvZsmULL7/8Mk8//TTr169n1apV3W6zbt06Nm/e7JIghaHDfngTcs0ZjDlXI2naV0FJkoRx6nUgy9TtEncBguBuPSYAq9XKmjVrWLRoEVlZWUybNo377ruPN954o8uyWrPZzLJlyxg/frzLAxa8l+KwYdu5Dm1sGroRkzpdRxMcjTYxi/o9n6I47YMcoSD4lh4TQF5eHhaLhezs7NZl2dnZmM1miouLO93mT3/6E9dddx2jRo1yXaSC13MU7UWx1GGY9BMkSepyPcO4uTgba3Ec3z6I0QmC7+kxAZSVlREYGEhAwA9jTUdFRbU+d66PPvqIEydOcOutt7owTGEosB/9BikgDO3wMd2up40biz4iDtvBzwYpMkHwTT32A7BYLBiNxnbLDAYDADabrd3y6upq/vKXv/Dcc8+1rtNfERGBA9q+RVRUkEv246m85fycjbXUn9pPyNSfEBET0uP6tdmXU/nJKoLlSowxSe4PcJB5y/vWH+LcvEePCcBkMnX4oW95bDKZ2i1//PHHueyyy5gwYcKAA6usbBjw5AtRUUGYzfUDjsVTedP52fZ/BoqMPW5Kr2IOHz2Nyk9ewbz7G4yTIwYhwsHjTe9bX4lzU4dGI/XrornHBBAbG0t9fT0WiwU/Pz+guZIXICYmpt26//nPfzCZTLz77rvAD4li4sSJ7N69u8/BCUOH/egWNJFJaMPjerW+NiAETcwoHEV7ME7+qZujEwTf1GMCyMjIwM/Pj507dzJjxgwAcnNziY6OJi6u/Zd548b2TfeWLl2KzWZj0aJFLgxZ8DbO6tPIlScxTruhT9vpEidg2/EOcmM1moAwN0UnCL6rx0pgk8nE/PnzWbx4Mbt27WLbtm0sW7aMm2++GYCamhrq65tvi0aMGNHuX0BAAP7+/owYMcK9ZyF4NGfRXgB0SZP7tJ1uxESgufWQIAiu16vB4O6//36sVisLFy7EaDQyf/58FixYAMBvf/tb4uLiePLJJ90aqOC9HKf2oQmPRxMY3qftNGHDkYKicJzcjWH0bPcEJwg+TFI8dOQtUQncM284P8VmoeGNuzBkXoox55peb9dybk1b38Se9xWBNz2HpDf2vKEX8Ib3rb/Euamjv5XAYjA4wa0cpw+B7ESb0L9e4brECeC04zx9yMWRCYIgEoDgVs6ifaD3Qxub0q/ttbGpoNHhOHPExZEJgiASgOA2iqLgOLUfXdyYDgO/9ZakM6CNTsZ55qiLoxMEQSQAwW3k6hKUxiq0iQMbFFAbm4ZccRLF3uSiyARBAJEABDdyljSX2+vixw1oP9ph6aA4cZYdd0VYgiB8TyQAwW2cZ44gBUWhCRzYUA7amBSQJJyloh5AEFxJJADBLRRFwXnmSPPV+wBJBj80ESNwiopgQXApkQAEt5BrTqNYG9C5IAFAczGQs7xATBIjCC4kEoDgFi1X6664A2jeT1pzfwBzoUv2JwiCSACCmzjPHEEKCEMKinLJ/rSxac37FfUAguAyIgEILtda/h+b3u3Uj32hMQUhBUUhl59wyf4EQRAJQHADpa4c5WxNc7GNC2mjRuI0iwQgCK4iEoDgcj+U/2e4dL/a6JEojVXIZ2tdul9B8FUiAQgu5yzPB2MAmtBhLt2vJioZALlC3AV4IllRqKix0Nhkx0MHGRbO0b8BWgShG86yArTRo1xW/t9CGzmiuUNY+YnmUUIF1VmsDr49VMaW/WcoLm/A5pAB0Gk1jE0K47zRMUwZHY1OK641PZFIAIJLKTYLcnUJhuQpLt+3pDehCR0u6gE8gKIofHuwjLc+O0pjk4OE6EBmT4wjJTGMquqzmGua2J1vZu/xSj767iQ3zE0jY4SY1tPTiAQguFTzj7OCNmaUW/aviRqJs2gviqK4/A5D6B2r3cmq/xxi5xEzKXEhXDMnhVHDg5Ekqd2kKddfnMruYxX86/Nj/G3Nbn4+K5l5U0eI982DiAQguJSzLB9obrHjDtqokTiOfoPSUIkUFOmWYwhds1gdPPP2Xo4V1zJ/9iguOy8RjabzH3RJkpiUFsW4keH846M83v26gKo6KzdckoZGJAGPIBKA4FLO8uNoQocjGQPcsv+WxOI0n0AjEsCgstmdLFu7h5Ol9fzPT8dy3uiYXm1n0Gu59SdjCQ828fF3RQT46bjqAvfcIQp9I2pmBJdRFAW5vABNtPu+3JqIBNBokUU9wKBSFIXXPj7MidN13NaHH/8WGkni6tmjmDl+GP/ZepKtB864KVKhL0QCEFxGqTejNNW7rfwfQNLq0YTH46wsctsxhI427jjFtwfL+NkFyUxOj+7XPiRJ4peXppORGMprGw5TUtHo4iiFvhIJQHCZ1vJ/N94BAGjCE5tnCBNtzQdFUVk973x1nMlpUfzo/BED2pdOq+G2n47DZNDx6kd5yLJ4D9UkEoDgMs7y46AzogmLc+txtJGJKE31KBbRI9jdnLLMqx8dJsCk4+bLM1zSgic4wMD1F6dScLqOjTtOuSBKob9EAhBcxmk+gTYqCUnj3o+VJiIRALlCFAO528ffFXGyrJ4bL0kn0E/vsv3mjI5hYmok6zYXUFUn5npWi0gAgksoTgdyZREaNzX/bEsbkQAg6gHcrLreyodbCpmYGkl2Rv/K/bsiSRLXXZSKLCt8sEVU6KtFJADBJeTqYnA63Nb+vy3J4N88NLRIAG61fnMBTlnh2otS3bL/qFA/LpwUx+Z9ZzhTKSqE1SASgOASLTN1DUYCANBGJIo7ADcqNjfwzf4zzJkUT3Son9uO86NpSRj1Wt77usBtxxC6JhKA4BKyuQCMAS6bAawnmohElNoyFLsoP3aH974uwGTQ8ePpSW49TrC/gUumJLDzqJnTolnooBMJQHAJp7kQbdTIQRvnRRuRCCjIVcWDcjxfUlRWz578Ci49L8GlFb9dmTM5Hr1Ow8Yd4o5usIkEIAyY4rAhVxWjjUwatGNqIptbAoliINf76NuTmAxaLpocPyjHC/Y3MH1cLFsPlFHbaBuUYwrNRAIQBkyuLAJFRhM9OOX/AFJAOBj8RUWwi5VVnWXH4XIunBRHgMn9V/8tLp6SgMMp8+UucUc3mEQCEAasZXx+beQgJgBJEhXBbrDhuyJ0Wg2XTEkc1OMOiwhgQkokX+4uweGUB/XYvkwkAGHAnOZCJL8QpIDBnfBDEx6PXFWCoogfDFdosNjZdrCUaeNiCQkwDPrxZ00YTv1ZO3uOVQz6sX2VSADCgMkVJ9BEJQ36RB+a8HhwWFHqxQ+GK2zeexq7Q+aiSYNT9n+uccnhhAUZ2bxPjBQ6WEQCEAZEsTchV58ZtPb/bWnDm3+onKIl0IDJssIXu0rISAwlPjpQlRi0Gg3TM2M5cKJSDA8xSEQCEAakuQxeGdQWQC1aBp0TTUEHbm9+BZV1TcxR6eq/xYzxw1EU+Ga/uAsYDCIBCAMif98DWBOVNOjHlgx+zUNCiAQwYF/uLiEsyMjENHVnWYsO9WP0iDC27i8Vw30PApEAhAFxVhQi+Yei8Q9V5fja8HiRAAaootbCwRNVzBw/DK2bR3LtjZwxMZTXWCgsrVc7lCGvV++23W5n8eLF5OTkkJOTw9KlS5Hlzlte5OXlceONNzJhwgTmzJnDypUrXRqw4FlkcyEaFYp/WmjC45FrS1GcdtVi8HZb9pcCMCNzmMqRNJuUFoVWI7E9r0ztUIa8XiWA5cuXs2XLFl5++WWefvpp1q9fz6pVqzqs19DQwIIFC0hLS+ODDz7gkUce4cUXX2TdunUuD1xQn2JvQq45g1aF4p8WmvB4UGTk6tOqxeDNZEXhm31nGJ0URqQbB33ri0A/PeNGhrPjcDmyKAZyqx4TgNVqZc2aNSxatIisrCymTZvGfffdxxtvvNGhjO706dOcf/75PPzwwyQmJjJnzhymTZvG9u3b3XYCgnrUrABuoQlvnhtAFAP1T15hNZV1TcwcP1ztUNo5b0wMVXVWjpeIWd/cqccEkJeXh8ViITs7u3VZdnY2ZrOZ4uL2X7q0tDSWLVuGVqtFURRyc3PZsWMH06ZNc33kgurUrABuoQmJAa1ONAXtp2/2nyHApGOSypW/55qQEolep2H7oXK1QxnSekwAZWVlBAYGEhAQ0LosKiqq9bmuTJkyhRtuuIGJEycyb948F4QqeBqn+YSqFcAAkkaLJnR484Q0Qp9YrA52HzUzZXQMep1W7XDa8TPqGJ8cwa5jZtEayI10Pa1gsVgwGo3tlhkMzd3EbbbOR+6TZZlXX32V0tJS/vjHP/LEE0/wyCOP9CmwiAjXdEaJigpyyX48lZrnd6rmFH5xKW6Lobf7VYaPxFK436vea0+I9YvcImwOmXnTk10aj6v2dcHkeHau2U2t1UlqwuAOM9IVT3jfXKnHBGAymTr80Lc8NplMnW6j0WjIzMwkMzOTxsZGHnnkER544IHWxNEblZUNyPLAMn9UVBBm89BtSqbm+Sn2JuwVJUgjst0SQ1/OzeYXg7P+a8qLS5GMAT1voDJP+Vx+sq2QyBATEQE6l8XjynMbGR2IRpL4YvtJQk09/lS5nae8b53RaKR+XTT3WAQUGxtLfX09FouldZnZbAYgJiam3brFxcV8/fXX7ZalpqZit9tpaGjoc3CC53JWnAQUVVsAtdCGN/cIdlaXqByJ96iut5JXWM35Y2MHfQyn3gr005MaH8JuMTic2/SYADIyMvDz82Pnzp2ty3Jzc4mOjiYuLq7dunv37uXuu++msfGHqd0OHjxIREQE4eHhLgxbUJtcUQigah+AFprvxwQSLYF677tDZSjA1LExPa6rpompkZSYGymvsfS8stBnPSYAk8nE/PnzWbx4Mbt27WLbtm0sW7aMm2++GYCamhrq65tvi2bPnk14eDiLFi2ioKCAzz77jOXLl3PHHXe49yyEQec0FyIFhKlaAdxCCggHvQm5StwB9Nb2vDJGxAYxLMKzi8wmpDU3ONlz1KxyJENTrwrW7r//fqxWKwsXLsRoNDJ//nwWLFgAwG9/+1vi4uJ48sknCQgI4JVXXuHxxx/n5z//OUFBQSxYsIAbb7zRrSchDD65olDV9v9tSZLU3CNYtATqlbLqsxSW1nPNhSlqh9Kj6FA/4qIC2JNfwSXnDe4kNb6gVwnAaDSyZMkSlixZ0uG51atXt3s8cuRIXnnlFddEJ3gkxWZBrilFlzJV7VBaacPicJzYiaIoHlum7Sl25DW3rZ+SEa1yJL0zPjmCjTtOYbE68DOqXxk8lKg/8pPgdX7oATz4cwB0RRMej2JtQLGInqM92Z5XRkpcCBEhnbfi8zSZyRE4ZYW8k9VqhzLkiAQg9Jn8/RzAavYAPtcPcwOIeoDulFQ0Umxu5LzR3nH1D5ASH4LJoGV/QaXaoQw5IgEIfdZcARyOxi9Y7VDUQ9gbAAAgAElEQVRatbYEEvUA3dqRV4YEZHtJ8Q+ATqthTFI4+wsqRa9gFxMJQOgzZ0WhR7T/b0vjF4xkChJ3AD3YecRMWkIooYHGnlf2IJnJ4VTVWTld0djzykKviQQg9IliO4tSW+oR7f/PpQmPxynuALp0uqKRkopGr7r6b5GZHAHA/oIqlSMZWkQCEPrE+f0IoNroZHUD6YQmLA65qgRF6XyyIl+Xe6S59c+k79vWe5PwYBNxkQEcPCHqAVxJJAChT5zmAgCP6QPQliY8HhxWlAbxI9GZ3MNmUuJDCAvyruKfFmOSwjlaXIvd4VQ7lCFDJAChT+TyE0jBMUgm14zW6kpa0RKoS6VVZyk2NzAl3fuKf1qMSQrD7pA5Viya+rqKSABCnzjNJ9BGeU77/7Y0rYPCiXqAc+38vvhncrr3Ff+0SE8MRauROFgo6gFcRSQAodfkszUojVUemwAkgz9SQLi4A+jErqNmRg4LJjzYOzp/dcZk0DFqeDCHCkWHMFcRCUDotdYOYB5YAdxCEx4vRgU9R2VtEyfO1Hv11X+LMUnhFJXW02Cxqx3KkCASgNBrzvICkDRoIz13UC5NWBxyzRkU2aF2KB5j17HmkTS9sfXPucaMDEcBMSyEi4gEIPSa03wCTXgcks5zW5Fow+NBdiDXisnEW+w+aiYuMoDYcH+1QxmwkcOC8DNqOSTqAVxCJAChVxRF8egK4BYtFcFiSIhmdWdtHDlVMySu/gG0Gg1p8aEcFncALiESgNArSr0ZrI1oojy3/B9AEzocJElUBH9v77EKFGVoFP+0yBgRRlm1hep6q9qheD2RAIRecZblA6CNHqVyJN2TdAak4BhREfy9XUfNRASbSIzxvH4b/ZWRGAYg7gJcQCQAoVecZcdBb2oddtmTacPixATxQJPNwcHCaiamRQ6pSXISYgIJMOnIKxIJYKBEAhB6xVmejzZqJJLG8z8ymvB4lLoyFIdN7VBUdaCgCodTZlLq0Cn+AdBIEumJYeIOwAU8/9ssqE6xW5ErT3l88U8LTXgcKApyzWm1Q1HV7mNmAkw6UhNC1A7F5TISQ6mobaKixqJ2KF5NJAChR86KQlBktDGeP4k4tJkcxofrARxOmb35lUxIiUTrBXdtfZUxorkeQBQDDczQ+2QILtdSAayJ8ZI7gOAY0OpxVp5SOxTVHD1Vw1mrg4lDqPVPW3GRAQT56zlSVKN2KF5NJAChR3JZPlJwDBpTkNqh9Iqk0aIJG+7TdwC7j1Wg12kYmxSudihuIUkSaQmhIgEMkEgAQrcURcFZfhytl1z9t9CEJyBX+eYdgKIo7DlWwdikcIwGrdrhuE16QiiVdU1U1Ip6gP4SCUDollJvRrHUeU35fwtteDyKpQ7ZUqd2KIPuVHkDlXVNTEiNVDsUt0pLCAWai7uE/hEJQOiWt3QAO5cmPAHwzYrgPfkVSEBWytBOAPHRzf0BRDFQ/4kEIHTLeeYo6P1af1C9xQ8tgXyvGGj3sQqShwcTEmBQOxS30kgSqfGhHBF3AP0mEoDQLWfpEbSxqV7RAawtjX8Ikl+wz90BVNU1cbK0fsgX/7RITwylXIwL1G/e9a0WBpVsqUOuOYN2WLraofSLJjwep48lgL3HKwGYMMSLf1qkJzbXAxw5JfoD9IdIAEKXnGeOAKDz2gSQgFxVgiLLaocyaPbmVxAVamJ4ZIDaoQyKxOjm+QGOnhITxfeHSABCl5ylR0FrQBOZpHYo/aINjwenDaXeNyaHsdqcHCqsJitlaA3+1h2NRmJUXIhoCdRPIgEIXXKeOYo2ZhSSVqd2KP3SUnHtKz2CDxU2D/7mK8U/LdITQjld0Uj9Wd8e/K8/RAIQOqXYziJXFnlt+T+AJmw4SBrkyiK1QxkUe/Ir8DNqW9vH+4qW8z1WLIqB+kokAKFTztJjgOLVCUDSGdCEDsfpAwlAVhT2Hq8kMzkCnda3vtZJscHotBpRDNQPvvVJEXrNcToPNFq00Z49BWRPNBEJyD5QBFR4pp66RhtZo3yr+AdAr9MwaniwSAD9IBKA0ClnyUG0MalIOqPaoQyINjIRpbEKuale7VDcam9+BZIEmaMi1A5FFakJoRSVNWCxOtQOxauIBCB0IJ+tbZ4AJn6c2qEMmCZiBMCQvwvYe7yClLgQAv30aoeiivSEUGRF4fhpUQ/QFyIBCB04Tx8CQBc/VuVIBk4T8f2YQJUnVY7EfarrrRSVNQz5sX+6MyouGI0kif4AfSQSgNCBo/ggGANar569mcYUhBQQjrNi6FYE7zteAUCWjxb/AJgMOhJjAjkm6gH6pFcJwG63s3jxYnJycsjJyWHp0qXIXfSuzM/PZ8GCBWRnZzNr1iz++te/YrWKcTq8haIoOEsOoosb43Xj/3RlqFcE782vJDLEd3r/diUtIZSCM3XYHb7T83ugevUNX758OVu2bOHll1/m6aefZv369axatarDeo2Njdxyyy1ERESwdu1a/vrXv7Jx40aWLVvm8sAF95BrzqA0VqON8/7inxbaiETkmtMojqHXUcjucHLoZBVZo3yn929XUuNDsTtkTpYO7Qp/V+oxAVitVtasWcOiRYvIyspi2rRp3HfffbzxxhsoitJu3W3btlFfX8/jjz/OqFGjmDp1KnfffTcffPCB205AcC1nyUFgaJT/t9BEJIIiI1efVjsUlztcVIPNLvts65+2UhNCADhaLIqBeqvHBJCXl4fFYiE7O7t1WXZ2NmazmeLi9iMtZmZm8vzzz2Mw/DAOuSRJNDQ0dEgWgmdyFO1FColBEzR0JhPXRiQC4ByCFcH7jldi0GnISPSt3r+dCfY3MCzCX9QD9EGPCaCsrIzAwEACAn4oX4yKimp9rq2YmBhycnJaHzudTlavXk1OTo7P3556A8VmwXk6D92IiWqH4lJScBTo/ZArhlYCUBSFfccrGD0iDIN+6M792xep8SHkl9QiiwvOXulxlC+LxYLR2L4zUMsVvs3WfZnqkiVLOHz4MP/+97/7HFhERGCft+lMVFSQS/bjqVx5fg2H9tEgO4nMmo6fB7xurjw3x/BRKNVFHvN5cEUcxeX1mGuamH9RmsecF6j7nZs8JpZNe89gcULSMNfH4Umvsyv0mABMJlOHH/qWxyaTqdNtnE4nf/rTn3j33Xd55plnyMjI6HNglZUNyPLAsnhUVBBm89CtEHL1+Vn2b0UyBVFvjKNB5dfN1efmDEnAfuBTykurVR/d1FXn9tWO5qatI6MDPOZzrvZ3Ljak+Tfpu30lBOhcW+qg9rl1R6OR+nXR3GMRUGxsLPX19VgsltZlZrMZaC7yOZfdbufee+9l3bp1rFixgrlz5/Y5KGHwKbIDR9E+tIlZQ6b5Z1vaqJEgO5Crh84MYfuOVxIXGUBkiJ/aoXiMyBATYUFGMTJoL/X4Tc/IyMDPz4+dO3e2LsvNzSU6Opq4uLgO6z/66KN8/fXXvPTSS1x00UWujVZwG+eZo2A7iy5paJX/t9BGJQHgNBeqGoerWKwOjp6qEa1/ziFJEqnxzRPEiIYnPesxAZhMJubPn8/ixYvZtWsX27ZtY9myZdx8880A1NTUUF/ffFv09ddf895773H//feTmpqK2Wxu/Sd4NkfhLtDq0cV5//g/nZGCosDgjzxEEkDeyWqcssL4ZJEAzpUaH0p1vZXK2ia1Q/F4vSoMvf/++7FarSxcuBCj0cj8+fNZsGABAL/97W+Ji4vjySefZMOGDQAsXryYxYsXt9vHvn37OlQmC55BkWUcJ3LRxY9D0g/N90iSJLRRI3FWnFA7FJfYd7wSk0FLSnyI2qF4nJYJYo4W1xAZKorHutOrBGA0GlmyZAlLlizp8Nzq1atb/37yySd58sknXRedMCicZw6jnK1BlzpN7VDcShuVhG3fxygOG5LO0PMGHkpRFPYXVDI2KdznJn/pjbioAPyMOo4V1zJt3DC1w/Fo4tMjYD+2FfR+6BKz1A7FrTSRSSA7kau8uyK4pKKR6nqrKP/vgqZNPYDQPZEAfJzisOI4kYs+Odurr4p7o7UiuKJQ1TgGav/xSgAyRfl/l1LjQzhTeVZMFN8DkQB8nOPkHrA3DfniHwApMBLJFISz3LvrAfYXVJIQHUhY0NCsr3EFMVF874gE4OPsx7YhBYR59eTvvSVJEproZOTyfLVD6TeL1cGx4lpx9d8DMVF874gE4MPkhkqcp/ahTzkfSfKNj4I2JqV5yOumBrVD6ZdDhc3NPzOTw9UOxaPpdRqShwVxTIwM2i3f+NYLnbLnfQWKgn7MhWqHMmi0MSkAOMsLVI6kf/YXVOJn1DIqTjT/7ElqQignSxtosomJ4rsiEoCPUpx27HlfoRsxYUgN/dwTbdRIkCScXlgM1NL8c8wI0fyzN9JaJ4qvUzsUjyU+RT7KUbADpake/VjfGq5D0pvQhCfgLPO+BHBaNP/sk5S4ECQJMT9AN0QC8FG2g58hhcSijRujdiiDThuTgrO8AKWLea091f6CKgDGjRTl/73hZ9SREB0oKoK7IRKAD3KcPoxcXoBh7EU+U/nbljZ6FNibkKtL1A6lT/YXVBIXFUB4cOfDsAsdpcWHUnC6DofTu5L9YPG9b7+Abdf7SH4h6DNmqR2KKlorgr2oGKjJ9v3on6L5Z5+kJYRiExPFd0kkAB/jOHME5+k8DFnzhnzP365IwdHfdwjzngTQMvqnSAB9k9pmYDihI5EAfEzz1X8w+jGz1Q5FNZIkoY1NbZ4DwUscKKjCaNCSKkb/7JOQAAOx4f4cLRIJoDMiAfgQx+k8nCWHMGRdjqTz7WEEtMMyUOrNyPUVaofSox+af4aJ5p/9kJYQytHi2gFPMTsUiU+Tj1BkGevWt5ACI9CP8a2mn53RDh8NNA+F7elKq85SUdvEOFH80y/pCaFYrA6Kzd7Z+9udRALwEfbDXyNXncI49VqfLftvSxMeh2QMxHHa8xNAS/PPTNH8s19aJ4gRzUE7EAnAByjWRmy576GNTUM3cora4XgESdKgHZbuFXcABwoqGRbhL2a36qeIEBMRwUaOipFBOxAJwAdYt7+DYm3AOO16JElSOxyPoR0+GqW+Arnec+esttqdHC6qYdxIUfwzEGkJoWKi+E6IBDDEOc4cwZ73Jfpxl6CNTFI7HI+iHZ4BgNODi4EOn6zG4ZQZL4Z/GJC0hFDqGm2UVp1VOxSPIhLAEKY4bFg3vYoUFIkx+yq1w/E4mrA4JFOQR9cD7C+oxKDXtJZjC/3T8vodEfUA7YgEMITZdn2AXFuKacbNSHrfbvbZGUmS0A7PwHn6kEcWDbQ0/xydGIZeJ76qAxEb7k9IgEH0BziH+FQNUc6yfGx7/4subSa6hEy1w/FYuvhMlMZqj5wovrTqLOaaJjH6pwtIkkR6YihHRD1AOyIBDEGKw4rlq5VIAeGYpl2vdjgeTZs4HgDHqb0qR9JRa/NP0f7fJdITQqmut2KusagdiscQCWAIsn73NkptGabZC5EMoulgdzT+oWgiR+As2qd2KB3s/775Z5Ro/ukSaYlhABwRxUCtRAIYYhwlh7Af/Az9uIvRfd/bVeieLmE8zrJjKNZGtUNp1WRzcKSoWlz9u9DwCH+C/PWiIrgNkQCGEMXaSNNXq9CExGI8b77a4XgNXWIWKAqOU/vVDqVVXmE1Dqcimn+6kCRJpCWEijuANkQCGEKatr2FcrYG04W3+vxgb32hiUpubg5a5Dn1AHuPV2IyaEXzTxfLSAyjsq6JClEPAIgEMGQ4ivbgOLoFw4Qr0EYnqx2OV5E0GrTx43Ce2o8iO9UOp7X559gkMfm7q6UnNifUvKJqlSPxDOLTNQQo1kaaNr+OJiwOw6SfqB2OV9IlTUKxNnhEr+BT5Q1U11tF8Y8bxEUGEOSv5/BJUQwEIgEMCdZv1zYX/cxagKTVqx2OV9IlZoHehOP4t2qHwt7jlQAiAbhBc3+AMA4XVYv+AIgE4PUcpw9jP7IJw/jLRdHPAEg6A7qkSdhP5KI47arGsu94BSNigwgJFPU47jA6sbk/QLmoBxAJwJspDjvWb15HCorEMPmnaofj9fSjpoLNgvPUAdViqG20UVBSx8SUSNViGOoyRjT3Bzh8UtQDiATgxWq+fR+55gym6b8UrX5cQBs/BskYiF3FYqB9+RUowIRUkQDcpWVcoMOiOahIAN5KrjdTs+VddCOzm8uvhQGTNDp0ydk4Tu5GsTepEsPuYxVEBBtJiA5U5fi+QJIkMkaEcfikqAcQCcBLWb9dC5KE8Xwx1o8r6VOng8OGPX/w7wKsdieHCquYkBIlJu5xs9EjwqhttHG60rfnBxAJwAs5TufhOJFL6PlXogkU88S6kiYmBU1EIvYDnw361WFeYTU2hyyKfwbBmO/rAQ4VVqkcibpEAvAyiixj3fYWUmAEIVNFm39XkyQJ/diLkKuLcZ45MqjH3n3MjJ9R29pZSXCfyFA/okP9yCv07YrgXiUAu93O4sWLycnJIScnh6VLlyLLcrfbVFdXM23aNI4fP+6SQIVmjmNbkCtPYcy5Bo2Y5MUt9CnngzEA+8HPBu2YTllm97EKxo+KFL1/B8mYpOb+AM4efsuGsl590pYvX86WLVt4+eWXefrpp1m/fj2rVq3qcv2qqipuvfVWKisrXRao8P0Uj7nr0EQlo0s+T+1whixJZ0CffgGOwl3IDYPzGT5aVEODxc7ktKhBOZ4Ao5PCabI5OXGmXu1QVNNjArBaraxZs4ZFixaRlZXFtGnTuO+++3jjjTc6LSPdsmULV155JXa7up1phiLbgU9RGqsw5lwjKgndzDD2IpAkbLs+HJTj5R41Y9BrxOxfgygjMRQJyPPheoAeE0BeXh4Wi4Xs7OzWZdnZ2ZjNZoqLO06j99VXX3HjjTeyYsUK10bq45SmBmx7/oM2MQvd8Ay1wxnyNEGR6EfPxn5kE3LNGbceS5YVdh0xk5kcgVGvdeuxhB8E+RtIjAnikA/XA/SYAMrKyggMDCQgIKB1WVRUVOtz53r44Ye55ZZb0GhEOaYr2fZ+BLYmMc7/IDJM/Alo9Vh3vOvW4+SX1FLbaCM7PdqtxxE6GpMURn5JLU02h9qhqKLHX2mLxYLR2L6y0WAwAGCz2dwTldCOfLYG24HP0KVMRRueoHY4PkPjH4Jh/GU4TuTiLHdfY4bcw+XotBox+JsKxo0MxykrPtsrWNfTCiaTqcMPfctjk8nknqiAiAjX9ISMigpyyX7UVLHx3yA7GHbx9ejD25/PUDi/rnjCuclz5nPq6Cbs37xK9G/+5rKWVy3n5nTK5B41M2VMDInxYS7Zt9o84X3rrdAwf4zv7ef4mXouPn9kj+t707n1Ro8JIDY2lvr6eiwWC35+zZNTm81mAGJiYtwWWGVlA7I8sI44UVFBmM3eXcMvN1TSuHMj+rQZ1DiDoM35DIXz64onnZth5m+wbFhGyUevYZo28J7Xbc/twIlKauqtTBwV4THnOxCe9L71VnpCKDsOlfLzmd0nAE8+N41G6tdFc49FQBkZGfj5+bFz587WZbm5uURHRxMXF9fnAwp9Y9v9IaCIiV5UpEvIRD/2IuwHNuIodu1Iod8eLMPPqCMrRRT/qGXcyHDKqy2UV/vesBA9JgCTycT8+fNZvHgxu3btYtu2bSxbtoybb74ZgJqaGurrPTMreju5vgL7kc3o0y9AEySGB1CTMecaNGFxWD59DmdFoUv2abU72XnUTHZ6FHqdaP2jlszk5uR74ITvNQftVVOd+++/n5ycHBYuXMi9997LlVdeyYIFCwD47W9/y5///Ge3Bumrmq/+JQwTf6R2KD5P0hnxu/w+JGMAlo+WuaRp6J5jFVhtTqaOjXVBhEJ/RYf5ERli4kCB7yWAHusAAIxGI0uWLGHJkiUdnlu9enWn28THx3PkyOCOpTKUyPVm7Ee+QT96FppAUTzgCTSB4fjPu5+zH/6Fsx/8BdPcO9ANH93v/W3Zf4awIKMY+0dlkiSRmRzB1gOl2B0yep3vNGH3nTP1MrbdH4IkYZggrv49iSY0Fv8fL0IyBWL571PY9n2CovR9LBlzjYWDJ6qYOX4YGtGrW3XjR0VgtTs5UuRbncJEAvBAcp0Z+5Et6EfPFsM9eyBNaCz+P3sUXWIW1m/XcPaDv+CsLunTPjbtPQ0SXJA13E1RCn0xekQYBp2Gvfm+NX6ZSAAeyLb7Q9BIGCZcoXYoQhckgx+mS36HafYtyDVnOPvOozRtfROlqaHHbR1OmW/2nWF8cgThwe7rSyP0nkGvZUxSOHuPV/jULGEiAXgYua4c+9Fv0I++EE3A0OgYNFRJkoQ+bToB1zyBPn0m9oOf0bD2/2E7+DlKN0MMbz9YSm2jjVkTRTNqTzI+JYKK2iZOVzSqHcqgEQnAwzRf/WsxZM1TOxShlzR+wZgu+BX+Vy1GG56Adctqzr73GM7ygk7Xf3/TcSKCjWQmi+I9T5I1qrmp9Z78CpUjGTwiAXgQubYU+9Hvy/7F1b/X0UYk4Pej/4dp7p0o1gbOvv841h3vojh/GGgsv7iWQyequGRKIloxYKJHCQsyMiImyKfqAcQn0INYd74PWp0o+/dikiShT55CwPzH0aVOw7b7Qyz//RuypQ6ADd+dJMhfLyp/PdTE1EiOl9RS22BVO5RBIRKAh3BWl+DI/xbD2Llo/EW7cG8nGQPwm70Q05zbcJpPcHbdnyg9fpTdxyqYN30kRoPo+euJJqVHoQC7jprVDmVQiATgIWw714PeiD7rcrVDEVxInzIV/588DIqM7ovlJBmq+fGMZLXDEroQFxlATLg/O0UCEAaL01yIo2AHhsxL0JiG1nCzAmijkjCfdxeNDi13BX+Kqb7jTHqCZ5Akiez0KA6fbJ6jeagTCcADWLe/jWQMxDBeXP0PRYqi8K/vanjV8SP0/oGc+dfjyHXlaocldGFSWhSyorD72NC/CxAJQGWOkkM4Sw5imPhjJIOf2uEIbpB7xMyx4lpmzxxPwLz/BUXm7Iblveo0Jgy+pNggIoJN7DwiEoDgRooiN1/9B0agH3Oh2uEIblDXaOOfG48wIjaoedyf0Fhir34Qpb4Cy6fPosi+ORetJ5MkiSkZ0Rw8UUX92aE97a1IACpy5H+LbD6BMftKJJ1B7XAEF1MUhTc+OYLF6mThFaNb2/2bEkZjmvUbnGeOYN3+jspRCp2ZOjYGp6yQe3hoF9WJBKASxWHFuv0dNJFJ6FKnqR2O4AZf7S5h11EzV12QTFxU++n69KnT0I+Zg33fx9hP5KoUodCVhOhAhkcGsO1QmdqhuJVIACqx7f0YpbEK4/m/QJLE2zDU7C+o5M1PjzF+VASXTEnodB3j+b9AE5VM01evINcN/fJmbyJJElPHxJBfXEtFjUXtcNxG/PKoQG6owrb3v+hGZqMblq52OIKL5ZfU8sL6A8RHBXDbT8ei0XQ+3r+k1eM393YALF++hCI7BzNMoQdTx8QA8O0QvgsQCUAF1m1vgQLGnGvVDkVwsT35FSxds5sQfwN3X52FydD9pHuaoChMM29GLsvHtuuDQYpS6I3IUD/S4kPYsv/MkB0iWiSAQeYo2ofjRC6GST9GExyldjiCizicMus3F/Dcu/sZHhnAol9OJizI2Ktt9SlTvx836AOcZflujlToi5lZwymrtnCkqEbtUNxCJIBBpDhsNG1ZjSZ0mOj0NYQUnK7jiX/u5IMtheSMieaB6ycSHNC3Vl2m6b9ECgjH8uVKFLtvDETmDaZkRONv1PH13tNqh+IWIgEMImvueyj1ZowzbkLSdl80IHi+0xWN/P39Azz+Ri6VtU3c/rNx3PLjsT0W+3RGMvhhmn0LSl051u/WuiFaoT8Mei3TxsWy80j5kBwhVPwKDRJnWT72fZ+gHz0b3fDRaocjDMDJ0nr+++1Jdh4uR6/X8KNpSVyek4ifcWBfJ93wDPSZl2Df/wm6ERPQJYx3UcTCQMyaMJzPdhbz+Y5TzBwXo3Y4LiUSwCBQHDYsX61CCgwXFb9eLL+klg++OcGBE1X4GbXMO38El0xJIMjfdZ34jFN+jrP4AE1f/4OA+Y8jmQJ73khwq7ioQFLjQ/jv1hNMGxM1pCbyGTpn4sGs365FqS3FNGuBGO/HCxWbG3j633v5y+qdnCyrZ/7sUTx1+3R+PmuUS3/8ASSdAdOFt6I01dP0zRtDtvWJt7n0vETKq86Se3ho9dcQdwBuZi/Ygf3Q5+jHX4Yuboza4Qh9YLE6ePfr43y5uwR/o46rLxzFnInxbp/MRRs5AsPkn2Hb8S6OERPQi57iqpuQGklcVAAff1fEeaOjkaTO+3Z4G5EA3EiuM9O06R9oopIxTpmvdjhCH+wvqOS1DYepabBy4cQ4fjYzmUA//aAd35B1Bc6ifTR98wbamBQ0wdGDdmyhI40kceXsFJ57ey+HT1YzOilc7ZBcQhQBuYlib8Ly6QoA/C66XbT68RJ2h8y/Pj/G0//ei59Rx6JfTubGS9IH9ccfQNJoMM35H5AkLF++LHoJe4ALJycQHGDgw62FQ6ZoTiQAN1AUmaYvVyJXFeM353bR4ctLVNU18eSbO9m44xQXTYrnsV9lM2p4iGrxaIIiMc34vpdw7jrV4hCaGfRafnT+CA4X1XDgRJXa4biESABuYMtdh6NwJ8aca9EliqZ83uDoqRoWv7aDM5VnufPKTG64JA29Tv2J2/UpU9GnX4Btz39wFO1ROxyfN3tiHJEhJt756jjyELgLEAnAxWwHPsW2+0P06Regz7xU7XCEXti87zRPrdmNn0nPH27OZnK6Z92xGaffiCYiEcuXK5Hrh1YrFG+j02q48oJkTpU38N0QGCROJAAXsh/dgnXrm+iSJmGcefOQaSkwVMmKwttf5fPqR4dJTwzlDzdNZlhEgNphdSDpDPhdfBcoMpaNK1BsQ3d4Ym+QMyaGEbFBrP0in8Ym7544Xgz1IHwAAA8dSURBVCQAF7Ef2UzTV6vQDh+Nac5tSBr1iw+ErtnsTl5cf4AN3xYxe8Jwfn91Fv6mwa3o7QtNcDR+F92BXFWC5Yu/i0phFWkkiV9dlkHDWTtvf+ndg/eJBOACtgOf0fT1K2jjx+J32e+7nN7ROsCKvLbbn7uvlseNHzzR6/30Jh5r7rp263e1Tdvn2sZgzV1H4wdPtD5X/+Z9ncbe9h9Ayeo/tO6nZftz99XZOXX2+Ny/axusPLnyc3YdMXPNhSlcE5GHTtv1V2H95oJeLWurqzg707JeV/tsWa5LyMQ4/UacRXvZ9tr/dWiJ8uSbOzvdtrP9tl237d/rNxd02E/L45b/z12n7f7f+uRw6+O2/7f8/b8vbOmwzbn7Onf7tsfv7XvR2X7O1dnr1VsjYoO45LwENu09Q16h91YIiwQwAIos07RtDdat/0Q3YiJ+l96NpOt6CGDbrvcHdLy225+7r5bHcumRXu+nN/HYdr3fbv2utmn7XNsYbLveRy498sN2jZWdxt72H4C16FDrflq277CvTs6ps8dt/y7Yvoklb+Ryuk7hzqsyuSwnEfvu7l+HD7YU9mpZW13F2ZmW9braZ9vlhjFz+NwyhkzHfqzfrW2XBI6equ10287223bdtn9/sKWww35aHrf8f+46bfe/ZuOR1sdt/2/5u6rO2mGbc/d17vZtj9/b96Kz/Zyrs9erL346YyQx4f68/OEharx0oDiRAPpJsTZi2fgM9v2foB87F9PFdyFpPbcIQYBvD5byf3WXoyhwd/AnTErzrMre3vrAMplvmtKw7/sY2453h0ybdG9j1Gu588pxWGwOXlx/AIdTVjukPhMJoB8cpcdofPdRnKcOYJxxE6bpN4oyfw9mVzS8ufEoL394iERdJY/+agrxOu+9bQeJd87moM+YhW3Pf7Bu+geK7FA7KJ8UHxXIry7P4FhxLa9tOOx1TUNF99Q+UOxWrDvXY9//CVJgBP4/fRhtdLLaYQndKDE38Pe6eZTsKubi7AQuP76akIBfUK92YAOkIGGc+Ssk/xBsuz5AbqgiQBqrdlg+aeqYWMqrLazffAKdVsNNl6Wj8ZIWgCIB9IKiKDhO5DaXudZXoE+/AOP51yEZ/NUOTeiCQ9HwwZYTfLilEJPiz+/mj2dCSiT1Bd51hdYdSZIwZl+FJiiKps2v8/9CCnAUJ6KLH6d2aD7nx9OSsDtk/rvtJBarg9/MG+32QQNdoVdFQHa7ncWLF5OTk0NOTg5Lly5Fljsv76qvr+e+++5j8uTJzJw5k1dffdWlAQ8mRZFxFO3h7PrFNH32fHN77B8/hGnWb8SPv4dSFIXcw+U8UfsT1m8+weT0KB4K+YAJKZFqh+Y2+vSZ+F/5KGdlA5aPlmL58mXkxmq1w/IpkiRx1QXJXD17FLmHy/nz6p2UmBvUDqtHvboDWL58OVu2bOHll1+msbGRBx54gODgYG699dYO6z788MOUl5fz5ptvUlhYyEMPPUR0dDRXXHGFy4N3F9lSh+P4d9gPfo5cW4oUGIFp1gJ0qdORhtBkEEOJQ9Gw7WApG74totjcQKzWyT3XZJGZHEH9y01qh+d22ohEltVdwYo5tdj2fYLjRC760RcSpgmgWhaTygwGSZK4fOoI4qMDefmDg/zx1R1ccl4CV0wd4bF9THpMAFarlTVr1vDMM8+QlZUFwH333ceyZcu45ZZb2vV2LSkpYePGjfznP/8hJSWFjIwM8vPzef311z06ASiKglJbhqN4P45T+3AWHwLFiSZqJKY5/4MueQqSRpSWeRqHU+bYqRp2Ha1gW83VnP3wEMMjA1hwxWjGbl1EaPLP1Q5xUNnRYTzvavQZs7Dmvof9wKf8IUThsH049mMmdIlZSEbP6+k81GQmR/DnW6fy9pf5bPi2iC92lTArazjTxsWSEB3oUSME9PirlpeXh8ViITs7u3VZdnY2ZrOZ4uJiEhISWpfv2bOH0NBQUlJS2q374osvYrfb0evVzYKKoqD8//buPKaJPo0D+LcHl6IFdqkg+iIxIBgUENQgaMKKGg9YjHggsroBAuFQY1SwakNUXGIkrxKqsYoXmrxRAxhEVzS7EBMIrwVNECIGlXgth3hVhZbS2T+Q0nIf1aHt8/lH5zfDzHcK9sGZ6fNTfAXzpQXqL81Qf3yLztbXUDc/B9PedVuQI3CA+dwV4LsuAs9uGqt5SQ8Fw8eHzoloUgvwTmWLN388wvN3X9Cu7ASfx8Ucs3cICgvDbBc7cDkcyMuN51r/SHEnC2H1t3ioF6zHzZzz8DV/ifb/SgEOB9y/OIPn4IqF5l/R2ewMziR7AKb7Wv0skyeYI3r1bCzzm45///kK92VvUPzwNYS2VnD/zRau0wSYZm8NB7sJrN4rGLIANDU1wdraGhMn9vzmYG9vr1mnXQCampogFOpOXGFvbw+VSoX379/D0dFx2MG43JFXSVVDJTrfvwI6O8B0dqCFo4bi+1egow2M4nvXm7z2R+g5XPAm2YPvHgjej38Y3Ek/71oxX2A/qvPq7+v5gq7vgfYyl8sZ1jFGum33cXofc7DteuccbB/dYy9U9qhRTof6P/X4zg9FO2OG9muPIcc/0QZzgAeAB3DMgCl8HlYs+A1u023gOs0GHTeKMFHrOn/vDIO9TkO9DkJbqz7r+xsb6PXoT+/z53I5A+6z97jQ1qrffffejjv5r6iwDMSfCMS/wp2gaqxDZ1M91O9q8Q8nJZQ/2hicmMFDe2EFOJbW2DVdgXbGDMqKFnDMLGDmvkSzX+0/tY+vfVztdf1t33tfw/36/o470HkPtJ+hXtfhGOn2MxwnI/7vntiyvANPXn5AbcMHNPxPjqeveu7RTLDgY7K1OSZa8DF7hh0We00d0TFGk6sbhxniUyQFBQU4duwYysrKNGNqtRoeHh64cOECFi3qma5OIpGgpKQE169f14y9fv0awcHBKC4uhrOz86hCEkII0b8h72haWlpCqVTqjHUvW1pajnpbQggh7BqyADg4OEAul6OtracFbUtLV0/yKVOm9Nm2e1235uZmmJmZwdbWVh95CSGE6MmQBcDd3R1WVlaorOzpnCeTySAUCuHk5KSzrbe3N1pbW/Hy5UvNWGVlJTw9PWFu3n+HTEIIIewY1iWg8PBwHDp0CFVVVSgvL0dmZia2bt0KAPj06RPk8q4naJycnBAUFISUlBTU1tbi7t27yMnJ0WxLCCFk/BjyJjDQ9VmAI0eOoKioCBYWFggPD8euXbvA4XAQFRUFJycnZGRkAOgqCGKxGKWlpRAIBIiOjqYCQAgh49CwCgAhhBDjQ30NCCHERFEBIIQQE0UFgBBCTJRJFIDS0lKEhIRg7ty5WLt2LWQyGduR9K6qqgoeHh54/vw521H0gmEYSKVSBAcHw8fHB5s2bcKjR4/YjjVqI2mpbmgaGxuxfft2LFy4EAEBARCJRPjy5QvbsfROLBZjw4YNbMfQK6MvADU1NUhMTERISAgKCwsxf/58JCQk4ONH4+mXrlQqsX//fqN5QwGA3Nxc5OTkIDU1Ffn5+Zg3bx6io6PR2NjIdrRR0W6p/vvvv6OgoADnzp1jO9aYqdVqJCYm4tu3b7h06RJOnz6Np0+fYt++fWxH06uKigpcu3aN7Rj6xxi5HTt2MElJSZpllUrFrF69mnnw4AGLqfTr+PHjzObNmxk3Nzemvr6e7Th6ERoaypw5c0ZnbPny5Uxubi5LiUavvb2d8fLyYkpKSjRjeXl5TEBAAKNWq1lMNna1tbWMm5sb09zcrBmTyWTMrFmzGLlczmIy/fn+/TuzbNkyJiIiglm/fj3bcfTK6P8HUF5ejpUrV2qWeTwebt26hcDAQBZT6U9NTQ3y8/ORmprKdhS9OnDgAEJDQ/uMG+KlhaFaqhsyR0dHnD17VtMhGOiaGIVhGHz9Ov5nxBqOEydOwNfXF/7+/mxH0TujLgByuRyfPn0Cj8dDQkIC/P39ERkZierqaraj6UVHRwdEIhH27t1rdL2W5s+fDwcHB81ySUkJGhoadLrPGoqhWqobMhsbGyxZskRn7OLFi3BxcdH5/hmqx48fo6ioCCkpKWxH+SkMfporhUIx4HXh7gloDh8+jOTkZCQlJeHGjRvYunUr7ty506eZ3Xgz2LkJBAJcvXoVQqEQoaGhBveb5FDnZmNjo1mur69HSkoK1qxZA29v718VUW/a2tpgYWGhM9bdG6t391xDJ5VKce/ePUilUrajjJlSqYRIJIJIJNL5eTQmBl8AqqurERkZ2e+6jRs3AgDCwsI0fz948CDKyspw8+bNfuc0Hk8GO7cVK1agoqIC+fn5vziVfgx2bklJSUhOTgYAPHnyBLGxsXB1dUV6evqvjKg3ptImXSKRICsrC2KxGIsXL2Y7zphJJBI4Oztj1apVbEf5aQy+APj5+aGurq7fdSqVCnl5eZg5c6ZmjMPhwMXFBW/fvv1VEUdtsHPLzs7G/fv3NXMtMz86eqxbtw7x8fGIj4//ZTlHY7Bz6/bw4UPExcXBy8sLp06dMtg3S+2W6lZWXTNVDdRS3VAdPXoUly9fRlpaGiIiItiOoxeFhYVoaWmBj48PgK5Lrp2dnfDx8UFRURGmTh35zF3jjcEXgMHw+XzMmTMHNTU1WLt2LYCux9ZevHgBX19fltONzZYtWxASEqJZbmpqQlRUFLKzs+Hp6cliMv149uwZ4uLisGDBAmRlZRl0O3HtlurdDx8M1FLdEGVnZ+PKlSvIyMhAWFgY23H0Jjc3FyqVSmdZJpPh5MmTfaa+NVRGXQAAICYmBjt37oSHhwf8/Pxw5coVtLa29vuEiSGxsbHRuS7J43VNLO3o6GgU1ytFIhGEQiHEYjE+f/6sGZ8wYYLOzVRDoN1SPSMjAwqFApmZmdi2bRvb0casrq4OEokEMTExCAgI0JkQys7OTvNzaYh6F2eBQABzc3OjmtrW6AvA0qVLkZ6eDolEgrS0NHh4eOD8+fNGU8GN0Zs3bzRPagUFBemsi42Nxe7du9mINSZ79uyBQqFATEyMpqV6dHQ027HGrLi4GGq1GlKptM+N39u3b+tcfiXjD7WDJoQQE2XUnwMghBAyMCoAhBBioqgAEEKIiaICQAghJooKACGEmCgqAIQQYqKoABBCiImiAkAIISaKCgAhhJio/wNkQdgLwJOeEQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# NOTE: we are storing y's as 1's or 0's but still doing the mean stuff correctly\n",
    "def get_sample(n, mean, var):\n",
    "    x_sample = []\n",
    "    y_sample = []\n",
    "    for i in range(n):\n",
    "        y = np.random.binomial(1, .5)\n",
    "        if (y == 0):\n",
    "            x = np.random.normal(-mean, var)\n",
    "        else:\n",
    "            x = np.random.normal(mean, var)\n",
    "        x_sample.append(x)\n",
    "        y_sample.append(y)\n",
    "    return x_sample, y_sample\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set(font_scale=1.3)\n",
    "\n",
    "def split_by_class(sample):\n",
    "    y_one = [pair for pair in sample if pair[1] == 1]\n",
    "    y_neg_one = [pair for pair in sample if pair[1] == 0]\n",
    "    return y_one, y_neg_one\n",
    "\n",
    "X, y = get_sample(100, 2, 1)\n",
    "sample = list(zip(X, y))\n",
    "class_one, class_neg_one = split_by_class(sample)\n",
    "sns.distplot(list(zip(*class_one))[0], hist=False, rug=True)\n",
    "sns.distplot(list(zip(*class_neg_one))[0], hist=False, rug=True)\n",
    "plt.title(\"Distribution Plot with Mean = 10\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithm Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3405430602325196"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.ensemble.forest import _generate_unsampled_indices\n",
    "from sklearn.ensemble import BaggingClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "X, y = get_sample(1000, 1, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "\n",
    "def rf_entropy_estimate_bagging(X, y, n_estimators = 200, max_samples = .63, bootstrap = True, depth = 4):\n",
    "    model = BaggingClassifier(DecisionTreeClassifier(max_depth = depth), \n",
    "                              n_estimators = n_estimators, \n",
    "                              max_samples= max_samples, \n",
    "                              bootstrap = bootstrap)\n",
    "    model.fit(X, y)\n",
    "    class_counts = np.zeros((X.shape[0], model.n_classes_))\n",
    "    for tree in model:\n",
    "        # get out of bag indicies\n",
    "        unsampled_indices = _generate_unsampled_indices(tree.random_state, len(X))\n",
    "        # get all node counts\n",
    "        node_counts = tree.tree_.n_node_samples\n",
    "        # get counts for out of bag samples\n",
    "        partition_counts = np.asarray([node_counts[x] for x in tree.apply(X[unsampled_indices])])\n",
    "        # get probability for out of bag samples\n",
    "        class_probs = tree.predict_proba(X[unsampled_indices])\n",
    "        # find total elements for out of bag samples\n",
    "        elems = np.multiply(class_probs, partition_counts[:, np.newaxis])\n",
    "        # store counts for each x (repeat fhis for each tree)\n",
    "        class_counts[unsampled_indices] += elems\n",
    "    # calculate p(y|X = x) for all x's\n",
    "    probs = class_counts/class_counts.sum(axis = 1, keepdims = True)\n",
    "    entropies = -np.sum(np.log(probs)*probs, axis = 1)\n",
    "    # convert nan to 0\n",
    "    entropies = np.nan_to_num(entropies)\n",
    "    return np.mean(entropies)\n",
    "\n",
    "np.warnings.filterwarnings('ignore')\n",
    "\n",
    "rf_entropy_estimate_bagging(X, y, 500, .63)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Mutual Information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import entropy\n",
    "\n",
    "def get_entropy(X):\n",
    "    probs = np.bincount(X)/len(X)\n",
    "    return entropy(probs)\n",
    "\n",
    "def get_mutual_information(X, y, n_estimators = 200, max_samples = .63, bootstrap = True, depth = 4):\n",
    "    cond_entropy = rf_entropy_estimate_bagging(X, y, n_estimators, max_samples, bootstrap, depth)\n",
    "    h_y = get_entropy(y)\n",
    "    return h_y - cond_entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean: 0\n",
      "MI 0.0020832183025618622\n",
      "Mean: 1\n",
      "MI 0.3798822472966329\n",
      "Mean: 5\n",
      "MI 0.6929851718108113\n"
     ]
    }
   ],
   "source": [
    "### TESTS ###\n",
    "print(\"Mean:\", 0)\n",
    "X, y = get_sample(1000, 0, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "print(\"MI\", get_mutual_information(X, y))\n",
    "\n",
    "print(\"Mean:\", 1)\n",
    "X, y = get_sample(1000, 1, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "print(\"MI\", get_mutual_information(X, y))\n",
    "\n",
    "print(\"Mean:\", 5)\n",
    "X, y = get_sample(1000, 5, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "print(\"MI\", get_mutual_information(X, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Permutation Test\n",
    "For now, we will just return a list of mutual information values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import shuffle\n",
    "from tqdm import tqdm_notebook\n",
    "def permutation_test(X, y, trials = 1000):\n",
    "    mutual_infos = []\n",
    "    for _ in tqdm_notebook(range(trials)):\n",
    "        shuffle(y)\n",
    "        mutual_infos.append(get_mutual_information(X, y))\n",
    "    return mutual_infos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "94df825cdf284169b1009af48e743355",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X, y = get_sample(1000, 1, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "mutual_informations = permutation_test(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAERCAYAAABowZDXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8U1X6+PFP0i3dW2jasgsUWqFASitFoCAMiMrMd3QGZ0CGwfmBiiiKwsimOIAoMriwuKGIIIzigsiIKLjhCIhSLCiLUCwiW/fSpEua5f7+CA0NdElptsLzfr2q5Jx7c889TfPce5Z7VIqiKAghhBANUHu7AEIIIZoHCRhCCCGcIgFDCCGEUyRgCCGEcIoEDCGEEE6RgNHMNPdBbc29/L5E6lJ42lUZMMaOHUtiYqL9Jykpid69ezNq1Ci+/vprbxevTnv27OGf//xnk99nxowZ/OUvf6k3PzExkbFjx9aan5ubS1JSEomJiRiNRqePu23bNp588slGl7c+J0+eJDExsd7f29ixY3nooYca9b5ZWVncdNNN9OjRgwULFjS1mC53cV029Dv1NRf/3s6cOcOdd95p/zzt3r2bxMREjh071uj3TExMZPfu3bVuM2fOHBITE1m8eHHTT8JLrFYrzz33HBkZGaSkpHD//feTm5vrkWP7e+QoPqhfv348+OCDgO1KzWAw8OabbzJx4kTeffddunfv7uUSXuq9997jxIkTHjmWSqUiMzOT4uJioqOjHfK2bdt2WVe3a9asISYmxlVFdKvly5ejKAorVqygTZs23i7OJS6uy0mTJjUqeHtbbGws69evp3PnzgDs2rWLXbt2ueS9VSoVn332Genp6Q7pVquVzz77zCXH8KYlS5awZs0apk+fTnR0NM899xwTJ07k/fffR6127z3AVXmHARAVFYVOp0On05GSkkJGRgZLly4lLCyM9evXe7t4Xte1a1eCgoL48ssvL8nbunUrXbt29UKpPEev16PT6bj++utp3769t4vToPbt29OlSxdvF8NpgYGB6HQ6wsPDXf7eOp2Ozz///JL0zMxMLBYLcXFxLj+mpxgMBtasWcO0adMYNWoUw4cPZ8WKFRw+fNgjrSNXbcCojUaj4ZprruH06dP2tB9++IHRo0fTs2dPMjIyWLZsGVar1Z4/ZMgQnnvuOW677TZ69erFhx9+yIwZM5g2bRovvvgi/fv3JyUlhdmzZ2MwGJg7dy6pqakMHDiQt99+2/4+tTUpvPXWWyQmJtrzP/jgA/bt20diYiInT54E4NixY0yYMAGdTkd6ejrz58+nsrLS/h6KovDKK6/Yb18XLFiAxWJxqi4yMjIuuSIrKioiMzOTG2+80SG9tmafxYsXM2TIEHv+d999x8cff+xwTvWdM0BVVRXPPvssQ4cOJTk5mb59+zJjxgwMBkOD51CXIUOG8MYbbzBz5kxSU1Pp06cP8+fPx2QyAZCYmEhWVhYbN250qOuPP/7Y/nv+3e9+x6uvvupwp5WYmMiKFSsYPnw4KSkp7N69m7Fjx7Jo0SIWLFjAddddR58+fXj22WcpLS3l4YcfRqfTMXToUId6VhSF1157jVtuuYXk5GTS0tIcmh2cqUuDwcCCBQsYPHgwPXv25I477uCHH36w52/YsIGBAweyfft2fv/735OcnMytt97Knj177NuUlZUxc+ZM+vfvT8+ePRk1alSddwEHDx4kMTGRffv22dOefPJJEhMTycvLs6fdf//9TJs2zaFJasOGDcycOROAnj17smHDBvv2+/fv5/bbb6dHjx7cfPPNbNu2rcHf77Bhwzh16hSHDh1ySP/0008ZMmQI/v6ODSt6vZ7HHnuM9PR0dDodEydO5NSpUw7bfPnll4waNQqdTmevi8zMTHv+2LFjWbx4MU8//TR9+/YlJSWFqVOn1vs5vbhpvObPjBkzat1n3759lJeX2/+uwHaxkJCQwDfffNNg3TSVBIwazGYzp06dsjdBHD58mHHjxhEVFcWyZcu46667WLlyJf/+978d9lu5ciV//OMfWbx4MX379gXgiy++YNeuXSxcuJCJEyfy3nvv8ac//YnKykqWLVuGTqdj/vz59i+jhkyaNIlBgwbRpUsX1q9fT2xsLLm5uYwZM4aysjKeffZZHnnkETZv3swjjzxi3+/VV19l2bJljB07lueff55jx47x8ccfO3XMoUOHsmPHDioqKuxpn3/+OYmJiY1upnn88cfp1q0b/fr1a9Qd3IIFC3jvvfe47777WLlyJRMmTOCjjz7i1VdfbdTxL7Zs2TJUKhXLli1j/PjxrF27lnfffReA9evX06VLFwYNGmSv67Vr1/Lwww/Tp08fXnjhBW677TaWLFlyyWfhpZdeYuLEicydO5eePXsC8Pbbb5Obm8vSpUv54x//yCuvvMLtt99O69ateemll4iLi2PGjBmUl5cDtt/Z8uXLGTNmDK+//joPP/ww3377LYsWLXKqLi0WC//4xz/YsmUL9913H0uXLiUkJIS///3vHDhwwL5dSUkJ8+bNY/z48bz44ouoVCoeeughzGazve6///575syZwyuvvEJ0dDT33nsvRUVFlxzz2muvpWXLlnz33Xf2tO+//x6AvXv32su1e/duMjIyHPa94YYbuPfeewFYu3YtN9xwgz3vySef5NZbb+Xll18mJiaGhx9+mPz8/Hp/t126dOGaa665JAh/9tlnl1zoWK1W7rnnHrZv387MmTNZvHgx+fn5/P3vf6esrAywXTROmjQJnU7Hyy+/zNNPP41er2fatGkOF19vvfUWOTk5LFq0iEceeYStW7fy0ksv1VnOxx9/nPXr19f6M2nSpFr3OX78OIGBgcTHxzukt23b1iPN1VdtH4aiKPY/DKvVytmzZ3n55ZcpLCxk5MiRgO2Pv127dixfvhw/Pz8AgoODmTt3LhMmTKBly5YAJCcnc+eddzq8f1VVFcuWLSMqKoqMjAzWr1+PSqViwYIFqNVqOnfuzKeffsqhQ4do27Ztg+Vt3749LVq0oKSkBJ1OB8Dq1atRqVS8+uqrhIWFAdC6dWvuvPNODh06RGJiIq+//jpjx47l7rvvBiA9PZ3Bgwc7VUeDBw/GYrGwY8cOhg4dCtiao4YPH+7U/jUlJCQQFhZmbwp0VnFxMTNnzuQPf/iDvfyZmZkOV3eXo3PnzvZO4379+vHll1/y9ddfc8cdd6DT6QgJCaFFixbodDosFgvLli1j5MiR9ivhAQMGoFKpeOmll5gwYQItWrQAbHcvt912m8OxQkJCeOaZZwgICCAtLY233nqLVq1aMW3aNMDWPHPHHXdw/PhxunXrRl5eHg888ABjxowBoE+fPvzyyy/2ZpaG6vLLL79k//79/Oc//yE1NRWAjIwMbrnlFl544QVefPFFAIxGI48//jgDBw4EQK1WM378eA4fPkxycjJ79+6lX79+9t93t27dePnllx0uIKqpVCr69+/P999/z1133YXBYODw4cMkJSWxd+9ebrrpJn788Uf0ej0ZGRn24AjQokULe7Nfz549CQoKsuc9+OCD9nqIi4tjxIgR7Nu3z/55rEv1XdvkyZMB+PHHHyktLaV///4O2/3vf/8jMzOT9evX2+uyT58+3HDDDbz77rvceeedHDt2jN///vcOV/3+/v7cf//9nD59mnbt2gEQGhrK8uXL7Xcwe/bs4euvv65zoEpCQkK951Abg8FASEgIKpXKIT00NJQzZ840+v0a66q9w9iyZQvdu3ene/fu9OjRg2HDhvHll18yb948evToAdiukPr3728PLmazmYyMDEwmk/2qCaBjx46XvH+HDh2Iioqyv46JiSEpKcneKVWdp9frL/scvv/+e1JTU9FoNPbyVb/evXs3OTk5FBcX278Q4EJTkzPCw8NJT0+3X6np9Xp27drFsGHDLrvMjbV06VL+8Ic/cObMGXbs2MEbb7zBsWPH7M1Hl6v66r9aXFycw5dYTb/88gslJSXcdNNNDum33HILJpPJoRmmts9Ct27dCAgIACAgIIDIyEi6detmz7/4s/Doo4/y//7f/6OgoIDdu3ezbt069u7d6/Q5Z2ZmotVq7cECwM/Pj+HDh18SaHv16mX/d/VVa3VASEtL45133mHSpEm89957WCwWpk+fXufdZf/+/e39BJmZmcTHx3PjjTfa/1Z27txJt27d7MHVGSkpKfZ/Vx/XmebIYcOGcfjwYfsd/NatW7nhhhsIDAx02O67774jKiqK5ORk+99QSEgIOp3OPtJq5MiR/Pvf/8ZgMLB//342btzIpk2bABx+J927d3do7qrvMwW2O67qY178U7PZuyar1XpJsKhWV7orXbV3GAMGDGDKlCmA7coqIiKCtm3bOlR6SUkJq1evZvXq1ZfsX7NdtvpOo6bQ0NBL0oKDg11RdIfy7d+/v9YRXXl5eZw7dw7gklFOLVu25Pjx404dY9iwYTz33HNYLBa+/PJLrrnmGjp16kRWVlaTy++MPXv28Pjjj5OdnW3/w9ZoNE2eg6DRaBxeq9XqOt+zuh4vHuFV/Xuv+QXmis9CdnY2s2fPJisri9DQULp3705QUJDT51xaWlrraLQWLVrYm1mq1ayH6s9+9ZfVo48+SkxMDB9++CGff/45/v7+3HzzzcyfP7/W8g8YMICysjIOHjzId999R2pqKqmpqbz44otUVFSwc+dOpy9Waitf9cVWXV+mNfXq1YvY2Fg+++wz7rzzTrZt21br0OqSkhJKSkpq/RuqvnAsKyvj0Ucf5ZNPPsHPz4+EhAR7q0DN30ljPlMAd955p0MTXk233XYbCxcuvCQ9PDy81iBUVlZmb2Vwp6s2YERERNg/EHUJDw/n97//Pbfeeuslea1bt3ZpeVQq1SWd0fVdnVSXb/jw4dx1112X5MXExNj3v7jNufoL0BlDhw5l7ty5ZGZm1toGXFNjy9/QOev1eu69914GDBjgMLx1ypQpDgMT3C0yMhKAgoICh/Tq19X5rmC1Wrn33nuJj49ny5YtdOzYEZVKxb///W+n26gjIiIuKStAYWFho8qq0WiYMmUKU6ZM4ciRI2zevJkVK1bQpUsX7rnnnku2r76L/u6778jMzOSPf/wjOp0OlUrFt99+S1ZWlv0izd1UKhVDhw7l888/p1+/fpw9e9bhTrtaREQEbdq0YcmSJZfkVQfFJ554gszMTNasWYNOpyMgIIDt27c71QFfn7lz514SwKtdfJFXrUOHDhiNRvLy8oiNjbWnnzx5kgEDBjSpPM64apuknJGSksLx48fp0aOH/cff35/nn3+ewsJClx4rJCSE3NxchyuSi5sPLh5j3bt3b3755Re6detmL19cXBzPPvssx48fp1OnTmi1WrZu3Wrfx2w28+233zpdrpiYGHQ6HVu2bOF///tfnQEjNDT0kslDNZvtait/Q+f8yy+/UFpayj/+8Q97sKisrGTv3r1OXWW6SqdOnYiKiuKTTz5xSN+yZQt+fn6XNG81RVFRESdOnGD06NF06tQJlUqF1Wpl165dDudc33j71NRU8vPzHerSYrGwdetWp/uPFEXh1ltv5Y033gBsw6wfeughrrnmGs6ePVvnfv379+ebb77hwIEDpKWlodFoSE5O5uWXX0aj0dR5fHfMHxg2bBh79+5lw4YNZGRkEBIScsk2KSkp5OXl0bJlS/vfUHJyMqtWrWLHjh2AbRLnkCFDuO666+xNizt37gSaNtu+U6dODt8tNX/q6tdMSUkhMDCQL774wp524sQJsrOzL5l34g5X7R2GMyZOnMiYMWOYOXMmt9xyC+fOneO5554jODi41rbqphgwYABr165l4cKFDB48mK+++uqSgBEREcGJEyfYtWsXKSkpjBs3jg8++ID77ruPUaNGYTabWb58OQUFBSQlJaFSqbjvvvuYN28e0dHRpKSk8Pbbb1NQUNCoCXQ33ngjzzzzDK1btyYpKanO8j/xxBOsWLGCHj168MEHH3DmzBmH5piIiAh+/vlndu/eTZ8+fRo8544dOxISEsKSJUsYP348er2e119/nby8PLeM36+Ln58fkyZN4qmnniI0NJSBAweSlZXFSy+9xNixYx36qpqqZcuWtGrVipUrVxIaGorVauWtt97i4MGDDk0eF9dlTTfccAPJyck8+OCDPPTQQ8TExLBu3TpOnTp1yaiuuqhUKnr27MmLL75ISEgIHTp0YNeuXeTk5PDYY4/Vud+AAQN47bXXiIqKsnfqpqam8tprrzF8+PBLhrTWPB+wBeF+/fo5VcaG9OnTh9DQUNauXVvnEwYGDx5Mly5dmDBhApMmTaJFixa88847bNu2jdGjRwO2QS2ffPIJqampxMTE8MUXX7Bu3Tqg4btoVwsLC+OOO+7gqaeewmKxEBsbyzPPPENSUlKtd1CuJncY9dDpdKxcuZKcnBzuu+8+FixYQO/evXn99dftVxquMnjwYB544AE2b97MxIkTOXPmDI8//rjDNn/5y18ICwvj7rvv5uDBg7Rr1461a9dSVVXFlClTmDVrFq1atWLNmjX2W9rRo0cze/ZsPvzwQyZPnkxISEijHyExdOhQTCZTvZ3df/3rX/nb3/7Ga6+9xv33349Go+H+++932GbcuHGUlpZy9913k5ub2+A5R0REsGTJEvLy8pg4cSILFy4kMTGRxx57jOPHj1NSUtKo82iKcePG8a9//YuvvvqKe+65h40bN/LQQw/VOV7+cqlUKpYuXYpareaBBx5gzpw5hIWF8eyzz1JRUcHPP/9sL0/NuqzJ39+flStXMnDgQBYtWsSDDz5IRUUFa9asadTdUPXotOXLlzN+/Hi2bdvGwoUL6/1CT01NJSQkhJSUFHufyHXXXQdQb//F9ddfT9++fXn00Uf58MMPnS5jffz9/e0jAusaGRgQEMDKlSvp1asX8+fPZ9KkSZw+fZqXXnrJXu4ZM2aQlpbGv/71L6ZMmcLhw4d54403CA4Odhjw4CnTpk1j9OjRLFu2jBkzZpCQkMDLL7/s9lneACpZolUIIYQz5A5DCCGEUyRgCCGEcIoEDCGEEE6RgCGEEMIpEjCEEEI4RQKGEEIIp1wxE/eKi8uwWn1rhHDLlmEUFl7+ug1XEqmLC6QuHEl9XOCpulCrVURHX/qMs4ZcMQHDalV8LmAAPlkmb5G6uEDqwpHUxwW+XBfSJCWEEMIpEjCEEEI4RQKGEEIIp0jAEEII4RQJGEIIIZwiAUMIIYRTJGAIIYRwyhUzD6M5MlvBaDI7pAUF+OMvYVwI4YMkYHiR0WTm+0OOq6Vdd20c/kHyaxFC+B65lhVCCOGURgUMo9HIiBEj+Prrr+1p2dnZjB8/nrS0NAYNGsTTTz+N0Wi05589e5a7776blJQUhg4dyqZNmxzes6F8IYQQvsHpto+KigoeeughsrOz7WllZWXcddddXHfddaxfv578/Hxmz56NxWJh1qxZANx33320atWKd999lz179jBr1izatm1L7969ncoXQgjhG5wKGD/99BPTp08nICDAIX3Xrl3o9XqeeOIJAgMD6dy5Mw8++CBPPvkks2bN4vvvv+fnn39m1apVREREkJCQwL59+3jzzTfp3bt3g/lCCCF8h1NNUrt27WLQoEG8/fbbDuk9evTghRdeIDAw0J6mUqkwGAwoisIPP/xAYmIiERER9vy0tDSysrIAGswXQgjhO5y6w7jrrrtqTY+LiyMuLs7+2mKx8Oabb5Keno5KpSI3N5fY2FiHfWJiYsjNtY0MaihfCCGE73Dp+M358+dz+PBh3nnnHcDW7xEUFOSwTWBgIBaLBbPZ3GC+v7/zxWvZMqzpJ+AGWm14nXlKUTnhYRqHtJCQILQtQtxdLK+ory6uNlIXjqQ+LvDlunBJwLBYLMydO5f333+fJUuWkJSUBIBGo6GkpMRh26qqKgICAvD3928wvzEKCw0+t/CIVhtOfr6+zvxyoxm9odIxrdxIvsXi7qJ5XEN1cTWRunAk9XGBp+pCrVZd1kV2k+dhmEwmHn74YT744AOWLl3K0KFD7Xnx8fHk5+c7bJ+fn29vhmooXwghhO9ocsCYM2cO27dv55VXXuF3v/udQ55Op+Pnn3/GYLiwRm1mZiYpKSlO5QshhPAdTQoY27dvZ8OGDfzzn/+kS5cu5Ofn23/ANuKpU6dOTJs2jSNHjrB+/Xo++ugjxo4d61S+EEII39GkgLFlyxYA5s2bx4ABAxx+jEYjarWa5cuXYzQaGTlyJK+99hpPPfUUOp3OdvAG8oUQQvgOlaIovtVTfJmaY6d3mbH2hw+GXoEPH5SOzQukLhxJfVxwxXd6CyGEuDpIwBBCCOGUK6/tw0fVtliSj7WgCSFEvSRgeEhtiyX16qr1UmmEEKLxpElKCCGEUyRgCCGEcIoEDCGEEE6RgCGEEMIpEjCEEEI4RQKGEEIIp0jAEEII4RQJGEIIIZwiAUMIIYRTJGAIIYRwigQMIYQQTpGAIYQQwikSMIQQQjhFAoYQQginSMAQQgjhFAkYQgghnCIBQwghhFMaFTCMRiMjRozg66+/tqfp9XqmTp1KamoqGRkZrFq1ymGfpuYLIYTwDU4v0VpRUcFDDz1Edna2Q/rs2bPJy8tj3bp1HD9+nJkzZxIbG8uIESNcki+EEMI3OBUwfvrpJ6ZPn05AQIBD+qlTp9i6dSsfffQRCQkJJCUlkZ2dzerVqxkxYkST84UQQvgOp5qkdu3axaBBg3j77bcd0rOysoiKiiIhIcGelpaWxoEDBzCZTE3OF0II4TucusO46667ak3Pzc0lNjbWIU2r1WI2mykoKGhyfqtWrRpzLkIIIdzI6T6M2lRUVBAUFOSQFhgYCEBVVVWT8xujZcuwRm3vKVptOABKUTnhYRqHvIAA/0vSQkKC0LYI8Vj5PKm6LoTUxcWkPi7w5bpoUsDQaDSXfLFXv9ZoNE3Ob4zCQgNWq9KofdxNqw0nP18PQLnRjN5Q6ZBvMl2aVl5uJN9i8VgZPaVmXVztpC4cSX1c4Km6UKtVl3WR3aR5GPHx8eTn5zuk5eXlERAQQHR0dJPzhRBC+I4mBQydTkdhYSE5OTn2tMzMTJKTkwkMDGxyvhBCCN/RpIDRpk0bBg8ezPTp0zl48CCffvopK1euZNy4cS7JF0II4Tua1IcBsHDhQubMmcPo0aOJjIxkypQp3HzzzS7LF0II4RtUiqL4Vk/xZfL1Tu8yo5nvD+U65PfqqmXfEcc+nOuujSM0qMlx3OdIx+YFUheOpD4uuKI7vYUQQlw9JGAIIYRwigQMIYQQTpGA4WWKolBYWskV0pUkhLiCScDwsp9/K2Hzzl859Guxt4sihBD1koDhRcfPlLLnUB4AR387J3cZQgifJgHDSyqMZlZ/fJgQTQCpiVrOlVWRX1LZ8I5CCOElEjC8wGpV+DrrNBVGMzektKZruyj8/VRknzzn7aIJIUSdJGB4QV5xBbnFFfxxYCdaRGgI8FdzTXwEx8+WUlll9nbxhBCiVhIwvKBIb2t6Su7Uwp6W0DYSs0Xhh4tmfgshhK+QgOEFxXojmkA/wkMuPJFXG6UhMjSQnT+d9WLJhBCibhIwvKBEbyQ63HGlQZVKRULbSI6f0XOqoMxLJRNCiLpJwPAwq1WhxFB1ScAA6NQ6AoC90iwlhPBBEjA8TF9ehcWq1BowgoP80UZp+PWsPLlTCOF7JGB4WLHeCEBULQEDoF1sOL+eLfVkkYQQwikSMDys2FCFSgVRobUvQdsuLozCUiP68ioPl0wIIeonAcPDivVGIkIC8fOrverbxdoWNfk1V5qlhBC+RQKGh9U2Qqome8CQfgwhhI+RgOFBVSYLhgpTvQEjRBOANkrDcQkYQggfIwHDg0oMtg7v+gIGQIf4CLnDEEL4HAkYHtTQCKlq18SHU3CuEkOFyRPFEkIIp7gkYJSWljJjxgzS09Pp168fc+bMoazMNlvZZDIxb9480tPTSU9PZ/HixVitVvu+DeVfSYr1RgL81YRq/OvdrkN8OCAd30II3+KSgDF37lyOHz/O6tWrefHFF/nuu+94+umnAXj22WfZsWMHK1as4LnnnmPjxo289tpr9n0byr+SFJ/v8FapVPVu1yHOFjBOSLOUEMKHuCRgfPXVV9x5550kJSWh0+kYM2YM3377LUajkbfeeotZs2bRq1cv+vXrx9SpU1mzZg2KojSYfyVRFIUSfe2PBLlYWHAAMZHS8S2E8C0uCRjR0dFs2rQJg8FASUkJn376Kd27d+fQoUNUVFSQlpZm3zYtLY38/HxOnjzZYP6VpKjUiMliJTqs4YABtmYp6fgWQvgSlwSMefPmsW/fPq677jr69u1LaWkp8+bNIzc3l7CwMEJDQ+3barVaAHJzcxvMv5KcKbT16USF1z7D+2LXxIeTV1JBeaV0fAshfEP9va9OysnJoXPnzjz//POYzWaeeuopZsyYwbBhwwgKcryiDgy0fWFWVVVRUVFRb35jtGwZ1oQzcB+t1tYfUZZ1GoC4mHDCggMACAjwJzxM47B9SEgQ2hYh9EyM4/3tv3DOaKFDuxZcCarrQkhdXEzq4wJfrosmB4wTJ07wxBNPsHXrVtq1awfAkiVLuOmmm+jdu/clX/zVrzUaDRqNpt78xigsNGC1+la/h1YbTn6+rVkpt7AMlQosZjN6gwUAk8mM3lDpsE95uZF8i4XoYNuvZt/hPFpFNq4ufFHNurjaSV04kvq4wFN1oVarLusiu8kB46effiI0NNQeLAA6duxIaGgoFRUV6PV6KioqCA4OBiA/37bWQ1xcHGq1ut78K0mJwUhIkD/qBkZIVQsLDiA6PIiT+QY3l0wIIZzT5D6MuLg4DAYDZ89eWFr0zJkzlJWVcf311xMcHExmZqY9b8+ePcTGxtKmTRuSkpLqzb+SlOiNhDQw/+JirWNCZfU9IYTPaHLA6NWrF9deey1Tp07lwIEDHDhwgKlTp5Kenk5qaiojR45k3rx57N27l127dvHMM88wbtw4wNbsVF/+laTEYCRUE9CofdrEhHKmoAzrFTbEWAjRPDW5Scrf358VK1awcOFCJkyYgEqlYtCgQcyYMQOAf/7znxiNRiZMmEBB1UwOAAAgAElEQVRQUBAjR45k/Pjx9v0byr8SVM/B0EYFN7itSq2izGgGICZKQ5XZyok8A2214fjLg1yEEF6kUq6QGXK+3OldWl7FlKXfcF1SLNdeE23P79VVy76L1u+umZZfUsGWb08wuHcb/jSoM6FBLhnU5hXSsXmB1IUjqY8LfL3TW65ZPaC41PbQwcb2YUSG2YYYl5x/aKEQQniTBAwPKNLbhs429NDBiwX6+xGq8bc/Fl0IIbxJAoYHVD/WPKSRnd4AUWFBlBhkfW8hhPdJwPCAolIjarWK4CC/Ru8bGRbIubIqn+ufEUJcfSRgeECRvpKosMAGH2tem6iwIKxWhYJzFW4omRBCOE8ChgcUlxqJcvIptRerfljhmcJyVxZJCCEaTQKGBxTpKxtclrUukaG2/c7IjG8hhJdJwHAzRVFsK+1d5h1GgL+asOAAucMQQnidBAw305ebMFsUp1baq0tUWKB9PQ0hhPAWCRhuVj0H43KbpAAiw4LIK67AbLG6qlhCCNFoEjDcrHqW9+U2SYHtDsNiVcgtlpFSQgjvkYDhZkXnJ+015Q6jeoTVaen4FkJ4kQQMNyvSV+KnVhEW0vhZ3tUiwwJRAadkMSUhhBdJwHCz4lIj0eFBTq+0Vxt/PzUxURpO5csdhhDCeyRguFlRaSUtmtAcVa11TJgs1yqE8CoJGG5WpDfSIkLT5PdpHRNCXnEFRpPFBaUSQojGk4DhRlbr+Ul7Ea64wwhFQTq+hRDeIwHDjc4ZjFisCi3CXXGHYVsd62SeNEsJIbxDAoYbVT9h1hV9GDGRGgL91ZyUjm8hhJdIwHCjgpLzAcMFfRhqtYrWMaHS8S2E8BoJGG5UdO78Y0HOr83dVG21YTIXQwjhNS4JGGazmX//+9/069ePtLQ0pk6dSmlpKQB6vZ6pU6eSmppKRkYGq1atcti3ofzmrEhvRKWC8BAXBYzYMErLTZwrkyVbhRCe55KAsXjxYjZt2sTzzz/PG2+8wZEjR1iwYAEAs2fP5tSpU6xbt47Zs2ezdOlSNm/ebN+3ofzmrLi0kojQQNTqy5+0V1NbbSiANEsJIbzCv6lvoNfrWbt2LUuXLqVPnz4ATJs2jWeeeYZTp06xdetWPvroIxISEkhKSiI7O5vVq1czYsSIBvObu6LSSqJCm97hXa2t1jZS6lSege7XtHDZ+wohhDOafIexZ88eAgICGDhwoD1t0KBBbNq0iaysLKKiokhISLDnpaWlceDAAUwmU4P5zV1xqZFIF/VfAESEBhIRGigjpYQQXtHkgHHixAnatGnDZ599xv/93/8xcOBA5syZg8FgIDc3l9jYWIfttVotZrOZgoKCBvObuyJ9pcs6vKu11YbymzRJCSG8oMkBo6ysjNOnT/Pqq68ye/ZsFi1axN69e5k5cyYVFRUEBTk2yQQG2r5Aq6qqGsxvzixWK+cMRvua3K7SVhvG6YIyrFbFpe8rhBANaXIfhr+/P2VlZSxatIjOnTsDMH/+fEaNGsW11157yRd/9WuNRoNGo6k3vzFatgy73FNwi8JzFSgKtG0VgVYbjlJUTniY4zkFBPg7lQYQEhKEtkUI13aKYev3v2FWq2mj9a1zbohWG+7tIvgMqQtHUh8X+HJdNDlgxMbGolar6dSpkz2t+t9Wq5X8/HyH7fPy8ggICCA6Opr4+Ph68xujsNDgU1fdv57VA+CnKOTn6yk3mtEbKh22MZmcSwMoLzeSb7EQGewHwP7DuQTiO+fbEK02nPx8vbeL4ROkLhxJfVzgqbpQq1WXdZHd5CaplJQUrFYrBw8etKcdPXoUtVrNbbfdRmFhITk5Ofa8zMxMkpOTCQwMRKfT1ZvfnJUYbCvtubLTG6B1y1BUKhlaK4TwvCYHjA4dOnDjjTcya9Ys9u/fz/79+5k7dy7Dhw+nTZs2DB48mOnTp3Pw4EE+/fRTVq5cybhx4wAazG/OqifXuWpYrUqtosxoxmRV0EYFk3NWj9nqkrcWQginNLlJCuDpp59m4cKFjB8/HkVRGD58OLNmzQJg4cKFzJkzh9GjRxMZGcmUKVO4+eab7fs2lN9cVd9hRIS65g7DaLKw74it+S5U48+xk+cwmsz4B7nkVyiEEA1yybdNSEgI8+bNY968eZfkRUVFsXTp0jr3bSi/uTpnqCI8JIAAf9c/rqtlpIacM3rOGYyESsAQQniIPHzQTUoMRqJd8JTa2sRE2t73RK70YwghPEcChpucK6tyycJJtYkO16ACTuTKyBIhhOdIwHCTcwbXLM1amwB/NZFhgRIwhBAeJQHDDRRFsd1huKlJCmz9GCdyDShK85mLIYRo3iRguEFZpRmzRXFbHwbY+jEMFSYKz106yU8IIdxBAoaLma1wtqgcAL/zcyfKjGZcPQm9ZWQwAMfPSrOUEMIzJGC4mNFk5rtDuQDkFZXz/aFcvj+Ui9nq2ll20eGB+KlV5Jwpden7CiFEXSRguEGF0QxAiCbAbcfwU6tpow2VgCGE8BgJGG5wIWC4d1Jd+7hwfs3VY5WObyGEB0jAcIMKowV/PxWBAX5uPU77uHAqjBZyz/eZCCGEO0nAcIMKo5lgDzyyo0O87bn5x89Ix7cQwv0kYLiBpwJGXIsQAv3V5JyVfgwhhPtJwHCDcg8FDD+1ig7x4fxyWgKGEML9JGC4QYXRTIiHniLbtV0Ux8/oqawye+R4QoirlwQMFzNWWTBbFIKD3NvhXS2xfRRWRSH75DmPHE8IcfWSgOFipeW2lfY80SQFkNAmEj+1ip9/K/HI8YQQVy8JGC52zuDZgKEJ9OeaVuEcPlHskeMJIa5eEjBcrLTMswEDILFdNMfP6DFWWTx2TCHE1UcChotdCBie6cMASGofhcWqkH1K+jGEEO4jAcPFSgxG1GoVQW6e5V1TQttI1CqVNEsJIdxKAoaLnSszEhLkj0qlcvuxVOcfn25RoH1cGId+Lcbs2ofiCiGEneca2q8SJfoqtz90sJrRZGHfkXwAQoMDOHS8CH25kegw9ywNK4S4urn0DmPOnDn85S9/sb/W6/VMnTqV1NRUMjIyWLVqlcP2DeU3RyUGo8cm7dUU3yIYq4I87lwI4TYu+2bbvXs377zzDj179rSnzZ49m7y8PNatW8fx48eZOXMmsbGxjBgxwqn85kZRFM4ZqoiNDvb4sbXRwahUcPTkOXp30Xr8+EKIK59LAkZFRQWPPfYYvXv3xmy2PaLi1KlTbN26lY8++oiEhASSkpLIzs5m9erVjBgxosH85qis0ozJYvVYk1RNgf5+tIzQcFQm8Akh3MQlTVLPP/88qampXH/99fa0rKwsoqKiSEhIsKelpaVx4MABTCZTg/nNUbHeCLh3pb36tIqxrcBnqGie9SeE8G1NDhhZWVls3ryZ6dOnO6Tn5uYSGxvrkKbVajGbzRQUFDSY3xzZA4YX+jAA2mlDURT48VihV44vhLiyNembraqqilmzZjFr1iyioqIc8ioqKggKchytExgYaN+vofzGatkyrNH7uJr5/Bd1bMtQwkNs5xIepgEgIMDf/u9qzqY5u21YaBCRoWc59FsJ/ze4S9NPyMW02nBvF8FnSF04kvq4wJfrokkB44UXXqBDhw7ccsstl+RpNJpLvvirX2s0mgbzG6uw0IDV6t21rU+cPocKsJot6A2VhIdp0BsqATCZzPZ/V3M2rTHbduvYgszDuZw5ew5/P9+ZZqPVhpOfLysDgtTFxaQ+LvBUXajVqsu6yG5SwPjvf/9Lfn4+KSkpAJhMJiwWCykpKTz++OPk5+c7bJ+Xl0dAQADR0dHEx8fXm98clRiMhIUEoFa7f9JeXXp0asmun85y5LcSul3TwmvlEEJceZoUMN588037qKjq13v27GHJkiX4+/tTWFhITk4OHTt2BCAzM5Pk5GQCAwPR6XT15jdHxfoqorw8aS6xfRQB/mqysgskYAghXKpJbRZt2rShQ4cO9p/IyEgCAwPp0KEDbdq0YfDgwUyfPp2DBw/y6aefsnLlSsaNG2fft7785qhYX+n1gBEY4Me1HaLJOlqAoni3iU4IcWVxayP3woULiY+PZ/To0SxYsIApU6Zw8803O53f3BTrjUSGef/uSJcQQ8G5Sk4XlHm7KEKIK4hLx39OnjyZyZMn219HRUWxdOnSOrdvKL85qTJZKKs0e/0OA6BXQgx8+jNZ2QW00Xp/9JgQ4srgO8Nomrlig20ORpQP3GFEhwfRIS6crOzmOZ9FCOGbJGC4SMn5SXuRPnCHAdA7UcuxU6UUlFR4uyhCiCuEBAwXqZ7l7QtNUgDXd48DYNeBs14uiRDiSiEBw0XsASPc+01SADGRwSS2i2LngVwZLSWEcAkJGC5SrDeiCfRDE+g7a1L1S44nt6icX2SNDCGEC0jAcJFig5HocN9ojqqWlhRLgL+aXT9Js5QQoukkYLhIid73AkZwkD8pXWLYfTAXs0UW+xZCNI0EDBcp0vvmWtr9kuMpqzTLI8+FEE0mAcMFrFbb0qxRPnaHAdC9YwsiQgLYKc1SQogmkoDhAqXlVVgVhRY+EDBUahVlRrP9p9JkpU+3eLKyCyg5P7lQCCEuh+8M6WnGLgyp9X7AMJos7Dvi+Nj4fj3i+XzPb3yx9yR/GtjZSyUTQjR3cofhAtUBw9c6vatpo4JJ6arly72nMJos3i6OEKKZkoDhAvaA4YOd3mBrphqoa01ZpZmvfjhlb64yy8ApIUQjSJOUCxSVVuKnVhEeGkhFle9dwRtNFgrPVdAyUsOWb38lMECNSqXiumvj8A+Sj4AQwjlyh+ECecUVaKOCUau8tzRrQ1QqFd2uiaa03MSpfFknQwjReBIwXCC3uIK46GBvF6NBHeLCCdH4c+B4kbeLIoRohiRgNJGiKOSVlBPXIsTbRWmQWq3i2g7R5BZVkFtU7u3iCCGaGQkYTVRiqKLKZCW2GdxhACS2jyIkyJ/Mn/PlKbZCiEaRgNFEecW2K/W4aN+/wwDw91PTq4ttze99siKfEKIRJGA0UW6xbUW75tCHUa1z6wgiwwL5747j8lBCIYTTJGA0UW5ROf5+KlpEaLxdFKep1Sp6d9WSV1zB//af8XZxhBDNhEsCxtmzZ3nggQdIT0+nf//+zJo1i9JS26I9er2eqVOnkpqaSkZGBqtWrXLYt6F8X5dbPaRW7btDamvTVhtK5zYRfPhNDhVGs7eLI4RoBpocMKxWK/fddx9lZWWsXr2al156icOHDzNz5kwAZs+ezalTp1i3bh2zZ89m6dKlbN682b5/Q/m+Lq+4vNn0X9SkUqm4bWBn9GVVvLf9mLeLI4RoBpo8zffnn3/mp59+4ptvvkGr1QK2IDBmzBhOnTrF1q1b+eijj0hISCApKYns7GxWr17NiBEjGsz3dVZFIa+4gm7XtPB2US5Lh/hwhqa1Y9ue3+jbLY4ubaO8XSQhhA9r8h1Gq1atePXVV+3BAmxXr4qisGfPHqKiokhISLDnpaWlceDAAUwmE1lZWfXm+7oSvZEqs7VZdXjXpFKrGJ7enhYRQbz+8SFKyqrk+VJCiDo1OWBERUUxcOBAh7Q33niDjh07UlhYSGxsrEOeVqvFbDZTUFBAbm5uvfm+Lu/8CKnYZtgkBbZnTO0/VkBKFy25RRWs+vgQRpP0ZwghaufyJ8+tWLGCbdu2sWLFCvbv309QkOMTXAMDAwGoqqqioqKi3vzGaNkyrAmlvjx7zy97em2CFu35md5KUTnhYRdGTFX/OyDA3yG9MWlN3b+h90wK03Ayv4wDvxSSf87INd2j6z3vy6XVhrvlfZsjqQtHUh8X+HJduDRgvPDCCyxdupQ5c+aQkZHBkSNHLvnir36t0WjQaDT15jdGYaEBq9WzM5ezTxTj76cCk5n8fD0A5UYzekMlYAsW1f82mS6kV3M2ran7O/OevRJa8luunpc/2E98VBCaQNdeS2i14fY6utpJXTiS+rjAU3WhVqsu6yLbZfMwnnzySZYtW8a//vUvxowZA0B8fDz5+Y6rv+Xl5REQEEB0dHSD+b4ur5kOqa2NJtCPAT1bkV9cwX8+O+rt4gghfJBLAsby5ctZu3YtCxcuZPTo0fZ0nU5HYWEhOTk59rTMzEySk5MJDAxsMN/X5TbTIbV1iW8Zwo192vHN/jN8dyjX28URQviYJgeMn3/+mRdeeIHx48fTv39/8vPz7T/x8fEMHjyY6dOnc/DgQT799FNWrlzJuHHjAGjTpk29+b6sekhtc3nooLNu6XcNHVtF8MaWw/yaq5eV+YQQdk1uqN66dStWq5UVK1awYsUKh7yPP/6YhQsXMmfOHEaPHk1kZCRTpkzh5ptvtm/TUL6vKtEbMTXjIbV1MVsVUrrEcKrAwNL39nNz3/Zcn9xKVuYTQjQ9YEyePJnJkyfXu83SpUvrzIuKiqo331dVP3Qwthmsg9FYYSEBZPRsxeeZp/j2QC59u8d7u0hCCB8gDx+8TLn2x5pfWXcY1dpow+iV0JJfTpfKAwqFEIAEjMuWV1Rhe0ptePN5Sm1j9ezckrbaUDZ8dYwjv5V4uzhCCC+TgHGZcs6U0kYbdkUMqa2LSqViQM9WtIzUsHzDj/a7KiHE1UkCxmUwW6zknCmlS9tIbxfF7QID/Jh4azIAz7+zD0OF7z/jSwjhHhIwLsOJXANVZutV83RXbVQwk//cg8LSSpZv+BGTjLMV4qokAeMyHD1pa89PaHPl32FU69I2iv834lqO/FbCyx/+JEu7CnEVkoBxGbJPniMmUkN0eFDDG19B+naLZ8ywrvxwtIBXPjwgQUOIq4wEjEZSFIWjJ0uumuYosK2bUWY0U2Y00zc5nj/f0JnMI/ms2CRBQ4iriUzfbaS8kgpKy01XRYd3NaPJwr4jFx4SGarx57ZBnfhg+y9Uvr+fif+XTIhGPkpCXOnkDqORsk+eAyDhKgoYtRnSuy1/vymRQ8eLeWLNHnKLZMitEFc6CRiNdPRkCSFB/rSOCfV2UbzuBl0bpo3SYagwMX/1Hn44kt/wTkKIZksCRiMdPXmOhLaRqFVX7oS9xkhsH81j49KIidSwbMOPvPDBjxTrjd4ulhDCDSRgNIKhwsSZwvKrajhtXWp2hIcEB/DwKB1/6H8N+48V8uhr37Jl969UGGV9cCGuJNJT2QjV/RdXU4d3XS7uCAeIDg9i5thUNmw/xrtfHmPzzl8Z3LsNQ3q39el1ioUQzpGA0QhHT5Xgp1bRsVWEt4vis7RRwTz8Fx05Z0r5+Ntf+XiX7ad755akdG6Jrov2qpu/IsSVQqUoiuLtQrhCYaEBq9V9p2KxWpnx8i60UcE8ckfvOrd796tstnx7wm3laCx/PxVmS+31ogIUoLo7RlEgKEBNldmKn9q2X4jGn/JKMyoVBAf5k9AmkuxT56ioNKM9/2j3m9Lb88nuE+QVV9hGjykKhaWVtIzQUGW2ogn052xROaVlVQAE+KsJCfKjvNJMa20oj/39OjbtyOGbH88QE6lhxphUFq7LJKl9NN/8eIYBPVoBcPhEMQAzxqSy8X+/2M/j8Ili+z7VZoxJBWDj/37h1oxO9u1vzehUb3rNvIu3q7lNQ+rafuP/fuGuP/UiP19/yXYL12Xay32x+vIac3xX7FPfuV3O/lptuL0+fNnl1Glj9/dUXajVKlq2DGv8fm4oyxVpz+F8CkuNDLuuXb3b+VKwAOoMFmALFmALFNWXDUaTFUW5sF95pdm+TXmlmf3HCimvNKMAecUV5J1fSKr6/9knz5F9qpRifRXZp0o5kWvgyG8l9mABYDJbOVdmwmRR+PWsgUnPbWfTjuMUlRo58ts5cs6UcuS3c/a0TTuOs2nHcY78do4jv9maBavTqtMBe3716+rtam7fUHrNvIu3a4y6tq/vWDXLfbH68hpzfFfs4+y5NfU4vqap5W6u512TBAwnKIrCJ9+dIK5FCL0SYrxdnCvO9Ret6Dd/9R4vlUQIUR8JGE448lsJv57VM7xPOxlO6wYd4h07xAfpWte7/X+2HXFncYQQdZCA4YRPdp8gPCSAfrK2tUdcHEAutn3f6UvSdvwoy8gK4W4SMBpwprCMfccKGdK7LYEBft4ujgD+OiThkrSVmw85vN554Cz55yo9VSQhrgo+ETBMJhPz5s0jPT2d9PR0Fi9ejNXq/aegWhWFD77+hQB/NYN7t/F2ccR5tS2LO+NvjiPXXvvvQWa8tBOwjVwTQjSdT8zDePbZZ9mxYwcrVqygrKyMRx55hIiICO6++26vlclitfL65kPs+Tmf2zI6EhES6LWyiIadLihzeD28TzvOFJaz/1ghW7/7zZ4+Z+VuOp+fqZ99qnEjj4S42nn9DsNoNPLWW28xa9YsevXqRb9+/Zg6dSpr1qzBW1NETGYrL37wE7sO5PLnQZ34Q/+OXimHuHxxLULQdbGNaFs2JcOeHhUWxHeH8gB48s0L8zaeWpfJqi2H2fLtrwDknCl1GArsTlZFsS97e4VMixJXKK/fYRw6dIiKigrS0tLsaWlpaeTn53Py5Enatat/3oMrFZVW8v3hPHb8eJaT+QbuGNqFoWmeO75wj5pTUe65NRlFUXjg+f8xuHcbvtx7CoDiUiOn88soOz/vpObQ3ukv25q2XvjgRyJCA4kKCyJM409wkD9BgX6oVSp7M5nJbKXKZAFsgyXKjWbKK02UV5oxVJgAmLjwcwwVVVQabdtNePpL+7HGP/0larWKoAA1mkB/NIH+BAfZ+s5W/PcAYcEBRIQGERkSQGhwACFBtnIEB/kRGOBHUIAfAf5Nvw5UFAVFsd1pW6yK/ad6cmxBSQVWLgQ41fnRg0WllajVtvrwV6vw81Pj76dCrVLZtxHNl9cDRm5uLmFhYYSGXnhcuFartec5GzBqa9duyA9H8zl6soTSMhPnyqrsazq0iQnj9sGdL2vORez52c9XkxBNQJPO++L9G3q/2vJrew+w/T4O5BTZ8w7kFNnTE9tH21//aVBnenXVUllp4qk3Mxk1tAvnDFVs+fZX2mrDUKlUFOmNnC4sw1LPZMhqsdHBfJV1CpUKNAH+aIL8CA6ylfHaji3wV0NQgB+7fjrLDSltUavhi8yTJHWIRlEULBYFk8WKJsifonOVxEYHc6agHKPZ4tQTDWKjg5n92rf4qdW2gKYCVCpqfmcriu0/VsX2xV+9j6JQ7zFio4NZvD6r1vRFb/1Q6z4qwM9fjb9aTWx0MAvX7SXAX42/n5rAQD8C1CoCA/xsr/1V+Pup8fNT46e2/Vulsv2Nq1Wg4vx5qFSoahyg+pVSPSVVOT859XxQUxSwYgt6VgWsFgWLomCxWDFbFCxWKybz+f9bFEzn69pqtd0FxkYHs2zDj/Zj2urSVpaa3z4K2N7faqXG4YmNDubFjT+hUqnwU3MhoKrVtrrxg9CQIMwmiy3NT1X/eaugfWwYnVo3/tl2l/N9CT7waJCNGzeyaNEidu7caU+zWq1ce+21rFq1in79+nmxdEIIIap5vQ9Do9FQVeXYVlz9WqPReKNIQgghauH1gBEfH49er6eiosKelp9ve2x2XFyct4olhBDiIl4PGElJSQQHB5OZeWHEyp49e4iNjaVNG5n7IIQQvsLrAUOj0TBy5EjmzZvH3r172bVrF8888wzjxo3zdtGEEELU4PVOb7DNxXjiiSfYvHkzQUFBjBw5kocffliG4QkhhA/xiYAhhBDC93m9SUoIIUTzIAFDCCGEUyRgCCGEcIoEjHo05rHrer2eqVOnkpqaSkZGBqtWrXJpvrd5si6ys7MZP348aWlpDBo0iKeffhqj0ei2c2ssT9ZFTRMmTOChhx5y6bk0lSfroqysjEcffZTrrruOvn37Mn/+/Esm/XqTJ+siNzeX++67j7S0NAYOHMgzzzyDxWJx27nZKaJOCxcuVG688UYlKytL2bFjh9K/f3/llVdeqXXbyZMnK3/961+VQ4cOKVu2bFF0Op3y0UcfuSzf2zxVFwaDQbnhhhuUf/7zn0p2draya9cuZciQIcqCBQs8cp7O8OTnotqGDRuUrl27KlOmTHHbeV0OT9bF/fffr4wYMULZv3+/snv3biUjI0NZsmSJ28/RWZ6si7/+9a/Kvffeq2RnZys7duxQ+vXrV+exXEkCRh0qKyuVXr16KV999ZU9bcOGDUr//v0Vq9XqsO3JkyeVxMRE5ejRo/a0ZcuWKbfffrtL8r3Nk3Wxbds2JTU1VTEajfb8Dz/8UElPT3fLuTWWJ+uiWl5entK/f39l5MiRPhUwPFkX2dnZSteuXZVDhw7Z89966y1lwoQJbjm3xvJkXZSUlChdu3ZV9u/fb89/6qmnlDvvvNMt51aTNEnVoaHHrteUlZVFVFQUCQkJDtseOHAAk8nU5Hxv82Rd9OjRgxdeeIHAwAsLVqlUKgwGg0+sFeHJuqg2d+5cRo0aRefOnd14Zo3nybrYtWsXnTp1IikpyZ4/atQoXn31VTeeofM8WRcajYaQkBDef/99qqqqOHv2LNu3b6d79+5uP08JGHVo6LHrF28bGxvrkKbVajGbzRQUFDQ539s8WRdxcXGkp6fb8ywWC2+++Sbp6ek+MZHTk3UB8PHHH5OTk+PV1Sfr4sm6OHHiBO3bt+c///kPw4cPZ8iQISxevNgnLqjAs3URFBTE3Llz+e9//4tOp2PQoEHExcXx4IMPuunsLvD6ehi+qqKigqCgIIe06qveizvaGtq2qfne5sm6uNj8+fM5fPgw77zzTtNOwkU8WRfFxcU8+eSTLF++3OGOy1d4si7KysrIzMzEZDLx9NNPU1RUxNy5czGbzcyYMcOl53U5PP03cuzYMdLT07nnnnsoKCjgiSeeYNGiRcyePdt1J1ULCRh1aMxj1xvatqn53ubJuqhmsViYO3cu77//PkuWLHFoivAmT9bFE088wU033YROp3PpObiKJ+vCz88Po9HI888/T0REBADl5SWx6XEAAALVSURBVOXMnj2bRx55BLXau40lnqyLb7/9ljfffJOvv/6asLAwAIKCgrj77ru55557iIlp/MJvzpKAUYeaj10PDrat1lbXY9fj4+PtedXy8vIICAggOjq6yfne5sm6ANvwxGnTpvHFF1+wdOlSfve737nr1BrNk3Xx0UcfodFoeP/994ELXxopKSn88EPtK9t5kifrIjY2Fq1Waw8WAJ06daKyspKioiK3fkk6w5N1sWnTJtq1a2cPFgDJyclYLBZOnz7t1rqQPow6NOax6zqdjsLCQnJycuxpmZmZJCcnExgY2OR8b/NkXQDMmTOH7du388orr/hUsADP1sXWrVvZtGkTGzduZOPGjQwZMoQBAwawceNG95+oEzxZF7179yY3N5fCwkJ7/tGjRwkLCyMqKsqNZ+kcT9ZFXFwcJ0+epLKy0p6fnZ0NQNu2bd11ijZuH4fVjM2fP18ZNmyYkpmZqezcuVPp37+/8uqrryqKoijFxcVKaWmpfdt77rlHuf3225UDBw4on3zyiaLT6ZSPP/7YZfne5qm6+Oqrr5SuXbsqa9euVfLy8hx+fIUnPxc1TZ8+3aeG1SqK5+rCYrEof/7zn5Vx48Yphw8fVnbu3KkMHDhQWbRokWdPuB6eqgu9Xq8MHDhQuf/++5WjR48q33//vXLzzTcrM2bMcPs5SsCoR2VlpfLoo48qKSkpSt++fZXFixfbx1T/7W9/U6ZPn27ftri4WJk8ebLSs2dPJSMjQ3njjTcc3qup+d7mqbqYPn260rVr11p/KisrPXOyDfDk56ImXwwYnqyLgoIC5cEHH1R0Op2Snp6uLFq0SKmqqnL/STrJk3WRk5Oj3H333UpaWpoycOBA5cknn/TI34c83lwIIYRTpA9DCCGEUyRgCCGEcIoEDCGEEE6RgCGEEMIpEjCEEEI4RQKGEEIIp0jAEEII4RQJGEIIIZwiAUMIIYRT/j+VYKgV1VFHqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(mutual_informations, hist=True, rug=True)\n",
    "plt.title(\"Permuted Mutual Informations with Mean = 0\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final Bit\n",
    "Return the p value!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "67b2465c7feb417095e9b716c3d4c6e2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-value when Mean is 0 0.421\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "11d9110d38294d05b7ef144b0f3396e8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-value when Mean is 1 0.0\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "db3de9fd4b7d4fc08aa8c76ba37aaaf7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-value when Mean is 1 0.0\n"
     ]
    }
   ],
   "source": [
    "def rf_independence_test(X, y, n_estimators = 200, max_samples = .63, bootstrap = True, depth = 4):\n",
    "    test_stat = get_mutual_information(X, y, n_estimators, max_samples, bootstrap, depth)\n",
    "    null_mutual_info = permutation_test(X, y)\n",
    "    p_value = sum(i >= test_stat for i in null_mutual_info)/len(null_mutual_info)\n",
    "    return p_value\n",
    "\n",
    "X, y = get_sample(1000, 0, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "print(\"P-value when Mean is 0\", rf_independence_test(X, y))\n",
    "\n",
    "X, y = get_sample(1000, 1, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "print(\"P-value when Mean is 1\", rf_independence_test(X, y))\n",
    "\n",
    "X, y = get_sample(1000, 5, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "print(\"P-value when Mean is 1\", rf_independence_test(X, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Compute the cdf of H(Y) evaluated at H(Y|X; under alternative)\n",
    "\n",
    "#### Reasoning behind this (?):\n",
    "$MI_0$(X,Y) - $MI_A$(X,Y; under alt) = H(Y|X; under alt) -H(Y|X; under null) \n",
    "\n",
    "Under the null, $X$ and $Y$ are independent. So, this means the conditional entropy $H(Y|X)$ should equal $H(Y)$. If we can construct the distribution of $H(Y)$, our p-value is just the probability that $H(Y|X)$ is close to $H(Y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To construct distribution of H(Y)\n",
    "\n",
    "We focus on 2 class case for now (unsure if this test will even work for more than 2 sample testing, my guess is no since mutual information involves just 2).\n",
    "\n",
    "We know that $Y \\sim Bernoulli(p)$. This is always the case since $Y$ is cateogrical and there are two classes. So $\\hat p$ is a binomial. We can use our sample to estimate p. Then, for each probability of $\\hat p$ which is easy with the binomial, we can contruct the entropy distribution as well by just plugging in that partiular value of $\\hat p$ into the entropy equation.\n",
    "\n",
    "### Entropy Equation:\n",
    "We are constructing the cdf of H(Y), where Y is Bernoulli. So the entropy is just $-plogp - (1-p)log(1-p)$. With our estimate, we get $\\hat p$ and the probabilty of seeing that $\\hat p$ is binomial. So, we can construct the distribution of H(Y) by assigning a probability to each entropy calculated at $\\hat p$ with the probability of seeing $\\hat p$\n",
    "### Code:\n",
    "Note in our case, our p parameter for the binomial should just be 1/2 since we are drawing evenly from $Y = 1$ and $Y = -1$.\n",
    "\n",
    "Note that in our code, we are also just considering the sum of y's instead of different $\\hat p$'s but these are identically distributed, $\\hat p$ just has a divided by n."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[nan, 0.007907255112232087, 0.014427214862176116, 0.020422924464179912, 0.02607783298339765]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/Cellar/jupyter/1.0.0_5/libexec/lib/python3.7/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log\n",
      "  if sys.path[0] == '':\n",
      "/usr/local/Cellar/jupyter/1.0.0_5/libexec/lib/python3.7/site-packages/ipykernel_launcher.py:12: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  if sys.path[0] == '':\n"
     ]
    }
   ],
   "source": [
    "#TODO: have to change to work with other categorical data not just 1 or 0\n",
    "from scipy.stats import binom\n",
    "import math\n",
    "\n",
    "def get_entropy_dist(y):\n",
    "    n = len(y)\n",
    "    p = y.count(1)/n\n",
    "    dist = {}\n",
    "    for i in range(0, n + 1):\n",
    "        p_hat = i/n\n",
    "        prob = binom.pmf(i, n, p)\n",
    "        dist[-p_hat*np.log(p_hat) - (1 - p_hat)*np.log(1 - p_hat)] = prob\n",
    "    return dist\n",
    "\n",
    "X, y = get_sample(1000, 0, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "\n",
    "\n",
    "dist = get_entropy_dist(y)\n",
    "print(list(dist.keys())[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's hard to visualize the distribution so we are going to monte carlo just to view the distribution. Also, problems with analytical is negative entropies from approximation (?). Better to just do monte carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "509\n",
      "1000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEkCAYAAACCIF6UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XdcU9f7B/BPwt57u0AELMpQUBHQuletlWorVeu3WrSKxYF1a9W6qmKrQK1QtGqto9afVpx1b6oILgRluCUsUXbW/f1BSYlJIAGSMJ7365uXcM695z73fGmeO849l8UwDANCCCFEjdjqDoAQQgihZEQIIUTtKBkRQghRO0pGhBBC1I6SESGEELWjZETUggZx1h31HWmOKBm1UH379oWrqysiIyOl1sfHx8PV1RWffPJJg273zZs3mDlzJrKyshqkvSNHjmDcuHHw9fVFt27dMH78eFy4cKFB2h4/fjxmzZrVIG1V5+rqij179gAAIiMj4e/vL/e6e/bswdatW2tc5t02q2+vrl69eoX//e9/qKioAAAkJCTA1dUVGRkZ9Wq3ISQnJ2Pw4MHo3LkzVq1aJVH//PlzuLq64uLFi1LX/+STTzB//nyJ8v/973+4cOECIiMja1w/PDwcXbt2BYfDwfXr1/HZZ5/RAUMdUDJqwVgsFk6fPi217tSpU0rZZmpqKo4fP17v/1gZhsHcuXOxZMkSeHp6IiIiAt9//z2sra0xefJk7N27t4EiVq7Ro0cjJiZG7uVjY2NRXFzcoG3K49q1a7h27Zrod3d3d+zbtw+tWrVq0O3URVRUFBiGQUxMDMaPH98gbR46dAgCgQC9e/fGlClT0K5dO6xcuRJcLldsuWvXriE+Ph7ffPMNbGxs0KNHDxgYGDSZv7/GhJJRC+bl5YXU1FQ8f/5crLy8vByXLl2Ci4uLmiKr3b59+3DkyBH8/PPP+Oabb9CrVy/06dMHERERCAoKwurVq1FQUKDuMGtla2sLd3f3Rt/muwwNDeHl5QUdHR2lbkceRUVF8PLygp+fH9q0aVPv9vh8Pn788UdMmjQJAKCtrY1ly5bhyZMniI2NFS3H5XKxfPly+Pr64tNPPxWVT5w4EVFRURKJi9SMklEL1qVLF1hYWEicHV26dAlmZmYSX2hcLhdRUVEYOHAgPDw8MHLkSJw7d05UX3Xp5tatWxg1ahQ6d+6MQYMGidpPSEjA559/DgAYOnSo6BIhl8vFunXrEBAQAA8PD4wfPx6pqak1xr5jxw70798fPXr0kKibPn06xowZg5KSEgCVZ1G//PILhg4dik6dOsHHxwfTp08Hh8MRrePq6oqYmBgMGjQI3t7eSEhIkGg3Ly8P8+fPR0BAALy8vBASElLrZSoOh4PQ0FB4e3ujX79+EpcQ372klpSUhDFjxsDb2xs9evTAvHnz8Pr1awCVl1ZfvHiB2NhY9O3bF0DlpcTly5dj/Pjx8PDwwNatW6Ve+svJycGkSZPg4eGBwYMH4+jRo6K6gwcPwtXVVXQJDgAyMjLg6uqKhIQEHDx4EAsWLAAAeHh44ODBg1Iv0x07dgwjR46Ep6cn+vXrh9jYWLEzYFdXVxw+fBjTp0+Hl5cX/P39ERUVVWP/CYVC7N69G8OGDYOHhweGDBmCP//8U6zN5ORkHDp0CK6urhIHVnVx6tQplJSUiPWhn58fRowYgZiYGNE2tm3bhpcvX2LlypVgsViiZbt16wahUIgjR47UO5aWhJJRC8Zms9G3b1+cOXNGrPzUqVMYNGiQxPJz5szBtm3bMG7cOERFRcHZ2RlTp04VS0hVy3300UfYunUr7O3tMXv2bBQWFsLd3R1Lly4FAPzwww8YPXo0AGDBggXYv38/pk6dis2bN0NbWxuff/45srOzpcbN4XCQmZmJwMBAqfUODg5YuHAhWrduDaDy0lZUVBTGjh2Lbdu2Yfbs2bh+/TrWrVsntt6WLVvw1VdfYfny5fDw8BCrKy4uxpgxY5CUlIQFCxZgw4YNeP36NT777DO8evVKahwCgQBffvkl0tPT8f333yM0NBSLFi2SuixQeYT/1VdfwcHBAT/99BOWLFmCK1euYPny5QAqL0dZWVlhxIgRYl/i+/fvh6enJyIjI9GvXz+pbW/duhUODg6IiopC165dER4ejsuXL8uMpbr3338fU6dOBQD89ttveP/99yWW+e233zB79mx069YN0dHRGDlyJDZt2oT169eLLbdy5Uq0adMGW7ZsER2Q1HSPb/369VizZg2GDh2Kn376CQEBAVi4cCF2794NoPIMuUOHDujduzf27dsHa2trmW0xDAM+ny/xedfx48fRq1cvaGlpiZXPnz8furq6iIiIQE5ODrZu3Yrp06ejXbt2YstpaGjg/fffx7Fjx2TGQiRpqjsAol4DBgzA1KlT8fr1a5iZmYHH4+HcuXOIjY3Fvn37RMulpqbi5MmTWL9+PT788EMAQK9evZCTk4Mff/wRffr0ES0bEhKC4OBgAICdnR0GDx6MhIQEDBo0CM7OzgAqj2htbW3x6NEjxMfHY+PGjRg2bBgAICAgAIMGDUJcXJzUL++qMxo7Ozu59jEnJwdhYWEYO3YsgMoj18zMTIkk3LdvX4wcOVJqG3/++Seys7Nx4sQJ0X2S7t27o1+/fti2bZvUOM+fP4+HDx/iyJEjokueBgYGCAsLk7qNjIwMFBYWYty4cfD29gZQeTksMzMTAPDee+9BW1sb1tbWeO+990TrWVhYYM6cOTX2gZ+fH1asWAGg8v+39PR0bN++HQEBATWuBwDm5uaiy18eHh4Sl+YEAgEiIyMxatQo0RlUQEAAWCwWtmzZgi+//BLm5uai8rlz5wIAevTogRMnTuDixYvo3bu3xHYLCgqwa9cuhIaGipJhQEAASkpKsHnzZnz66afw8vKCvr4+zM3N4eXlVeN+TJ48WWadk5OT6Od//vkHU6ZMkdoPc+bMweLFi/Hq1Su0a9cOEydOlNree++9hxMnTkAoFILNpmN+eVAvtXB+fn7Q09MTnd1cu3YNenp6Ev9hJyYmgsViYfDgwWLlQ4cORWpqqthNdU9PT9HPtra2AICysjKp279x4waAyi+Z6keqPXv2lHqpDKg88gQqL+HIY/HixZg4cSLy8vKQkJCA3bt349atW+DxeGLLOTo6ymwjMTERnTp1Erthb2RkhMDAQCQmJkpd59atW7CzsxO799a3b1+ZX07Ozs4wNTXF1KlTsXLlSly9ehX+/v744osvaty/d4/MpRkwYIDY7++//z5u3bpV63ryyMzMRGFhodS/DR6Ph9u3b4vKqv9tsFgs2NjYoLS0VGq7d+7cAY/Hk9puYWGhKEnLa/HixThw4IDEp/r/P6WlpSgsLBT93b5r1KhR6Nq1K+7cuYNVq1ZBU1P68by9vT1KS0vx5s0bhWJsyejMqIXT1tZG7969cfr0aQQFBeHUqVMYMGCA2DVwoHJItpGREbS1tcXKLSwsAEB0fwaA2JFz1RevrMRRWFgIoPJs5V1VR9PvqjojknV5DACys7NFXyjp6elYtGgRkpOTYWBgAHd3d+jo6EiM6KvaF2nevn0rtd7c3Bz37t2TuY6ZmZlYmZaWFkxMTKQub2hoiF27diEqKgoHDhzArl27YGFhgfDwcHz88ccyY6spblnLmJmZobS0FAKBoNZ1a1P1hWtpaSl1m9UPVN49q2Kz2TJHVirSrjzatm2Lzp07S5Tr6emJfi4qKpIoq47FYsHPzw9PnjwROzuV1WZxcbHE3wCRjpIRwYABAzB//nyUlJTg7Nmz2Lhxo8QyJiYmKCoqApfLFUtIeXl5ovq6MDY2hra2Nn7//XeJOllHnebm5nBzc8PVq1cxZswYifoXL16gb9++WLp0KYKDgzF16lTY2tri+PHjcHR0BIvFwvr16/H06VOF4szJyZEoz8/Pl7nvJiYmyM/PFysTCoU1fom6uLhg8+bN4HK5uHbtGuLi4rB48WL4+/vLPFqXx9u3b6XGraGhITrwqJ6YZJ2tSFO1/1V/C1Xq+7dRvV0jI6MGa7cmpqamABRPdO+q6m9KRPKjy3QEvXr1AsMw+OmnnyAUCuHr6yuxTJcuXcAwDE6cOCFWfvz4cXTs2BG6urpybavqElsVb29vcLlcCAQCdO7cWfQ5fPiwxLaqGzduHP7++2/RZb7qIiMjoaWlhf79+6OgoABPnz5FcHAwnJycwGKxIBQKce3aNbkv8wGV+3/v3j2x0VrFxcW4fPmyzHsVvr6+4HA4uHPnjqjs2rVrEpcHq9f5+fmhoKBAdMYaHh4OoVAoSoR1vf9w9epVsd9Pnz6Nrl27AgD09fUBQGzAyLuXHmvarpOTE0xNTaX+bWhoaEgMBpGXh4cHtLS0pLZrZmYm1+VJReno6MDMzExspGVdcDgcGBoawtDQsIEia/7ozIjAwMAAPXv2xI4dO/DRRx9JJAwA6NixI/r3749ly5ahsLAQjo6OiI+PR0JCQq3Dc6urOsI9c+YMdHR04O7ujj59+iAsLAxff/01WrdujVOnTmH37t3YsGGDzHY+/vhjnDlzBiEhIfj888/RvXt3lJSU4P/+7/9w9uxZrFy5EjY2NmAYBnZ2doiLi4OBgQGEQiH27NmDlJQUuRNo1fZ27tyJSZMmYcaMGdDW1kZMTAyEQiEmTJggdZ3AwEB4e3tjxowZ+OabbyAQCBARESExSqtKp06dwGKxEBYWhpCQEABATEwMWrVqhY4dOwKoPEO7ffs2kpOTa71hX93p06exadMm+Pr64s8//8SjR49Eo/S6desGLS0trFy5EiEhIcjMzMTOnTvF1jc2NgZQmQh69uwpVqehoYFp06ZhzZo1MDAwQK9evZCcnIwtW7Zg/PjxorMNRZmbm+Ozzz5DdHQ0hEIhvLy8cPHiRRw8eBALFy6U+nfaEHr06CF2n6subt++DT8/vwaKqGWgMyMCoPJSHY/Hk7jRXV1ERAQ++eQTxMTEIDQ0FBkZGdiyZQv69+8v93Y6dOiAYcOGYfPmzdi2bRuAymHegwcPxqZNmzB58mTcvHkT69evx/Dhw2W2w2azERUVhZkzZ+LixYsICwvDkiVLUFJSgm3btomGjbNYLGzevBlsNhthYWFYunQpDA0NsXHjRpSVlSEtLU2uuI2MjLBr1y64urpi8eLFmDdvHszMzLBv3z44ODjIjPHnn39Gly5dsHjxYqxbtw5hYWEyLy8ZGRkhNjYWmpqaCA8Px8yZM2FgYIC4uDhRAvvyyy+RlpaGkJAQqcOSZQkPD8fNmzcxZcoUPHz4ED///LPojMXMzAw//vgjXr16hZCQEPz111/YtGmT2Pp+fn7o0aMHFi9ejMOHD0u0P2HCBCxbtgznz5/HlClTcOjQIcyaNUvqNDuKmD9/PkJDQ3HgwAF89dVXuHr1KlavXi16Xk0Z+vfvj3/++UfmGWxtBAIB/vnnH4X+uyAAi147Tggh/6k6KFu4cCEGDhyo8Prnzp3Dt99+i9OnT0sM+CGy0ZkRIYRUo6WlhdDQUPz66691Wn/Hjh2YMmUKJSIFUTIihJB3jBo1Cpqamjh//rxC6129ehXl5eWih76J/OgyHSGEELWjMyNCCCFqp/JklJ2djbCwMHTv3h3+/v5YuHCh6AGxqtFK1T/V54jKzs7G5MmT4e3tjf79++Ovv/6SaLumekIIIY2TSp8zEgqFCA0NhampKXbs2AEul4tly5ZhwYIFiI6OxqNHjxAUFITZs2eL1qk+fUhoaCjs7Ozwxx9/4ObNm1i4cCFatWqFLl26yFVPCCGkcVJpMkpLS8O9e/dw+fJlWFlZAQAWLVqEsWPHori4GBkZGRgwYICorrobN24gLS0N27dvh7GxMZydnXH79m3s2rULXbp0qbVeEa9fl0AorNutNAsLQ+Tn128qkZaI+k1x1Gd1Q/2muJr6jM1mwczMoN7bUGkysrOzQ2xsrFiyYbFYYBhGlIxkTfOelJQEV1dX0ZPgAODj4yN6+r+2ekUIhUydk1HV+kRx1G+Koz6rG+o3xSm7z1R6z8jU1BS9evUSK/v111/h6OgILS0tvH79GsePH0e/fv0wYMAAREREiF7dy+FwJF6cZWlpKZpDqrZ6QgghjZda56aLiYnB33//jZiYGNHri42MjBAdHY3Hjx9j1apVePv2LZYvX46ysjKJ6ee1tbUhEAjA5/NrrZc1A7Q0Fhb1m9zQysqo9oWIBOo3xVGf1Q31m+KU3WdqS0bR0dHYvHkzli5dKnp99PXr10VTrru5uQGonFNr0aJF0NXVFb37pgqXy4WWlhY0NTVrrVdEfn5xnU9JrayMkJtbVKd1WzLqN8VRn9UN9ZviauozNptV7wN4QE3PGa1evRqRkZFYtmyZ6FXQgOS7P5ydncHn81FQUABbW1vk5uaK1efm5oouzdVWTwghpPFSeTKKiorCb7/9hrVr14pNmbF3717069dP7B0zKSkpMDQ0hLW1Nby8vJCWlib20qvExER4e3sDQK31hBBCGi+VJqO0tDRER0dj0qRJ8Pf3R25urujTs2dPFBQUYMWKFXj8+DHOnj2LdevWISQkBGw2Gz4+PnBycsKcOXPw8OFD7Nu3D/Hx8Rg/fjwA1FpPCCGk8VLp3HSRkZEyh1ofO3YMhYWFWL9+PVJTU2FsbIxPP/0U06ZNE70W+fnz51iyZAkSExNhY2ODGTNm4IMPPhC1UVu9vOiekepRvymO+qxuqN8Up4p7RjRRqhSUjFSP+k1x1Gd101j7jS8EKniSL0zU0dKEpppnEVVFMqLXjhNCSCNQwePjxgPJ5yJ9O9pAU6f5f1XTrN2EEELUjpIRIYQQtaNkRAghRO0oGRFCCFE7SkaEEELUjpIRIYQQtaNkRAghRO0oGRFCCFE7SkaEEELUjpIRIYQQtaNkRAghRO0oGRFCCFE7SkaEEELUjpIRIYQQtaNkRAghRO0oGRFCCFE7SkaEEELUjpIRIYQQtaNkRAghRO0oGRFCCFE7SkaEEELUjpIRIYQQtaNkRAghRO0oGRFCCFE7SkaEEELUjpIRIYQQtaNkRAghRO0oGRFCCFE7SkaEEELUjpIRIYQQtVN5MsrOzkZYWBi6d+8Of39/LFy4EG/fvgUAFBUVITw8HF27dkVgYCC2b98utm596wkhhDROmvIsVFFRAR0dnXpvTCgUIjQ0FKamptixYwe4XC6WLVuGBQsWIDo6GosWLUJOTg52796Nx48fY8GCBbC2tsawYcMAoN71hBBCGie5kpGfnx8GDhyI4cOHw8/PD2x23U6o0tLScO/ePVy+fBlWVlYAKhPI2LFj8eLFC5w6dQrx8fFwdnaGm5sb0tPTsWPHDgwbNqze9YQQQhovubLK4sWLkZubi8mTJ6NXr15YuXIl7ty5o/DG7OzsEBsbK0pEAMBiscAwDG7evAlTU1M4OzuL6nx8fHD//n3weDwkJyfXq54QQkjjJdeZUVBQEIKCglBQUIDjx4/j6NGj2L17N1q3bo3hw4fjgw8+gKOjY63tmJqaolevXmJlv/76KxwdHZGfnw9ra2uxOisrK/D5fOTl5YHD4dSr3s7OTp5dJYQQogZyJaMq5ubmGDt2LMaOHYvU1FR89913iI6Oxk8//QQPDw98/vnnCl0Si4mJwd9//42YmBjcuXNH4r6UtrY2AIDL5aKsrKxe9YqwsDBUaPl3WVkZ1Wv9lor6TXHUZ3XTGPuNKSiFkaGuRLm+vg6szPXVEJE4ZfeZQskoLy8Px48fx7Fjx3D79m1YWVlh0qRJCAgIwJUrV7BgwQLcvHkT3377ba1tRUdHY/PmzVi6dCkCAwPx8OFDiaRR9buuri50dXXrVa+I/PxiCIWMQutUsbIyQm5uUZ3Wbcmo3xRHfVY3jbXfSiv4KCoulywvrUCuQKCGiP5TU5+x2ax6H8ADciaj/fv34+jRo7h58yZ0dXUxcOBAhIWFoUePHmCxWAAqBznw+Xzs27ev1mS0evVq7Ny5E8uWLUNwcDAAwNbWFrm5uWLL5eTkQEtLC2ZmZvWuJ4QQ0njJlYxWrFgBf39/rFu3Dv3795c5zNvDwwOGhjVnyKioKPz2229Yu3YtPvroI1G5l5cX8vPzkZWVJbr/lJiYiE6dOkFbW7ve9YQQQhovuZLR77//Djc3N6lf6hUVFXjw4AG8vLwwdOjQGttJS0tDdHQ0vvzyS/j7+4udydja2qJPnz6YN28eli1bhmfPniEuLg6rV68GADg4ONSrnhBCSOPFYhim1psjHTt2xL59++Dh4SFRd+3aNXz11Ve4fft2rRuLjIxEVFSU1Lpjx47BwsICS5cuxYULF2BiYoJJkyZhwoQJomUKCwvrVS8vumeketRviqM+q5vG2m8lFXzceMCRKPftaAMDHYVu7zc4VdwzkpmMvvnmGzx79gwAkJycjA4dOsDAwEBiuadPn8LQ0BCnTp2qdzCNBSUj1aN+Uxz1Wd001n5r6clI5h5++umn+PPPPwFUJiN7e3uYm5u/EwQb3t7e+OSTT+odCCGEkJZLZjLy8fGBj4+P6Pdp06ahdevWKgmKEEJIyyIzGXG5XNGAheXLl4vKZKERa4QQQupKZjLy9PQUDVrw8PAQPU8ky4MHDxo8OEIIIS2DzGS0evVq0WW51atX15qMCCGEkLqSmYxGjhwp+jkoKEglwRBCCGmZ5H4x0d69e3H27FkAwJ07dzBkyBB4e3tj3rx5KC0tVVqAhBBCmj+5klFUVBRWrlyJly9fAqh8IZ5AIMC0adOQkJCADRs2KDVIQgghzZtcyejgwYOYM2cOxo0bh9TUVDx69AihoaEICQnBnDlzcPLkSWXHSQghpBmTKxnl5eWhU6dOAIDTp09DU1MTffr0AQBYWlrSZTpCCCH1IlcyatWqFZKSksDj8XDs2DF4e3vD2NgYAHD8+HG0bdtWqUESQghp3uRKRpMmTcIPP/wAPz8/ZGVlYeLEiQCAMWPG4MCBA5gyZYpSgySEENK8yTX73scff4w2bdogKSkJXbt2RdeuXQEAffv2xbx58+Dt7a3UIAkhhDRvck8F6+vrC19fX7GyyZMnN3hAhBBCWh65k9GlS5dw4cIFlJWVQSgUStSvWbOmQQMjhBDScsiVjKKjoxEZGQlra2vY2NiAzRa/1URTBRFCCKkPuZLRnj17MGHCBCxYsEDZ8RBCCGmB5BpNV1JSInquiBBCCGlociWjgIAAXL58WdmxEEIIaaHkukw3YMAArFq1Cq9evYKHhwd0dXXF6lksFr16nBBCSJ3JlYzmzp0LADh69CiOHj0qUU/JiBBCSH3IlYxSU1OVHQchhJAWTO73GQGAQCBAamoqLl68iOLiYuTn5ysrLkIIIS2I3A+9/v7774iKikJBQQFYLBYOHDiAzZs3QyAQYNOmTTAwMFBmnIQQQpoxuc6Mdu/ejVWrVuHTTz/Fb7/9BoZhAACfffYZ7t27hx9++EGpQRJCCGne5EpG27dvR2hoKGbMmCE2KWrv3r0RHh6Ov//+W2kBEkIIaf7kSkY5OTnw9PSUWteqVSu8fv26QYMihBDSssiVjNq3b48TJ05Irbt48SLat2/foEERQghpWeQawBAWFobQ0FBwOBwEBgaCxWLh6tWr+OOPP3DgwAFs3LhR2XESQghpxuQ6M+rTpw9+/vlnvH79GqtXrwbDMIiIiMCtW7ewceNGDBw4UNlxEkIIacbkHtrdq1cv9OrVC2VlZXj79i0MDQ1pODchhJAGIfPMKCsrS+onOzsbpaWlyMnJEStXVEVFBYYNG4aLFy+Kynbt2gVXV1exz5QpU0T12dnZmDx5Mry9vdG/f3/89ddfYm3WVk8IIaRxknlmNGTIEIVemvfgwQO5ly0rK8OsWbOQnp4uVv7o0SMEBQVh9uzZojIdHR3Rz6GhobCzs8Mff/yBmzdvYuHChWjVqhW6dOkiVz0hhJDGSWYy2rlzp+jnZ8+eYcWKFRg9ejT69+8PS0tLFBYW4tKlS/j999+xZMkSuTd47949zJs3D1paWhJ1GRkZGDBgAKysrCTqbty4gbS0NGzfvh3GxsZwdnbG7du3sWvXLnTp0qXWekIIIY2XzGTUrVs30c8bN27E5MmTERoaKraMj48PDA0NERcXhw8//FCuDV67dg29e/dGWFiYxLNLGRkZmDx5stT1kpKS4OrqCmNjY7HtR0VFyVVPCCGk8ZJ71u53E1EVFxcXPHnyRO4NhoSESC3Pz8/H69evcfz4caxYsQJsNhuDBw/G119/DW1tbXA4HFhbW4utY2lpCQ6HAwC11hNCCGm85EpGLi4u2LNnD/z8/KCp+d8qFRUV2LZtGzp37lzvQDIyMgAARkZGiI6OxuPHj7Fq1Sq8ffsWy5cvR1lZmdj9IwDQ1taGQCAAn8+vtb563LWxsDCs175YWRnVa/2WivpNcdRnddMY+40pKIWRoa5Eub6+DqzM9dUQkThl95lc39Dz5s3DpEmT0K9fP/j5+cHMzAz5+fm4fPkyBAIBduzYUe9AunXrhuvXr8PMzAwA4ObmBgAIDw/HokWLoKuri8LCQrF1uFwutLS0oKmpWWu9IvLziyEUMnXaDysrI+TmFtVp3ZaM+k1x1Gd101j7rbSCj6Licsny0grkCgRqiOg/NfUZm82q9wE8IOdDr127dkV8fDwGDx6MzMxMnD17Fk+ePEFQUBDi4+NFiaO+qhJRFWdnZ/D5fBQUFMDW1ha5ubli9bm5uaJLc7XVE0IIabzkPmVo1aoVFixYoLRA9u7di9jYWPz9999gsytzZEpKCgwNDWFtbQ0vLy9ERUWhuLgYhoaVWTgxMVE0i3ht9YQQQhovhd70qkwBAQEoKCjAihUr8PjxY5w9exbr1q1DSEgI2Gw2fHx84OTkhDlz5uDhw4fYt28f4uPjMX78eACotZ4QQkjj1WiSUatWrfDLL78gNTUVH330EZYtW4bg4GDRDAxsNhtRUVGoqKjAqFGj8Msvv2DNmjXw8vKSq54QQkjjxWKqXttKRGgAg+pRvymO+qxuGmu/lVTwceOB5KMovh1tYKCj2CCshtZoBjAQQgghyiQz3So6+amjo2O9gyGEENIyqWWiVEIIIaQ6lU/EoKpCAAAgAElEQVSUSgghhLxL5ROlEkIIIe+SawBDamoqPDw8pNYpOlEqIYQQ8i65klHVRKl8Pl+svCEnSiWEENJyNZqJUgkhhLRcciWjqolSd+3ahaSkJLx58wampqYICgrChAkTpL6ZlRBCCJFXo5kolRBCSMsldzLicDjYtm0bEhISkJeXh19++QUnTpyAi4sLhg4dqswYCSGENHNyj6b78MMPcebMGXh4eCA/Px8CgQDFxcUIDw/H0aNHlR0nIYSQZkyuM6NVq1ahc+fOiImJgVAoxP79+wEAixcvhlAoRGxsLIYNG6bUQAkhhDRfcp0Z3blzB+PHjwebzZaYImjQoEF4/PixMmIjhBDSQsiVjIyNjfHq1Supdc+ePYORkVGDBkUIIaRlkesy3fDhw7Fx40aYmJggICAAAMBisXDnzh1ERkbSJTpCCCH1IlcymjVrFjgcDmbNmiW6TBccHAwul4vevXtj1qxZSg2SEEJI8yZXMtLS0kJERASmTp2KGzduoLCwEMbGxvD29sZ7772n7BgJIYQ0c3Iloxs3bqBz585wdnaGs7OzWB2Hw8Eff/yB6dOnKyVAQgghzZ9cAxjGjx+P0aNHIzMzU6IuOzsb0dHRDR4YIYSQlkOuZAQAfD4fo0aNQnx8vDLjIYQQ0gLJnYzWrFmDL774AnPnzsW3334LLperzLgIIYS0IHInIzabja+//hpbt27FyZMnMWbMGDx79gwaGhrKjI8QQkgLIHcyqhIYGIiDBw9CU1MTQUFBuH79ujLiIoQQ0oIonIwAwN7eHrt378bw4cOxYcOGho6JEEJICyPX0O41a9agdevWYmVaWlpYunQpfHx8cOnSJaUERwghpGWQmYy4XC60tbUBQDTdj7RBC/3790f//v2VFB4hhJCWQGYy8vT0xL59++Dh4QEPDw+J2bqrY7FYSElJUUqAhBBCmj+ZyWj16tWiS3OrV6+uMRkRQggh9SEzGY0cOVL0c1BQkEqCIYQQ0jLJTEYbN26UuxEWi6XwzN0VFRUICgrCvHnz0KtXLwBAUVERli1bhvPnz0NfXx8TJ07EF198IVqnvvWEEEIaJ5nJSJFpfxRNRmVlZZg1axbS09PFyhctWoScnBzs3r0bjx8/xoIFC2BtbS0aQFHfekIIIY2TzGR09uxZpWzw3r17mDdvHrS0tMTKX7x4gVOnTiE+Ph7Ozs5wc3NDeno6duzYgWHDhtW7nhBCSOOl8EOvPB4PXC4XXC4XFRUVeP36Nc6fPy/3+teuXUPv3r2xd+9esfLk5GSYmpqKvaLCx8cH9+/fB4/Hq3c9IYSQxkuuh14zMjKwZMkSJCcng2EYqcs8ePBArg2GhIRILedwOLC2thYrs7KyAp/PR15eXr3r7ezs5IqPEEKI6smVjFatWoWXL19i7ty5OH/+PDQ1NdGvXz9cuXIFZ86cwbZt2+odSFlZGXR0dMTKqh665XK59a5XhIWFoULLv8vKyqhe67dU1G+Koz6rm8bYb0xBKYwMdSXK9fV1YGWur4aIxCm7z+RKRklJSVi2bBlGjBgBAwMD/PnnnwgODkZwcDDmzZuHnTt3ws/Pr16B6OrqSiSNqt91dXXrXa+I/PxiCIXSzwBrY2VlhNzcojqt25JRvymO+qxuGmu/lVbwUVRcLlleWoFcgUANEf2npj5js1n1PoAH5LxnJBQK4eDgAABo37490tLSRJfrPvzwQ9y9e7fegdja2iI3N1esLCcnB1paWjAzM6t3PSGEkMZLrmTk5OSEpKQkAJXJqKysDGlpaQAqL6+VlpbWOxAvLy/k5+cjKytLVJaYmIhOnTpBW1u73vWEEEIaL7mS0fjx4/HDDz9g/fr1MDExgb+/P+bOnYu4uDisX78enp6e9Q7EwcEBffr0wbx585CSkoKTJ08iLi4OEyZMaJB6QghpzB49L8SZm8/xprhC3aGohVzJKCgoCBs2bICJiQmAyldKmJqaYvPmzbC0tMTy5csbJJi1a9fC1tYWwcHBWLVqFWbOnIkhQ4Y0WD0hhDRG5249R9Sfd/EirwTnbr0Al6fee0TqwGJkjdVuwWgAg+pRvymO+qxuGlO/8QVC/P73Q5xPfon3HM1ha6aHc0kvYG9pgD5dHMBmseDb0QYGOnKNNVMaVQxgkHsPX758iZSUFBQXF0ut/+ijj+odDCGEtCQXkl/ifPJLDOneBoN7tEViWg66dbRBQgoHtx/lwdvFSt0hqoxcyWjfvn347rvvwOfzpdazWCxKRoQQogAhw+B04nM42hljdB9nlFRUfr+6tDZBwdty3M0sgL2VgZqjVB25ktFPP/2Evn37YtGiRTRMmhBCGsD9rAJwCkoRMvw9sXIWiwXfjtbIfPkWj181jsuJqiBXMiopKUFwcDBsbGyUHQ8hhLQIZxKfw8RAG75u1hJ1mhps2Fsa4FlOscwp2JobuUbTDR8+XKFXShBCCJGNU1CKOxn56O1lD00N6V/Dra0NUVrOx/PcEhVHpx5ynRnNmTMHQUFBGDx4MDp27Ch1ep01a9Y0eHCEENIcnbn1HBpsFt73dpC5TCvryvtFdzPy4NbaVFWhqY1cyWjt2rV48uQJHBwc8PLlS4l6FovV4IERQkhzVM7l48rdV/B1s4apoY7M5XS1NWFlqoe7mQUY/b7q4lMXuZLR0aNHMXv2bEyePFnZ8RBCSLOWkMJBWYUA/bq2qnXZ1jaGuJWWi/w35bAwUWzC56ZGrntG+vr6cHd3V3YshBDS7N1IzYGNuT6c7I1rXba1VeXDpMnpecoOS+3kSkZjx45FTEwMiopazjBDQghpaEWlXKQ+KYSPq5VctzdMDLVhbaaH5Ee5tS7b1Ml1mS4zMxP37t2Dv78/2rZtCwMDyQex3n2NOCGEEHFJj/IgZBj4uEoO55als5MFzie9QGk5H/q66p0WSJnk2jNNTU0MHDhQ2bEQQkizdjMtB1amumhjI/9cbh7tLXAm8TnuZeWjW8fm+6ynXMmoZ8+e6N69O6yt5c/mhBBC/lNSzsODx6/R36c1SrmSs3LLmpu5rZ0xtLXYSH/+plknI7nuGS1duhR37txRdiyEENJsJT/Kg0DIoFN7c9x4wJH48IVCqetpsFloZ2OEzFdvVRyxasmVjBwdHfHixQtlx0IIIc1WYlouzI110NbGSOF1nRxM8JRTBB5fesJqDuS6TBcYGIiIiAicPHkSjo6OsLCwkFhm9uzZDR4cIYQ0B2UVfNzLykffLq3qNEmAk50x+AIGT3OK0N7eRAkRqp9cyejIkSOwtLREdnY2srOzJepZLBYlI0IIkeF2eh74AsVG0VXX3qEyAWW+fNuyk9HZs2eVHQchhDRbN9NyYWqoDScHY5RJGbxQGzMjHZgZ6SDzZfO9b6TQoPXU1FQkJiaipKQEpqam8PLygouLi7JiI4SQJq+cy8fdzHz08rQHux7zeDrZGSPz5ZsGjKxxkSsZ8Xg8zJ07F8ePHweLxYKhoSGKiorAYrEwYMAAREREQEtLS9mxEkJIk3MnIx88vlDqe4sU4eRgjMSHuXhbyoWxvnYDRdd4yDWabtOmTbhw4QK+//573L59Gzdu3MCdO3ewZs0aXLlyBdHR0cqOkxBCmqQbqTkwMdCGs0P97vU42VXOZddcL9XJlYz++usvzJo1CyNGjIC2dmVG1tbWxkcffYQZM2bg8OHDSg2SEEKaogquAHcz8tHV1Qpsdv1etdPO1hhsFqtlJ6OioiI4OTlJrXNyckJ+fn6DBkUIIc3Bncx8cPnCOo+iq05HWwOtrAyQ1UzvG8mVjNzc3GS+djw+Ph7t27dv0KAIIaQ5uJmaA2N9Lbg00JtaneyNkfnqLYSMjLmDmjC5BjCEhYVh4sSJyMnJweDBg2FhYYH8/HycOHECV69eRWRkpLLjJISQJqWCJ8CdjHz4dbKt9yW6Kk72Jjif/BLZ+aWwt5R8e0JTJlcy8vPzQ3R0NCIjI7F06VIwDAMWiwVXV1dERkaif//+yo6TEEKalHuZ+ajgCeDratVgbVa9kC/j5ZuWmYwAoG/fvujbty9KS0tRVFQEQ0NDqe81IoQQAly/z4GRvhZc2jTMJToAsLXQh56OBh5nFyHQo8GabRRkJiMulyt9BU1NmJmZSSxTNcqOEEJaurclXCSn52GAT2tosOW6NS8XNouFNtZGePyq+b11W2Yy8vDwkHtCPxaLhZSUlAYLihBCmrKr97IhEDII9LRr8Lbb2RnhTOIL8AVCaGo0XKJTN5nJaPXq1TUmo5ycHMTGxqK4uBienp5KCY4QQpoahmFw8fZLOLcygZ1Fw9/KaGtrBL5AiJd5JWhTh9dRNFYyk1FQUJDMlXbv3o3Y2Fiw2WwsX74cn376qVKCI4SQpubR8zfILijFxB4dldK+o23lIIbH2UXNKhkpdI734MEDjB49GitXrkTv3r1x/PjxBk9EZ8+ehaurq9jngw8+AFD58G14eDi6du2KwMBAbN++XWzd2uoJIUTZLt1+CV1tjXrPRSeLlZke9HQ08CS7ed03kms0XVlZGX788Ufs3r0b9vb2iIuLQ8+ePZUSUHp6Ovz8/LB+/fr/gtSsDHPRokXIycnB7t278fjxYyxYsADW1tYYNmyYXPWEEKJMpeV83EjNQc9OttDR1lDKNtgsFtraGOFxdvOaFqjWZHT69GmsWrUKeXl5+PLLLzFt2jSljpxLT09Hhw4dYGUlPjb/xYsXOHXqFOLj4+Hs7Aw3Nzekp6djx44dGDZsWK31hBCibAkPOODyhQj0tFfqdtrZGuN04vNmNYhB5l5kZ2dj2rRp+Prrr9GqVSscPnwYM2fOVPoQ7vT0dDg6OkqUJycnw9TUFM7OzqIyHx8f3L9/Hzwer9Z6QghRJoFQiFP/PEUbG0O0s1XuvZyqQQwvckuUuh1VknlmNGTIEJSXl8PIyAhWVla1viYiIiKi3sEwDIOsrCzcuHEDu3btQnl5OQIDA/HNN9+Aw+HA2lr8GqyVlRX4fD7y8vJqrbeza/ghloQQUuX6fQ44r8swPaiz3I/F1FVVsnvCKUJbJSc+VZGZjDp16iT6OTc3VyXBvHz5EqWlpWCxWIiIiEBeXh7Wrl2L2bNnw8vLCzo6OmLLV52lcblclJWV1VivCAsLw3rsBWBl1Tz+OFSN+k1x1Gd109D9xhcIcfT6E7RvZYKBPR1RXMZDWTlf6rIaWoCRoa5EuZaWptRyfX0dWJnri5VZWBjCQFcTnMJylf0NKHs7MpPRrl27lLphaRwcHJCQkAATExPRkYWZmRlGjRqFbt26SSSVqt91dXWhq6tbY70i8vOLIRTWbVZcKysj5OY2r1EuqkD9pjjqs7ppyH7jC4EKHh9X775Cdn4ppoxwx5MXhRAyQGIqR+o6ni5WKCoulyjn8fhSy0tLK5ArEEiUt7Y2xIOsfJX8DdTUZ2w2q94H8ICCQ7tVwdTUVOwUt+oeEI/HkzhDy8nJgZaWFszMzGBra1tjPSGENLQKHh/X72fjr8tZsDTRRWk5DzcecMAXCpW+7Xa2xnieWwy+QPnbUoVGlYwuXLgAX19fFBcXi8pSUlLAZrMxYsQI5OfnIysrS1SXmJiITp06QVtbG15eXjXWE0KIMqQ/f4OScj48nS2Vfq+ounZ2RuALmGYziKFRJaMuXbpAT08P8+fPR0ZGBhISErBo0SKMGjUKDg4O6NOnD+bNm4eUlBScPHkScXFxmDBhAgDUWk8IIQ3tTXEFkh7lwtpMD/aW+rWv0ICqBi40l+eNGlUyMjIyQlxcHMrKyvDJJ59gxowZCAgIwJIlSwAAa9euha2tLYKDg7Fq1SrMnDkTQ4YMEa1fWz0hhDQUhmGw5/QjCAQM/NxtVXpWBADWpnrQ19FE1qvmkYzkfp+RqnTo0AFxcXFS60xNTbF582aZ69ZWTwghDeXSnVe4n1UAXzdrmBiq/lYAi8VCewcTZLxoHsmoUZ0ZEUJIU5BXWIY9Zx6hQ2sTuLVtuJfnKcrZwRgv8kpQWt70H+ynZEQIIQrg8QWIOZICFoBxA11VfnmuOmcHEwBAxsumf3ZEyYgQQuQkZBjEHX2A9Bdv8L8hbjA3VuwZxobmaG8MFqtyRF9TR8mIEELk9Of5DPzzIAej32+Pbh1t1B0OdLU10draEOkvKBkRQkizxxcCJ/55iuMJTxHgYYdAL3uUVPBRx4laGpSzgwkyX72FQAUP2ioTJSNCCKnF+aRn2H82HQ5WBmhna4SbqTkqm2mhNs4OJqjgCpr8w6+UjAghpAYXkl/g978fwd5SH7297MFmq2/AgjRVgxia+qU6SkaEECLD+eQX2HEiDR3bmaGPt0OjfJGdhYkuTAy1KRkRQkhzdPKfp9h5Ig0e7S0QMtwdGo0wEQGVD786O5g0+RF1jbN3CSFExfhCoKSCj+JyHvadS8e+s+nw7mCJL4Z1bLSJqIqzgwny3pSjsLhC3aHUWePuYUIIUZEKHh//pGQj+uBdnEx4CudWJnB3MkfSw9xGMVChJqKHX5vwpTpKRoQQAoDHF+Li7VdIe1qI99qZwc/dBmw1zq6giDY2RtDUYONRE75UR8mIENLilVXw8fOhe3iSXYQurlbo6mql1ml+FKWlyUaHVia4m5mv7lDqjJIRIaRFK3hbju9330L6izfw72yLTo7mTSoRVfFytsSr/FJwXpeqO5Q6oWRECGmxMl68wXc7boJTWIYpI9zR/t97L02RZwdLAMDtR3lqjqRuKBkRQpqsqhFw0j78WsYcXL33Ct//ngRtLTYWj++K99qZqyZoJbE21YODpQGS05tmMmp0L9cjhBB5VfD4uPGAI7XOt6MNNHUkv+JKyniIi0/BlXvZcGtjimkjO8NQTwslFXxlh6t0ns6WOJHwFCXlPBjoaqk7HIXQmREhpMV48OQ1vo44h2v3ORjesx1mf+oFQ72m9aVdEy9nSwgZpkkOZKAzI0JIs/citxiHL2fhZlou7C0NsGB8F7S3b7r3h2RxsjeGkb4Wbqfno8d7tuoORyGUjAghzZJQyCD1yWtcuP0S/6RwoKOtgeE92+HzD9xR9LZM3eEpBZvNgkd7CyQ9zANfIGyUc+nJQsmIENJslJbzkFtYjuyCUhy6mIW3pVzoaGlgcI82GNK9LQz1tKCro4kidQeqRF7OlrhyNxvpz9/Ara2ZusORGyUjQkiTwxcI8Sq/FI9evEHSw1wUFleg4G0Fyv4dhKDBZqGzkwV6uNvAs70ldLQ11Byx6rg7mkNTg4WkR3mUjAghpCFVcAVIe1aIlMcFePT8DZ7nFoP379htNosFYwMt2JrrwdJED5amujA30kGPTnYwkDKarrnT1daER3tLXL33CiN7OUJXu2n0QdOIkhDS4hSX8XDrYS5upOYg9clrCIQMNDVYcLI3QR9vB7SzNYKVmT6eZL+V+sI7Fpsldbi2bilXFeGr1ZDubXDrYS4u3n6Fgb6t1R2OXJrO3S0icuhSZqNoozZrdyfWK4ba6mtqvz5t1xZ3XdqctvGCwtua89MVhdqatPasQu1MlLF8yLpzcrctbVlp7cpTNnHtWfD4QtxMzcGPf9zGrMjL+PV4KnJel2KAT2uEf+oFvoDB/LFdMKZfB/Rwt4WthT5+O/VQrJ2dJ9IAABU8AW484ODrHy7ixgOO6PPZkuMyY6ha992fwzdflrtcWhs1/Vv1c/jmy6Lfq3++/uGiRN/Vpr2DCVxam+LkP0/BFzTuGcerUDJqgv668rhRtFGbh89qnkG4thhqq6+p/fq0XVvcdWmznCtQeFsFb6W/m0ZWW4yC7cgiEEq2JKttacvKo2p2hKr3B6U8eQ0AmBV1GT8duoennCIM+PeIfu0UP3zS1xnujk17hgRVG9qjLV4XVSAhRfpDwY0NXaYjhKhcBa/y8tm2+BRkvnyLNyWVl87c2pqhx3s2cG1jBjabhRMJT1HKFUBLUxM8fuU61S+91TEXirVR9XNd22qsOjuZo5WVIY5dfwK/TraN/nUYlIwIISohEArx+FUR7mcV4EZaDgAg6VEerEz10MPdBtfvczB2kCtuP8xF4r/1AHDjAQeeLla4/TBX9HsVTxerOsVSvY2qn+vaVmPFYrEwtEcbxBxJwe30PHh3aNz7R8mIEKIU5Vw+Hr8qQuartwCAsE2XUVbBBwuAo70xACCot5NoOp7r95vG5aSmxLejNQ5ezMSfFzLh1sYMeo14dGHjjYwQ0iSUlvPAeV2G7IJSZOdXvktn/s/XkFtYJnavydfNCu+1M0fHtmZga7Dx9Q8Xm9W8cI2RBpuNCUPc8MO+24j56z6+/thD6sjDxqDZJSMej4c1a9bg6NGjAIDRo0dj9uzZYLNprAYhtak+sIDHF6K0nIeiMh4A4PLdbLwtqcDrospP3pvKKXWm/3hJtH7VbQl7KwP0cLeFk70RHO2MMWPzZfxvSEfRcs1hhuymwr2dOT4b0AG/nXqIAxcy8EkfZ3WHJFWzS0YbN27ElStXEBMTg5KSEsydOxfGxsaYPHmyukMjROkEQiEquAKU//up4FWOvEt6mItyXlU5HxX/jsjbfuwByrkClHH5KK8QoKS8MvHM2nwJfIH4Hf1tR1MAAHo6GjAz0oWliS6ecorxUaAjCosrYGygDSN9Lew+9Qge7S3Qzd0WTLVRAQ0x8IDUTd8urfAyrwQnEp7CylQP73vZN7q32TarZFRRUYE9e/Zg06ZN8PT0BACEh4cjIiICISEhja7zCanCMAy4PCHKuXxRIgGAWw9zUVbBR2kFH6XllR8AiDp4F3whgzfFFSgr56OMW1kesu681PYjD96VWn47Ix+62hr/fjRhZaaPV/mlcGltCm0tDehosaGrrYkLyS+xeIIP7C30xZ7on7j2LPr5tJb6TqEKnkA06ABomIEHpO6C+3cA53UZdp1MQ9LDXAT37wA7CwN1hyXSrJLRgwcPUFZWBh8fH1GZj48PcnNz8fz5c7Ru3TSeRCaVX85A5ZE+w+DfD1N5D6LaUXXFu8/csICqQw4eXyi6bCT6F/8dkAj/3QYLqPVAhWEYCIQMBAIGAqEQPL4Q3KoPT4Dyf4/6r93PRnkFX3S2UVbtbGPDvmRweQJweQLw+Ax4/MrYv4o4Dy5P+oOJUe8kEW2tysvNr/JLYKCnBQM9LVgY66KolIui0jfwdLaAliYbWhpsaGqy8V47c2w9fB/f/s8XOtoa0NGqTDx8hsGMHy/ho0BHsfY9XawQvvkyfNysJWKxtTSAQMhIXGKjs5ymQYPNxoxRHjh76wUOX87E0rh/4N/ZFp2dLODS2hRG+tpqja9ZJSMOhwNDQ0MYGPyX7a2srER18iajutzge5z9Fk+yi6Gvn4+S0n8fMmSq/fPvFx/DAAwYMEzll+F/X7SVv1d9UFUvZCCs/kUMBtZmeth/Ll3hGKuzNtPDgfPp/35zV+5v5Zfyf79XxsuI4hcylQ85CoX/fjELGQgEQgiEDPgCBnyBEHyBEAKhEAJh5Ta+23lTYh+q/rU208OiX67/10//9s2//6usj02ocR++3f6PzLolcTWvuzDmukQ569+6BTHX/g1K9gOf0to8fDnrv7ZYgI6WBlhgwdpMDwK+ENqabOjpaMLCRBdaGmzcTs9Dz8520NJgQ0uLDW0tDejraEJDg4V9px9hygh36GhrQk9bA1pamkjJysfBCxkYEegEQwMdFJdUoKOjOR5kFeBNCRd+nezEYvJ0sYK1WSZsLcWPgDWZynj133kbqKYGW2q5tZkeBEIGD7IKJMrBgtjyVetrarBF5e+2KW071mZ6YutIa0vedmpap/rv1X9+t626/lx9P6qXVa+r/jsAsZ+laahBB9psDQzu3gYBHnY4kfAUSem5SH1aCADo7GiB8YNdFY6hoWJjMVWHoM3AoUOHsG7dOly9elVUJhQK0bFjR2zfvh09e/ZUY3SEEEJkaVZDzHR1dcHlik+CWPW7rq6uOkIihBAih2aVjGxtbVFUVISysv/e4pibW3kD1cbGRl1hEUIIqUWzSkZubm7Q09NDYuJ/MyHfvHkT1tbWcHBwUGNkhBBCatKskpGuri5GjRqFFStW4NatW7h27RoiIiIwYcIEdYdGCCGkBs1qAANQ+azRypUrcfToUejo6GDUqFGYPXs2PWNECCGNWLNLRoQQQpqeZnWZjhBCSNNEyYgQQojaUTIihBCidpSM3sHj8bBixQp0794d3bt3x4YNGyAUSp83rKSkBIsXL4avry969OiB7777Tuyh24yMDEyYMAHe3t7o168ffv31V7H18/LyEBYWBm9vbwQEBCAqKkqZu6ZUquy32uqbiobss7S0NIwfPx7e3t4YOnQojh07JrZ+dnY2Jk+eDG9vb/Tv3x9//fWXUvdNmVTZb7du3UJwcLCo32JiYmRuqzFTZZ9V3+bw4cOxYcMG+YJkiJi1a9cyAwcOZJKTk5krV64w/v7+zNatW6UuO336dGbYsGHMnTt3mISEBCYwMJDZtGkTwzAMU1FRwfTt25dZsmQJ8/jxY+bkyZOMp6cnEx8fzzAMwwgEAubjjz9mxo0bxzx8+JA5ffo006VLF+bAgQMq29eGpKp+q62+KWmoPsvPz2e6devGzJgxg3n48CFz/PhxpmvXrsy5c+dE6wcFBTGhoaHMo0ePmD179jDu7u5MYmKiKnazwamq3169esV06dKF+f7775msrCzm9OnTTLdu3Zhff/1VVbvaYFT5t1YlMjKScXFxYdavXy9XjJSMqikvL2c8PT2Z8+fPi8oOHjzI+Pv7M0KhUGzZ9PR0xsXFhXnw4IGobM+ePcyXX37JMAzD3L9/n3FxcWHy8/NF9dOnT2cWLVrEMAzDnD9/nvHw8GByc3NF9T/88AOzePFipeybMqmy32qrbyoass+2bdvGBAQEMBUVFaL6LVu2MGPGjGEYhmH++ecfxt3dnXnz5o2ofv78+czMmTOVsm/KpMp+2759O3VUiKAAAAlHSURBVDN48GCxNn/66Sdm5MiRDb5fyqTKPquSlpbG+Pv7M0OHDpU7GdFlumpqewVFddeuXYOTkxPc3NxEZWPGjEFsbCwAwNTUFCwWC/v374dAIEB6ejoSExPh7u4OALh+/Tp8fX1haWkpWn/mzJn47rvvlLmLSqHKfqutvqloyD57+vQp3N3doa393ysAOnbsiLt374LP5yMpKQmurq4wNjYW21ZycrKydk9pVNlv/fr1w5o1a8TaZLFYKCoqUsauKY0q+wwABAIBFixYgPDwcJibm8sdJyWjamp7BUV1T58+RZs2bfD7779j0KBB6Nu3LzZs2AAer/LdNfb29pgzZw6io6PRuXNnDBs2DIGBgQgODhat7+DggKioKPTt2xeDBg1CXFyc6D0+TYkq+622+qaiIfvM0tIS2dnZYuu8ePECPB4Pb9++BYfDgbW1+PuJLC0tJbbTFKiy31q3bg0vLy9RXWlpKfbv3w8/Pz9l7Z5SqLLPAGDbtm0wMTHByJEjFYqTklE1ZWVl0NHRESurOgJ4dzbwkpISJCYm4vTp0/j++++xePFiHDlyBBEREQAAPp+PzMxMfPDBB9i/fz/Wr1+Pc+fOiW62l5SU4OjRo3j8+DE2b96M6dOnY8uWLfjtt9+Uv6MNTJX9Vlt9U9GQfTZ48GCkp6fjl19+AZfLxcOHD7F9+3YAlTeRZW1LIBCIjmabClX2W3U8Hg/h4eEoLS1FaGiosnZPKVTZZ1lZWfjll1+wfPlyxQNV9Ppjc1Z1M666srIyxsXFReJm75IlS5hOnTqJXYc/cuQI4+HhwQgEAubAgQNMv379GD6fL6rfu3cv07VrV4bP5zMTJ05kAgMDGR6PJ6r/+eefmUGDBilp75RHlf1WW31T0ZB9xjAMc+zYMaZbt26Mm5sb06tXL2bnzp2Mi4sLU1xczCxfvpyZOnWqWJvnz59n3N3dlbR3yqPKfqve/pQpUxhvb28mOTlZiXunHKrqs6KiIiY4OJjZvn27aN1x48bJfc+oWb3ptb6qv4JCT6/yrYuyXkFhbW0NKysrsevwTk5OKC8vR0FBAe7evQsXFxdoaGiI6jt16oSioiIUFhbC2toaXC4XmpqaYuu/fPlSmbuoFKrst9rqLSwslLmrDaYh+8zS0hJDhgzB4MGDkZOTA0tLS1y+fBkWFhYwMDCAra0t7t4Vf315bm6uxKW7pkCV/QYAxcXFmDJliugMwNPTU0V72nBU1Wdv3rxBYmIiUlJSsGnTJgBAeXk5kpKScO7cORw9erTGOOkyXTWKvIKiS5cu4HA4yM/PF5U9evQIhoaGMDU1hY2NDTIyMsTuAVXVm5ubo0uXLkhPTxc7TX706JHcr0ZvTFTZb7XVNxUN2WcJCQmYNm0aWCwWbGxsoKGhgTNnzqBHjx4AAC8vL6SlpaG4uFi0fmJiIry9vZW8lw1Plf0mEAgwffp0ZGZmYteuXU0yEQGq6zMbGxucOnUKhw8fxqFDh3Do0CF06tQJQUFBiImJqT3Qupz2NWffffcdM2DAACYxMZG5evUq4+/vz8TGxjIMwzCvX79m3r59yzDMf88JTZgwgUlNTWWuXr3K9OrVi1m3bh3DMAzz4sULxsvLi/n222+ZrKws5vz580xAQAATGRnJMAzDFBcXM++//z4ze/ZsJiMjgzl16hTTtWtXZvfu3erZ8XpSVb/VVt+UNFSf5eTkMF5eXszWrVuZZ8+eMb/++ivj4eEhGp4rEAiYESNGMFOmTGHS0tKYvXv3/n979xMKzR/HAfwdrbRtNnvCwd+yJSursOtP1HIQObhaRWojB0siRWmzByJH7MWBwk1JeyAHNA57ILmplYSIHBZtu+r7Ozz9Nvvwe34tuzOzPe9XTWnm+535zLR6NzPLR5SWloqTkxNlTvyH5Lpua2trwmg0ir29PfHw8BBZPv5ZQbKQ65r9LpbHdAyj3wSDQTExMSHMZrOwWCxibm4u8l18u90uxsbGImMfHx/F4OCgKC8vF9XV1WJ2dlaEQqHI9rOzM2G324XZbBY2m00sLi5Gvde4uroSvb29wmQyifr6+siHIxnJed3+b3uyiOc1Ozw8FG1tbaKsrEx0dHQISZKijnV9fS26u7uFyWQSTU1NYnt7W56TTAC5rpvdbhfFxcWflpqaGvlONk7k/Kx9FEsYsYUEEREpju+MiIhIcQwjIiJSHMOIiIgUxzAiIiLFMYyIiEhxDCMiIlIcw4hIRl1dXRgaGvpy2/r6OoxGY9S6ra0tOJ1O+P1+lJaWwul0fjn36OgIRqMRq6urCIfDaG1txcXFRdzrJ0oUhhGRSj0/P2Nubg7Dw8MoLCyEw+GA1+vF8fFx1LhQKASXywWz2YzOzk5oNBoMDAxgcnIyKVtk09+JYUSkUisrK6ioqEBubi4AoK+vD/n5+XC5XFEtDjweD25vb+F2u5GS8utX+t9/ZLm/v69I7USxYhgRqVAwGMTm5iZaWloi69LS0jA1NQW/3x/pIXN9fQ2Px4P+/n4UFRVFxqakpKC5uTnp+jzR34thRCQzIQTe398/LR8fqUmShNfXVzQ0NETNtVqtaG9vx9LSEp6enjAzM4O8vDw4HI5Px7HZbPD5fJF2AURqxn5GRDLzer3wer1/HOPz+VBQUACtVvtp2/j4OA4ODjAwMICzszNsbGxAo9F8GldSUgLgV7uAj3dYRGrEMCKSWV1d3Zffitvd3cXy8jIA4ObmBllZWV/ONxgMGBkZwcTEBLq7u1FWVvblOJ1OB71ej7u7u/gVT5QgDCMimWVkZMBkMn1af35+Hvn55eUFOp3uP/dhtVoBALW1tX88Vnp6OgKBwDcrJZIP3xkRqZBer4/qzPpdgUAAmZmZcaiIKLEYRkQqlJWVhfv7+x/tIxAI4O3tDdnZ2XGqiihxGEZEKmSxWOD3+3/0iO309BSpqamoqqqKY2VEicEwIlIhi8UCrVYLSZK+vQ9JklBZWQm9Xh/HyogSg23HiVRqfn4eFxcXWFpainnu+/s7Ghoa4Ha70djYGP/iiOKMd0ZEKtXT04PT01NcXl7GPHdnZwc5OTkMIkoaDCMilTIYDBgdHcXCwkJM88LhMBYXFzE9PZ2gyojij4/piIhIcbwzIiIixTGMiIhIcQwjIiJSHMOIiIgUxzAiIiLFMYyIiEhx/wARh+exufaHCwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X, y = get_sample(1000, 0, 1)\n",
    "X = np.array(X).reshape(-1, 1)\n",
    "\n",
    "print(y.count(1))\n",
    "probs = np.random.binomial(len(y), y.count(1)/len(y), 1000)/len(y)\n",
    "entropies = [-p_hat*np.log(p_hat) - (1 - p_hat)*np.log(1 - p_hat) for p_hat in probs]\n",
    "print(len(entropies))\n",
    "sns.distplot(entropies, hist = True, rug = True)\n",
    "plt.title(\"Monte Carlo distribution of H(Y)\")\n",
    "plt.xlabel(\"H(Y)\")\n",
    "plt.ylabel(\"Normalized density\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}