jhrcook · January 27, 2020 01:51 · jhrcook · Jan 27, 2020
diff --git a/EXAMPLE_CODE_FOR_BEN.ipynb b/EXAMPLE_CODE_FOR_BEN.ipynb
 {
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boostrapping confidence intervals\n",
    "\n",
    "**resource**: http://www2.stat.duke.edu/~ar182/rr/examples-gallery/BootstrapConfidenceIntervals.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.stats.api as sms\n",
    "import plotly.express as px\n",
    "\n",
    "np.random.seed(1)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create mock data\n",
    "\n",
    "I created mock data to replicate what is on one Excel spread sheet.\n",
    "I used 4 WT and 5 *KRAS* mutant mice and created data for 4 proteins.\n",
    "The first two proteins and second two proteins have the same real average (`real_wt_means` and `real_kras_means`), but the first of each pair has a smaller standard deviation (`std_devs`).\n",
    "These values are not known to us in reality, but we will use bootrapping to estimate the distributions and use that to establish confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>WT0</th>\n",
       "      <th>WT1</th>\n",
       "      <th>WT2</th>\n",
       "      <th>WT3</th>\n",
       "      <th>KRAS0</th>\n",
       "      <th>KRAS1</th>\n",
       "      <th>KRAS2</th>\n",
       "      <th>KRAS3</th>\n",
       "      <th>KRAS4</th>\n",
       "      <th>protein_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140.608634</td>\n",
       "      <td>84.706090</td>\n",
       "      <td>86.795706</td>\n",
       "      <td>73.175784</td>\n",
       "      <td>121.635191</td>\n",
       "      <td>42.461533</td>\n",
       "      <td>143.620294</td>\n",
       "      <td>80.969827</td>\n",
       "      <td>107.975977</td>\n",
       "      <td>proteinA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>82.544074</td>\n",
       "      <td>202.347556</td>\n",
       "      <td>44.209850</td>\n",
       "      <td>77.430796</td>\n",
       "      <td>73.116195</td>\n",
       "      <td>179.363861</td>\n",
       "      <td>23.007611</td>\n",
       "      <td>87.930025</td>\n",
       "      <td>38.549911</td>\n",
       "      <td>proteinB</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>300.844275</td>\n",
       "      <td>311.656304</td>\n",
       "      <td>277.987616</td>\n",
       "      <td>322.894474</td>\n",
       "      <td>218.031814</td>\n",
       "      <td>210.049887</td>\n",
       "      <td>218.017119</td>\n",
       "      <td>186.325443</td>\n",
       "      <td>197.542195</td>\n",
       "      <td>proteinC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>234.496140</td>\n",
       "      <td>281.247834</td>\n",
       "      <td>337.124883</td>\n",
       "      <td>251.583747</td>\n",
       "      <td>172.227253</td>\n",
       "      <td>151.897911</td>\n",
       "      <td>140.835605</td>\n",
       "      <td>153.012771</td>\n",
       "      <td>199.113478</td>\n",
       "      <td>proteinD</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          WT0         WT1         WT2         WT3       KRAS0       KRAS1  \\\n",
       "0  140.608634   84.706090   86.795706   73.175784  121.635191   42.461533   \n",
       "1   82.544074  202.347556   44.209850   77.430796   73.116195  179.363861   \n",
       "2  300.844275  311.656304  277.987616  322.894474  218.031814  210.049887   \n",
       "3  234.496140  281.247834  337.124883  251.583747  172.227253  151.897911   \n",
       "\n",
       "        KRAS2       KRAS3       KRAS4 protein_name  \n",
       "0  143.620294   80.969827  107.975977     proteinA  \n",
       "1   23.007611   87.930025   38.549911     proteinB  \n",
       "2  218.017119  186.325443  197.542195     proteinC  \n",
       "3  140.835605  153.012771  199.113478     proteinD  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_wt = 4\n",
    "n_kras = 5\n",
    "real_wt_means =   [100, 100, 300, 300]\n",
    "real_kras_means = [100, 100, 200, 200]\n",
    "std_devs = [25, 70, 20, 70]\n",
    "n_proteins = len(real_wt_means)\n",
    "\n",
    "protein_data = []\n",
    "\n",
    "for avg_wt, avg_kras, sd in zip(real_wt_means, real_kras_means, std_devs):\n",
    "    df = pd.DataFrame()\n",
    "    for i in range(n_wt):\n",
    "        df['WT' + str(i)] = np.abs(np.random.normal(avg_wt, sd, 1))\n",
    "    \n",
    "    for i in range(n_kras):\n",
    "        df['KRAS' + str(i)] = np.abs(np.random.normal(avg_kras, sd, 1))\n",
    "    \n",
    "    protein_data.append(df)\n",
    "\n",
    "protein_data = pd.concat(protein_data).reset_index(drop=True)\n",
    "protein_data['protein_name'] = ['protein' + x for x in ['A', 'B', 'C', 'D']]\n",
    "protein_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To help you understand what these data mean, I show below the real normal distributions that the measurements we have were pulled from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAKTCAYAAAC3hhQdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdfbhVdZn/8fdHPIIImgfQeBAPKaH4hMn4AOlY2ViNSaNWYpM6meg11ijplDZT0a/RaRrTfjP2RGlQoyhlGpFTImrkQCL6IxXRkVHCI6iAOEAJCdy/P9b34Oaw9zn7POyz9j7787qufZ21v+vpXotzbu71XU+KCMzMzMwsP3vkHYCZmZlZvXNBZmZmZpYzF2RmZmZmOXNBZmZmZpYzF2RmZmZmOXNBZmZmZpYzF2TWZZKWSTq1h9Y1TdJ/pOGRkjZL6tNNy/6OpC+k4VMlNXfHctPyTpb0THctz8y6j3NYWct2DqswF2S2G0kh6dByp4+IIyLiwQ6u40FJGyT17XCAb653VUQMiIjt7azrQkkPlbG8SyPiK52Np9U6d9mHEfGbiBjTHcs2s7ZVModJWinp9VRIbZD0C0kHdSZO5zAr5IKszkjaswpiaAJOBgI4M9dgku46QjWzyqqGHAZ8MCIGAEOBl4F/zzke57BewAVZL5CO2K6R9FQ6YvuBpH5p3KmSmiV9TtJLwA9S+8WSVkh6VdIcScNS+4K02N+lI8CPpvYzJC2V9JqkhZKObrX+09LwNEmzJf1Q0qZ0KmB8q5DPB34LzAAuaGfbRkn6dVrWPGBwwbimdBS3Z/p+oaTn0rTPS/qYpMOB7wAnpe15LU07Q9K3Jd0j6Q/Au1LbP7Va/+clrUvb+LGC9gclfbLg+84j2GL7sPXpA0mHp2W8lvbRmQXjZkj6Zjry3iTpYUmHtLWfzGpZDeYwACJiC/ATYGwb2+Yc5hxWFhdkvcfHgNOBQ4C3A/9YMO6tQCNwMDBF0ruBfwY+QnaE93vgdoCIOCXNc0zqSr9D0juAW4BLgEHAd4E5Kn268cy0vLcAc4CbWo0/H7g1fU6XdGAb23Ub8ChZEvsKJQo4SfsA/wa8PyIGAhOApRGxHLgUWJS25y0Fs50HXAsMBIqdDnhrWu/wtN7pktrtsi+2D1vF2gD8HLgXOAD4NHBrq2VPBr4M7A+sSHGa9Wa1lMMAkNQf+CjZAWYpzmHOYWVxQdZ73BQRL0TEq2S/+JMLxu0AvhQRWyPidbLEd0tEPBYRW4FryI6+mkos+2LguxHxcERsj4iZwFbgxBLTPxQR96TrIn4EHNMyQtI7yZLq7Ih4FPgfsqSyG0kjgT8DvpBiX0CWBErZARwpae+IWBMRy9qYFuBnEfFfEbEjHekW07LuXwO/IPsPoKtOBAYAX42IP0XE/cBcdv03+2lELI6IbWSF67huWK9ZNauJHJbcnXqqNgLvBf612EKcw5zDOsIFWe/xQsHw74FhBd/XtvpjHZamASAiNgPryY6iijkYuDJ1Tb+WEtFBrdZR6KWC4T8C/fTmdR8XAPdGxLr0/TZKn7YcBmyIiD+02rbdpGk+SnYkuSZ1lR9WYrktXmhnfLF1l9rmjhgGvBARO1otu3D/t96HA7phvWbVrFZyGMCHUk9VX+BTwK8lvbXIcpzDMs5hZXBB1nsU3uUzElhd8D1aTbuaLEEBO7vKBwEvllj2C8C1EfGWgk//iJjVkQAl7U12dPbnkl5K14NMBY6R1PoIFGANsH+Kr3DbioqIX0XEe8lOYTwNfK9lVKlZ2gm52Lpb9usfgP4F44ol41JWAwdJKvz7G0np/W9WD6o+h7WWett+CmwH3llkEucwK5sLst7jMkkjJDUCnwfuaGPa24C/kTQuXUNxHfBwRKxM418G3lYw/feASyWdoMw+kv5S0sAOxvghssQ1lqz7ehxwOPAbsuvKdhERvweWAF+WtFc63fnBYguWdKCkM1Py2QpsTutq2Z4RkvbqYLwUrPtk4Azgx6l9KXCWpP7Kbg2/qNV8rfdhoYfJkuFnJTUoe/7RB0nXwJjVqVrIYbtIy5pEdp3U8tbjncOsI1yQ9R63kV1g+Vz6/FOpCSNiPvAF4E6yI7hDgHMLJpkGzExd+x+JiCVk12DcBGwgu0Dzwk7EeAHwg/TsnZdaPmm5H1Px29nPA04AXgW+BPywxLL3AK4kO3J7Ffhz4G/TuPuBZcBLktYVn72ol8i2dzXZNRCXRsTTadyNwJ/IktbMNL7QNAr2YeGIiPgT2UXD7wfWAd8Czi9Ytlk9qoUc1uLnkjaTXUN2LXBBG9d7OYdZWRTRXo+nVTtJK4FPRsR9ecdiZtZRzmFm7iEzMzMzy50LMjMzM7OcVeyUpbKnLC8guy14T+AnEfGldMHmHUATsBL4SERsSPNcQ3Zh4Xbg7yLiVxUJzszMzKyKVLIgE7BPRGxOT/V9CLgcOAt4NSK+KulqYP+I+JykscAs4HiyZ5zcB7w92nnpqpmZmVmtq9gpy8hsTl8b0ieASWR3dJB+figNTwJuT08Ufp7sLpjjKxWfmZmZWbUo9piBbqPs7fOPAocC34yIhyUdGBFrACJijaQD0uTD2fV9YM2UfuoyAIMHD46mpqbuD9zMqtajjz66LiKG5B1Hd3AOM6svbeWvihZk6XTjOElvAe6SdGQbk6vYInabSJoCTAEYOXIkS5Ys6ZZYzaw2SCr66pla1NTU5BxmVkfayl89cpdlRLwGPAi8D3hZ0tAU2FDglTRZM7u+OmMEu746o2VZ0yNifESMHzKkVxwkm5mZWZ2rWEEmaUjqGWt5h+FpZO/mmsObL5O+APhZGp4DnCupr6RRwGhgcaXiMzMzM6sWlTxlOZTstQt9yAq/2RExV9IiYLaki4BVwIcBImKZpNnAU8A24DLfYWlmZmb1oGIFWUQ8DhxbpH098J4S81xL9l4wM0veeOMNmpub2bJlS96h9Kh+/foxYsQIGhoa8g7FzLqgHnNYZ/JXRS/qN7Oua25uZuDAgTQ1NZE93q/3iwjWr19Pc3Mzo0aNyjscM+uCesthnc1ffnWSWZXbsmULgwYNqotE1kISgwYNqqsjarPeqt5yWGfzlwsysxpQL4msUD1us1lvVW9/z53ZXhdkZtauqVOn8o1vfGPn99NPP51PfvKTO7+fffbZ7LvvvowbN47GxkZGjRrFuHHjOO200/II18xsp1rJX76GzKzGTJ/evcubMqX9aSZMmMCPf/xjrrjiCnbs2MG6devYuHHjzvGrV69m3rx5nHDCCVx44YWcccYZnHPOOd0bqJn1Cj2dw2olf7mHzMzaNXHiRBYuXAjAsmXLOPLIIxk4cCAbNmxg69atLF++nGOP3e2majOz3NVK/nIPmZm1a9iwYey5556sWrWKhQsXctJJJ/Hiiy+yaNEi9ttvP44++mj22muvvMM0M9tNreQvF2RWcd3dPV2onNNt1j1ajjIXLlzIZz7zGV588UUWLlzIfvvtx4QJE/IOz6xbVTJvgXNXT6uF/OVTlmZWlgkTJrBw4UKeeOIJjjzySE488UQWLVrEwoULmThxYt7hmZmVVAv5ywWZmZVl4sSJzJ07l8bGRvr06UNjYyOvvfYaixYt4qSTTso7PDOzkmohf7kgM7OyHHXUUaxbt44TTzxxl7b99tuPwYMH5xiZmVnbaiF/+RoysxqT17Unffr02eVWcYAZM2bsNl2xNjOzFnnksFrIX+4hMzMzM8uZCzIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMzM8uZCzIzMzOznLkgM7N2DRgwYOfwPffcw+jRo1m1ahXTpk1j+PDhjBs3jrFjxzJr1qxd5tu2bRuDBw/mmmuu2aV97ty5HHvssRxzzDGMHTuW7373uz2yHWZWn2ohh/k5ZGa1ZkU3v2Tv0PIfCjR//nw+/elPc++99zJy5EgApk6dylVXXcWzzz7LcccdxznnnENDQwMA9957L2PGjGH27Nlcd911SOKNN95gypQpLF68mBEjRrB161ZWrlzZvdtkZtXLOawo95CZWVl+85vfcPHFF/OLX/yCQw45ZLfxo0ePpn///mzYsGFn26xZs7j88ssZOXIkv/3tbwHYtGkT27ZtY9CgQQD07duXMWPG9MxGlEnSQZIekLRc0jJJl6f2aZJelLQ0fT5QMM81klZIekbS6flFb2bFVHsOc0FmZu3aunUrkyZN4u677+awww4rOs1jjz3G6NGjOeCAAwB4/fXXmT9/PmeccQaTJ0/eeSqgsbGRM888k4MPPpjJkydz6623smPHjh7bljJtA66MiMOBE4HLJI1N426MiHHpcw9AGncucATwPuBbkvrkEbiZ7a4WcpgLMjNrV0NDAxMmTODmm2/ebdyNN97ImDFjOOGEE5g2bdrO9rlz5/Kud72L/v37c/bZZ3PXXXexfft2AL7//e8zf/58jj/+eK6//no+8YlP9NSmlCUi1kTEY2l4E7AcGN7GLJOA2yNia0Q8D6wAjq98pGZWjlrIYS7IzKxde+yxB7Nnz+aRRx7huuuu22Xc1KlTeeaZZ7jjjjs4//zz2bJlC5B19d933300NTVx3HHHsX79eh544IGd8x111FFMnTqVefPmceedd/bo9nSEpCbgWODh1PQpSY9LukXS/qltOPBCwWzNlCjgJE2RtETSkrVr11YoajMrVAs5zAWZmZWlf//+zJ07l1tvvbXoUeZZZ53F+PHjmTlzJhs3buShhx5i1apVrFy5kpUrV/LNb36TWbNmsXnzZh588MGd8y1dupSDDz64B7ekfJIGAHcCV0TERuDbwCHAOGAN8PWWSYvMHsWWGRHTI2J8RIwfMmRIBaI2s2KqPYf5LkszK1tjYyO//OUvOeWUUxg8ePBu47/4xS9y3nnn0dDQwLvf/W769u27c9ykSZP47Gc/yw033MDXvvY1LrnkEvbee2/22WcfZsyY0YNbUR5JDWTF2K0R8VOAiHi5YPz3gLnpazNwUMHsI4DVPRSqmZWpmnOYIooexNWE8ePHx5IlS/IOw9oxvZvvcC40pfy7nWvW8uXLOfzww/MOIxfFtl3SoxExvpLrlSRgJvBqRFxR0D40Itak4anACRFxrqQjgNvIrhsbBswHRkfE9rbW4xxWvSqZt6A+cleLes1hHc1f7iEzM9vdRODjwBOSlqa2zwOTJY0jOx25ErgEICKWSZoNPEV2h+Zl7RVjZmaFXJCZmbUSEQ9R/Lqwe9qY51rg2ooFZWa9mi/qNzMzM8uZCzKzGlDL13p2Vj1us1lvVW9/z53ZXhdkZlWuX79+rF+/vq4SWkSwfv16+vXrl3coZtZF9ZbDOpu/fA2ZWZUbMWIEzc3N1NtDRPv168eIESPyDsPMuqgec1hn8pcLMrMq19DQwKhRo/IOw8ysU5zDyuNTlmZmZmY5c0FmZmZmljMXZGZmZmY5q1hBJukgSQ9IWi5pmaTLU/s0SS9KWpo+HyiY5xpJKyQ9I+n0SsVmZmZmVk0qeVH/NuDKiHhM0kDgUUnz0rgbI+L6wokljQXOBY4gexfcfZLe7tePmJmZWW9XsR6yiFgTEY+l4U3AcmB4G7NMAm6PiK0R8TywguxFvWZmZma9Wo9cQyapCTgWeDg1fUrS45JukbR/ahsOvFAwWzNFCjhJUyQtkbSknp5pYmZmZr1XxQsySQOAO4ErImIj8G3gEGAcsAb4esukRWbf7bG+ETE9IsZHxPghQ4ZUKGozMzOznlPRB8NKaiArxm6NiJ8CRMTLBeO/B8xNX5uBgwpmHwGsrmR8ZmZmu1gxHYDDKv3Y9BVdmPfQKd0WhlWPSt5lKeBmYHlE3FDQPrRgsr8CnkzDc4BzJfWVNAoYDSyuVHxmZmZm1aKSxwATgY8DT0hamto+D0yWNI7sdORK4BKAiFgmaTbwFNkdmpf5DkszMzOrBxUryCLiIYpfF3ZPG/NcC1xbqZjMzMzMqpFfLm5dk663aEtnrsV4epuvkTCz3mvBgs7P+/T9bY+f4vRZk/zqJDMzM7OcuSAzMzMzy5kLMjMzM7OcuSAzMzMzy5kLMjMzM7OcuSAzMzMzy5kLMjMzM7OcuSAzMzMzy5kLMjMzM7OcuSAzMzMzy5kLMjMzM7OcuSAzMzMzy5kLMjMzM7Oc7Zl3AGbFHLbn9PImXNGBhR46pVOxmJmZVZp7yMzMzMxy5oLMzMzMLGcuyMzMzMxy5oLMzMzMLGe+qN9q2oIF5U/79P3lTzvF1/+bmVkPcg+ZmZmZWc5ckJmZmZnlzAWZmVkrkg6S9ICk5ZKWSbo8tTdKmifp2fRz/4J5rpG0QtIzkk7PL3ozq0W+hswAmF7mc1hbO8y/QdY7bQOujIjHJA0EHpU0D7gQmB8RX5V0NXA18DlJY4FzgSOAYcB9kt4eEdtzit/Maoz/OzUzayUi1gBr0vAmScuB4cAk4NQ02UzgQeBzqf32iNgKPC9pBXA8sKhnI68vnT2QbIsPMi0vPmVpZtYGSU3AscDDwIGpWGsp2g5Ikw0HXiiYrTm1mZmVxQWZmVkJkgYAdwJXRMTGtiYt0hYlljlF0hJJS9auXdsdYZpZL+CCzMysCEkNZMXYrRHx09T8sqShafxQ4JXU3gwcVDD7CGB1seVGxPSIGB8R44cMGVKZ4M2s5rggMzNrRZKAm4HlEXFDwag5wAVp+ALgZwXt50rqK2kUMBpY3FPxmlnt8+WLZma7mwh8HHhC0tLU9nngq8BsSRcBq4APA0TEMkmzgafI7tC8zHdYmllHuCAzM2slIh6i+HVhAO8pMc+1wLUVC8rMejWfsjQzMzPLmQsyMzMzs5y5IDMzMzPLma8hMzOzmnTYnhV4VL9ZTtxDZmZmZpazihVkkg6S9ICk5ZKWSbo8tTdKmifp2fRz/4J5rpG0QtIzkk6vVGxmZmZm1aSSPWTbgCsj4nDgROAySWOBq4H5ETEamJ++k8adCxwBvA/4lqQ+FYzPzMzMrCpUrCCLiDUR8Vga3gQsJ3vZ7iRgZppsJvChNDwJuD0itkbE88AK4PhKxWdmZmZWLXrkGjJJTcCxwMPAgRGxBrKiDTggTTYceKFgtubU1npZfjGvmZmZ9SoVL8gkDSB7Qe8VEbGxrUmLtMVuDX4xr5mZmfUyFS3IJDWQFWO3RsRPU/PLkoam8UOBV1J7M3BQwewjgNWVjM/MzMysGlTyLksBNwPLI+KGglFzgAvS8AXAzwraz5XUV9IoYDSwuFLxmZmZmVWLSj4YdiLwceAJSUtT2+eBrwKzJV0ErAI+DBARyyTNBp4iu0PzsojYXsH4zMzMzKpCxQqyiHiI4teFAbynxDzXAtdWKiYzMzOzauQn9ZuZmZnlzAWZmZmZWc5ckJmZmZnlzAWZmZmZWc5ckJmZmZnlzAWZmZmZWc5ckJmZmZnlzAWZmZmZWc5ckJmZmZnlzAWZmZmZWc5ckJmZmZnlzAWZmZmZWc5ckJmZmZnlzAWZmZmZWc5ckJmZmZnlzAWZmZmZWc72zDsA60ErppccdZh/E8zMzHLjHjIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMrQtItkl6R9GRB2zRJL0pamj4fKBh3jaQVkp6RdHo+UZtZrSqrIJM0v5w2M7Nq04X8NQN4X5H2GyNiXPrck5Y3FjgXOCLN8y1JfToftZnVmzbfYCipH9AfGCxpf0Bp1L7AsArHZmbWaV3NXxGxQFJTmaubBNweEVuB5yWtAI4HFnU0bjOrT+29UvoS4Aqy5PUobya0jcA3KxiXmVlXVSp/fUrS+cAS4MqI2AAMB35bME1zatuNpCnAFICRI0d2IQwz603aPGUZEf83IkYBV0XE2yJiVPocExE39VCMZmYdVqH89W3gEGAcsAb4empXkWmjRFzTI2J8RIwfMmRIJ8Mws96mvR4yACLi3yVNAJoK54mIH1YoLjOzbtGd+SsiXm4ZlvQ9YG762gwcVDDpCGB1Z+I1s/pUVkEm6UdkR4VLge2pOQAXZGZW1bozf0kaGhFr0te/AlruwJwD3CbpBrJTpKOBxV2J28zqS1kFGTAeGBsRRbvgzcyqWKfyl6RZwKlkNwU0A18CTpU0jqygW0l2nRoRsUzSbOApYBtwWURsL7ZcM7Niyi3IngTeSnbNhJlZLelU/oqIyUWab25j+muBazsWmplZptyCbDDwlKTFwNaWxog4syJRmZl1H+cvM6t65RZk0zq6YEm3AGcAr0TEkaltGnAxsDZN9vmCByteA1xEdo3H30XErzq6TjOzIqblHYCZWXvKvcvy151Y9gzgJna/cPbGiLi+sKHVU66HAfdJeruvwTCzrupk/jIz61Hlvjppk6SN6bNF0nZJG9uaJyIWAK+WGcfOp1xHxPNAy1Ouzcy6pDP5y8ysp5XbQzaw8LukD9H5gslPuTazHtPN+cvMrCLK6iFrLSLuBt7diVn9lGszy1UX8peZWcWU+2DYswq+7kH2XJ8OP5PMT7k2s57WXfnLzKySyr3L8oMFw9vIHog4qaMr81OuzSwH3ZK/zMwqqdxryP6mowv2U67NrBp0Jn+ZmfW0ck9ZjgD+HZhIVkw9BFweEc2l5vFTrs2sGnQmf5mZ9bRyL+r/AdlpxWFkdz/+PLWZmVU75y8zq3rlFmRDIuIHEbEtfWYAvsXRzGqB85eZVb1yC7J1kv5aUp/0+WtgfSUDMzPrJs5fZlb1yi3IPgF8BHiJ7Plh5wC+UNbMaoHzl5lVvXIfe/EV4IL0VH0kNQLXkyU6M7Nq5vxlZlWv3B6yo1uSGUBEvAocW5mQzMy6lfOXmVW9cguyPSTt3/IlHWGW27tmZpYn5y8zq3rlJqWvAwsl/YTsOT4fwc8MM7Pa4PxlZlWv3Cf1/1DSErIX8go4KyKeqmhkZmbdwPnLzGpB2d32KYE5iZlZzXH+MrNq5+sorG4ctuf08ideUeZ0h07pVCxmZmaFyr2o38zMzMwqxAWZmZmZWc5ckJmZmZnlzNeQmZlZxUzvwKWbHXWY/wezXsQ9ZGZmZmY5c0FmZmZmljMXZGZmZmY5c0FmZmZmljMXZGZmZmY58z0qZmZmNaTdt46U+6aRQn7rSO7cQ2ZmZmaWMxdkZmZmZjlzQWZmZmaWMxdkZmZmZjnzRf1mZma9yIIFHZ/n6fvLm26Kr/2vGPeQmZmZmeXMBZmZmZlZzlyQmZmZmeXMBZmZWRGSbpH0iqQnC9oaJc2T9Gz6uX/BuGskrZD0jKTT84nazGqVCzIzs+JmAO9r1XY1MD8iRgPz03ckjQXOBY5I83xLUp+eC9XMap0LMjOzIiJiAfBqq+ZJwMw0PBP4UEH77RGxNSKeJ3t5zfE9EqiZ9QouyMzMyndgRKwBSD8PSO3DgRcKpmtObWZmZfFzyMzMuk5F2qLohNIUYArAyJEjKxlTVWn3hdhmdc49ZGZm5XtZ0lCA9POV1N4MHFQw3QhgdbEFRMT0iBgfEeOHDBlS0WDNrHa4h6zGTO/CQeZh/tc266o5wAXAV9PPnxW03ybpBmAYMBpYnEuEZlaTKtZD5lvGzayWSZoFLALGSGqWdBFZIfZeSc8C703fiYhlwGzgKeCXwGURsT2fyM2sFlWyz2QGcBPww4K2llvGvyrp6vT9c61uGR8G3Cfp7U5oZpaXiJhcYtR7Skx/LXBt5SIys96sYj1kvmXczMzMrDw9fVF/l28ZlzRF0hJJS9auXVvRYM3MzMx6QrXcZVn2LeO+Q8nMzMx6m54uyLp8y7iZmZlZb9PTBVnLLeOw+y3j50rqK2kUvmXczMzM6kjF7rJMt4yfCgyW1Ax8iewW8dnp9vFVwIchu2VcUsst49vwLeNmZmZWRypWkPmWcTMzM7PyVMtF/WZmZmZ1ywWZmZmZWc78dkOzIhYsKG+6p+/v+LKnTOn4PGZm1ru5h8zMzMwsZ+4hq1YrphdtPsz/YmZmZr2Oe8jMzMzMcuaCzMzMzCxnLsjMzMzMcuaCzMzMzCxnLsjMzMzMcuaCzMzMzCxnLsjMzMzMcuaCzMzMzCxnLsjMzMzMcubnvpt1wWF7Fn+jQptWlDHNoX7hpZlZPXEPmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnO/C5LMzOzOlf2e3nLeRcv+H28neAeMjMzM7OcuSAzMzMzy5kLMjMzM7OcuSAzMzMzy5kLMjMzM7OcuSAzMzMzy5kfe2Fm1kGSVgKbgO3AtogYL6kRuANoAlYCH4mIDXnFaGa1xT1kZmad866IGBcR49P3q4H5ETEamJ++m5mVxQWZmVn3mATMTMMzgQ/lGIuZ1ZhcCjJJKyU9IWmppCWprVHSPEnPpp/75xGbmVkZArhX0qOSWh5JfmBErAFIPw8oNqOkKZKWSFqydu3aHgrXzKpdnj1k7u43s1o1MSLeAbwfuEzSKeXOGBHTI2J8RIwfMmRI5SI0s5pSTacs3d1vZjUhIlann68AdwHHAy9LGgqQfr6SX4RmVmvyKsg63d1vZpYnSftIGtgyDPwF8CQwB7ggTXYB8LN8IjSzWpTXYy8mRsRqSQcA8yQ9Xe6MqYCbAjBy5MhKxWdmVsqBwF2SIMuht0XELyU9AsyWdBGwCvhwjjGaWY3JpSAr7O6XtEt3f0Ssaau7PyKmA9MBxo8fHz0Vs5kZQEQ8BxxTpH098J6ej8jMeoMeP2Xp7n4zMzOzXeXRQ+bufjMzM7MCPV6QubvfzMzMbFfV9NgLMzMzs7rkl4ubmdmuVkzvtkUd5v9lzMriHjIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMzM8uZCzIzMzOznPn+lwqY3g03KPnOJDPrKa1zlvOPWc9zD5mZmZlZzlyQmZmZmeXMHdNmZmZWlgULypvu6fs7vuwpUzo+T2/iHjIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMzM8uZCzIzMzOznLkgMzMzM8uZH3vRVSt2fyy/n3JtbSnntvHO3DIOvm3czKxWuYfMzMzMLKmyeLoAACAASURBVGcuyMzMzMxy5pNrZmZm1q0O23P3y3nataKd8Yf27msy3ENmZmZmljMXZGZmZmY5c0FmZmZmljMXZGZmZmY5c0FmZmZmljPfZWlWhTp1hxK0fZdSL79DycyslrmHzMzMzCxnLsjMzMzMcuaCzMzMzCxndXsN2fROXqLTml8kbmZmZl3lHjIzMzOznNVX/86KN7vF3LNlZr3Giq51+TsfmuXPPWRmZmZmOfNxkZmZmVW/LvYEF1VFz2d0QWZmVgPauhHJpxytN1iwoHLLPuWUyi27u1TdKUtJ75P0jKQVkq7OOx4zs3I5f5lZZ1XVcZWkPsA3gfcCzcAjkuZExFP5RmbWC/Ty7v68OX+ZWVdUVUEGHA+siIjnACTdDkwCnNDMclbsdMLT93d9uVN6T03n/GVWpUqdDq2mHFZtBdlw4IWC783ACTnFYlZzKnkNhrXL+cvMOq3aCjIVaYtdJpCmAC316GZJz3RxnYOBdV1cRi2vvxpiyHv91RBDja7/ki6v+JI3F1FuDAd3eaWV0W7+gorkMMj/96c71Po2OP78dWIbujWHlaNk/qq2gqwZOKjg+whgdeEEETEd6LaLYSQtiYjx3bW8Wlt/NcSQ9/qrIYZ6X3+1xNBF7eYv6P4cBr1i39X8Njj+/NX6NlTbXZaPAKMljZK0F3AuMCfnmMzMyuH8ZWadVlU9ZBGxTdKngF8BfYBbImJZzmGZmbXL+cvMuqKqCjKAiLgHuKcHV1mBZwHU1Poh/xjyXj/kH0O9rx+qI4YuySF/taj5fUftb4Pjz19Nb4Midrvm1MzMzMx6ULVdQ2ZmZmZWd+qmIJN0kKQHJC2XtEzS5al9mqQXJS1Nnw9UOI6Vkp5I61qS2holzZP0bPq5f4XWPaZgO5dK2ijpikrvA0m3SHpF0pMFbSW3WdI16dUzz0g6vULr/1dJT0t6XNJdkt6S2pskvV6wL75TofWX3Ofdvf1txHBHwfpXSlqa2iuxD0r9/fXY70FvIenDaR/ukDS+1bii+0zScSnvrJD0b5KKPaIjF6qR103lnce6qtb/BiX1k7RY0u9S/F9O7TURf1kioi4+wFDgHWl4IPDfwFhgGnBVD8axEhjcqu1rwNVp+GrgX3ogjj7AS2TPRKnoPgBOAd4BPNneNqd/k98BfYFRwP8AfSqw/r8A9kzD/1Kw/qbC6Sq4/UX3eSW2v1QMrcZ/HfhiBfdBqb+/Hvs96C0f4HBgDPAgML6c3x1gMXAS2bPS/hN4f97bkeLqk+J8G7BXin9s3nGViDXXPNYN8df032D63R2QhhuAh4ETayX+cj5100MWEWsi4rE0vAlYTvZk7WowCZiZhmcCH+qBdb4H+J+I+H2lVxQRC4BXWzWX2uZJwO0RsTUingdWkL2SplvXHxH3RsS29PW3ZM+MqogS219Kt29/ezGk3pKPALO6up421l/q76/Hfg96i4hYHhHFHiZbdJ9JGgrsGxGLIvuf6of0TI4px87XTUXEn4CW101VnbzzWFfV+t9gZDanrw3pE9RI/OWom4KskKQm4FiyChvgU+nU1S2q0OnCAgHcK+lRZU/sBjgwItZA9kcDHFDhGCB7RlLhf8A9uQ+g9DYXe/1MpQvnT5D1GrQYJen/Sfq1pJMruN5i+zyP7T8ZeDkini1oq9g+aPX3V02/B7Wu1D4bnoZbt1eDWv93rsnf31r9G5TUJ11a8QowLyJqKv721F1BJmkAcCdwRURsBL4NHAKMA9aQnbqppIkR8Q7g/cBlkk6p8Pp2o+yhlWcCP05NPb0P2lLW62e6bWXSPwDbgFtT0xpgZEQcC3wGuE3SvhVYdal93qPbn0xm1+K8YvugyN9fyUmLtNXNLeGS7pP0ZJFPW71HpfZZNe/Lao6tK6p2u2r5bzAitkfEOLIzGsdLOrKNyasu/vZU3XPIKklSA9kv4q0R8VOAiHi5YPz3gLmVjCEiVqefr0i6i6wL9WVJQyNiTTq98EolYyArBh9r2fae3gdJqW0u6/Uz3UHSBcAZwHvSqRwiYiuwNQ0/Kul/gLcDS7pz3W3s8x7b/rTuPYGzgOMKYqvIPij290cV/B5Uo4g4rROzldpnzex6Sr6a9mWt/zvX1O9vb/kbjIjXJD0IvI8ajL+UuukhS9fJ3Awsj4gbCtqHFkz2V8CTreftxhj2kTSwZZjswvInyV6vckGa7ALgZ5WKIdmlR6Qn90GBUts8BzhXUl9Jo4DRZBckdytJ7wM+B5wZEX8saB8iqU8aflta/3MVWH+pfd4j21/gNODpiNh5SqsS+6DU3x85/x70MkX3WTqNs0nSienf4Xwqn2PKVeuvm6qZ399a/xtMeanlbvi9SbmLGom/LHnfVdBTH+CdZN2VjwNL0+cDwI+AJ1L7HGBoBWN4G9ldH78DlgH/kNoHAfOBZ9PPxgrG0B9YD+xX0FbRfUBW/K0B3iA7armorW0G/oHsjphn6Ia7wUqsfwXZ9QUtvwvfSdOenf5tfgc8BnywQusvuc+7e/tLxZDaZwCXtpq2Evug1N9fj/0e9JYPWQHfTNaL+TLwq/b2GTCerOj/H+Am0kPBq+GTfg/+O8X2D3nH00acueaxboi/pv8GgaOB/5fif5I37wqvifjL+fhJ/WZmZmY5q5tTlmZmZmbVygWZmZmZWc5ckJmZmZnlzAWZmZmZWc5ckJmZmZnlzAWZVRVJD0oaX4HlnippQhnTnSnp6g4ue3P7U3VMJZZpZpXnHFa5ZfZ2dfWkfsuHpD3jzRd557WeU4HNwMK2lhERc6itB1OaWYU5h1lPcA9ZLyOpSdLTkr6f3n13q6TTJP2XpGclHZ+m2ye91PqR9BLpSQXz/0bSY+kzIbUPlbRA0tK03JNT++aCdZ8jaUYaniHpBkkPAP/Sxvr2lnS7spds3wHsXWK7Vkr6F0mL0+fQEutplHR3Wt5vJR2t7EW6lwJTU/wnp6c+35nieUTSxLS8CyXdVLDsf5O0UNJzks4pY///fVre45K+nNr+RdLfFkwzTdKVpaY3q2fOYc5hdSvvJ9P6070foInsZdlHkRXcjwK3kL1odRJwd5ruOuCv0/BbyJ6UvQ/Zk/z7pfbRwJI0fCVvvlmgDzAwDW8uWPc5wIw0PIPs/Yx92lnfZ4BbUvvRKfbxRbZrZcH6zwfmlljPvwNfSsPvBpam4WnAVQXLuw14ZxoeSfY6EYALgZsKlv3jtB/HAitK7PPN6edfANPTvt4jxXUKcCzw64Lpn0rrLDp96/3qjz/19HEOcw6r149PWfZOz0fEEwCSlgHzIyIkPUGW7CD7QzpT0lXpez+yP7DVwE2SxgHbyV4qDdk7525R9nLauyNiaRlx/DgitrezvlOAfwOIiMclPd7G8mYV/LyxxHreSfbqHyLifkmDJO1XZFmnAWMltXzfV+k9o63cHRE7gKckHdhGbJBt41+Qvd4DYAAwOiJulnSApGHAEGBDRKyS9HfFpgcWtLMes97OOQznsHrjgqx32lowvKPg+w7e/DcXcHZEPFM4o6RpZO/HO4bsiGcLQEQskHQK8JfAjyT9a0T8kOzdaC36tYrjD4WLLrE+Wi2jLVFiuPV62pqvxR7ASRHxepF4ChXuy2LLptX4f46I7xYZ9xOyo++3AreXMb1ZPXMOKz1fC+ewXsbXkNWvXwGfVvrrlXRsat8PWJOOqD5O1rWPpIOBVyLie8DNwDvS9C9LOlzSHmQvPe7o+hYAH0ttR5J1+Zfy0YKfi0pMU7i8U4F1EbER2AQUHj3eC3yq5Us6mu6qXwGfkDQgLXO4pAPSuNuBc8kS2k/KmN7M2uYc5hzWq7iHrH59BfgG8HhKMCuBM4BvAXdK+jDwAG8euZ0K/L2kN8ju9Dk/tV9Ndt3AC8CTZF3WHVnft4EfpG7+pcDiNmLuK+lhsgOJySWmmVawvD8CF6T2nwM/UXYh7qeBvwO+mabbkywJXtrGutsVEfdKOhxYlHL2ZuCvyf4TWJZOJ7wYEWvam74rcZjVCecw57BeRRHl9rSa5UfSSrILZdflHYuZWUc5h1l7fMrSzMzMLGfuITMzMzPLmXvIzMzMzHLmgszMzMwsZy7IzMzMzHLmgszMzMwsZy7IzMzMzHLmgszMzMwsZy7IrMskLUuv+OiJdV0o6aGC75slva2blv15Sd9Pw02SQlK3vM1C0sgUa5/uWJ6ZdY+ezF9pfdMk/Uca7ta8IOk7kr6Qhk+V1Nwdy03LO1nSM+1PaZ3lgsx2kwqRQ8udPiKOiIgHO7D88yQtSYlojaT/lPTOzsQaEQMi4rl21ldWYoqI6yLik52Jo8g6V0o6rWDZq1Ks27tj+WZWXCXzV/q7fl3SJkmvSVoo6dL0HswOKzcvtD4QbWN5l0bEVzoTS5F17rIfI+I3ETGmO5ZtxbkgqzPd1ePThfV/hux9cNcBBwIjyd49NynnuPxeV7MqVyV/px+MiIHAwcBXgc+Rvaw8V+59r30uyHqBdNR2jaSnJG2Q9ANJ/dK4UyU1S/qcpJeAH6T2iyWtkPSqpDmShqX2BWmxv0s9WB9N7WdIWlpwVHh0q/WfloanSZot6YfpKHKZpPFp3H7A/wEui4ifRsQfIuKNiPh5RPx9iW0blOLbKGkxcEir8TuP4iR9IO2DTZJelHSVpH2A/wSGpe3ZLGlYivMnkv5D0kbgwsJTCQU+IWl16sm7smC9MyT9U8H3nb1wkn5EVmj+PK3vs61PgaYY5qT9v0LSxQXLKrkPzXqbWslfrUXE/0bEHOCjwAWSjiyxfaMk/Totbx4wuGBc67xwoaTn0rTPS/qYspd3fwc4KW3Ta2naGZK+LekeSX8A3tU6L6XpPi9pXdrOjxW0PyjpkwXfd/bCFduPanWmQdLhaRmvpf10ZsG4GZK+KekXaVselrRL7rbduSDrPT4GnE5WsLwd+MeCcW8FGsmO6KZIejfwz8BHgKHA74HbASLilDTPMakr/Q5J7wBuAS4BBgHfBeZI6lsiljPT8t4CzAFuSu0nAf2AuzqwXd8EtqQ4P5E+pdwMXJKOXo8E7o+IPwDvB1an7RkQEavT9JOAn6Q4by2xzHcBo4G/AK5WwWnIUiLi48AqsiPpARHxtSKTzQKagWHAOcB1kt5TML7UPjTrjWohfxUVEYvJ/pZPLjHJbcCjZIXYV4ALik2UDh7/DXh/ymETgKURsRy4FFiUtuktBbOdB1wLDASKndJ8a1rv8LTe6ZLaPe1YbD+2irUB+DlwL3AA8Gng1lbLngx8GdgfWJHitDa4IOs9boqIFyLiVbJf/MkF43YAX4qIrRHxOlnyuyUiHouIrcA1ZEdfTSWWfTHw3Yh4OCK2R8RMYCtwYonpH4qIe9J1ET8Cjkntg4B1EbGtnA1S1gV/NvDF1Jv2JDCzjVneAMZK2jciNkTEY+2sYlFE3B0RO9J+KebLad1PkB2dTy4xXdkkHQS8E/hcRGyJiKXA94GPF0xWah+a9Ua1kL/aspqsaNyFpJHAnwFfSPEvICtkStkBHClp74hYExHL2lnvzyLiv1IO21JimpZ1/xr4BVkh21UnAgOAr0bEnyLifmAuu/67/TQiFqd8fyswrhvW26u5IOs9XigY/j1Zz0uLta3+WIelaQCIiM3AerKjqGIOBq5MXdOvpS7zg1qto9BLBcN/BPqlLvn1wGCVfx3IEGBPdt+2Us4GPgD8Pp0iOKmd5b/QzvjW07Ter501DHg1Ija1Wnbh/i+1D816o1rIX20ZDrxapH0YsCH11LcomsPSNB8l6w1bk073HdbOetvLYcXW3V057IWI2NFq2W3lsAHdsN5ezQVZ73FQwfBIsiO2FtFq2tVkSQrY2VU+CHixxLJfAK6NiLcUfPpHxKwOxriI7PTjh8qcfi2wjd23raiIeCQiJpF1od8NzG4ZVWqWMmIotV//APQvGPfWDix7NdAoaWCrZZfa/2a9XS3kr6Ik/RlZIVLslOEaYP8UY4u2ctivIuK9ZKdinwa+1zKq1CzthFds3eXmsLasBg7SrneXOod1kQuy3uMySSMkNQKfB+5oY9rbgL+RNC5dR3Ed8HBErEzjXwYKn+31PeBSSScos4+kv2xVULQrIv4X+CLwTUkfktRfUoOk90va7TqrdMrgp8C0NO1YSl9/sVe6AHa/iHgD2Ai03Er+MjBI2U0FHfWFtO4jgL/hzf26FPiApEZJbwWuaDVf631YuF0vAAuBf5bUT9kFxhdR+jo2s96u6vNXa5L2lXQG2fVm/5Eua9hFRPweWAJ8OeWodwIfLLG8AyWdmQqorcBmds1hIyTt1YlQW9Z9MnAG8OPUvhQ4K+W3Q8lyUKGSOQx4mKyg+2zK4aem7bq9E/FZ4oKs97iN7ALL59Lnn0pNGBHzgS8Ad5IdwR0CnFswyTRgZure/0hELCG7DuMmYAPZBZoXdibIiLgB+AzZRbtryY5eP0XWo1XMp8i6ul8CZpDusirh48BKZXdNXgr8dVrn02QX0T+XtqkjXfa/Jtve+cD1EXFvav8R8DtgJdl+b/0fyD8D/5jWd1WR5U4GmsiONO8iu0ZmXgfiMutNaiJ/JT+XtIksd/0DcAPZwVop5wEnkJ3S/BLwwxLT7QFcSZYTXgX+HPjbNO5+YBnwkqR1HYj1JbJtXk12wHdpyocANwJ/Iiu8ZrL7AeE0CvZj4YiI+BPZzQ/vB9aRPbro/IJlWycoopyzNlbNJK0EPhkR9+Udi5lZRzh/mWXcQ2ZmZmaWMxdkZmZmZjnzKUszMzOznLmHzMzMzCxnLsjMzMzMclbTT/4ePHhwNDU15R2GmfWgRx99dF1EDMk7ju7gHGZWX9rKXzVdkDU1NbFkyZK8wzCzHiSprddn1RTnMLP60lb+8ilLMzMzs5y5IDMzMzPLmQsyMzMzs5zV9DVkZvXgjTfeoLm5mS1btuQdSo/q168fI0aMoKGhIe9QzKwL6jGHdSZ/uSAzq3LNzc0MHDiQpqYmJOUdTo+ICNavX09zczOjRo3KOxwz64J6y2GdzV8+ZWlW5bZs2cKgQYPqIpG1kMSgQYPq6ojarLeqtxzW2fzlgsysBtRLIiuU5zZL6idpsaTfSVom6cupvVHSPEnPpp/7F8xzjaQVkp6RdHpuwZtVoXrLYZ3ZXhdkZtauqVOn8o1vfGPn99NPP51PfvKTO7+fffbZ7LvvvowbN47GxkZGjRrFuHHjOO200/IItztsBd4dEccA44D3SToRuBqYHxGjgfnpO5LGAucCRwDvA74lqU8ukZvZLmolf/kaMrMaM3169y5vypT2p5kwYQI//vGPueKKK9ixYwfr1q1j48aNO8evXr2aefPmccIJJ3DhhRdyxhlncM4553RvoD0oIgLYnL42pE8Ak4BTU/tM4EHgc6n99ojYCjwvaQVwPLCo56I2qw09ncNqJX9VrIdM0kGSHpC0PHX5X57ap0l6UdLS9PlAwTzu8jerQhMnTmThwoUALFu2jCOPPJKBAweyYcMGtm7dyvLlyzn22GNzjrJ7SeojaSnwCjAvIh4GDoyINQDp5wFp8uHACwWzN6e2YsudImmJpCVr166t3AaYGVA7+auSPWTbgCsj4jFJA4FHJc1L426MiOsLJ27V5T8MuE/S2yNiewVjNLMyDBs2jD333JNVq1axcOFCTjrpJF588UUWLVrEfvvtx9FHH81ee+2Vd5jdKuWecZLeAtwl6cg2Ji92wUiUWO50YDrA+PHji05jZt2nVvJXxQqydPTYciS5SdJyShwxJu7y74wVnej7PbSMc1RmrbQcZS5cuJDPfOYzvPjiiyxcuJD99tuPCRMm5B1exUTEa5IeJLs27GVJQyNijaShZL1nkPWIHVQw2whgdc9GamXrTN5sj/NqVauF/NUjF/VLagKOBR5OTZ+S9LikWwruUiq7y9/Met6ECRNYuHAhTzzxBEceeSQnnngiixYtYuHChUycODHv8LqVpCGpZwxJewOnAU8Dc4AL0mQXAD9Lw3OAcyX1lTQKGA0s7tmozayUWshfFS/IJA0A7gSuiIiNwLeBQ8juXFoDfL1l0iKz79ad7+svzPIxceJE5s6dS2NjI3369KGxsZHXXnuNRYsWcdJJJ+UdXncbCjwg6XHgEbJryOYCXwXeK+lZ4L3pOxGxDJgNPAX8ErjMl1uYVY9ayF8VvctSUgNZMXZrRPwUICJeLhj/PWBu+lpWl7+vvzDLx1FHHcW6des477zzdmnbvHkzgwcPzjGy7hcRj5P16rduXw+8p8Q81wLXVjg0M+uEWshfFSvIlD0V7WZgeUTcUNA+tOUuJeCvgCfT8BzgNkk3kF3U7y5/syLKeUxFJfTp02eXW8UBZsyYsdt0xdrMzFrkkcNqIX9VsodsIvBx4Il06zjA54HJksaRnY5cCVwCWZe/pJYu/224y9/MzMzqRCXvsnyI4teF3dPGPO7yNzMzs7rjVyeZmZmZ5cyvTjKrAm29SmTcOCh2Q/GQIZWLx8zMepZ7yMzMzMxy5oLMzMzMLGcuyMysXQMGDNg5fM899zB69GhWrVrFtGnTGD58OOPGjWPs2LHMmjVrl/m2bdvG4MGDueaaa3Zpnzt3LsceeyzHHHMMY8eO5bvf/W6PbIeZ1adayGG+hsysxvRrThec/W83LbAD7+CbP38+n/70p7n33nsZOXIkAFOnTuWqq67i2Wef5bjjjuOcc86hoaEBgHvvvZcxY8Ywe/ZsrrvuOiTxxhtvMGXKFBYvXsyIESPYunUrK1eu7KaNMbOq193vEu0lOcw9ZGZWlt/85jdcfPHF/OIXv+CQQw7Zbfzo0aPp378/GzZs2Nk2a9YsLr/8ckaOHMlvf/tbADZt2sS2bdsYNGgQAH379mXMmDE9sxFmVreqPYe5IDOzdm3dupVJkyZx9913c9hhhxWd5rHHHmP06NEccMABALz++uvMnz+fM844g8mTJ+88FdDY2MiZZ57JwQcfzOTJk7n11lvZsWNHj22LmdWfWshhLsjMrF0NDQ1MmDCBm2++ebdxN954I2PGjOGEE05g2rRpO9vnzp3Lu971Lvr378/ZZ5/NXXfdxfbt2cs3vv/97zN//nyOP/54rr/+ej7xiU/01KaYWR2qhRzmgszM2rXHHnswe/ZsHnnkEa677rpdxk2dOpVnnnmGO+64g/PPP58tW7YAWVf/fffdR1NTE8cddxzr16/ngQce2DnfUUcdxdSpU5k3bx533nlnj26PmdWXWshhvqi/HnXmgsoOXDRpvVP//v2ZO3cuJ598MgceeCAXXXTRLuPPOussZs6cycyZM5k8eTIPPfQQL7zwAn379gXgBz/4AbNmzeLEE09kyZIlnHrqqQAsXbqUgw8+uKc3x8zqTLXnMBdkZla2xsZGfvnLX3LKKacwePDg3cZ/8Ytf5LzzzqOhoYF3v/vdOxMZwKRJk/jsZz/LDTfcwNe+9jUuueQS9t57b/bZZx9mzJjRg1thZvWqmnOYIqLLC8nL+PHjY8mSJXmHka/uvn24FPeQVVTbr05azqhRh+/WXg+vTlq+fDmHH77rtkt6NCLG5xRSt3IOy0kl8qZzZEnF/o7rQUfzl68hMzMzM8uZCzIzMzOznLkgMzMzM8uZCzKzGlDL13p2Vj1us1lvVW9/z53ZXhdkZlXuj3/sx6ZN6+sqoUUE69evp1+/fnmHYmZd1K9fP9avr58c1tn85cdemFW5554bATTTv//aXdrXrcsnnp7Sr18/RowYkXcYZtZFI0aMoLm5mbVr17Y/cS/Rmfzlgsysym3b1sB///eo3dqn+C57M6sBDQ0NjBq1ew6zXfmUpZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmZmZ5cwFmZmZmVnOXJCZmbUi6SBJD0haLmmZpMtT+zRJL0pamj4fKJjnGkkrJD0j6fT8ojezWuSXi5uZ7W4bcGVEPCZpIPCopHlp3I0RcX3hxJLGAucCRwDDgPskvT0itvdo1GZWs9xDZmbWSkSsiYjH0vAmYDkwvI1ZJgG3R8TWiHgeWAEcX/lIzay3cEFmZtYGSU3/v737D5brrO87/v4gO4ZgSOxYNsKyYjeIAUGJzWg8FBNGBQYMocgkmMpTiFI8vc7UhtBAWpnMFGcYtaYNEMKvicDCgoKNggErDFNwBImaQhEChC3JdqxgjS2kkUSgA05bN5K//WPPrRZ579Ve6e49Z+99v2Z29uyzzzn7PXuuHn33Oec8D3AZ8M2m6IYkdyfZmOScpuxC4OG+1fYzfQInST/DU5aSNIUkZwN3AG+tqp8k+QjwLqCa5/cAbwIyYPWaYpsTwATAsmXLRhG2xtyGDbOznYmJ2dmO5oY9ZJI0QJIz6SVjn6qqzwFU1aGqOlZVjwEf5fhpyf3ARX2rLwUODNpuVW2oqpVVtXLx4sWj2wFJY8WETJJOkCTALcC9VfXevvIlfdVeC+xqlrcAa5KcleQSYDmwfa7ilTT+PGUpSY93BfBG4J4kO5uydwDXJLmU3unIfcB1AFW1O8lmYA+9OzSv9w5LSTNhQiZJJ6iqv2bwdWFfmmad9cD6kQUlaV7zlKUkSVLLTMgkSZJaZkImSZLUspFdQ5bkIuATwNOAx4ANVfX+JOcCnwEupndR7Our6sfNOjcC1wLHgLdU1ZdHFZ807mY6VpFjEklSd42yh2xyLrhnAy8Arm/me1sHbK2q5cDW5vWJc8FdCXw4yaIRxidJktQJI0vIppkLbjWwqam2CbiqWXYuOEmStCDNyTVkJ8wFd0FVHYRe0gac31Qbai64JBNJdiTZceTIkVGGLUmSNCdGnpCdOBfcdFUHlD1uLjinHZEkSfPNSBOyQXPBAYcmpx9png835UPPBSdJkjSfjPIuy4FzwdGb820tcHPzfGdf+aeTvBd4Os4FpzE107sf2DEABQAAF8xJREFUJUka5dRJU80FdzOwOcm1wEPA1eBccJIkaeEaWUI2zVxwAC+dYh3ngpMkSQuOI/VLkiS1zIRMkiSpZSZkkiRJLRvlRf2SJC0Me2fv9upnNf8z33fUCWgXEnvIJEmSWmZCJkmS1DITMkmSpJaZkEmSJLXMhEySJKll3mUpSdI8NFvz6k54s+ecsIdMkiSpZSZkkiRJLTMhkyRJapkJmSRJUstMyCRJklpmQiZJktQyEzJJkqSWmZBJkiS1zIRMkiSpZY7UL0nSLNi2re0INM5MyCRJ3bZ3luYAkjrMU5aSJEktMyGTpBMkuSjJ15Lcm2R3kt9tys9NcleSB5rnc/rWuTHJ3iT3J3lFe9FLGkcmZJL0eEeBt1XVs4EXANcnWQGsA7ZW1XJga/Oa5r01wHOAK4EPJ1nUSuSSxpIJmSSdoKoOVtV3muWfAvcCFwKrgU1NtU3AVc3yauD2qnq0qh4E9gKXz23UksaZCZkkTSPJxcBlwDeBC6rqIPSSNuD8ptqFwMN9q+1vyiRpKCZkkjSFJGcDdwBvraqfTFd1QFlNsc2JJDuS7Dhy5MhshClpHjAhk6QBkpxJLxn7VFV9rik+lGRJ8/4S4HBTvh+4qG/1pcCBQdutqg1VtbKqVi5evHg0wUsaOyZkknSCJAFuAe6tqvf2vbUFWNssrwXu7Ctfk+SsJJcAy4HtcxWvpPHnwLCS9HhXAG8E7kmysyl7B3AzsDnJtcBDwNUAVbU7yWZgD707NK+vqmNzH7akcWVCJkknqKq/ZvB1YQAvnWKd9cD6kQUlaV7zlKUkSVLLTMgkSZJaZkImSZLUMhMySZKklpmQSZIktcyETJIkqWUmZJIkSS0zIZMkSWqZCZkkSVLLHKlfw9m74dTWe8bE7MYhSdI8ZA+ZJElSy0aWkCXZmORwkl19ZTcl+UGSnc3jVX3v3Zhkb5L7k7xiVHFJkiR1zSh7yG4FrhxQ/r6qurR5fAkgyQpgDfCcZp0PJ1k0wtgkSZI6Y2QJWVVtA340ZPXVwO1V9WhVPQjsBS4fVWySJEld0sY1ZDckubs5pXlOU3Yh8HBfnf1NmSRJ0rw313dZfgR4F1DN83uANwEZULcGbSDJBDABsGzZstFE2ZZTvZNRkiSNtTntIauqQ1V1rKoeAz7K8dOS+4GL+qouBQ5MsY0NVbWyqlYuXrx4tAFLkiTNgTlNyJIs6Xv5WmDyDswtwJokZyW5BFgObJ/L2CRJktoyslOWSW4DVgHnJdkPvBNYleRSeqcj9wHXAVTV7iSbgT3AUeD6qjo2qtgkSZK6ZGQJWVVdM6D4lmnqrwfWjyoeSZKkrnKkfkmSpJYNlZAl2TpMmSR1je2XpHEw7SnLJE8Efp7edWDncHx4iqcCTx9xbJJ0ymy/JI2Tk11Ddh3wVnqN17c53qD9BPjQCOOSpNNl+yVpbEybkFXV+4H3J3lzVX1gjmKSpNNm+yVpnAx1l2VVfSDJC4GL+9epqk+MKC5JmhW2X5LGwVAJWZJPAr8C7AQmxwcrwAZNUqfZfkkaB8OOQ7YSWFFVA+eXlKQOs/2S1HnDjkO2C3jaKAORpBGx/ZLUecP2kJ0H7EmyHXh0srCqXjOSqCRp9th+Seq8YROym0YZhCSN0E1tByBJJzPsXZZ/NepAJGkUbL8kjYNh77L8Kb27kgB+DjgT+PuqeuqoApOk2WD7JWkcDNtD9pT+10muAi4fSUSSNItsvySNg2HvsvwZVfUF4CWzHIskjZztl6QuGvaU5W/0vXwCvXF9HNNHUuedavuVZCPwauBwVT23KbsJ+FfAkabaO6rqS817NwLX0ht89i1V9eXZ2gdJ89+wd1n+s77lo8A+YPWsRyNJs+9U269bgQ/y+BH931dVf9RfkGQFsAZ4Dr3JzP8iyTOr6hiSNIRhryH7l6MORJJG4VTbr6raluTiIauvBm6vqkeBB5PspXed2jdO5bMlLTxDXUOWZGmSzyc5nORQkjuSLB11cJJ0ukbQft2Q5O4kG5Oc05RdCDzcV2d/UyZJQxn2ov6PA1vodcVfCPx5UyZJXTeb7ddH6E1UfilwEHhPU54BdQdep5ZkIsmOJDuOHDkyqIqkBWjYhGxxVX28qo42j1uBxSOMS5Jmy6y1X1V1qKqOVdVjwEc5PnzGfuCivqpLgQNTbGNDVa2sqpWLF9uMSuoZNiH7YZI3JFnUPN4A/N0oA5OkWTJr7VeSJX0vX0tv4nLo9cCtSXJWkkuA5cD204pa0oIy7F2Wb6J3t9H76HXDfx3wQn9J4+CU2q8ktwGrgPOS7AfeCaxKcmmznX3AdQBVtTvJZmAPvTs5r/cOS0kzMWxC9i5gbVX9GCDJucAf0WvoJKnLTqn9qqprBhTfMk399cD604hT0gI27CnL5002ZgBV9SPgstGEJEmzyvZLUucNm5A9oe/27slfmMP2rklSm2y/JHXesI3Se4CvJ/ksvWsnXo9d85LGg+2XpM4bdqT+TyTZQW9C3gC/UVV7RhqZJM0C2y9J42DobvumAbMRkzR2bL8kdd2w15BJkiRpREzIJEmSWmZCJkmS1DITMkmSpJaZkEmSJLXMhEySJKllJmSSJEktMyGTJElqmfO5SSexYUPbEUiS5jt7yCRJklpmQiZJktQyEzJJkqSWmZBJkiS1bGQJWZKNSQ4n2dVXdm6Su5I80Dyf0/fejUn2Jrk/yStGFZckSVLXjLKH7FbgyhPK1gFbq2o5sLV5TZIVwBrgOc06H06yaISxSZIkdcbIErKq2gb86ITi1cCmZnkTcFVf+e1V9WhVPQjsBS4fVWySJEldMtfXkF1QVQcBmufzm/ILgYf76u1vyiRJkua9rlzUnwFlNbBiMpFkR5IdR44cGXFYkiRJozfXI/UfSrKkqg4mWQIcbsr3Axf11VsKHBi0garaAGwAWLly5cCkTdLjncqMAxMTsx+HJOnx5rqHbAuwtlleC9zZV74myVlJLgGWA9vnODZJkqRWjKyHLMltwCrgvCT7gXcCNwObk1wLPARcDVBVu5NsBvYAR4Hrq+rYqGKTJEnqkpElZFV1zRRvvXSK+uuB9aOKR5Ikqau6clG/JEnSgjXXF/VLkqQxcio3BA3iTULTs4dMkiSpZSZkkiRJLTMhkyRJapnXkEmSZs/eWbrgSFpg7CGTJElqmQmZJA2QZGOSw0l29ZWdm+SuJA80z+f0vXdjkr1J7k/yinailjSuPGUpSYPdCnwQ+ERf2Tpga1XdnGRd8/rfJVkBrAGeAzwd+Iskz3TGkfGwbVvbEUj2kEnSQFW1DfjRCcWrgU3N8ibgqr7y26vq0ap6ENgLXD4ngUqaF0zIJGl4F1TVQYDm+fym/ELg4b56+5sySRqKpywl6fRlQFkNrJhMABMAy5YtG2VMGnPPOmP271i976jD5XeVPWSSNLxDSZYANM+Hm/L9wEV99ZYCBwZtoKo2VNXKqlq5ePHikQYraXyYkEnS8LYAa5vltcCdfeVrkpyV5BJgObC9hfgkjSlPWUrSAEluA1YB5yXZD7wTuBnYnORa4CHgaoCq2p1kM7AHOApc7x2WkmbChEySBqiqa6Z466VT1F8PrB9dRJLmM09ZSpIktcweslFwLjdJkjQD9pBJkiS1zIRMkiSpZSZkkiRJLfMaMo3WqVxP9wxHkpYkLSz2kEmSJLXMhEySJKllJmSSJEktMyGTJElqmQmZJElSy0zIJEmSWmZCJkmS1DITMkmSpJaZkEmSJLXMhEySJKllJmSSJEktMyGTJElqmQmZJElSy0zIJEmSWmZCJkmS1DITMkmSpJaZkEmSJLXMhEySJKllJmSSJEktMyGTJElqmQmZJElSy85o40OT7AN+ChwDjlbVyiTnAp8BLgb2Aa+vqh+3EZ/mrw0b2o5gvJzK9zUxMftxSNJ812YP2T+tqkuramXzeh2wtaqWA1ub15IkSfNel05ZrgY2NcubgKtajEWSJGnOtJWQFfCVJN9OMnmC44KqOgjQPJ8/aMUkE0l2JNlx5MiROQpXkiRpdFq5hgy4oqoOJDkfuCvJfcOuWFUbgA0AK1eurFEFKEmSNFda6SGrqgPN82Hg88DlwKEkSwCa58NtxCZJkjTX5jwhS/LkJE+ZXAZeDuwCtgBrm2prgTvnOjZJkqQ2tHHK8gLg80kmP//TVfVfk3wL2JzkWuAh4OoWYpMkSZpzc56QVdX3gV8dUP53wEvnOh5JminHUpQ027o07IUkjRPHUpQ0a0zIJGl2OJaipFPW1rAXkjTOJsdSLOBPm+F4fmYsxWZYn8dpxl6cAFi2bNlcxSu1bramrpuv07OZkEnSzDmWoqRZ5SlLSZohx1KUNNtMyCRpBhxLUdIoeMpSkmbGsRQlzToTMkmaAcdSlDQKnrKUJElqmQmZJElSy0zIJEmSWmZCJkmS1DITMkmSpJaZkEmSJLXMhEySJKlljkMmSRpL27a1HYE0e+whkyRJapkJmSRJUstMyCRJklpmQiZJktQyEzJJkqSWmZBJkiS1zIRMkiSpZY5Dpu7Zu2Hm6zxjYvbjkCRpjpiQSdJCdSo/fiSNhAmZxlb/KN33fbW9OCRJOl1eQyZJktQyEzJJkqSWecpSkqQF4llnzO51g/cd9Yaq2WJCdjJe9CrNyIYZ/pOZsD2XJE9ZSpIktc2ETJIkqWUmZJIkSS0zIZMkSWqZCZkkSVLLTMgkSZJa5rAXkiRpbMx0aJ2pdG3IHXvIJEmSWmYPmeaFUx192lGmNTYcpFqa1xZWQmaDJkmt27at7Qik7llYCZk6zUZakrRQmZBJkqRTMtuTlcPCvZSkcwlZkiuB9wOLgI9V1c0thyRJQ7H9ksZH1+7W7NRdlkkWAR8CXgmsAK5JsqLdqCTp5Gy/JJ2OrvWQXQ7srarvAyS5HVgN7Gk1Kkk6OdsvaRYs1NOgneohAy4EHu57vb8pk6Sus/2SdMq61kOWAWX1MxWSCWAy1X0kyf0jj+q484AfzuHnjdp82p9T3JfrZj2QWTKfjg1Msz/XzfwQ/PLpBjMiJ22/YM7bsK7/HXU9Puh+jF2PDzoR47QNzWnFN8M2bMr2q2sJ2X7gor7XS4ED/RWqagPQyoBiSXZU1co2PnsU5tP+zKd9AfdnTJ20/YK5bcO6/r13PT7ofoxdjw+6H2NX4uvaKctvAcuTXJLk54A1wJaWY5KkYdh+STplneohq6qjSW4AvkzvtvGNVbW75bAk6aRsvySdjk4lZABV9SXgS23HMYX5NvfSfNqf+bQv4P6MpQ62X13/3rseH3Q/xq7HB92PsRPxpepx15xKkiRpDnXtGjJJkqQFx4SskWRjksNJdvWV3ZTkB0l2No9X9b13Y5K9Se5P8op2op5akouSfC3JvUl2J/ndpvzcJHcleaB5Pqdvnc7u0zT7M3bHKMkTk2xP8r1mX/6wKR/XYzPV/ozdsZkvklzZfLd7k6xrMY5B7Wpn/s673k6OU1uRZFGS7yb5YhdjTLIvyT1NW7SjizFSVT56p21fDDwf2NVXdhPw9gF1VwDfA84CLgH+FljU9j6cEOMS4PnN8lOAv2ni/k/AuqZ8HfDucdinafZn7I4RvfGqzm6WzwS+CbxgjI/NVPszdsdmPjzo3VDwt8A/An6u+a5XtBTLoHa1M3/nXW8nx6mtAH4P+DTwxa4d5+Zz9wHnnVDWqRjtIWtU1TbgR0NWXw3cXlWPVtWDwF5606Z0RlUdrKrvNMs/Be6lN2r4amBTU20TcFWz3Ol9mmZ/ptLZ/ameR5qXZzaPYnyPzVT7M5VO78888P+ncKqq/wtMTuE056ZoVzvzd971dnJc2ookS4FfBz7WV9ypGKfQqRhNyE7uhiR3N13vk92ZYzVFSpKLgcvo/bq6oKoOQq8xAs5vqo3NPp2wPzCGx6jp3t8JHAbuqqqxPjZT7A+M4bGZB7r+/Xby77yr7eSYtBV/DPxb4LG+sq7FWMBXknw7vdkyOhejCdn0PgL8CnApcBB4T1M+1BQpXZDkbOAO4K1V9ZPpqg4o69w+DdifsTxGVXWsqi6lN5r75UmeO031Tu8LTLk/Y3ls5oFx/X5bi7vL7WTX24okrwYOV9W3h11lQNlcHOcrqur5wCuB65O8eJq6rcRoQjaNqjrU/GN4DPgox7ssh5oipW1JzqTXyHyqqj7XFB9KsqR5fwm9X10wBvs0aH/G/RhV1f8E/hK4kjE+NpP692fcj80Y6/r326m/83FpJzvcVlwBvCbJPnqnx1+S5L90LEaq6kDzfBj4PL32qFMxmpBNY/JANV4LTN4ptAVYk+SsJJcAy4Htcx3fdJIEuAW4t6re2/fWFmBts7wWuLOvvLP7NNX+jOMxSrI4yS82y08CXgbcx/gem4H7M47HZp7o+hROnfk773o7OQ5tRVXdWFVLq+pien9rX62qN3QpxiRPTvKUyWXg5fTao87ECHiX5eQDuI3eaZV/oJcdXwt8ErgHuLs5QEv66v8BvTsv7gde2Xb8A/bnRfS6WO8GdjaPVwG/BGwFHmiezx2HfZpmf8buGAHPA77bxLwL+PdN+bgem6n2Z+yOzXx5NP82/qb5jv+gxTgGtaud+Tvvejs5bm0FsIrjd1l2JkZ6dxx/r3nsnvw30aUYq8qR+iVJktrmKUtJkqSWmZBJkiS1zIRMkiSpZSZkkiRJLTMhkyRJapkJmTolyV8mWTmC7a5K8sIh6r0myboZbvuRk9eamVFsU9Lo2YaNbpvz3RltB6D5L8kZVXW05c9ZBTwCfH26bVTVFro1iKakltmGaS7YQzbPJLk4yX1JPpZkV5JPJXlZkv+e5IEklzf1ntxM+PytJN9Nsrpv/f+W5DvN44VN+ZIk25LsbLb7a035I32f/boktzbLtyZ5b5KvAe+e5vOelOT29Cag/gzwpCn2a1+SdyfZ3jyeMcXnnJvkC832/keS56U3afDvAP+mif/XmhGw72ji+VaSK5rt/XaSD/Zt+0+SfD3J95O8bojv//eb7d2d5A+bsncn+dd9dW5K8rap6ksLmW2YbdiCNZej+PoY/QO4GDgK/GN6Cfe3gY30JktdDXyhqfcfgDc0y79Ib1TvJwM/DzyxKV8O7GiW38bx0Y0XAU9plh/p++zXAbc2y7cCXwQWneTzfg/Y2JQ/r4l95YD92tf3+b/F8dGgT/ycDwDvbJZfAuxslm8C3t63vU8DL2qWl9GbOgXgt4EP9m37z5rvcQWwd4rv/JHm+eXAhua7fkIT14uBy4C/6qu/p/nMgfVP/F59+FhID9sw27CF+vCU5fz0YFXdA5BkN7C1qirJPfQaO+j9Q3pNkrc3r59I7x/YAeCDSS4FjgHPbN7/FrAxvYl4v1BVO4eI48+q6thJPu/FwJ8AVNXdSe6eZnu39T2/b4rPeRHwm832vprkl5L8woBtvQxYkWTy9VPTzHV2gi9Ub3LsPUkumCY26O3jy+lNdQJwNrC8qm5Jcn6SpwOLgR9X1UNJ3jKoPrDtJJ8jzXe2YdiGLTQmZPPTo33Lj/W9fozjxzzAb1bV/f0rJrkJOAT8Kr1fPP8HoKq2JXkx8OvAJ5P856r6BL154CY98YQ4/r5/01N8HidsYzo1xfKJnzPdepOeAPyTqvrfA+Lp1/9dDto2J7z/H6vqTwe891l6v76fBtw+RH1pIbMNm3q9SbZh84zXkC1cXwbenOZfb5LLmvJfAA42v6jeSK9rnyS/DByuqo8CtwDPb+ofSvLsJE8AXnsKn7cN+BdN2XPpdflP5Z/3PX9jijr921sF/LCqfgL8FOj/9fgV4IbJF82v6dP1ZeBNSc5utnlhkvOb924H1tBr0D47RH1J07MNsw2bV+whW7jeBfwxcHfTwOwDXg18GLgjydXA1zj+y20V8PtJ/oHenT6/1ZSvo3fdwMPALnpd1jP5vI8AH2+6+XcC26eJ+awk36T3Q+KaKerc1Le9/wWsbcr/HPhsehfivhl4C/Chpt4Z9BrB35nms0+qqr6S5NnAN5o2+xHgDfT+E9jdnE74QVUdPFn904lDWiBsw2zD5pVUDdvTKrUnyT56F8r+sO1YJGmmbMN0Mp6ylCRJapk9ZJIkSS2zh0ySJKllJmSSJEktMyGTJElqmQmZJElSy0zIJEmSWmZCJkmS1LL/B1D152vgH16HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x792 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(10, 11))\n",
    "for ax, avg_wt, avg_kras, sd, pname in zip(axes.flatten(), real_wt_means, real_kras_means, std_devs, protein_data.protein_name):\n",
    "    wt_sample = np.random.normal(avg_wt, sd, 1000)\n",
    "    kras_sample = np.random.normal(avg_kras, sd, 1000)\n",
    "    ax.hist(wt_sample, color='blue', alpha=0.4, label='WT')\n",
    "    ax.hist(kras_sample, color='orange', alpha=0.4, label='KRAS')\n",
    "    ax.set_title(pname + ' distribution')\n",
    "    ax.legend(loc='upper right')\n",
    "    ax.set_ylabel('count')\n",
    "    ax.set_xlabel('measured protein level')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate observed ratios\n",
    "\n",
    "First, we must calculate the *observed* ratio, the ratio of the average protein levels that we actually got from our data.\n",
    "I used the ratio of average G12D levels to average WT levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>WT0</th>\n",
       "      <th>WT1</th>\n",
       "      <th>WT2</th>\n",
       "      <th>WT3</th>\n",
       "      <th>KRAS0</th>\n",
       "      <th>KRAS1</th>\n",
       "      <th>KRAS2</th>\n",
       "      <th>KRAS3</th>\n",
       "      <th>KRAS4</th>\n",
       "      <th>protein_name</th>\n",
       "      <th>ratio</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140.608634</td>\n",
       "      <td>84.706090</td>\n",
       "      <td>86.795706</td>\n",
       "      <td>73.175784</td>\n",
       "      <td>121.635191</td>\n",
       "      <td>42.461533</td>\n",
       "      <td>143.620294</td>\n",
       "      <td>80.969827</td>\n",
       "      <td>107.975977</td>\n",
       "      <td>proteinA</td>\n",
       "      <td>1.031260</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>82.544074</td>\n",
       "      <td>202.347556</td>\n",
       "      <td>44.209850</td>\n",
       "      <td>77.430796</td>\n",
       "      <td>73.116195</td>\n",
       "      <td>179.363861</td>\n",
       "      <td>23.007611</td>\n",
       "      <td>87.930025</td>\n",
       "      <td>38.549911</td>\n",
       "      <td>proteinB</td>\n",
       "      <td>0.791017</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>300.844275</td>\n",
       "      <td>311.656304</td>\n",
       "      <td>277.987616</td>\n",
       "      <td>322.894474</td>\n",
       "      <td>218.031814</td>\n",
       "      <td>210.049887</td>\n",
       "      <td>218.017119</td>\n",
       "      <td>186.325443</td>\n",
       "      <td>197.542195</td>\n",
       "      <td>proteinC</td>\n",
       "      <td>0.679071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>234.496140</td>\n",
       "      <td>281.247834</td>\n",
       "      <td>337.124883</td>\n",
       "      <td>251.583747</td>\n",
       "      <td>172.227253</td>\n",
       "      <td>151.897911</td>\n",
       "      <td>140.835605</td>\n",
       "      <td>153.012771</td>\n",
       "      <td>199.113478</td>\n",
       "      <td>proteinD</td>\n",
       "      <td>0.591849</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          WT0         WT1         WT2         WT3       KRAS0       KRAS1  \\\n",
       "0  140.608634   84.706090   86.795706   73.175784  121.635191   42.461533   \n",
       "1   82.544074  202.347556   44.209850   77.430796   73.116195  179.363861   \n",
       "2  300.844275  311.656304  277.987616  322.894474  218.031814  210.049887   \n",
       "3  234.496140  281.247834  337.124883  251.583747  172.227253  151.897911   \n",
       "\n",
       "        KRAS2       KRAS3       KRAS4 protein_name     ratio  \n",
       "0  143.620294   80.969827  107.975977     proteinA  1.031260  \n",
       "1   23.007611   87.930025   38.549911     proteinB  0.791017  \n",
       "2  218.017119  186.325443  197.542195     proteinC  0.679071  \n",
       "3  140.835605  153.012771  199.113478     proteinD  0.591849  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_wt_data(df, row_i):\n",
    "    return df.loc[row_i, ['WT0', 'WT1', 'WT2', 'WT3']]\n",
    "\n",
    "\n",
    "def get_kras_data(df, row_i):\n",
    "    return df.loc[row_i, ['KRAS0', 'KRAS1', 'KRAS2', 'KRAS3', 'KRAS4']]\n",
    "\n",
    "\n",
    "protein_data['ratio'] = np.empty(len(protein_data))\n",
    "\n",
    "for i in range(len(protein_data)):\n",
    "    wt_data = get_wt_data(protein_data, i)\n",
    "    kras_data = get_kras_data(protein_data, i)\n",
    "    \n",
    "    wt_mean = np.mean(wt_data)\n",
    "    kras_mean = np.mean(kras_data)\n",
    "    \n",
    "    protein_data.loc[i, ['ratio']] = kras_mean / wt_mean\n",
    "\n",
    "protein_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boostrap confidence intervals of ratios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_array(a):\n",
    "    \"\"\"Take a sample from the array `a` of the same length.\"\"\"\n",
    "    return np.random.choice(a, size=len(a), replace=True)\n",
    "\n",
    "\n",
    "def mean_sample(a):\n",
    "    \"\"\"Calculate the mean from a random sample of the array.\"\"\"\n",
    "    s = sample_array(a)\n",
    "    return np.mean(s)\n",
    "\n",
    "\n",
    "def confidence_intervals(a):\n",
    "    \"\"\"Calculate the upper and lower 95% CIs of an array.\"\"\"\n",
    "    # Just a wrapper around a `statsmodels` object + method.\n",
    "    u, l = sms.DescrStatsW(a).tconfint_mean()\n",
    "    return u, l\n",
    "    \n",
    "    \n",
    "    \n",
    "def bootstrap_ratio_CI(g1, g2, nboots=1000):\n",
    "    \"\"\"\n",
    "    Bootsrap the mean and confidence interval of the ratio of `g1` to `g2`.\n",
    "    \"\"\"\n",
    "    \n",
    "    ratios = []\n",
    "    for i in range(nboots):\n",
    "        ratios.append(mean_sample(g2) / mean_sample(g1))\n",
    "    \n",
    "    m = np.mean(ratios)\n",
    "    l, u = confidence_intervals(ratios)\n",
    "    \n",
    "    return m, l, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>WT0</th>\n",
       "      <th>WT1</th>\n",
       "      <th>WT2</th>\n",
       "      <th>WT3</th>\n",
       "      <th>KRAS0</th>\n",
       "      <th>KRAS1</th>\n",
       "      <th>KRAS2</th>\n",
       "      <th>KRAS3</th>\n",
       "      <th>KRAS4</th>\n",
       "      <th>protein_name</th>\n",
       "      <th>ratio</th>\n",
       "      <th>boot_ratio</th>\n",
       "      <th>lower_ci</th>\n",
       "      <th>upper_ci</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140.608634</td>\n",
       "      <td>84.706090</td>\n",
       "      <td>86.795706</td>\n",
       "      <td>73.175784</td>\n",
       "      <td>121.635191</td>\n",
       "      <td>42.461533</td>\n",
       "      <td>143.620294</td>\n",
       "      <td>80.969827</td>\n",
       "      <td>107.975977</td>\n",
       "      <td>proteinA</td>\n",
       "      <td>1.031260</td>\n",
       "      <td>1.059094</td>\n",
       "      <td>1.007954</td>\n",
       "      <td>1.110235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>82.544074</td>\n",
       "      <td>202.347556</td>\n",
       "      <td>44.209850</td>\n",
       "      <td>77.430796</td>\n",
       "      <td>73.116195</td>\n",
       "      <td>179.363861</td>\n",
       "      <td>23.007611</td>\n",
       "      <td>87.930025</td>\n",
       "      <td>38.549911</td>\n",
       "      <td>proteinB</td>\n",
       "      <td>0.791017</td>\n",
       "      <td>0.817061</td>\n",
       "      <td>0.746572</td>\n",
       "      <td>0.887550</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>300.844275</td>\n",
       "      <td>311.656304</td>\n",
       "      <td>277.987616</td>\n",
       "      <td>322.894474</td>\n",
       "      <td>218.031814</td>\n",
       "      <td>210.049887</td>\n",
       "      <td>218.017119</td>\n",
       "      <td>186.325443</td>\n",
       "      <td>197.542195</td>\n",
       "      <td>proteinC</td>\n",
       "      <td>0.679071</td>\n",
       "      <td>0.680348</td>\n",
       "      <td>0.674519</td>\n",
       "      <td>0.686178</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>234.496140</td>\n",
       "      <td>281.247834</td>\n",
       "      <td>337.124883</td>\n",
       "      <td>251.583747</td>\n",
       "      <td>172.227253</td>\n",
       "      <td>151.897911</td>\n",
       "      <td>140.835605</td>\n",
       "      <td>153.012771</td>\n",
       "      <td>199.113478</td>\n",
       "      <td>proteinD</td>\n",
       "      <td>0.591849</td>\n",
       "      <td>0.595329</td>\n",
       "      <td>0.583709</td>\n",
       "      <td>0.606950</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          WT0         WT1         WT2         WT3       KRAS0       KRAS1  \\\n",
       "0  140.608634   84.706090   86.795706   73.175784  121.635191   42.461533   \n",
       "1   82.544074  202.347556   44.209850   77.430796   73.116195  179.363861   \n",
       "2  300.844275  311.656304  277.987616  322.894474  218.031814  210.049887   \n",
       "3  234.496140  281.247834  337.124883  251.583747  172.227253  151.897911   \n",
       "\n",
       "        KRAS2       KRAS3       KRAS4 protein_name     ratio  boot_ratio  \\\n",
       "0  143.620294   80.969827  107.975977     proteinA  1.031260    1.059094   \n",
       "1   23.007611   87.930025   38.549911     proteinB  0.791017    0.817061   \n",
       "2  218.017119  186.325443  197.542195     proteinC  0.679071    0.680348   \n",
       "3  140.835605  153.012771  199.113478     proteinD  0.591849    0.595329   \n",
       "\n",
       "   lower_ci  upper_ci  \n",
       "0  1.007954  1.110235  \n",
       "1  0.746572  0.887550  \n",
       "2  0.674519  0.686178  \n",
       "3  0.583709  0.606950  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add new columns to be filled by bootstrapping.\n",
    "for c in ['boot_ratio', 'lower_ci', 'upper_ci']:\n",
    "    protein_data[c] = np.empty(len(protein_data))\n",
    "\n",
    "\n",
    "# Perform boostrapping on each protein.\n",
    "for i in range(len(protein_data)):\n",
    "    \n",
    "    wt_data = get_wt_data(protein_data, i)\n",
    "    kras_data = get_kras_data(protein_data, i)\n",
    "    \n",
    "    mean, lower_ci, upper_ci = bootstrap_ratio_CI(wt_data, kras_data, nboots=100)\n",
    "    \n",
    "    protein_data.loc[i, ['boot_ratio']] = mean\n",
    "    protein_data.loc[i, ['lower_ci']] = lower_ci\n",
    "    protein_data.loc[i, ['upper_ci']] = upper_ci\n",
    "\n",
    "\n",
    "protein_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "\n",
    "The scatter plot function in `plotly` requires the error values be represented as distances from the data point.\n",
    "Therefore, we must subtractthe lower interval from the observed ratio and subtract the observed ratio from the upper interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>protein_name</th>\n",
       "      <th>ratio</th>\n",
       "      <th>boot_ratio</th>\n",
       "      <th>log_ratio</th>\n",
       "      <th>lower_ci</th>\n",
       "      <th>upper_ci</th>\n",
       "      <th>log_lower_ci_diff</th>\n",
       "      <th>log_upper_ci_diff</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>proteinA</td>\n",
       "      <td>1.031260</td>\n",
       "      <td>1.059094</td>\n",
       "      <td>0.044408</td>\n",
       "      <td>1.007954</td>\n",
       "      <td>1.110235</td>\n",
       "      <td>0.032979</td>\n",
       "      <td>0.106457</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>proteinB</td>\n",
       "      <td>0.791017</td>\n",
       "      <td>0.817061</td>\n",
       "      <td>-0.338219</td>\n",
       "      <td>0.746572</td>\n",
       "      <td>0.887550</td>\n",
       "      <td>0.083428</td>\n",
       "      <td>0.166119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>proteinC</td>\n",
       "      <td>0.679071</td>\n",
       "      <td>0.680348</td>\n",
       "      <td>-0.558365</td>\n",
       "      <td>0.674519</td>\n",
       "      <td>0.686178</td>\n",
       "      <td>0.009704</td>\n",
       "      <td>0.015019</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>proteinD</td>\n",
       "      <td>0.591849</td>\n",
       "      <td>0.595329</td>\n",
       "      <td>-0.756698</td>\n",
       "      <td>0.583709</td>\n",
       "      <td>0.606950</td>\n",
       "      <td>0.019981</td>\n",
       "      <td>0.036346</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  protein_name     ratio  boot_ratio  log_ratio  lower_ci  upper_ci  \\\n",
       "0     proteinA  1.031260    1.059094   0.044408  1.007954  1.110235   \n",
       "1     proteinB  0.791017    0.817061  -0.338219  0.746572  0.887550   \n",
       "2     proteinC  0.679071    0.680348  -0.558365  0.674519  0.686178   \n",
       "3     proteinD  0.591849    0.595329  -0.756698  0.583709  0.606950   \n",
       "\n",
       "   log_lower_ci_diff  log_upper_ci_diff  \n",
       "0           0.032979           0.106457  \n",
       "1           0.083428           0.166119  \n",
       "2           0.009704           0.015019  \n",
       "3           0.019981           0.036346  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "protein_data['log_ratio'] = np.log2(protein_data['ratio'])\n",
    "\n",
    "protein_data['log_lower_ci_diff'] = protein_data.log_ratio - np.log2(protein_data.lower_ci)\n",
    "protein_data['log_upper_ci_diff'] = np.log2(protein_data.upper_ci) - protein_data.log_ratio\n",
    "\n",
    "protein_data[['protein_name', 'ratio', 'boot_ratio', 'log_ratio',\n",
    "             'lower_ci', 'upper_ci', 'log_lower_ci_diff', 'log_upper_ci_diff']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "        <script type=\"text/javascript\">\n",
       "        window.PlotlyConfig = {MathJaxConfig: 'local'};\n",
       "        if (window.MathJax) {MathJax.Hub.Config({SVG: {font: \"STIX-Web\"}});}\n",
       "        if (typeof require !== 'undefined') {\n",
       "        require.undef(\"plotly\");\n",
       "        define('plotly', function(require, exports, module) {\n",
       "            /**\n",
       "* plotly.js v1.51.2\n",
       "* Copyright 2012-2019, Plotly, Inc.\n",
       "* All rights reserved.\n",
       "* Licensed under the MIT license\n",
       "*/\n",
No results found