PheMedDemo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import phemed as phe\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from importlib import reload\n",
    "import sys\n",
    "import warnings\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import chi2\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm_notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bootstrap as boots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'bootstrap' from '/Users/davidburstein/Documents/GitHub/PheMED/bootstrap.py'>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reload(phe)\n",
    "reload(boots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Generate Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(20)\n",
    "n_snps = 5000\n",
    "gold_means = .03*np.ones(n_snps)\n",
    "diluted_means = .5*gold_means\n",
    "sample_gold = 1000\n",
    "sample_diluted = 1000\n",
    "se_gold = 1/np.sqrt(500)\n",
    "se_diluted = 1/np.sqrt(10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TWAS_N*alpha^2 = SC_RNASEQ\n",
    "#alpha = sqrt(SC_RNASEQ/TWAS_N) = sqrt(1000/20000) = sqrt(1/20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.22360679774997896"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(1/20)\n",
    "#alpha = [1, .223]\n",
    "#alpha = [1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_results = pd.DataFrame(columns = ['SNP','CHR','POS','BETA_Gold','BETA_Diluted','SE_Gold','SE_Diluted'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_results['SNP'] = np.arange(1,n_snps + 1)\n",
    "df_results['CHR'] = 1\n",
    "df_results['POS'] = df_results['SNP']\n",
    "df_results['BETA_Gold'] = np.random.normal(gold_means, se_gold)\n",
    "df_results['BETA_Diluted'] = np.random.normal(diluted_means, se_diluted)\n",
    "df_results['SE_Gold'] = se_gold\n",
    "df_results['SE_Diluted'] = se_diluted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_results['SNP'] = df_results['SNP'].apply(lambda x: 'rs' + str(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SNP</th>\n",
       "      <th>CHR</th>\n",
       "      <th>POS</th>\n",
       "      <th>BETA_Gold</th>\n",
       "      <th>BETA_Diluted</th>\n",
       "      <th>SE_Gold</th>\n",
       "      <th>SE_Diluted</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>rs1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.069529</td>\n",
       "      <td>0.010957</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>rs2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0.038759</td>\n",
       "      <td>0.015121</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>rs3</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>0.045990</td>\n",
       "      <td>-0.000488</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>rs4</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>-0.074794</td>\n",
       "      <td>-0.010719</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>rs5</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>-0.018515</td>\n",
       "      <td>0.011009</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4995</td>\n",
       "      <td>rs4996</td>\n",
       "      <td>1</td>\n",
       "      <td>4996</td>\n",
       "      <td>0.078014</td>\n",
       "      <td>0.004748</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4996</td>\n",
       "      <td>rs4997</td>\n",
       "      <td>1</td>\n",
       "      <td>4997</td>\n",
       "      <td>0.040007</td>\n",
       "      <td>0.025823</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4997</td>\n",
       "      <td>rs4998</td>\n",
       "      <td>1</td>\n",
       "      <td>4998</td>\n",
       "      <td>-0.057730</td>\n",
       "      <td>0.006685</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4998</td>\n",
       "      <td>rs4999</td>\n",
       "      <td>1</td>\n",
       "      <td>4999</td>\n",
       "      <td>-0.003057</td>\n",
       "      <td>0.001067</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4999</td>\n",
       "      <td>rs5000</td>\n",
       "      <td>1</td>\n",
       "      <td>5000</td>\n",
       "      <td>0.038027</td>\n",
       "      <td>0.004365</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5000 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         SNP  CHR   POS  BETA_Gold  BETA_Diluted   SE_Gold  SE_Diluted\n",
       "0        rs1    1     1   0.069529      0.010957  0.044721        0.01\n",
       "1        rs2    1     2   0.038759      0.015121  0.044721        0.01\n",
       "2        rs3    1     3   0.045990     -0.000488  0.044721        0.01\n",
       "3        rs4    1     4  -0.074794     -0.010719  0.044721        0.01\n",
       "4        rs5    1     5  -0.018515      0.011009  0.044721        0.01\n",
       "...      ...  ...   ...        ...           ...       ...         ...\n",
       "4995  rs4996    1  4996   0.078014      0.004748  0.044721        0.01\n",
       "4996  rs4997    1  4997   0.040007      0.025823  0.044721        0.01\n",
       "4997  rs4998    1  4998  -0.057730      0.006685  0.044721        0.01\n",
       "4998  rs4999    1  4999  -0.003057      0.001067  0.044721        0.01\n",
       "4999  rs5000    1  5000   0.038027      0.004365  0.044721        0.01\n",
       "\n",
       "[5000 rows x 7 columns]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_results.to_csv('data/sim_data.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Run PheMed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['BETA_Gold', 'BETA_Diluted']"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_vars = [var for var in df_results.columns if var.startswith('BETA')]\n",
    "beta_vars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['SE_Gold', 'SE_Diluted']"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "se_vars = [var for var in df_results.columns if var.startswith('SE')]\n",
    "se_vars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "betas = df_results[beta_vars]\n",
    "ses = df_results[se_vars]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = minimize(lambda x: phe.nll_data(betas,ses, x), np.ones(2), \n",
    "                        options={'maxiter': 300}, method = 'Nelder-Mead')\n",
    "#code normalizes first entry to be 1\n",
    "optimizer.x[0] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " final_simplex: (array([[1.        , 1.92982354],\n",
       "       [0.56018589, 1.92982354],\n",
       "       [0.56013706, 1.92982354]]), array([2449.79496822, 2449.79496822, 2449.79496822]))\n",
       "           fun: 2449.794968222299\n",
       "       message: 'Optimization terminated successfully.'\n",
       "          nfev: 115\n",
       "           nit: 56\n",
       "        status: 0\n",
       "       success: True\n",
       "             x: array([1.        , 1.92982354])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2449.794968222299"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#negative log likelihood of alpha\n",
    "phe.nll_data(betas, ses, optimizer.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 1.92982354])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optimizer.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Alpha</th>\n",
       "      <th>NLL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>3451.180461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>0.109820</td>\n",
       "      <td>3440.359554</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>0.119639</td>\n",
       "      <td>3429.589099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>0.129459</td>\n",
       "      <td>3418.869390</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>0.139279</td>\n",
       "      <td>3408.200718</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>495</td>\n",
       "      <td>4.960721</td>\n",
       "      <td>3682.197143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>496</td>\n",
       "      <td>4.970541</td>\n",
       "      <td>3687.489877</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>497</td>\n",
       "      <td>4.980361</td>\n",
       "      <td>3692.779979</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>498</td>\n",
       "      <td>4.990180</td>\n",
       "      <td>3698.067399</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>499</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>3703.352086</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>500 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        Alpha          NLL\n",
       "0    0.100000  3451.180461\n",
       "1    0.109820  3440.359554\n",
       "2    0.119639  3429.589099\n",
       "3    0.129459  3418.869390\n",
       "4    0.139279  3408.200718\n",
       "..        ...          ...\n",
       "495  4.960721  3682.197143\n",
       "496  4.970541  3687.489877\n",
       "497  4.980361  3692.779979\n",
       "498  4.990180  3698.067399\n",
       "499  5.000000  3703.352086\n",
       "\n",
       "[500 rows x 2 columns]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_results = phe.nll_search2D(betas, ses)\n",
    "plot_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Log Likelihood')"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEjCAYAAAD31uwUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUZfbA8e9JJyGFEloIvfcSAcWGKIKiYAOs6KLorvVnd+1tV9fee8GCgIKCHWyoq5TQO4QeQgklkNCTnN8f92YdwqSQdpPJ+TzPPMm8t507mcyZe9/3niuqijHGGFMaQV4HYIwxpuqzZGKMMabULJkYY4wpNUsmxhhjSs2SiTHGmFKzZGKMMabULJkYAERkiYic6v7+kIh8VIJ1NBMRFZEQ9/m3IjLS/f1KEfm9TIMuOA4VkValXMc6ETm9rGIqS76xicg/ReTtctjGpSIytazX6667rYjME5FMEbmpPLbhZ5snicgKn+cl/vuKSBMRyRKR4LKLsOqzZFJJVMSHV2FJQlU7quovZbk9VR2kqmPKcp2BSkROdZPgnceynKr+S1WvLuW2j/gS4K73Y1UdUJr1FuJO4BdVjVbVF0u7Mvd9fdhNTpkislJEXhaRhnnzqOpvqtq2hOs/4n9TVTeoak1VzSlt7IHEkokxlcNIYKf7M9A1BZaUZEHfhJfPeFWNBmoD5wENgDm+CcWUL0smVYCIXCMiKSKyU0SmiEgjn2kDRGSFiOwWkVdFZLqIHPM31YKOjEQkVEQ+EZGJIhImIkEicreIrBaRHSIyQURqF7DOX/LHIiJPi8guEVkrIoN82hu5+7bT3ddrfKaFi8jzIpLmPp4XkXCf6XeIyGZ32t8K2cd+IrLI5/kPIjLL5/nvIjLUZ5FuIrLQfW3Hi0iEz7yDRWS+iGSIyB8i0iXfa3l7Qcv6iSsSuBC4HmgtIkn5pl8uIuvd1/vefNP+d7TpHt2k5pvue0qsl4gki8geEdkqIs+6s/3q/sxwT98cL/lOS7pHLteJyCr37/eKiIg7LVhEnhGR7e7f9Yb8Rzo+6/kJ6Ae87G6rjYjEisgHIpLu7ud9IhLkzn+liPxXRJ4TkZ3AQwW9jgCqelhVlwDDgXTgtoJeG5+Y3heRx3ye/29eEfkQaAJ86cZ7pxx9Orew9+5D7v/IB+IcNS3J//cNFJZMKjkROQ34NzAMaAisB8a50+oCnwH3AHWAFcAJZbjtGsAXwEFgmKoeAm4ChgKnAI2AXcArxVxlbzfGusB/gHfyPpCAT4BUd50XAv8Skf7utHuBPkA3oCvQC7jPjXEgcDtwBtAaKOxU4Z9AKxGp634QdAIai0i0u689gd985h8GDASaA12AK91t9gDeBa7Fed3fAKb4JriCli3ABUAW8CnwPXBF3gQR6QC8BlzuvjZ1gMaFrKswLwAvqGoM0BKY4Laf7P6Mc0/f/FnA8oOB43D+BsOAM932a4BBOH+fHjjvD79U9TSc1/gGd1srgZeAWKAFzvvqCuAqn8V6A2uAesDjxdlR9xTUZOCk4sxfyHouBzYA57jx/sfPbIW9dwHOxfmfjQOmAC+XJqbKypJJ5Xcp8K6qzlXVgziJ43gRaQacBSxR1Umqmg28CGwpo+3GAN8Bq4GrfM4PXwvcq6qpbjwPARf6+xbqx3pVfctd1xic5FhfRBKBE4G7VPWAqs4H3sb5AAXnNXhEVbepajrwsM+0YcB7qrpYVfdSyDdXVT0AJON8eCYBC4Hfgb44yWqVqu7wWeRFVU1T1Z3AlzgfluB8eL6hqjNVNcftFzrorqOoZf0ZiXOaJgcYC1wsIqHutAuBr1T1V/f1vh/ILWRdhTmMm0xVNUtVZxzj8k+oaoaqbgB+5q99GoaTpFJVdRfwRHFXKE4n9nDgHlXNVNV1wDP89fcFSFPVl1Q1W1X3H0O8aTinvcpNMd67AL+r6jfu3/dDnGQccCyZVH6NcI5GAFDVLGAHkOBO2+gzTXG+IZWFPjjfqJ/QI6uBNgU+d0/vZADLgBygfjHW+b9Ep6r73F9r4uzHTlXN9Jl3Pc4+Qr7XwP29kc+0jfmmFWY6cCpOQpkO/ILzbfgU97nfeIF9bqzgvAa35b0G7uuQ6BNTYcsewf0w6gd87DZNBiKAs93n+f/Ge3H+/iUxCmgDLBeR2SIy+BiXL2if8v8NfH8vSl0gjKP/vgk+z49lfb4ScPqhylNR7104+nWLKOaXryrFkknll4bz4QWAiEThnOrYBGzG55SHe8qopKdA8puKc3rtRxHxTRQbgUGqGufziFDVTaXYVhpQW0Sifdqa4Oxj3vSm+aalub9vxvkg951WmPzJZDoFJ5OCbAQez/caRKrqJ8Vc3tflOP+HX4rIFpzTORH8darriP1z+1fqFLCuvUCkz7zBQHzec1VdpaoX45wuehL4zH0/lbZ0+BHvQ478exRlO84RU/6/r+/76Zjjc/tczuHI05YFOeJ1w+m891XY9ot671Yblkwql1ARifB5hOCc9rhKRLq55+T/Bcx0Twd8DXQWkaHuvNdz9D9CfkH5thFe0Izu+eGxOAmlrtv8OvC4iDQFEJF4ERlSmp1W1Y3AH8C/3Zi64HyLzvu2/glwn7utusADQN4Q5wnAlSLSwf2gfbCIzf0BtMXpd5nldtY2xTkv/2thC/p4C7hORHqLI0pEzs73gVJcV+Cctuvm87gAOFtE6uD0iQ0WkRNFJAx4hIL/b1fifOs92z1Ndh/gO1DhMhGJV9VcIMNtzsHpqM7F6bMoiQnAzSKSICJxwF3FXdA99TMB5z0V7b6vbuWvv+8xEWfASHuc90wD4NkiFgGYD5wlIrVFpAFwS77pWyngtSnGe7fasGRSuXwD7Pd5PKSqP+KcJ5+I8w2wJTACQFW3AxfhdGbvADrg9AkcLGQbF+fbxurCAlLVR3E64X8QZ9TWCzidiFNFJBOYgfNBXFoXA81wvul9DjyoqtPcaY/h7NdCYBEw121DVb8Fngd+AlLcn4Xtz153+SXugAJwOubXq+q24gSqqsk4/SYv4wxASKHwDna/RKQPzj6/oqpbfB5T3HVe7Ca763GS+mZ3e35PZarqbuAfOOfsN+F84/addyCwRESycP6OI9zz/PtwOrb/656268OxeQvnSHYhMA/nfZyNk6iK40Y31jU4fVhjcQY4HIvh7n5l4Lw/dwA9VTWt8MUApx9jAbAOZz/G55v+b5wvMxkicruf5Qt771YbYjfHChzuoX0qcKmq/ux1PKb8icgjQGNVLXBIdEUTZ8j366ratMiZTcCwI5MqTkTOFJE493TVPwHBOVowAc7tI+sArPU4jhoicpaIhIhIAs6pxs+9jMlUPEsmVd/xOKeqtuN0OA49xuGTpuqai9Px/ZbHcQhOv88unNNcy3D6tUw1YsmkilPVh1S1jlvnqLeqzixqGRGpLyK/ulfkPlMWcYjHxe8kX2FCKUWxR8lXFLAsich5IrLRfa26l2ZdqtpdVfuoaplcWyQlLBqpqvtU9Tj3PVhPVa9S1T351n1UNYTyIiUsVGpKx5JJgJBjKxQ5GudIJkZVbyuL7Ws5Fr8Tp9zFIfmrkN9iEfm3iMT6bL/EhQnzJx4tRVHAYniav67+nldO2yiQiAwRpwzMHnHKn/wozgWwZVI0siwc43v5WNd9qojkusk8S0RSxSl3ctwxrMOSlR+WTKqnpsBSLcHoC/Ffb6kiLsD6jzqF/OJxSm30wRl9FFUB2y5LBRY5LO/X0U2YH+DUq4rFKfXyKiW/or6qSlPVmkA0zvtoOfCbHFkCxRwjSyYBSNwifeKnqKKIvI9TvuNO95vZ6VJI8Ub5q6jdKBHZAPxURFte8btfRORRcYr0ZYrIVPnrWhVE5Ar5q3jh/cX9NuoOZZ2NU++oDm4NJynkfin5T7H4zisiedeWLHBfj+GSryigiLR315EhTqG+c32mvS9O0cOv3f2cKSIt/cQQLs7Q1WB3W6vd9nUicpeILAT2up3YRW3vVXHuFZPlvr4NxCl+uUtElhdy+qwbsFZVf1RHpqpOdMuj5C8amff3vEqc03K7xCn0eJw4BSwzRORln7iO+Lae//2Q77VoKSI/uX/77SLysTjXp/gtrOi29xGnoGaGiCwQ99477rTm4hQ4zRSRaThX1RfJfQ1SVfUBnOHUT/qs8wV3v/eIyBwROcltH4gz0GW4G98Ct/0qEVnmxrBGRK4tTgyBxJJJ4PJbVFFVr8S5oOo/7qmWHyhe8cZTgPb8VdyvoDZfl+B82NfDKZlxO/yveOGrODW3GuJ8S04oYB1+ueUrplH6Qn55RQ67uq/HEdcYiHPx35c41x/Uw7km4mMR8T0NdjFOB3QtnOtDjipGqKoH3W/DedtqmW/5s3EKAUoxtjcM54LEujjXFP2J0xmfV/izoAv15gLtxKnA209E/JZ4yac3TgHN4TjX89yLU0yzIzBMRE4pxjryE5xrNxrhvH8ScWuqqZ/CiuKMEPsa59qi2jjvo4kiknd1/1hgDs7+P0rJyvhPAnrIX0e6s3GSb213/Z+KSISqfodz4fB4N768OlvbcAphxuC8558TpyBotWHJJHD5LapYwLzFKd74kKruzTdSzF+br/dUdaU7fQJ/FQa8EPhSVX93Lxx8gJKV9Cj3Qn44p0Fq4tQoO6SqPwFf4SSAPJNUdZY6xTY/pvCijv68qKob3depONv7XFXnqFO48nPggKp+4P6txwN+j0xUdQ1OKZkEnL/HdvdIp7Ck8qh7NDgV58LCT9QpuLkJp1TJMQ8iUNUUVZ3mJth0nORXWFK6DPhGnWKJue4Fgck4V603walkfL+7vl9xkvGxSsNJcnFujB+p6g51iks+g1NJoMB+NFX9WlVXu0c703G+DJTqi05VY8kkcBVUVNGf4hRv9Fdsr6gCfMUqDOjGV5LihRVVyG+jOiVI8hRVyK843/h9+b6OxdneVp/f9/t5XuD2VXWGqg5T1XicD7uTcY42ClLibRVEROqJyDgR2SQie3BKpxR2aqopcJEcWVjzRJwvSI2AXepUNshTVLFPfxJwvtBkuDHe5p622u1uL7awGEVkkIjMEOeeJhk4Fb2LdbotUFgyMVC84o3+jhxKWj4hf4HKGhRcvNAv99v06ZRNIb/CpAGJ4t6syVXWhfx8X8eK2J6zUafvaRLOfV1K61he43/j7HMXde6tchnOUcH/Qss3/0bgw3zvzyhVfQLnvVRLjhyIUVSxT3/OA+aq6l63f+QunNOJtVQ1DtjtE+MR8YlzwfBEnJF69d35v8m3TwHPkomBcijeWITPgHNE5ARxihc+TDH/8dyO7J449cJ2Ae8VY7H5wPkiEinOiKZR+aYXWMgPmInzQXmnOEUET8W5OHRcceItgXLbnjjFIq8RkXru83Y4AxnKomLCfOBkca43isW5705BonFuBpbh9ofckW96/r/HRzjvlzPFuatjhDiDJBqr6nqcU14Pi3Mn0BNxXq8iiSNBRB4ErsbpWM+LLxunAGaIiDyA0xfiG18zn4QfhnMaLB3IFmewS4mGqVdllkwMlF/xRr/UKV54I84H5GYgE6cDs7AClXe6se3EGd46Bzgh3+mNgjwHHML5EBjD0RVdHwLGuKdQhuWL9RDOB+4gnGtzXgWuUNXlxdjuMSvn7WW4614kzsiy73D6XPzdPfCYuP0Y43GKPc7B6ecpyMM4d2TcjdOxPinf9CMKK6pTmXcIzod9Os6Ryh389fl1Cc77dSdOKZcPigi3kbv/WTgd7Z2BU91+IXDudvktThXm9cABjjwV+an7c4eIzHUHg9yE0w+1y41nShExBBwr9Gg8556yygBaq6qndaaMMSVjRybGEyJyjnvaKQrnXPMinBLgxpgqyJKJ8coQnM7mNJzrGEaU5Ip8Y0zlYKe5jDHGlJonRyYi8pQ4ZR8WisjneaUUfKbnVaC93adtoIisEJEUEbnbp725OCUsVonIeHd0kDHGmArkyZGJiAwAflLVbBF5EkBV7/KZPhGn+NxMVX1anLLmK4EzcO4kOBvnlqZLRWQCzhXI40TkdWCBqr5WVAx169bVZs2alfm+VXU7djjXDtapc0yXfRhjqok5c+Zsdy96PUJFVHs9is8QPHCGoV6Y90REhuLcC9p3yGcvIMUtB4GIjAOGiMgy4DScoXjgDPt8CCgymTRr1ozk5ORS7EVgev/99wG48sorPY3DGFM5iYjfCgOVoQP+bzhjunFH9tyFMw7dVwJHjvNOddvqABluTSTfdr9EZLSIJItIcnp6ehmFb4wxptySiYj8IM5NjPI/hvjMcy/OlaZ5F5E9DDynqln5V+dnE1pIu1+q+qaqJqlqUnz8UUdpxhhjSqjcTnOpaqH3phCRkTglm/v7DAntjVOt9j841TtzReQAzhW1iT6LN8YZUrodiBOREPfoJK/dGGNMBfKkz0ScG8zcBZziU9EWVT3JZ56HgCxVfdkthd5aRJrjFLwbAVyiqioiP+P0uYzDuY/B5IrbE2OMMeBdn8nLOMXUpolzP+rXC5vZPeq4AadmzjJgglvfCZykdKuIpOD0obxTfmEbY4zxx6vRXK2KMc9D+Z5/g1PWOf98a3BGexljjPFIZRjNZYwxporz5MjEmIqgqmzPOsSGnfvYkXWQjH2Hydh/iMM5Sk6ukqtKRGgwNcNDqBkeQr3ocBrF1aBBbAQRocFeh29MlWLJxASEnFxl2eY9LEjNYOHG3SxO28267XvZeyinROurFx1Ou4YxdGwUQ4eGMXRLjCOxdmTRCxpTTVkyMVXWzr2H+GHZVn5dmc7vKdvJ2HcYgLjIUDonxHJcs9o0qxNJkzqRxNeMIC4ylLjIUMJCgggWIUiE/Ydz2Hswm8yD2Wzdc4C0jAOkZexn/Y59LN28h7d+XUN2rjNyvXGtGpzQsg4ntKzLyW3iqR1lZeCMyWPJxFQp+w5l8+2iLXy5MI3fV20nO1epFx1O/3b1Oal1XXo0qUVi7RqIFO/221HhIUSFh1APaBlf86jpB7NzWLU1i7kbdvFHyg6+X7KVCcmpBAn0bl6HQZ0bMLBjA+rFRJTxnhpTtVgyMVXCmvQsPpyxns/mpJJ5IJuEuBpcfVILBndpSMdGMcVOHscqPCSYTgmxdEqI5Yrjm5GTqyxJ283UJVv5dvFmHpi8hIemLOGk1vGMOC6R/u3rExZi41pM9WPJxFRqC1MzeOmnFKYt3UposDCoU0Mu69OU45rVKrcEUpjgIKFL4zi6NI7j9jPbsmprJlMWpPFpcip//3gudaLCuLBnY0ae0IxGcTUqPD5jvGLJxFRKi1J38/TUFUxfmU5sjVBu7t+ay/o0JT463OvQjtC6fjS3DWjLLae34ddV6YyftZG3f1/LO7+v5ewuDbnmpBZ0Soj1Okxjyp0lE1OppGXs5+nvVzBp3iZqR4Vx58C2XN6nKdERoV6HVqjgIKFf23r0a1uP1F37eO+/6xg3awOT56dxcpt4bj2jDd0S44pekTFVlCUTUykcys7lzV9X8/LPKeQq/P3Ulvz91JbEVPIk4k/jWpHcP7gDN/Vvzccz1/PWr2sY+sp/6d+uHv93Rhs7UjEByZKJ8dyc9bu4Z9JCVm7N4qzODbhnUPuAuKYjtkYo/zi1FVcc34wxf6zjzV/XMPil3zm/ewJ3DmxHg1gbAWYChyUT45kDh3P4z3creO+PtTSMieCdkUn0b1/f67DKXM3wEK7v14rLj2/Ka7+s5p3f1/Lt4i1cd0pLRp/cghphdrW9qfosmRhPrNiSyc3j5rF8SyaX92nKXYPaUTM8sN+OMRGh3DWwHZf0asK/v13Gcz+sZNzsDTwwuAMDOzXwZHSaMWXFBsSbCqWqfDxzPee8/Dvbsw7y3pXH8ejQTgGfSHwl1o7k1Ut7Mn50H+Iiw/j7x3MZ/eEcNu/e73VoxpSYJRNTYQ4czuGuiQu59/PF9GlRh29vPpl+7ep5HZZnereow5Qb+nLPoHb8tiqdM579lTF/rCMnt8A7TxtTaVkyMRUiLWM/w9/4kwnJqdx4Wiveu/K4SnfNiBdCg4O49pSWTL3lFLo3iePBKUsY/safbNixr+iFjalELJmYcjd3wy7Oeel3Vqfv5Y3Le3LbgLYEB1n/gK8mdSL54G+9eOairqzYmsnAF35l3KwNqNpRiqkaLJmYcjVt6VYueWsGNSNC+OL6vpzZsYHXIVVaIsIFPRvz/S0n0y0xjrsnLeKaD5JJzzzodWjGFMmSiSk3H81Yz7UfJtO2fjQT/34CreodXZXXHK1RXA0+GtWb+wd34NdV2xn0wq/8vmq712EZUyhLJqbMqSrPTF3BfV8s5tS29fhkdB/q1rT+kWMRFCSMOrE5X914IrUiw7j83Zk8N22ldc6bSsuSiSlTqsrjXy/jpZ9SGJ6UyJuX9yQyrPoM+y1rbepHM/mGvpzXPYEXflzFFe/OtNNeplKyZGLKjKry8JdLefv3tVx5QjOeuKAzIcH2FiutyLAQnrmoK/+5oAvJ63Zx1ou/MWPNDq/DMuYI9p9uykRurnL/5MW8/8c6Rp3YnAfP6WBXdJchEWHYcYl8cX1fosNDuPTtmXzw5zob7WUqDU+SiYg8JSLLRWShiHwuInE+07qIyJ8iskREFolIhNve032eIiIvivtJJSK1RWSaiKxyf9byYp+qM1XlwSlL+GjGBq49pQX3nd3eEkk5ad8whsk39KVf23gemLyEf36+iEPZuV6HZYxnRybTgE6q2gVYCdwDICIhwEfAdaraETgVOOwu8xowGmjtPga67XcDP6pqa+BH97mpQM9MXcmHM9Yz+uQW3D2wnSWSchYdEcqblydxfb+WfDJrI5e8NcP6UYznPEkmqjpVVbPdpzOAxu7vA4CFqrrAnW+HquaISEMgRlX/VOe4/gNgqLvMEGCM+/sYn3ZTAd7+bQ0v/5zCiOMSuWeQJZKKEhQk3HFmO166uDuL03Yz5OXfWbxpt9dhmWqsMvSZ/A341v29DaAi8r2IzBWRO932BCDVZ5lUtw2gvqpuBnB/FljsSURGi0iyiCSnp6eX6U5UR58mb+Sxr5dxdueGPH5eZ0skHjinayM+u+4EAIa98Sc/L9/mcUSmuiq3ZCIiP4jIYj+PIT7z3AtkAx+7TSHAicCl7s/zRKQ/4O9T6ph7HlX1TVVNUtWk+Pj4Y94n85efl2/jrokLOal1XZ4d3tXKo3ioU0IsX9zQlxbxUYwaM5uPZqz3OiRTDZXbBQCqenph00VkJDAY6K9/DUlJBaar6nZ3nm+AHjj9KI19Fm8MpLm/bxWRhqq62T0dZl/NytnStD3cMHYuHRrF8PplPQkPsZs7ea1edATjRx/PjZ/M474vFpO6az93ntmWIEvypoJ4NZprIHAXcK6q+pZH/R7oIiKRbmf8KcBS9/RVpoj0cUdxXQFMdpeZAox0fx/p027KwbY9Bxg1ZjbREaG8M/I4oqrRfUgqu6jwEN68vCeX9m7C69NXc/P4+RzMzvE6LFNNePVJ8DIQDkxzz7PPUNXrVHWXiDwLzMY5jfWNqn7tLvN34H2gBk4fS14/yxPABBEZBWwALqqwvahm9h3KZtSYZHbvP8yn1x1P/Ri7h3llExIcxGNDO5FYO5Invl1OeuYB3h55XLW6+ZjxhifvMFVtVci0j3BOa+VvTwY6+WnfAfQv0wDNUXJzlVvHL2BJ2m7euiKJjo1ivQ7JFEBEuO6UltSPCef2TxdyyVszeP+qXtSOCvM6NBPAKsNoLlMFvDZ9Nd8t2cI/z2pP//b1vQ7HFMN53RvzxmU9WbElk2Fv/Gm3BTblypKJKdL0lek8PXUFQ7o1YtSJzb0OxxyD0zvUZ8zferFl9wEufO1P1m7f63VIJkBZMjGF2rhzHzePm0fb+tH8+3y7lqQq6tOiDuNG92H/4Rwuev1Plqbt8TokE4AsmZgC7T+Uw7UfziE3V3nDSslXaZ0SYplw7fGEBQvD3/yT+RszvA7JBBhLJqZAD01ZwtLNe3h+RDea1onyOhxTSq3q1eTTv5/g3Gzr7ZnM3bDL65BMALFkYvyasiCN8ckbub5fS05rZx3ugSIhrgbjr+1DnZphXPHOLOas3+l1SCZAWDIxRzl4OJd/TlpEz6a1uOX0Nl6HY8pYw9gajL/2eOpFh3PFO7OYtdYSiik9SybmCKqwalsmQQIvjOhGqN0pMSDVj4lg3Og+NIiN4Mr3ZtmdG02p2SeFOcLGXfvIOpjNkxd0oXGtSK/DMeWoXkwE40YfT0JcDa58bxZ/rN7udUimCrNkYv7nvynbScvYT/2YCAZ1buh1OKYCxEeH88noPjSpHcnVY5JJXmenvEzJWDIxAOw5cJg7Pl1AjdBgmtaxI5LqpG7NcD66ujcNYiK46r3ZLEy1YcPm2FkyMQA88uVStuw5QMt6NQmyCxOrnXrREXx8TW/iokK5/J1ZLNtsFzaaY2PJxDBt6VY+m5PKP05tZdVlq7GGsTUYe3UfaoQGc9nbM0nZluV1SKYKsWRSze3IOsg9kxbSoWEMN/Vv7XU4xmOJtSMZe01vRIRL357B+h1Wy8sUjyWTau7+yYvZvf8wzw7vSliIvR0MtIivycdX9+ZQdi6XvDWTTRlWbdgUzT49qrHvl2zhm0VbuOX0NrRrEON1OKYSadsgmg9H9WbPgcNc/vZMdmQd9DokU8lZMqmm9hw4zAOTF9O+YQyjT27hdTimEuqUEMt7Vx7Hpoz9XPnebLIOZnsdkqnELJlUU//5bjnpmQd54vzOdpW7KVBSs9q8dlkPlm7ew+gPku2e8qZA9ilSDSWv28lHMzZw5QnN6ZoY53U4ppI7rV19nrqwC3+s3sH/jZ9PTq56HZKphGwcaDVzMDuHuyctIiGuBrcNsCKOpnjO79GYnXsP8djXy4iLXMzjQzvZjdLMESyZVDOv/bKalG1ZvHfVcUTZNSXmGFx9Ugt27D3Ea7+spm5UGLcOaOt1SKYSsU+TamR1ehav/ryac7s2ol/bel6HY6qgO89sy86sQ7z4Uwq1o8K4sm9zr0MylYQnfSYi8pSILBeRhSLyuYjEue2hIjJGRBaJyDIRucdnmYEiskJEUkTkbp/25iIyU0RWich4EQnzYp8qO1XloSlLCELIQ28AACAASURBVA8J4v7BHbwOx1RRIsLj53ViQIf6PPzVUr5bvNnrkEwl4VUH/DSgk6p2AVYCeUnjIiBcVTsDPYFrRaSZiAQDrwCDgA7AxSKS94n4JPCcqrYGdgGjKnA/qoypS7fy26rt/N8ZbYiPDvc6HFOFhQQH8eLF3emWGMfN4+bb3RoN4FEyUdWpqpo3aH0G0DhvEhAlIiFADeAQsAfoBaSo6hpVPQSMA4aI0wN4GvCZu/wYYGgF7UaVceBwDo98uZS29aO54vimXodjAkBEaDBvX5FEw9gIrh6TzNrtVnaluqsMQ4P/Bnzr/v4ZsBfYDGwAnlbVnUACsNFnmVS3rQ6Q4ZOY8tr9EpHRIpIsIsnp6elluxeV2Gu/rGZTxn4eOrcjIXZNiSkjdWqG8/5VvRARrnxvFtvtKvlqrdw+WUTkBxFZ7OcxxGeee4Fs4GO3qReQAzQCmgO3iUgLwN8YRC2k3S9VfVNVk1Q1KT4+voR7VrVs2LGP16av5pyujTi+ZR2vwzEBplndKN4emcTWPQe4ekwy+w/ZRY3VVbklE1U9XVU7+XlMBhCRkcBg4FJVzUsAlwDfqephVd0G/BdIwjniSPRZfWMgDdgOxLmnxXzbjevRr5cSEiT886x2XodiAlSPJrV4YUR3FqRmcNO4eXZRYzXl1WiugcBdwLmqus9n0gbgNHFEAX2A5cBsoLU7cisMGAFMcZPQz8CF7vIjgckVtR+V3S8rtjFt6VZuPK01DWNreB2OCWBndmzAg4M7MG3pVh75cgl/fT801YVXJ9BfBqKBaSIyX0Red9tfAWoCi3ESyHuqutDtE7kB+B5YBkxQ1SXuMncBt4pICk4fyjsVuB+VVnZOLo9/vYymdSL524nNvA7HVANX9m3ONSc1Z8yf63nn97Veh2MqmCcXLapqqwLas3CGB/ub9g3wjZ/2NTh9LcbHhORUVm3L4rVLexAeEux1OKaauGdQe1J37efxb5bRvG4U/dvX9zokU0FsaE8AyjqYzbPTVpLUtBYDOzXwOhxTjQQFCc8O60anRrHc9Mk8u5d8NWLJJAC9MX0127MOcu/Z7a0Yn6lwNcKCeXtkEtERoVw9JpltmQe8DslUAEsmAWbz7v289dsazunaiO5Nankdjqmm6sdE8PbIJHbuPcToD+Zw4LANGQ50lkwCzNPfryQ31ynIZ4yXOiXE8tzwbszfmMEdny20EV4BzpJJAFm8aTeT5qVyVd9mJNaO9DocYxjYqQF3DWzHlwvSeOHHVV6HY8pRgaO5RGQRhV9N3qVcIjIl9uR3y4mtEco/+vkdLGeMJ647pQUp27J4/odVtIivybldG3kdkikHhQ0NHuz+vN79+aH781Jg39GzGy/9sXo7v63azn1ntye2RqjX4RjzPyLCv87vxMad+7j90wU0rlWDHtafF3AKPM2lqutVdT3QV1XvVNVF7uNu4MyKC9EURVV5+vsVNIiJ4LI+VhXYVD7hIcG8fnlPGsREMPqDOWzevd/rkEwZK06fSZSInJj3REROAKLKLyRzrH5avo25GzK4qX9rIkLtAkVTOdWOCuPtkUnsP5TNtR/aCK9AU5xkMgp4RUTWicg64FWcsvGmEsjNVZ76fgVN60RyUVLjohcwxkNt6kfz3PBuLEzdzT2TFtkIrwBSZDJR1Tmq2hXoAnRV1W6qOrf8QzPF8dWizSzfksmtZ7Qh1O5VYqqAAR0bcOsZbfh83iar4RVAivz0EZFYEXkW+An4UUSeEZHY8g/NFOVwTi7PTl1BuwbRnNPFRsiYquOGfq0Y1KkB//pmGb+urD43qgtkxfkq+y6QCQxzH3uA98ozKFM8E+eksm7HPm4b0JagICubYqqOoCDh6Yu60qZ+NDeMncs6u+1vlVecZNJSVR9077++RlUfBlqUd2CmcAezc3jhx1V0S4zj9Pb1vA7HmGMWFR7CW1ckERQkXPNBMlkHs4teyFRaxUkm+/ON5uoL2Lg+j32anMrm3Qe49Yw2VszRVFmJtSN59ZIerNm+l/8bP59cu0tjlVWcZPJ3/hrNtR7nxlbXlm9YpjCHsnN57ZfVdG8Sx0mt63odjjGlckKrutx/dnumLd3K81Zypcoq8uZYqjof6CoiMe5zu0GBxz6bk8qmjP08fl4nOyoxAWHkCc1YunkPL/64ivYNohnUuaHXIZljdKyjuX6y0VzeOpyTyys/p9A1MY5T2sR7HY4xZUJEeHRoJ7o3ieO2Txewamum1yGZY2SjuaqYSXOdo5Kb+7eyoxITUMJDgnnt0p5EhgVz7YdzyDxw2OuQzDGw0VxVyOGcXF7+OYUujWPp19ZGcJnA0yA2gpcv6cF6tyikXSFfddhorirki3mb2LhzPzed1tqOSkzA6tOiDvcMasf3S7by2vTVXodjiqnIDnjgOuADt59EgJ3AleUZlDlatntU0rFRDP3tuhIT4Ead2Jz5GzN4+vsVdE6I5aTW1j9Y2RWnNtcCn9pcnVW1u6ouKP/QjK+vF21m/Y593NTfjkpM4BMRnrygC63q1eSmT+aRustuoVTZFWc0V7iIXALcANwiIg+IyAOl3bCIPCoiC0VkvohMFZFGbruIyIsikuJO7+GzzEgRWeU+Rvq09xSRRe4yL0qAfdqqKq/9sprW9WpyRvv6XodjTIWICg/hjcuTyM5R/v7RXCtZX8kVp89kMjAEyAb2+jxK6ylV7aKq3YCvgLwENQho7T5GA68BiEht4EGgN9ALeFBE8m7X9po7b95yA8sgvkrjl5XpLN+SybWntLQaXKZaaV43iueGd2PRpt3c/8Vi65CvxIrTZ9JYVcv8wznfxY9R/HW/+SHAB+q8a2aISJyINAROBaap6k4AEZkGDBSRX4AYVf3Tbf8AGAp8W9Yxe+W1X1bTKDbC7p1tqqXTO9TnxtNa8dJPKXRrEselve1uopVRcY5M/hCRzuWxcRF5XEQ24txXPu/IJAHY6DNbqttWWHuqn3Z/2xstIskikpyeXjXKXs9Zv4tZa3cy6qQWhIXY/UpM9XTL6W04uU08D01ZwrwNu7wOx/hR4KeT2wexEDgRmCsiK9w+jLz2IonIDyKy2M9jCICq3quqicDHOH0y4IwYy09L0H50o+qbqpqkqknx8VVjdMjr01cTFxnKiOMSvQ7FGM8EBwkvjuhG/ZgI/v7RXLZnHfQ6JJNPYae5Bpd25ap6ejFnHQt8jdMnkgr4fnI2BtLc9lPztf/itjf2M3+Vt2prJtOWbuXm/q2JCi/OGUljAldcZBivX9aTC177gxvGzuWjUb0JsbuLVhqF/SV2qep6nFIq/h6lIiKtfZ6eCyx3f58CXOGO6uoD7FbVzcD3wAARqeV2vA8AvnenZYpIH3cU1xU4gwaqvNenryEiNIiRJzTzOhRjKoVOCbE8NrQTM9bs5LkfVnodjvFR2NfdsThHJ3M4+nSSUvqSKk+ISFsgF1iPc3EkwDfAWUAKsA+4CkBVd4rIo8Bsd75H8jrjccrkvw/UwOl4r/Kd72kZ+5k8fxOX9WlK7agwr8MxptK4KCmR5HW7eOXn1fRsWovT2tlw+cqgwGSiqoPdn83LY8OqekEB7QpcX8C0d3EKT+ZvTwY6lWmAHnv7t7UocPVJ5fLyG1OlPTykIws37eb/xi/g65tOpHGtSK9DqvYK64DvUdijIoOsbjL2HeKTWRsY0rWR/ZMY40dEaDCvXtqD3Fzl+rHzOJSd63VI1V5hp7meKWSaAqeVcSzGNXbWBvYfzuGak604szEFaV43iqcu6sJ1H83lX98s46FzO3odUrVW2GmufhUZiHEczsnlgz/W07dVHdo3jPE6HGMqtYGdGvK3vs15979r6dm0FufYhb2eKU5trkgRuU9E3nSftxaRUg8bNv59s2gzW/YcYNSJ1ldiTHHcPagdPZrEcffEhaxOz/I6nGqrOIO03wMOASe4z1OBx8otompMVXnn97W0iI/i1DZWZt6Y4ggLCeLlS3oQFhLEPz6ay/5DVhDSC8W90+J/gMMAqrof/1edm1JKXr+Lham7+Vvf5lbQ0Zhj0CiuBs+P6M7KbZncZwUhPVGcZHJIRGrgligRkZaA1TIoB+/8tpa4yFAu6NG46JmNMUc4pU08N57WmolzU5mQvLHoBUyZKk4yeRD4DkgUkY+BH4E7yzWqamjDjn1MXbqFS3o1oUZYsNfhGFMl3dy/NSe2qssDk5ewJG231+FUK8W50+I04HycW/V+AiSp6i/lG1b1894fawkS4Yrjm3kdijFVVnCQ8PyIbsRFhvKPj+ey58Bhr0OqNoozmusRVd2hql+r6lfATvcIxZSRPQcOM2H2RgZ3aUiD2AivwzGmSqtbM5yXL+lB6q793PnpQus/qSDFOc3VRETuAecWvsAXwKpyjaqamTB7I3sP5TDqRLtI0ZiycFyz2tw1sC3fLdnCe/9d53U41UJxkslVQGc3oXwJ/KyqD5VrVNVITq7y/h/r6NW8Np0bx3odjjEB45qTWnB6+3r8+9tlLEzN8DqcgFdkbS6gO/ACMBzniGS61eYqOz8t30bqrv1cZWXmjSlTIsLTF3UlvmY4N4ydZ/0n5aywI5NnfB5PALuADu7zp8s/tOrhgz/X0TA2gjM6WBltY8paXGQYL13SnU0Z+7ln4iLrPylHVpvLQ6vTs/ht1XZuO6ON3THOmHLSs2ltbh/Qlie/W87xM+twWZ+mXocUkApMJiJymap+JCK3+puuqs+WX1jVw4d/ric0WBjRq4nXoRgT0K49uQUz1uzgka+W0r1JHB0bWf9kWSvs63CU+zPaz6NmOccV8PYezGbinFTO6tyQ+Ohwr8MxJqAFBQnPDutKrchQbhw7j6yD2V6HFHAKTCaq+ob78+H8D2BPhUUYoD6ft4nMg9lccbwdchtTEerUDOeFEd1Zt2Mv931u/SdlraQn6v2e+jLFo6p8+Od6OjaKoUeTWl6HY0y10adFHW45vQ1fzE+z+l1lrKTJxEralsKstTtZsTWTK45vioi9lMZUpOv7taJvqzo8OGUJK7Zkeh1OwChpMrHjw1L44M/1xNYI5dyuCV6HYky1ExwkPDe8GzXDQ7l+7Fz2HbL+k7JQ2EWLmSKyx88jE7B7Y5bQlt0H+H7JFoYlNbbqwMZ4pF50BM8P78bq9CwenLzE63ACQmEd8NGqGuPnEa2qBQ4pLg4ReVREForIfBGZKiKN3PZL3faFIvKHiHT1WWagiKwQkRQRudunvbmIzBSRVSIyXkTCShNbeRs7awM5qjbW3RiPndi6Ljf0a8Wnc1KZNDfV63CqPK+ulHtKVbuoajfgK+ABt30tcIqqdgEeBfLuOx8MvAIMwrkK/2IR6eAu8yTwnKq2xrlKf1TF7caxOZyTy7hZGzilTTxN60QVvYAxplzd3L81vZrX5r4vFpOyze4fXxqeJBNV9R1aHIXbB6Oqf6jqLrd9BpB3y8FeQIqqrlHVQ8A4YIg4vdenAZ+5840BhpZ3/CX147JtbMs8yKW97ajEmMogJDiIF0d0JyI0mBvGzuXAYbt/fEl5VsNDRB4XkY3Apfx1ZOJrFPCt+3sC4DuOL9VtqwNkqGp2vvZK6ZNZG2gQE0G/tvFeh2KMcTWIjeDZYV1ZviWTR75a6nU4VVa5JRMR+UFEFvt5DAFQ1XtVNRH4GLgh37L9cJLJXXlNfjahhbQXFNNoEUkWkeT09PSS7FaJbdy5j19XpTPsuESrw2VMJXNq23pcd0pLxs7cwJcL0rwOp0oqsiPdHb2V/wN6N5AM3Kaqa/wtp6qnFzOGscDXOPeaR0S6AG8Dg1R1hztPKpDos0xjIA3YDsSJSIh7dJLX7peqvonbD5OUlFShw5vHz3YOrIYfl1jEnMYYL9w2oA2z1u7gn5MW0S0xjsTakV6HVKUU5yvys8AdOKePGgO3A2/h9Fu8W5KNikhrn6fnAsvd9ibAJOByVV3pM89soLU7cisMGAFMUacews/Ahe58I4HJJYmpPB3OyWVC8kZObRNPQlwNr8MxxvgRGhzECyO6g8CNn8zjcE6u1yFVKcVJJgNV9Q1VzVTVPe63+7NUdTxQ0logT7invBYCA4Cb3fYHcPpBXnWHDScDuEcdNwDfA8uACaqaNzj8LuBWEUlxl32nhDGVm5+WOx3vl1jHuzGVWmLtSJ44vwvzN2bwzNSVRS9g/qc414vkisgw/hoxdaHPtBKdKlLVCwpovxq4uoBp3wDf+GlfgzPaq9IaO3MD9WPCrePdmCrg7C4N+T0lkdenr6Zvqzqc1Nr+b4ujOEcmlwKXA9vcx+XAZSJSg3wd5+ZoeR3vw5Os492YquKBwR1pXa8m/zd+AemZB70Op0oo8tPNvbbjHFWt6z7OUdUUVd2vqr9XRJBVWV5l0uF2AyxjqowaYcG8dEl3Mg8c5rZPF5Cba+UIi1JkMhGRxiLyuYhsE5GtIjJRRBoXtZyB7Jxcxs+2jndjqqJ2DWK4b3AHfl2Zztu/+x20anwU57zLe8AUnOKOCcCXbpspwo9ux/vFdlRiTJV0We8mDOzYgP98t4IFGzO8DqdSK04yiVfV91Q12328D1iPVDGMm+V0vJ/Wrp7XoRhjSkBEePKCLtSPieDGT+aReeCw1yFVWsVJJttF5DIRCXYflwE7ilyqmtuy+wDTV6ZzUU/reDemKouNDOWFEd3YlLGf+75YbLf7LUBxPuX+BgwDtgCbcYYGX1WeQQWCiXNTyVW4sKd1LxlT1SU1q80t/VszeX4an82xcvX+FGc01wZVPVdV41W1nqoOBc6vgNiqLFXl0+SN9Gpem2Z1rdS8MYHgH/1a0adFbR6YvITV6VauPr+Snn+5tUyjCDDJ63exbsc+hiVZHS5jAkVwkPD88O5EhAZx49h5Vq4+n5ImE3/Veo1rwuyNRIUFc1bnBl6HYowpQw1iI3j6oq4s3byHJ75d7nU4lUpJk4n1QBVg78Fsvl60mcFdGhEZVqq7GxtjKqH+7etzVd9mvP/HOqYt3ep1OJVGgclERDJFZI+fRybONSfGj68XbWbfoRyGHWcd78YEqrsHtaNjoxju+GwBm3fv9zqcSqHAZKKq0aoa4+cRrar2lbsAnyZvpEV8FD2alLSgsjGmsgsPCeali7tzKDuXW8bNJ8fKrXh3295AtCY9i9nrdnFRz0Sc29MbYwJVi/iaPDKkEzPX7uTln1K8DsdzlkzK0GdzUgkOEi7oUWlvQ2+MKUMX9EhgaLdGvPDjSmat3el1OJ6yZFJGsnNymTg3lVPaxFMvJsLrcIwxFUBEeOy8ziTWjuSWcfPI2HfI65A8Y8mkjPy2ajtb9xxkWJJ1vBtTndQMD+Gli7uTnnWQOz9bWG3LrVgyKSOfztlI7agwTmtX3+tQjDEVrEvjOO48sx1Tl27loxnrvQ7HE5ZMysDOvYeYtnQr53VPICzEXlJjqqNRJzbnlDbxPPb1MlZsyfQ6nApnn3xl4MsFaRzOUSvqaEw1FhQkPH1RV6IjQrnxk7nVrtyKJZMyMGneJto3jKF9wxivQzHGeCg+OpxnhnVl5dYsHvt6qdfhVChLJqW0Oj2LBRszOL+7DQc2xsApbeIZfXILPpqxge8Wb/E6nApjyaSUPp+7iSCBId2swowxxnH7gLZ0TojlrokLScuoHuVWPEsmIvKoiCwUkfkiMlVEGuWbfpyI5IjIhT5tI0VklfsY6dPeU0QWiUiKiLwoFXT5eW6u8vm8TfRtVdeuLTHG/E9YSBAvXtydwzm53DK+epRb8fLI5ClV7aKq3YCvgAfyJohIMPAk8L1PW23gQaA30At4UETyCmC9BowGWruPgRWxA7PX7WRTxn4u6GEd78aYIzWvG8WjQzoxa+1OXvk58MuteJZMVHWPz9MojixrfyMwEdjm03YmME1Vd6rqLmAaMFBEGgIxqvqnOlcLfQAMLd/oHZ/P20RkWDADOtq1JcaYo53fI4Eh3Rrx/A8rSV4X2OVWPO0zEZHHRWQjcCnukYmIJADnAa/nmz0B2OjzPNVtS3B/z9/ub3ujRSRZRJLT09NLFfuBwzl8vXAzAzs1sPuWGGP8EhEeG9qJhFo1uHncfHbvP+x1SOWmXJOJiPwgIov9PIYAqOq9qpoIfAzc4C72PHCXquYfpO2vH0QLaT+6UfVNVU1S1aT4+PiS7ZTrh2VbyTyYzfnd7RSXMaZg0RGhvDiiO1v3HOCfkxYFbLmVcv1KraqnF3PWscDXOH0iScA4tw+9LnCWiGTjHHGc6rNMY+AXt71xvva00sRdHJ/P3UT9mHCOb1mnvDdljKniujepxa0D2vCf71Zw0uy6jOjVxOuQypyXo7la+zw9F1gOoKrNVbWZqjYDPgP+oapf4HTGDxCRWm7H+wDge1XdDGSKSB93FNcVwOTyjH171kGmr0xnaPcEgoPsviXGmKJdd3JL+raqw8NfLiVlW+CVW/Gyz+QJ95TXQpzEcHNhM6vqTuBRYLb7eMRtA/g78DaQAqwGvi23qHHKp2Tnqp3iMsYUW1CQ8OywbtQIC+bGT+YHXLkVz3qOVfWCYsxzZb7n7wLv+pkvGehUZsEV4fN5m+jQMIa2DaIrapPGmABQPyaCpy7swqgxyTzx7XIeOrej1yGVGbsC/hilbMtkYepuzre7KRpjSqB/+/pceUIz3v9jHT8u2+p1OGXGkskxmuSWTznXyqcYY0ro7kHtaN8whjs+W8jWPQe8DqdMWDI5RmkZ+51b80Zb+RRjTMlEhAbz0sXd2X8oh1snzCc3AMqtWDI5Rs+P6M4blyd5HYYxpoprVa8mD57Tgf+m7OCNX9d4HU6pWTIpAbubojGmLAw/LpGzOzfkmakrmLdhl9fhlIp9KhpjjEdEhH+d35n6MRHcNG4eew5U3XIrlkyMMcZDsTVCeWFENzbt2s/9XyyusuVWLJkYY4zHkprV5pbT2zB5fhqT5m7yOpwSsWRijDGVwPX9WtGreW3un7yYNelZXodzzCyZGGNMJRAcJDw/vBuhwUHcPG4+h7JzvQ7pmFgyMcaYSqJRXA2evKALizbt5umpK7wO55hYMjHGmEpkYKcGXNanCW/+uobpK0t3E7+KZMnEGGMqmfvO7kCb+jW5bcJ80jMPeh1OsVgyMcaYSsYpt9KDzAPZ3P7pgipRbsWSiTHGVEJtG0Rz3+AOTF+Zzrv/Xet1OEWyZGKMMZXUZb2bMKBDfZ78bjmLUnd7HU6hLJkYY0wlJSI8eUEX6kSFc9O4eew9mO11SAWyZGKMMZVYragwnh/RjXU79vLglCVeh1MgSybGGFPJ9WlRhxv7teKzOalMnl85y61YMjHGmCrgpv6t6dm0Fvd9vpgNO/Z5Hc5RLJkYY0wVEBIcxPPDu4HATePmcTincpVbsWRijDFVRGLtSJ44vwvzN2bw3LSVXodzBE+SiYg8KiILRWS+iEwVkUY+005125eIyHSf9oEiskJEUkTkbp/25iIyU0RWich4EQmr6P0xxpiKcnaXhgxPSuS16av5I2W71+H8j1dHJk+pahdV7QZ8BTwAICJxwKvAuaraEbjIbQ8GXgEGAR2Ai0Wkg7uuJ4HnVLU1sAsYVaF7YowxFezBczvQvG4Ut4yfz869h7wOB/AomajqHp+nUUBerYBLgEmqusGdb5vb3gtIUdU1qnoIGAcMEREBTgM+c+cbAwwt7/iNMcZLkWEhvHRxdzL2HeaOTxdUirszetZnIiKPi8hG4FLcIxOgDVBLRH4RkTkicoXbngBs9Fk81W2rA2Soana+dmOMCWgdG8Vy96B2/Lh8G2P+WOd1OOWXTETkBxFZ7OcxBEBV71XVROBj4AZ3sRCgJ3A2cCZwv4i0AcTPJrSQ9oJiGi0iySKSnJ5edUo7G2OMP1f1bcZp7erxr2+XszRtT9ELlKNySyaqerqqdvLzmJxv1rHABe7vqcB3qrpXVbcDvwJd3fZEn2UaA2nAdiBORELytRcU05uqmqSqSfHx8aXfSWOM8ZCI8NSFXYitEcqNn8xl3yHvyq14NZqrtc/Tc4Hl7u+TgZNEJEREIoHewDJgNtDaHbkVBowApqhzovBn4EJ3+ZHuOowxplqoUzOc54Z1Y832vTz61VLP4vCqz+QJ95TXQmAAcDOAqi4DvgMWArOAt1V1sdsncgPwPU5ymaCqeUVq7gJuFZEUnD6Udyp2V4wxxlsntq7LtSe35JNZG/lm0WZPYggpepayp6oXFDLtKeApP+3fAN/4aV+DM9rLGGOqrdsGtOHPNTu4e+JCujSOpXGtyArdvl0Bb4wxASA0OIgXR3QjV+GWcfPJruByK5ZMjDEmQDStE8VjQzuRvH4XL/2UUqHbtmRijDEBZGj3BM7vkcBLP61i5podFbZdSybGGBNgHhnSiSa1I7ll/Hwy9lVMuRVLJsYYE2Bqhofw4sXd2Z51kLsnLqqQciuWTIwxJgB1aRzHHWe25bslWxg7a0O5b8+SiTHGBKirT2zBSa3r8siXS1m5NbNct2XJxBhjAlRQkPDMsK5ER4Rw49h5HDicU37bKrc1G2OM8Vy96AieuqgrK7Zm8q9vlpXbdiyZGGNMgOvXth5Xn9icD/5cz9QlW8plG5ZMjDGmGrhjYFs6NorhzokL2bx7f5mv35KJMcZUA+Ehwbx0cXe6No4rl/V7UujRGGNMxWsRX5Mxfyufurh2ZGKMMabULJkYY4wpNUsmxhhjSs2SiTHGmFKzZGKMMabULJkYY4wpNUsmxhhjSs2SiTHGmFKTirhpSmUkIunA+iJmqwtsr4BwKhvb7+rF9rt6Ke1+N1XV+PyN1TaZFIeIJKtqktdxVDTb7+rF9rt6Ka/9ttNcxhhjSs2SiTHGmFKzZFK4N70OwCO239WL7Xf1Ui77bX0mxhhjSs2OTIwxxpSaJRM/RGSgiKwQkRQRudvreCqKiLwrIttEZLHXsVQkEUkUkZ9FZJmILBGRm72OqSKISISIzBKRBe5+Woh7PwAACD5JREFUP+x1TBVJRIJFZJ6IfOV1LBVFRNaJyCIRmS8iyWW6bjvNdSQRCQZWAmcAqcBs4GJVXeppYBVARE4GsoAPVLWT1/FUFBFpCDRU1bkiEg3MAYYG+t9cRASIUtUsEQkFfgduVtUZHodWIUTkViAJiFHVwV7HUxFEZB2QpKplfn2NHZkcrReQoqprVPUQMA4Y4nFMFUJVfwV2eh1HRVPVzao61/09E1gGJHgbVflTR5b7NNR9VItvlyLSGDgbeNvrWAKFJZOjJQAbfZ6nUg0+WIxDRJoB3YGZ3kZSMdxTPfOBbcA0Va0W+w08D9wJ5HodSAVTYKqIzBGR0WW5YksmRxM/bdXi21p1JyI1gYnALaq6x+t4KoKq5qhqN6Ax0EtEAv70pogMBrap6hyvY/FAX1XtAQwCrndPbZcJSyZHSwUSfZ43BtI8isVUELfPYCLwsapO8jqeiqaqGcAvwECPQ6kIfYFz3f6DccBpIvKRtyFVDFVNc39uAz7HOa1fJiyZHG020FpEmotIGDACmOJxTKYcuR3R7wDLVPVZr+OpKCISLyJx7u81gNP/v70zjbVriuL474+YSgyVmONFaYUGrdkz1JiQ4NUQnrHS9JOSKq8hocQYisaYGD4UEbRqSA3lodI5VNUsFWkEIRFCqoSU5cNe7+X0vHvuPXd4vWm6fslNzt1nD+vu9d5eezhnLeCr9ko1+JjZ9Wa2h5l1kP6/3zWzi9ss1qAjaYg/YIKkIcCpQMue3AxjksPM1gITgTdJB7Ezzezz9kq1fpD0LLAEGCHpe0nj2y3TeqITuIQ0Q13hn9PbLdR6YFdgnqRPSJOoXjPbaB6T3QjZGVgo6WPgfeA1M5vbqsrj0eAgCIKgaWJlEgRBEDRNGJMgCIKgacKYBEEQBE0TxiQIgiBomjAmQRAEQdOEMQkaQtJYSSZpv0xaRy2Pw2XylGi7S9JUv75Z0g/+OO8XkrpLlL9Z0rUl2tg/8/0WSSc3I7fXM0PSKvfUu1LSU5J2z9x/PfPuxx/FNfXnnyRp60rlG5DtHkkn1pB7haTlko5qpI1Wk9dT0D7CmASN0k3yMntBG9qeAjyS+T7dXYKcBTzqb7M3SxfQP0iZ2VQze7sF9QL0mNlBwAjgI9K7Hpt7O6f72+hlmQT0G5MGymd5EKgWcqHH+/k64NGylUrarEF5yrCOnsowyPJstIQxCerGfVh1AuMpMCaSxkl6RdJcjw1zU+b2ppIe9xgab/nb10iaIOkDn7XPzs64M/UOB/6u5ELbzL4G/gR28LzDvP0PJS3IrqIy9Q1oU9LRwJnANJ+JD/OZ+bmSTpM0M1N+jKQ5fn2qpCU+c5/l/VSIe+2dDvxE8pXUF29ip5yMY5SJuSHpIe/fq4DdSMZoXr68pMmSPvPPJE/rUIrbMqD/zexbYKikXarJDcwH9inqP0+fIek+l+suSYdLWqwUP2SxpBGeb5yklyXN8ZXPRJf7I0lLJe1YpMsCPVXUeV6eGr8vaIAwJkEjdAFzzWwl8Kuk0QX5DgcuAg4GzpN0qKfvCzxsZgcAvwHnePqLZnaYz9q/JBmrPJ3A8kqNuRxfu98hSLGurzSzQ4BrWXc108eANs1sMcmFTo+ZHWxm32Ty9wJHKrmjADgfeN4H8BuAk92R3jJgckG/5FkODDB0tTCzB0h+404wsxOy9yQdAlwOHAEcCUyQNMpvF/V/nyydNZo+A/jUr6vpbDipP64huWk5zsxGAVOBOzL5RgIXkv5ebgf+9HxLgEs9zwBdFuipms6z8gQtJpZ7QSN0k1x4Q3KU103lAb7XzH4BkPQicAzwMrDKzFZ4ng+BDr8eKek2YHtgG5JLmzy7Aj/n0q6WNAHYG3dU6KuCo4FZUr8j6C0q1FemzX7MbK2kucAZkl4gxcSYAhxP2m5Z5O1tThoMy1DJU3WzHAO8ZGZroL//jyUNvkX9D8kV/W4FdU6TdAOp//uMRrX+m2Vm//r1dsCTkvYleeHObkXO8zgyqyX9Dszx9E+BA8vqskS+rDxBiwljEtSFpKHAiaRBxIBNAZM0pUL2vK+evu9/Z9L+Bbby6xmkCIcfSxoHjKlQ51+kgSnLdDO7R9LZwFOShpFW3b/5Hn81yrSZ53ngClIgsQ/MbLXS6NVrZjUfAKjAKOCdKvfXsu4uwpYl6qxmoIr6v6/uvwrK9ZjZC7m0GRT335rM9a0kozFWKWbMewXy/Jf5/h9pjCqry1r51hSkBy0gtrmCejmXFNZ3LzPrMLM9gVWkmXCeUyTt6HvyXcCiGnVvC/yodIB+UUGeL/H9+jzuOn4ZcJnHI1kl6TxInoElHVRHm6v9XiXeA0YDE0iGBWAp0Cmp7yxha6XznUJcpqtIq61qDve+BfaXtIWk7YCTSsg5H+hyOYYAY4EF1eRxhlOfJ9kyOoM0AfjBr8fVUT81dNn/++vQeTAIhDEJ6qWbFAchy2zSnneehcDTwApgtpktq1H3jaQIh70Uu0KfD4xSZh8jxy3AZEmbkAa38UpeUj+ncvjlojafA3r8IHhYtoBvlbxKOjR/1dN+Jg2Szyp54V1K8TnINJdpJXAY6czjn4K8mNl3wEzgE+AZ0hNgfTwGvOEHy9kyy0mrhvf99z1hZtlyA3CDsA/JIJeljM4A7gbulLSItJqtlyJd5vVURufBIBBeg4NBwbc8DjWziYNQ9/3AnBY+qhuQ3h0CRpvZje2WJdjwiJVJsCFyB5l3K4KWsRlwb7uFCDZMYmUSBEEQNE2sTIIgCIKmCWMSBEEQNE0YkyAIgqBpwpgEQRAETRPGJAiCIGiaMCZBEARB0/wP0Ev90jZIMx8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_results['LL'] = -1*plot_results['NLL']\n",
    "sns.lineplot(data = plot_results, x = 'Alpha', y = 'LL')\n",
    "plt.axvline(x = 1, color = 'grey')\n",
    "plt.suptitle('Log Likelihood when Adjusting for Dilution')\n",
    "plt.title('Inferring Dilution from Simulated Data')\n",
    "plt.xlabel('Alpha (Relative Dilution) Parameter')\n",
    "plt.ylabel('Log Likelihood')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Meta-Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SNP</th>\n",
       "      <th>CHR</th>\n",
       "      <th>POS</th>\n",
       "      <th>BETA_Gold</th>\n",
       "      <th>BETA_Diluted</th>\n",
       "      <th>SE_Gold</th>\n",
       "      <th>SE_Diluted</th>\n",
       "      <th>Z</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.069529</td>\n",
       "      <td>0.010957</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>1.622046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0.038759</td>\n",
       "      <td>0.015121</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>1.731773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>0.045990</td>\n",
       "      <td>-0.000488</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.362680</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>-0.074794</td>\n",
       "      <td>-0.010719</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>-1.646765</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>-0.018515</td>\n",
       "      <td>0.011009</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.846813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4995</td>\n",
       "      <td>4996</td>\n",
       "      <td>1</td>\n",
       "      <td>4996</td>\n",
       "      <td>0.078014</td>\n",
       "      <td>0.004748</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>1.127127</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4996</td>\n",
       "      <td>4997</td>\n",
       "      <td>1</td>\n",
       "      <td>4997</td>\n",
       "      <td>0.040007</td>\n",
       "      <td>0.025823</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>2.725419</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4997</td>\n",
       "      <td>4998</td>\n",
       "      <td>1</td>\n",
       "      <td>4998</td>\n",
       "      <td>-0.057730</td>\n",
       "      <td>0.006685</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.102326</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4998</td>\n",
       "      <td>4999</td>\n",
       "      <td>1</td>\n",
       "      <td>4999</td>\n",
       "      <td>-0.003057</td>\n",
       "      <td>0.001067</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.070844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4999</td>\n",
       "      <td>5000</td>\n",
       "      <td>1</td>\n",
       "      <td>5000</td>\n",
       "      <td>0.038027</td>\n",
       "      <td>0.004365</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.737671</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5000 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       SNP  CHR   POS  BETA_Gold  BETA_Diluted   SE_Gold  SE_Diluted         Z\n",
       "0        1    1     1   0.069529      0.010957  0.044721        0.01  1.622046\n",
       "1        2    1     2   0.038759      0.015121  0.044721        0.01  1.731773\n",
       "2        3    1     3   0.045990     -0.000488  0.044721        0.01  0.362680\n",
       "3        4    1     4  -0.074794     -0.010719  0.044721        0.01 -1.646765\n",
       "4        5    1     5  -0.018515      0.011009  0.044721        0.01  0.846813\n",
       "...    ...  ...   ...        ...           ...       ...         ...       ...\n",
       "4995  4996    1  4996   0.078014      0.004748  0.044721        0.01  1.127127\n",
       "4996  4997    1  4997   0.040007      0.025823  0.044721        0.01  2.725419\n",
       "4997  4998    1  4998  -0.057730      0.006685  0.044721        0.01  0.102326\n",
       "4998  4999    1  4999  -0.003057      0.001067  0.044721        0.01  0.070844\n",
       "4999  5000    1  5000   0.038027      0.004365  0.044721        0.01  0.737671\n",
       "\n",
       "[5000 rows x 8 columns]"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_results['Z'] = phe.compute_meta_scores(betas, ses, optimizer.x)\n",
    "df_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SNP</th>\n",
       "      <th>CHR</th>\n",
       "      <th>POS</th>\n",
       "      <th>BETA_Gold</th>\n",
       "      <th>BETA_Diluted</th>\n",
       "      <th>SE_Gold</th>\n",
       "      <th>SE_Diluted</th>\n",
       "      <th>Z</th>\n",
       "      <th>TEMA_P</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.069529</td>\n",
       "      <td>0.010957</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>1.622046</td>\n",
       "      <td>0.104794</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0.038759</td>\n",
       "      <td>0.015121</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>1.731773</td>\n",
       "      <td>0.083314</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>0.045990</td>\n",
       "      <td>-0.000488</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.362680</td>\n",
       "      <td>0.716844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>-0.074794</td>\n",
       "      <td>-0.010719</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>-1.646765</td>\n",
       "      <td>0.099606</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>-0.018515</td>\n",
       "      <td>0.011009</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.846813</td>\n",
       "      <td>0.397100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4995</td>\n",
       "      <td>4996</td>\n",
       "      <td>1</td>\n",
       "      <td>4996</td>\n",
       "      <td>0.078014</td>\n",
       "      <td>0.004748</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>1.127127</td>\n",
       "      <td>0.259689</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4996</td>\n",
       "      <td>4997</td>\n",
       "      <td>1</td>\n",
       "      <td>4997</td>\n",
       "      <td>0.040007</td>\n",
       "      <td>0.025823</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>2.725419</td>\n",
       "      <td>0.006422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4997</td>\n",
       "      <td>4998</td>\n",
       "      <td>1</td>\n",
       "      <td>4998</td>\n",
       "      <td>-0.057730</td>\n",
       "      <td>0.006685</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.102326</td>\n",
       "      <td>0.918498</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4998</td>\n",
       "      <td>4999</td>\n",
       "      <td>1</td>\n",
       "      <td>4999</td>\n",
       "      <td>-0.003057</td>\n",
       "      <td>0.001067</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.070844</td>\n",
       "      <td>0.943522</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4999</td>\n",
       "      <td>5000</td>\n",
       "      <td>1</td>\n",
       "      <td>5000</td>\n",
       "      <td>0.038027</td>\n",
       "      <td>0.004365</td>\n",
       "      <td>0.044721</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.737671</td>\n",
       "      <td>0.460714</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5000 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       SNP  CHR   POS  BETA_Gold  BETA_Diluted   SE_Gold  SE_Diluted  \\\n",
       "0        1    1     1   0.069529      0.010957  0.044721        0.01   \n",
       "1        2    1     2   0.038759      0.015121  0.044721        0.01   \n",
       "2        3    1     3   0.045990     -0.000488  0.044721        0.01   \n",
       "3        4    1     4  -0.074794     -0.010719  0.044721        0.01   \n",
       "4        5    1     5  -0.018515      0.011009  0.044721        0.01   \n",
       "...    ...  ...   ...        ...           ...       ...         ...   \n",
       "4995  4996    1  4996   0.078014      0.004748  0.044721        0.01   \n",
       "4996  4997    1  4997   0.040007      0.025823  0.044721        0.01   \n",
       "4997  4998    1  4998  -0.057730      0.006685  0.044721        0.01   \n",
       "4998  4999    1  4999  -0.003057      0.001067  0.044721        0.01   \n",
       "4999  5000    1  5000   0.038027      0.004365  0.044721        0.01   \n",
       "\n",
       "             Z    TEMA_P  \n",
       "0     1.622046  0.104794  \n",
       "1     1.731773  0.083314  \n",
       "2     0.362680  0.716844  \n",
       "3    -1.646765  0.099606  \n",
       "4     0.846813  0.397100  \n",
       "...        ...       ...  \n",
       "4995  1.127127  0.259689  \n",
       "4996  2.725419  0.006422  \n",
       "4997  0.102326  0.918498  \n",
       "4998  0.070844  0.943522  \n",
       "4999  0.737671  0.460714  \n",
       "\n",
       "[5000 rows x 9 columns]"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#create a 2 sided p-value from z-statistic\n",
    "df_results['TEMA_P'] = 2*norm.cdf(-1*np.abs(df_results['Z']))\n",
    "df_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Create Confidence Intervals for Effective Dilution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for real data\n",
    "#start with random clumping SNPs in 250kb region, r^2 threshold = .5\n",
    "#pvalues for each snp are produced through simulation of a uniform distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for confidence intervals- important that observations are sorted by CHR, POS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "628c900707314176a8a57507162dc524",
       "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": [
      "\n",
      "0\n",
      "    Shift  CovRatio   BaseCov  AbsShift  Lambda  LambdaG\n",
      "-1    0.0       1.0  0.000985       0.0       1      0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/davidburstein/Documents/GitHub/PheMED/bootstrap.py:78: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_corr_mini = df_corr.loc[df_corr['AbsShift'] <= block_size]\n",
      "/Users/davidburstein/Documents/GitHub/PheMED/bootstrap.py:79: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_corr_mini['Lambda'] = df_corr_mini['Shift'].apply(lambda x: trap_function(x, total_distance = block_size))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_results = df_results.sort_values(by = ['CHR','POS']).reset_index(drop = True)\n",
    "corr_stats = boots.compute_autocorr(df_results, 'BETA_Gold', shift_range = 500)\n",
    "block = boots.compute_block_size(corr_stats, df_results)\n",
    "print(block)\n",
    "g, d = boots.compute_g_and_d(corr_stats, 2*block)\n",
    "final_block_size_gold = boots.compute_final_block_size(df_results, g, d)\n",
    "final_block_size_gold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a13e9428fda940368f528773ab765a34",
       "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": [
      "\n",
      "0\n",
      "    Shift  CovRatio   BaseCov  AbsShift  Lambda  LambdaG\n",
      "-1    0.0       1.0  0.000084       0.0       1      0.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr_stats = boots.compute_autocorr(df_results, 'BETA_Diluted', shift_range = 500)\n",
    "block = boots.compute_block_size(corr_stats, df_results)\n",
    "print(block)\n",
    "g, d = boots.compute_g_and_d(corr_stats, 2*block)\n",
    "final_block_size_diluted = boots.compute_final_block_size(df_results, g, d)\n",
    "final_block_size_diluted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_block_size = np.max([final_block_size_gold, final_block_size_diluted,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fb58a14c66034790b8b0015cf3c98b70",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=2000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(18)\n",
    "df_results_sim = pd.DataFrame(columns = ['Sample', 'Alpha', 'Convergence'])\n",
    "for i in tqdm_notebook(range(2000)):\n",
    "    df_results_sample = boots.circular_block_bootstrap(df_results,final_block_size)\n",
    "    #df_meta_sample = df_meta.sample(n = df_meta.shape[0], replace = True)\n",
    "    betas = df_results_sample[beta_vars]\n",
    "    ses = df_results_sample[se_vars]\n",
    "    optimizer = minimize(lambda x: phe.nll_data(betas,ses, x), np.array((1,1)), \n",
    "                        options={'maxiter': 300}, method = 'Nelder-Mead')\n",
    "    df_results_sim.loc[i] = [i, optimizer.x[1], optimizer.message]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_alpha = df_results_sim.sort_values(by = 'Alpha', ascending = False)['Alpha']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Top Performing Alpha, 100 Bootstrap Samples')"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEjCAYAAAAlhuZMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwdVZ338c83CwRIWLKA7AFkURAD04CKcVg1orIoIhlFYJDAMzKDg48bIKAjjwsKwsiAQXgFRBbZJDKKIjKgPmwJhBBkC5sEsjT7lgCB3/xxTmPl5nb3re6+XfeS7/v1uq++dapu1a/rVtWvTp1zqxQRmJmZNWpI1QGYmVl7ceIwM7NSnDjMzKwUJw4zMyvFicPMzEpx4jAzs1KcOFYQkr4kaZGklySt1uRlrZyXs14zl9PNsidJmjvQ05p1R9IISSFpg6pjGSxOHCXkg2HX601JiwvDnx3gZd0iaUmed6ekyySt3cd5rQZ8H5gYESMj4uWBjLVWRLyal/Nks5Yh6ci8s+7drGX0h6TtJf1e0tOSltQZP07SryW9LOkRSfvXjD9Y0t/y93+5pDW6WU7XQevlPO1Tkn4uadQA/A9lkvBWkpb2d5ll5P/9DElP5P/9YUnfH8wYVlROHCXkg+HIiBgJ/A34RKHsF01Y5Bfyst4FvAP4QdkZSBoGrAsMjYj7+/D5IZJacTs5GHgm/21FrwIXA0d2M34q8CywNnAYcJ6kzQEkbQecAXyG9N0JOL2X5W2Zt5V3AusDx/X3HxhoeVscSCeS9o3tgVHAHsDsAV6G1RMRfvXhBTwK7FFTtgpwJjAfmAecAgzP4yYBc4FvkQ54DwOf7mH+twCfKwx/GZhRWM6PgceBBcB/AivXLOebwELgBuBlIICXgN/m6f4RuAN4Pi9rh5plfxu4FVgCbJDLTgJuy/O5EhgD/BJ4IY/fIH9+RF5e1/AlOd7fAS8CfwE2LizvY8CDwHN5umX+9zrrZgvgTWB/0gF6dGHcJGBuYXgB8FXgvrzep9ZZV8cCncATwGcLn90PuCv/f48Bx/ZhO9kGWFJTthawtGYdXAaclN+fCpxXGLc1sBgYUWf+y6zrXHYMML0wvBHwm/z/PwAc3Ns2m7/bxXk9v5RfY4CdgTvzOlkAfDfPZ1FhG3sJ2I6UNP+Y5/8scDywFfA/OZZO4HxgVCPfV53//Q/AkT2s+xOAR/I2Nwf4WGFcV2w/Ie0DDwIdwJS8HSwEDixMfwlpP7shz+96YP1utvee9s93ANeStvWngT9WfSzry6vyANr1Rf3E8QPgT8BYYB3gduC4PG5SPlh8F1iJdHb0CrBJN/N/6+BJOiv9M3BOHj4buBxYE1iDdEA+sWY5387LWSXvrEsL81477/gHAMOAQ/JOvEZh2Q8DW5IOIsNy2b3AeGB03tHuIyWgYcClwFn58/USxyLSmeHwHPu0PO4dpAPNx/O4rwKv03PiOBm4Kb9/APiXwrh6ieNOYD1gXP5Oji9M+zrp7Hw4KVG8CIzM43cnHbSH5NifASaV3E7qJY73A8/WlB0PXJbf/w44umb8a8DWdeZfu67HkA7MxxamuRU4DViZdHB8Bti5wW12bs3y7iSf8JDO8nfK75fZxnLZkaRt8XBgKH/fFncjbZvvyNvV9xr5vur8798hJYYju1k3XTW2IcBB+bsdW4jtdeCfSNvvKaSTg9NybHuTkt2Iwjb8XP7uRpD2wT908x30tH+eRqo9DsvL+VDVx7K+vCoPoF1f1E8cTwC7FYb3Ae7L7yeRzt5HFMZPB77SzfxvIdUUniOdCZ5POmAPIx1E1i9Muytwb2E5L5NrOrmsNnEcTj7wFsruJJ9h5WUfWyeeLxeGzwSuKgx/Grglv6+XOH5SmPaTwKz8fgpwQ2HcEFKSqZs48vi/kc80STW4Wwvj6yWOQ2qWfU9h2ueBIYXxLwATuln22eQz7BLbSb3EsSfwaE3ZvwLX5vd/Kcacy54G3ldn/l3r+vm8rbxJOrteJ4/fPG93qxQ+cxpwdoPbbG3iuI2UaMfUlHeXOB7oZf0cCNzcyPdV57PDgaOBm0k1z3nA5B6WdR/wkUJsdxfG7ZDX4xqFspeBrQrb8LTCuNF5+nHF7Z3e988fkGqXm5bZjlrt1YrXrtuSJJHOoB4rFD9Gut7cpTMiltSM76nn0RERsWZEbBARB0fEM3n64cA9kp6T9BzwK1ItosuCiHi9h/muVxNnvVgfr/O5hYX3i+sMj+xhmQsK718pTLtecVkR8SbpYNadXUnr+bI8/AtgR0lb9fCZ4v9Su8478zKXi03SzpJuzJ0TnifVzMb2sJxGvQSsXlO2OumMuLvxowrj69k6ItYkHcQuAG6StBLpf+2MiMWFaR8D1m9wm611MLAt8ICkWyV9pIdpoWY7krRe7ujxhKQXgJ+x/Drt6ft6S0S8HhGnR8T7SZf/TgUukLRZXtZhkmYX9pN31iyrdvt9NSKerykrbtPF7fQZ0vdUG1tv++fJwJPADZLmSjqm3v/W6pw4Bkik04kFwMaF4o1Y9iA4VtKImvFlex7NJ1X/N8tJZc2IWCMixhTD6WUeT9bEWS/W3uYxUOaTztSA1BhP7weuIaQdcwFwEynWz/fwmQ0L78us81+SLsFtGBFrANNIDdX9dR+wuqSNCmXvBe7J7+/JwwBIejfwBvBQbzOOiNeAc0ntQFuQ/tdxklYpTLYR8EQD2+xy20BE3BsRnyEdCM8ArswJqrvtpbb8FNKZ/DYRsTrwBZZfp6W/r4h4JSJOJdU8tpK0BaltYQqpDWxNUntWf76/t+KSNJqUVObXTNPj/hkRz0fE0RGxMfAp4HhJO/cjpko4cQysi4ETJY3JXWePAy4sjB8OfFPSSpJ2I12yuKLMAnJN4jzgdEljlWwoac8Ss5kObCdpf0nDJH2etINeWyaWATId2EnSXrnXzTGks8flSBpJunRxCDCh8Pq/wEE99P76N0nrShoLfJ2UDHqUz8ZHAk9HxBJJHyBdjitOs0DSgd19Pp8krJSHR+QDLBHxLHAN8G1Jq0rahXRZqKtn3oXApyS9L//P3wIuramtdhd3V5vVS6Sz9bmknkbfUfp9zfak5Nu1rJ622YXA2jmGrvl/XtKYiHiDdHksSJfHFgFDa5JhPaNybC/kaeudcTf0fUn6sqSJed0OlzSF1JZyF+m7e5PUdjdE0pGkGkd/7CNpJ0krk9pXboiIRcUJets/Je0taZO8fT1POiF4o59xDTonjoF1AvBX0hnjLNK16mIX2kdJZyMLSBvXoRHxcB+W8yXSWdgM0sZ3LSV2iohYSGr8O4507fwo4OMR8VwfYumXiJgPTCadvT5Fqn3cTTpzrLU/qWH34ohY0PUi9bxZndToWs8lpN4wD+Z599qtOZ+NHwn8UNKLpEb7rstjSFqVdBC8rZtZbEm61DGT1Ci9mGW7ih5OumzyFKkmc1hEPJiXfQfpO76cdPAeQrqW35P7Jb1EWj8HAPtExIv5/zgAeDdpu7uU1K72p/y5nrbZu0iJ/bF82WU0qRPD/XmdfBc4ICKW5mT4A2BmnnZCN3GeAHyQtN1eRf0Tp0a/r1dJ281CUuI6FNg3IubldXg2aR+ZD2yS3/fHhcD3SN/Zu+i+K3hP++e7SJ0XXiTVln8YEbf0M65Bp9xgY00maRKpgbi/Zz1va/mMeQHpNzI3D8D8FgD7R8Sf+x3csvPdg9R199CBnO+KrlnfV39JugSYExHfqTqWVuAah1VO0kclrZEv7ZxIaqCeWXFYPYqIPzhp2IrKicNawYdI/fEXkX47sV9u5DWzFuRLVWZmVoprHGZmVspA33RsUI0dOzbGjx9fdRhmZm1l5syZT0XEuL5+vq0Tx/jx45kxo7897MzMViySau8cUYovVZmZWSlOHGZmVooTh5mZleLEYWZmpThxmJlZKU4cZmZWihOHmZmV0rTEIek8SYskzSmUXSppVn49KmlWLh8vaXFh3NnNisvMzPqnmT8AnAb8hPQYSwDyk8MAkPQj0r3quzwUEd3dw9/MzFpE0xJHRNwkaXy9cfnpVwfQ/YN3zFZYs57q9UF/TBg7otdpzJqlqjaOicDCrieeZZtIulPSjZImdvdBSVMkzZA0o7Ozs/mRmpnZMqpKHJNJzzruMh/YKCK2Iz2D+CJJq9f7YERMjYiOiOgYN67P9+gyM7M+GvTEkR8N+kkKD6CPiFcj4un8fibwELDFYMdmZma9q6LGsQdwX0TM6yqQNE7S0Px+U2Bz4OEKYjMzs140szvuxcDNwJaS5kk6LI86kGUvU0F6dOhsSXcBlwNHRsQzzYrNzMz6rpm9qiZ3U35InbIrgCuaFYuZmQ0c/3LczMxKceIwM7NSnDjMzKwUJw4zMyvFicPMzEpx4jAzs1KcOMzMrBQnDjMzK8WJw8zMSnHiMDOzUpr5BEAza5JGHvbUCD8QyvrCNQ4zMyvFicPMzErxpSqzQTJQl5fMquYah5mZleIah9kAcG3CViSucZiZWSlOHGZmVooTh5mZleLEYWZmpTQtcUg6T9IiSXMKZSdJekLSrPzaqzDuG5LmSrpf0keaFZeZmfVPM2sc04BJdcpPi4gJ+fUbAEnvBg4Ets6f+S9JQ5sYm5mZ9VHTEkdE3AQ80+Dk+wCXRMSrEfEIMBfYsVmxmZlZ31XRxnGUpNn5UtZauWx94PHCNPNy2XIkTZE0Q9KMzs7OZsdqZmY1BjtxnAVsBkwA5gM/yuWqM23Um0FETI2IjojoGDduXHOiNDOzbg1q4oiIhRHxRkS8CZzD3y9HzQM2LEy6AfDkYMZmZmaNGdTEIWndwuB+QFePq+nAgZJWlrQJsDlw22DGZmZmjWnavaokXQzsAoyVNA84EdhF0gTSZahHgSMAIuIeSb8E/gosBb4YEW80KzYzM+u7piWOiJhcp/jcHqY/GTi5WfGYmdnA8C/HzcysFCcOMzMrxYnDzMxKceIwM7NS/ARAs1746X5my3KNw8zMSnHiMDOzUpw4zMysFCcOMzMrxYnDzMxKceIwM7NSnDjMzKwUJw4zMyvFicPMzEpx4jAzs1KcOMzMrBQnDjMzK8WJw8zMSnHiMDOzUpw4zMysFCcOMzMrpWmJQ9J5khZJmlMoO0XSfZJmS7pK0pq5fLykxZJm5dfZzYrLzMz6p5k1jmnApJqy64BtImJb4AHgG4VxD0XEhPw6solxmZlZPzQtcUTETcAzNWW/j4ilefAWYINmLd/MzJqjyjaOfwZ+WxjeRNKdkm6UNLG7D0maImmGpBmdnZ3Nj9LMzJZRSeKQdBywFPhFLpoPbBQR2wHHABdJWr3eZyNiakR0RETHuHHjBidgMzN7y6AnDkkHAx8HPhsRARARr0bE0/n9TOAhYIvBjs3MzHo3qIlD0iTga8DeEfFKoXycpKH5/abA5sDDgxmbmZk1ZlizZizpYmAXYKykecCJpF5UKwPXSQK4Jfeg+hDwbUlLgTeAIyPimbozNjOzSjUtcUTE5DrF53Yz7RXAFc2KxczMBo5/OW5mZqU4cZiZWSlOHGZmVkrT2jjM2sGsp5ZUHYJZ23GNw8zMSnHiMDOzUpw4zMysFCcOMzMrpaHEIWmbZgdiZmbtodFeVWdLWon0cKaLIuK55oVkZq2k0Z5nE8aOaHIk1ioaqnFExAeBzwIbAjMkXSRpz6ZGZmZmLanhNo6IeBA4nnR3238EzsjPD/9ks4IzM7PW02gbx7aSTgPuBXYDPhER78rvT2tifGZm1mIabeP4CXAOcGxELO4qjIgnJR3flMjMrOn8y3nri0YTx17A4oh4A0DSEGBERLwSET9vWnRmZtZyGm3j+AOwSmF41VxmZmYrmEYTx4iIeKlrIL9ftTkhmZlZK2s0cbwsafuuAUn/ACzuYXozM3ubarSN40vAZZKezMPrAp9pTkhmZtbKGkocEXG7pK2ALQEB90XE602NzMzMWlKZmxzuAGwLbAdMlvT53j4g6TxJiyTNKZSNlnSdpAfz37VyuSSdIWmupNnFS2NmZtY6Gv0B4M+BHwIfJCWQHYCOBj46DZhUU/Z14PqI2By4Pg8DfBTYPL+mAGc1EpuZmQ2uRts4OoB3R0SUmXlE3CRpfE3xPsAu+f35wP+QbmOyD3BBXsYtktaUtG5EzC+zTDMza65GL1XNAd4xQMtcpysZ5L9r5/L1gccL083LZcuQNEXSDEkzOjs7BygkMzNrVKM1jrHAXyXdBrzaVRgRew9gLKpTtlwNJyKmAlMBOjo6StWAzMys/xpNHCcN4DIXdl2CkrQusCiXzyPdtr3LBsCTy33azMwq1Wh33BslbQxsHhF/kLQqMLSPy5wOHAx8L/+9ulB+lKRLgJ2A592+YdY+Grlhoh/29PbQaK+qw4HLgZ/movWBXzXwuYuBm4EtJc2TdBgpYewp6UFgzzwM8BvgYWAu6U68/1Li/zAzs0HS6KWqLwI7ArdCeqiTpLV7/ghExORuRu1eZ9rIyzEzsxbWaK+qVyPita4BScOo03BtZmZvf40mjhslHQuskp81fhnw6+aFZWZmrarRxPF1oBO4GziC1B7hJ/+Zma2AGu1V9Sapwfqc5oZjZmatrqHEIekR6v8Yb9MBj8jMzFpamXtVdRkBfBoYPfDhmJlZq2uojSMini68noiIHwO7NTk2MzNrQY1eqio+G2MIqQYyqikRmZlZS2v0UtWPCu+XAo8CBwx4NGZm1vIa7VW1a7MDMTOz9tDopapjehofEacOTDhmZtbqyvSq2oF0B1uATwA3seyDl8zMbAVQ5kFO20fEiwCSTgIui4gvNCswMzNrTY3ecmQj4LXC8GvA+AGPxszMWl6jNY6fA7dJuor0C/L9gAuaFpWZmbWsRntVnSzpt8DEXHRoRNzZvLDMutfIk+bMrHkavVQFsCrwQkScDsyTtEmTYjIzsxbW6KNjTwS+BnwjFw0HLmxWUGZm1roarXHsB+wNvAwQEU/iW46Yma2QGm0cfy0iQlIASFqtiTGZ2dtUI+1TE8aOGIRIrD8aTRy/lPRTYE1JhwP/TB8f6iRpS+DSQtGmwAnAmsDhpCcNAhwbEb/pyzLMzKx5Gu1V9cP8rPEXgC2BEyLiur4sMCLuByYASBoKPAFcBRwKnBYRP+zLfM3MbHD0mjjywf13EbEH0Kdk0YPdgYci4jFJAzxrMzNrhl4bxyPiDeAVSWs0YfkHAhcXho+SNFvSeZLWasLyzMysnxrtVbUEuFvSuZLO6Hr1Z8GSViL11LosF50FbEa6jDWfZZ8BUvzcFEkzJM3o7OysN4mZmTVRo43j/51fA+mjwB0RsRCg6y+ApHOAa+p9KCKmAlMBOjo6YoBjMjOzXvSYOCRtFBF/i4jzm7DsyRQuU0laNyLm58H9gDlNWKaZmfVTb5eqftX1RtIVA7VQSasCewJXFop/IOluSbOBXYF/H6jlmZnZwOntUlWxq9OmA7XQiHgFGFNTdtBAzd/MzJqntxpHdPPezMxWUL3VON4r6QVSzWOV/J48HBGxelOjMzOzltNj4oiIoYMViJmZtYcyz+MwMzNr+HccZmaDwnfQbX2ucZiZWSlOHGZmVooTh5mZleI2DmspjVzfNrNqucZhZmalOHGYmVkpThxmZlaKE4eZmZXixGFmZqU4cZiZWSlOHGZmVooTh5mZleLEYWZmpfiX42bWdhq9w4DvotscrnGYmVkpThxmZlZKZZeqJD0KvAi8ASyNiA5Jo4FLgfHAo8ABEfFsVTGamdnyqq5x7BoREyKiIw9/Hbg+IjYHrs/DZmbWQqpOHLX2Ac7P788H9q0wFjMzq6PKxBHA7yXNlDQll60TEfMB8t+1az8kaYqkGZJmdHZ2DmK4ZmYG1XbH3TkinpS0NnCdpPsa+VBETAWmAnR0dEQzAzQzs+VVVuOIiCfz30XAVcCOwEJJ6wLkv4uqis/MzOqrJHFIWk3SqK73wIeBOcB04OA82cHA1VXEZ2Zm3avqUtU6wFWSumK4KCKulXQ78EtJhwF/Az5dUXxmZtaNShJHRDwMvLdO+dPA7oMfkZmZNarVuuOamVmLc+IwM7NSnDjMzKwUJw4zMyvFicPMzEpx4jAzs1KcOMzMrBQnDjMzK8WJw8zMSnHiMDOzUpw4zMysFCcOMzMrxYnDzMxKceIwM7NSnDjMzKwUJw4zMyvFicPMzEpx4jAzs1Kqeua4vc3MempJ1SGY2SBxjcPMzEoZ9MQhaUNJN0i6V9I9ko7O5SdJekLSrPzaa7BjMzOz3lVxqWop8OWIuEPSKGCmpOvyuNMi4ocVxGRmZg0a9MQREfOB+fn9i5LuBdYf7DjMzKxvKm3jkDQe2A64NRcdJWm2pPMkrdXNZ6ZImiFpRmdn5yBFamZmXSpLHJJGAlcAX4qIF4CzgM2ACaQayY/qfS4ipkZER0R0jBs3btDiNTOzpJLuuJKGk5LGLyLiSoCIWFgYfw5wTRWxmdnbRyPdxCeMHTEIkby9VNGrSsC5wL0RcWqhfN3CZPsBcwY7NjMz610VNY6dgYOAuyXNymXHApMlTQACeBQ4ooLYzMysF1X0qvozoDqjfjPYsZiZWXm+5Yj1yLcSMbNavuWImZmV4sRhZmalOHGYmVkpThxmZlaKG8fNbIXmHwmW5xqHmZmV4sRhZmalOHGYmVkpbuNYgfnHfWaNcTvIslzjMDOzUpw4zMysFF+qepvyZSgzaxbXOMzMrBQnDjMzK8WJw8zMSnHiMDOzUtw4bmY2AFak33q4xmFmZqW4xtGG3NXWzKrkxDEAVqQqqpn13UCe9FV5TGm5S1WSJkm6X9JcSV+vOh4zM1tWS9U4JA0FzgT2BOYBt0uaHhF/rTay/nOtxMzeLloqcQA7AnMj4mEASZcA+wBNSRytdjB324WZtYNWSxzrA48XhucBOxUnkDQFmJIHX5J0/yDFVjQWeKqC5faFY22edorXsTZPO8XbFevG/ZlJqyUO1SmLZQYipgJTByec+iTNiIiOKmNolGNtnnaK17E2TzvFO1Cxtlrj+Dxgw8LwBsCTFcViZmZ1tFriuB3YXNImklYCDgSmVxyTmZkVtNSlqohYKuko4HfAUOC8iLin4rDqqfRSWUmOtXnaKV7H2jztFO+AxKqI6H0qMzOzrNUuVZmZWYtz4jAzs1KcOMyspUmq102/ZbVTvH2N1YnDzFqSpOEA0SYNsZK2hvaIV9JI6HusThx9JGkvSZ+oOo5GSBonab2aspY9K5K0q6Rdqo6jEZI2lrRFTVlLrts222b3Bs6QdL6krSWNrjqmnkj6CPALSZtXHUtv8ro9V9IleZvYqOw8nDj6QNKewCnAy1XH0htJ+wP/DVwt6T8kTYR0ptGKBzhJHwZ+DLxWdSy9yev2V8AFkk6V9DlozXXbZtvsNsBZwGXAQuCLwOclrV9pYN3IB+ITgS9GxIOt9t0X5ZOcs4AzgJuBnYEvS9qy1HzaoFbVUvKZ8DRg/4iYkat8KwPPRsSbVcZWS9IY4GrSjrcQOAoYAfz/iLiyytjqkfRR4Bxgj4i4T9IqpG30lYpDW46k1Ug/Tv0a6Sack4HtSDfp/HGVsdVqp20WQNJuwOERMTkPfxzYBfgbcEFEPFdheMvISeJO4MWImChpHeAg0j2hLiVtDy9WGWORpPcAx0fEZ/Lw9sBewGjg1IiY18h8XOMob1R+PStpLeAS4BfAf0rap9LIljecdIBYEhELgNNIO9/7Jb2v0sjqWwNYD1iUb7F/AXCZpO9L2q7a0JYzhLR+h+bE9kvSD1c3lXRApZEtbzXaZ5sFuAPYWNInASLiGuAmYBvSAbll5DaCDwKjJV0OXEj6YfVawNHAthWGV899pG30/wBExB3Ab4GlwBbQ2KVWJ44GSZooaaOI+DXwr8CNpKre1cDhwCPApHwmWilJIyQNycniauAwSetFxNPARXmyvaqLcFk53qERcQnpzsfzgfuBPwBfB0YCn6swxOXks8jLga9K2iwP30jaMXfq8cODLCL+m9bfZneUtLOknXONYhqwU1dbV0RMB5YAX6ouyr8rxDsxIl4iPRJiW+DmiPheRBxB2o4/X2mgvNXGOQQgIl4Hjiet2wNy2UzgGfI+1kiDuRNHAyTtTtrpLpQ0PCIuIl3+OTcizomIx4GfkjJ2pWdEuVo/FZieGxQvI911+EBJ60fEU8DpwK6SKj97q4l3bET8DPgCcFFE/DQi7gZOBv5B0riKY50k6VhJ38rr9gJgFnC0pHdGxAvAz4H3SerXbasHONYxeZv9F1pzm/0I6bLfx0j72GHArcDzwF6SPpMnvQtY0nUQrEpNvOdLOiYiXgbeFREnFOKbS3r0w9AKY92XdAyYXIhjJnA98DFJ/5bLngCGSFq5oRlHhF89vEhn5jNI2fjHwM6FcSq8/yTwF2BMhbHuStq5JgLnARcWyk8hJYytgU8DfwZGVrxui/GemzfwIXnckMJ0BwB/BEZVGOsuwN3A3sD3SJdT3g/8A/BN0oFkZ1Jbx23A6BaJ9bvAbGBinekq3WZJJzQrk2oXB+Sy7Ug1zSNIl6YOJd389ArSs3q2rXC91ot3AnAd8FVglcK0XyAdoLeuMN6NgXtIbS0n521zWB43Btgtx3g58Bjw3obnXdU/1Q4vYEtSZv5gHj4fOLPOdEfmnXObiuP9DnBCfr8V6ez3GGATUlX6K6RrxdcD27fA+q2Ndxrp0tTmpLYDgMNycql63X4L+EZheDrwJ+A9wEqkGug1pB5sla7bOrFeDdzSdWAgXWk4ohW22RzP1/K2MDIPbwPcAPxzHh5Fuvy3XtWxdhPv1qQTmy/m4XcDvwbeU3GcI/N+Pzpvn6fn5LFSYZphwDuBtcvM25eqevYs8IWI+HMePg74oKS32geUfqT0IjA5IuZUEGNRV8PX0cC1wAJgM9IGszgiTgH2BfaO1ChWtdp4O0lPgTwdeHdet6OAf2qBdfsksHLuNQPpTO054GfA8Ig4E9gf2K8F1m1trHcAi4BzJK0aqSfVS7TGNgspgY0BNpM0LMf0FeA4SdtHxIsRcWtEtMqzeWrjvYcU7zG5K/GDpHV7d5VBRmp7mR0Rz5Bq9A8A7wM+BZDb5pZGxNyIWFRm3u6O26C8gUyGGC0AAAUqSURBVCyVdCzwekSckhugW6Y7o6RNgA+TztjXj793Z/w+6WztoCrjq9VDvN8D3hERh7TKOs7dGE8mJbc1SVX+fSRNBa6J1HjbEnqJdXqkXkotRdIppIPxGaQurC9JOhM4JyJmVRvd8nqI92cRcWe10S1LkiIicieIQ4B1SCeUu5IupT1bdp4t9TyOVpF7+LxRLIuIpfntDGCqpF9HxH2DH92yirFGxCPAT5V+CfpFSetExELS7wxG1Pu/WjTee4HVqo63a/l5x7tb0r8D6wJrA7/Jk7XEDxVLxFpZQy2ApHeSktmciFjSVR4RX5H0A9IltCWSHifVjk+pJtKkneLtIdbI28XLwJmSrib1APtoX5IGuMaxDElbRMQD+X23By1JZwCvA1+p6my4p1glrQucTaqargx8CDioyqpzO8VbYjuYQuqp9KmIeGgwYyzE0E6xfhz4f8DTpMuoJ0fEnNxT8fU8za6kg9oWpPbEv1YRa7vF20OsinyQz72q3klqm9sjImb3eYGD0UjTDi/g48ArpG6gXWVDa6bparD9AOlSSsvFyt9PBiaQenb8B7Blq67bVou3we1gGKmqfyUwwbE2FOsHSG1a2+Xh/yI94bNr/JCa6YdVFWu7xdtbrHWmH9ffZbrGwVu3j7iCtHN9gLQRfC6PGxb5MlXuD/90dZGWinVEFKqrVWmneEvEukZEPC9plYhY7FgbivcDwBYRMS0PjyPdXuYzEfFqLtsBWCcirimeKTveAYm1g3SyOyDtW04cmdLdY18g3cvpbNJtOj5XGP9e0u8NftYCB7jeYp1Aug3Cz4BXq9wBczxtE2+D28GHgKldO2VV2izWocBqEfFCfr8uqcvqhyOiU9IGpMba6yLd8aBS7RRvFbE6cdShdHPAqaQurJ+TtC2p58+fomS3tWZrp1ihveJ1rM0haRgp2V0dEbsr3VV4O+CkaKEbAnZpp3gHK1Ynjm4o3Y7jFNJlgCHAhyJifrVR1ddOsUJ7xetYm0fSNNL9nD4MHBr9aawdBO0Ub7NjdXfcbkTEU5JmAx8F9mzlHbCdYoX2itexDjxJIt1ZeGL+u3tEPFhtVN1rp3gHK1Ynjm4o3X56L9J1wkp/AdqbdooV2itexzrwchvWa5L+A7i9VQ/CXdop3sGK1ZeqetAKPX0a1U6xQnvF61ibo+qeU2W1U7zNjtWJw8zMSvFNDs3MrBQnDjMzK8WJw8zMSnHiMDOzUpw4zBokaT9JIWmrPDxeUo8PQmpkGrN248Rh1rjJpGe1H1h1IGZVcuIwa4CkkcDOpGegL5c4JB0i6WpJ10q6X9KJhdFDJZ0j6R5Jv5e0Sv7M4ZJul3SXpCskrTo4/41Z/zhxmDVmX+DaSA9NekbS9nWm2RH4LOnZIp/Ot7KGdLPBMyNia9Jzyj+Vy6+MiB0i4r2kpx4e1tT/wGyAOHGYNWYycEl+f0kernVdRDydn3txJelW8QCPxN+fmz0TGJ/fbyPpT5LuJiWcrZsSudkA872qzHqRb1m+G+lAH6TndgfpSWtFtbdh6BouPgvjDWCV/H4asG9E3CXpEGCXgYvarHlc4zDr3f7ABRGxcUSMj4gNgUeADWqm21PS6NyGsS/wl17mOwqYL2k4qcZh1hacOMx6Nxm4qqbsCuDYmrI/Az8HZgFXRMSMXub7TeBW4DrSM6PN2oJvcmg2APKlpo6IOKrqWMyazTUOMzMrxTUOMzMrxTUOMzMrxYnDzMxKceIwM7NSnDjMzKwUJw4zMyvlfwH7xm4tvQB/oAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(top_alpha, color = 'skyblue', kde = False)\n",
    "plt.xticks(rotation = 45)\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Top Performing Alpha, 100 Bootstrap Samples')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Alpha</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0.025</td>\n",
       "      <td>1.828874</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>0.500</td>\n",
       "      <td>1.929954</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>0.975</td>\n",
       "      <td>2.024713</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          Alpha\n",
       "0.025  1.828874\n",
       "0.500  1.929954\n",
       "0.975  2.024713"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(df_results_sim['Alpha'].quantile([.025,.50,.975]))"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}