figures/main-text/F4-predictors.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import pickle\n",
    "%matplotlib inline\n",
    "import scipy\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def getData(path):\n",
    "    df = pd.read_csv(path, sep=',', header=0, index_col=0)\n",
    "    return df\n",
    "\n",
    "names = pickle.load(open('../../data/hosts/common_names.p', 'rb')) # Convert between common and scientific names\n",
    "names['Chimp_pts'] = 'Pan_troglodytes_schweinfurthii'\n",
    "\n",
    "bmr = getData('../../data/hosts/body_mass_rawdata.csv') # body mass values\n",
    "hgd = getData('../../data/hosts/hostdistance_rawdata.csv') # host genetic distance\n",
    "hgd_names = { i: i.split('|')[2] for i in list(hgd.columns.values)+list(hgd.index.values) } # fix taxa -> host names\n",
    "hgd.rename(columns=hgd_names, index=hgd_names, inplace=True)\n",
    "hgd.rename(index={'P._t._schweinfurthii': 'Pan_troglodytes_schweinfurthii'}, \n",
    "           columns = {'P._t._schweinfurthii': 'Pan_troglodytes_schweinfurthii'},\n",
    "          inplace=True)\n",
    "\n",
    "bf = getData('../../beast/main/discreteTraits/results/discreteTraits_2016-10-25_bf.csv') # Bayes factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = [] # Make tidy data\n",
    "for from_host in bmr.index.values[1:]:\n",
    "    for to_host in bmr.index.values[1:]:\n",
    "        if from_host != to_host:\n",
    "            transmission = 1 if bf.at[names[from_host], names[to_host]] >= 36.0 else 0 # Code transmission as binary variable\n",
    "            bmr_val = bmr.at[from_host, 'mass'] / bmr.at[to_host, 'mass']\n",
    "            hgd_val = hgd.at[from_host, to_host]\n",
    "            data.append({ 'bmr': bmr_val, 'hgd': hgd_val, 'transmission': transmission, 'intercept': 1 } )\n",
    "            \n",
    "data = pd.DataFrame(data, columns=['transmission', 'intercept', 'bmr', 'hgd'])\n",
    "data['log_bmr'] = data['bmr'].map(np.log10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def make_quantiles(df, column):\n",
    "    df['%s_quantiles'%column] = pd.qcut(df[column], 5, labels=range(5)) # Cut the data into quantiles (makes new categorical series column)\n",
    "    x_means = [] # [x value means for each quantile]\n",
    "    y_means = [] # [(y mean, (lower, upper 95% binary y CI)) for each quantile]\n",
    "    y_ci = []\n",
    "    \n",
    "    for q in range(5):\n",
    "        obs = df.loc[df['%s_quantiles'%column] == q] # observations in this quantile\n",
    "        n_obs = len(obs) # number of observations in quantile\n",
    "        n_transmissions = obs['transmission'].sum() # number of 1 values in quantile\n",
    "        xmean = obs[column].mean() # mean predictor value in this quantile\n",
    "        ymean = float(n_transmissions) / float(n_obs)\n",
    "        ystddev = np.sqrt(float(n_obs)*ymean*(1-ymean))\n",
    "        y95ci = 1.96*(ystddev/np.sqrt(n_obs))\n",
    "        x_means.append(xmean)\n",
    "        y_means.append(ymean)\n",
    "        y_ci.append(y95ci)\n",
    "    return {'x': x_means, 'y': y_means, 'y_ci': y_ci}\n",
    "    \n",
    "log_bmr_quantiles = make_quantiles(data, 'log_bmr')\n",
    "hgd_quantiles = make_quantiles(data, 'hgd')\n",
    "# bmr_quantiles = make_quantiles(data, 'bmr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.202398\n",
      "         Iterations 7\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.184935\n",
      "         Iterations 8\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:           transmission   No. Observations:                  506\n",
      "Model:                          Logit   Df Residuals:                      504\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Mon, 09 Jan 2017   Pseudo R-squ.:               0.0008416\n",
      "Time:                        18:36:51   Log-Likelihood:                -102.41\n",
      "converged:                       True   LL-Null:                       -102.50\n",
      "                                        LLR p-value:                    0.6779\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "log_bmr       -0.1251      0.301     -0.415      0.678        -0.716     0.465\n",
      "intercept     -2.9188      0.202    -14.453      0.000        -3.315    -2.523\n",
      "==============================================================================\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:           transmission   No. Observations:                  506\n",
      "Model:                          Logit   Df Residuals:                      504\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Mon, 09 Jan 2017   Pseudo R-squ.:                 0.08705\n",
      "Time:                        18:36:51   Log-Likelihood:                -93.577\n",
      "converged:                       True   LL-Null:                       -102.50\n",
      "                                        LLR p-value:                 2.396e-05\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "hgd           -4.9234      1.383     -3.561      0.000        -7.633    -2.213\n",
      "intercept     -1.4423      0.380     -3.793      0.000        -2.188    -0.697\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "log_bmr_logit = sm.Logit(data['transmission'], data[['log_bmr', 'intercept']]).fit() # Logistic regression summary\n",
    "hgd_logit = sm.Logit(data['transmission'], data[['hgd', 'intercept']]).fit()\n",
    "print log_bmr_logit.summary()\n",
    "print hgd_logit.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhMAAADNCAYAAAACCPJSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmYVNWZ+PHvvbf2rm56ge5mbQFZDI4oqIhKBLXdM4qA\nQQk4mokmM8EYNSqPE5dfZEg0MXFUJiEzMWLyiCZqjBgXdh2igIRdQEVoaaD3tbprv+f3R3VVd/VW\nvVSvvp/n6aeq7nLOudX33PvWueeeqymlFEIIIYQQXaT3dQGEEEIIMbBJMCGEEEKIbpFgQgghhBDd\nIsGEEEIIIbpFggkhhBBCdIsEE0IIIYToFktfFyDqyJEjXHvttQBkZGTwwQcfYLH0m+IJIQagZcuW\n8frrr7eYbrPZGDp0KJdccgk//OEPSUtL64PSCTF49JuWiTfeeAMATdOoqqpi8+bNfVsgIcSgoWla\n3F8wGOTUqVO89NJL3HnnnchwO0J0T78JJtauXYumabFK/Ze//KWPSySEGCwyMzPZsmULmzdvZtOm\nTbz55pv88z//MwC7d++WHy9CdFO/CCY+/vhjTp48CcCYMWNQSrFlyxaqqqr6uGRCiMFA0zSys7PJ\nyckhNzeX008/nQcffDA2/5NPPunD0gkx8PWLYOKvf/0rAE6nk4ceegiAUCjEW2+91ZfFEkIMYpqm\nxd6np6f3YUmEGPj6PJgIBoO8++67aJrG7NmzmTVrFpmZmYBc6hBCJJ9SioqKCp544gkg8iPmiiuu\n6ONSCTGw9fntElu2bKG6uhpN07jiiivQdZ3LL7+cV155hf3793PkyBHGjx/f18UUQgxgZWVlTJ48\nucV0p9PJk08+ybBhw/qgVEIMHn3eMhG9xGG327nkkksAuPrqq2PzpXVCCNFdze/miP55vV7+8Ic/\nUFpa2tdFFGJA69NgwuPxsHnzZjRNY+zYsXz88cd88MEHBAIBbDYbSinefPPNviyiEGIQaH43x9q1\na3n00Uex2Wx8+OGHLF26tK+LKMSA1qeXOd5++20CgQCapnHw4EG+853vtFimuLiYv//971x44YV9\nUEIhxGAQvZujqfHjx3PgwAFeeeUV9uzZQ0FBAXl5eX1UQiEGtj5tmWja6tBWMyTIpQ4hRM+wWq2x\n9x6Ppw9LIsTA1mctE8XFxezYsQNN07juuut48sknWyxzzTXX8MUXX7Bu3Trq6+txuVx9UFIhxECn\nlKK4uDj2ORgMsnfv3tgPFbfbzYQJE/qqeEIMeH0WTLzxxhsopWLBRGtuuOEGnnrqKXw+H++++y5z\n587t5VIKIQaD8vLyWAfv5jRNY+nSpdhstl4ulRCDR59d5ogOn52RkcHFF1/c6jLXX389hmEAtPqw\nHiGESKS1y6cWi4WUlBTOPPNM/vM//5MlS5b0dTGFGNA0JU+4EUIIIUQ39Pk4E0IIIYQY2CSYEEII\nIUS3SDAhhBBCiG6RYEIIIYQQ3SLBhBBCCCG6RYIJIYQQQnSLBBNCCCGE6BYJJoQQQgjRLRJMCCGE\nEKJbJJgQQgghRLdIMCGEEEKIbpFgQgghhBDdIsGEEEIIIbrF0tcF6Ck+n4/9+/czbNiw2GPMhRAd\nk5ubi8UyaA8PQohWhEIhioqKulT/B+3RYv/+/SxatKiviyHEgLRhwwZGjRrV18Xosoceeojc3Ny+\nLoYQA0pNTQ2rV6/uUv0ftMHEsGHDAPjjH/8oBxUhOuGFF14Y8HUmNzeXpUuX9nUxhBhQCgsLWb16\ndZfWHbTBRPTSRm5u7oD+hSVEb0tLS5NLHEKITpEOmEIIIYToFgkmhBBCCNEtEkwIIYQQolskmBBC\nCCFEt/S7YGLPnj0sXry4xfSNGzcyf/58Fi5cyJ/+9Kc+KJkQoidJ3Rdi4OpXXbb/53/+hzfeeIOU\nlJS46aFQiJ/+9Ke89tpr2O12br75Zi677DIyMzP7qKRCiGSSui/EwNavgom8vDyee+457r///rjp\nR44cIS8vD7fbDcD06dPZsWMHV155ZV8Us99QSnH48GHKy8vJyspi0qRJaJrWYplDhw5x6NAhNE1j\n8uTJrS7XXppAwnwATNNk48aNfP755yilyM3NZdKkSRw/fpx9+/ZhmiZ5eXlkZGRQUlKCzWZjxIgR\n+P1+7HY7Q4cOZeLEiXz66accPHiQ6upqhgwZwhlnnBHL0zRNNm3aRFFREVarlREjRsStl6iMnf3+\nRO8YSHW/vf0mOq+srCy2X/v9fhwOB1lZWUyYMIHNmzfHRhmcM2cOuq63um50neb7ZXS5aB0BGDJk\nCJMnT2by5MmxZcPhMK+++irFxcVomsaQIUMYMWJEXJ6Jtqez38uhQ4fYtm1brO7OmDEjrky9rbPH\nv74W/V8cOnQIpVSL/2l/1q+Cifz8fE6cONFiusfjITU1NfY5JSWF2tra3ixav3T48GEOHDgAQFFR\nEQCTJ09uscxHH30UO+hUVla2ulx7aQIJ8wHYtGkTu3btwu/3Y5om1dXVHDt2DI/Hg2masfQNw8Bm\nsxEMBjlx4kTsQFdcXMyJEyf48ssvKSsrIxgMUlpaSlVVVSzPTZs2sW/fPoLBYGx+bm4uJ06ciC3X\nXhk7+/2J3jGQ6n57+010Xk1NDdXV1VitVoLBIEOGDCEtLY29e/dy8uRJAEpLSwG47LLLWl03uk7T\n9KPLffjhh7GgI5pWVVVV7IQJ8Oqrr/L5558TDodj9bG8vDwuz0Tb09nvZdOmTVRUVGCaJrW1tdTV\n1cWVqbd19vjX16L/22h5m/9P+7N+FUy0xe124/F4Yp/r6upilay5bdu2sX37dmpqanqreH0memBo\n63N0WvSAA+D3+1tdrrNptqaoqAjTNFFKAZFfRj6fL/YZIpG3aZqx4MLn8+F0OmNlLCoqigUjEGnt\naFrm6MGu6frR6Q6HI2EZO7utom91pe5DZCS/ntLefhN9H92ffT4fhmHEPpeUlMSt2zRgb75u9LW1\n/KJ1JFq3mteTpnk1XaZ5nom2pzPKy8vxer2x/JRS+Hy+Pq1XnT3+9bWBVt6m+l0HTCDu5AMwfvx4\nCgoKqKmpIRAIsGPHDs4+++xW150xYwZLly7l1ltv7Y2i9qmsrKx2P0en2e322Ge73d7qcu2l2ZF8\nIDLaqK7rsSY5wzBwOBxxTXSapqHreqyZNRoARMuYm5uL3W6Pzdd1Pa7M0WGem6/ffPjn9raxvW0V\nfSsZdX/p0qU9Oupte/tN9H10f26+f2dnZ8et23S/bb5u9LW1/KJ1RNO0WJ1qXrejeUXrX7TOJKor\nXa0HWVlZOJ3OWH6apsUu1fSVzh7/+tpAK29T/bJlIrozrl27Fq/Xy4IFC1i2bBm33347SikWLFjQ\nolJ+FUX7MzTv39B8meh1w6bXDLuSZnv5AMyZMwelVI/1mYjmAXSoz0Qyvj/RuwZC3W9vv4m+70yf\nibbWbdpnorX8W+sz0XTZefPmtdlnoqPb09nvRSnVos9EX9arzh7/+lq0bE37TPTn8jalqeY/BQaJ\nwsJCLrvssgH/9EMhetszzzwz4B+SNRi2QYje1p3zZr+8zCGEEEKIgUOCCSGEEEJ0iwQTQgghhOgW\nCSaEEEII0S0STAghhBCiWySYEEIIIUS3SDAhhBBCiG6RYEIIIYQQ3SLBhBBCCCG6RYIJIYQQQnSL\nBBNCCCGE6BYJJoQQQgjRLUl9aqjf7+eDDz6gqKgIXdfJzc3l/PPPx+12JzMbIYQQQvQjSWuZ2LVr\nF/n5+axevZo9e/awa9cufv/733PVVVfx4YcfJisbIYQQQvQzSWuZePjhh1m1ahWTJ0+Om37o0CHu\nv/9+/vrXv7a7vlKKRx99lMOHD2Oz2Vi+fDmjR4+OzV+3bh2//vWv0XWdG2+8kZtvvjlZRRdC9DGp\n/0IMbEkLJkzTbBFIAEyePBmlVML1169fTyAQYM2aNezZs4cVK1awcuXK2PwVK1bwxhtv4HA4uPba\na7nuuutITU1NVvGFEH1I6r8QA1vSLnOMGDGCVatWUVlZGZtWW1vLb3/7W0aOHJlw/Z07dzJr1iwA\npk6dyv79++PmW61Wqqur8fv9AGialqyiCyH6mNR/IQa2DrVMnDhxgurq6rgWhilTpsQt87Of/Yz/\n9//+H3PmzIlb7pJLLmH58uUJ8/B4PHG/NCwWC6ZpouuReOf2229n3rx5uFwu8vPzpVOnEIOI1H8h\nBraEwcSTTz7JH/7wB7KysmLTNE1jw4YNcctlZmbyq1/9inA4TGVlJaZpkpWVhWEYHSqI2+2mrq4u\n9rnpgeTUqVP84Q9/YOPGjbhcLu677z7effddrrzyyhbpbNu2je3bt1NTU9OhfIUQfS8Z9T9a9wEK\nCwt7r/BCiMTBxNtvv817771HTk5OhxI0DIOhQ4fGTdu0aRNz5sxpd71p06axadMmrrrqKnbv3s3E\niRNj8/x+P4ZhYLPZ0DSNzMzMNoOFGTNmMGPGDAoLC1m9enWHyiyE6FvJqP/Rug/wzDPP9FrZhRAd\nCCaGDx/e4UCiLS+99FLCYCI/P5+tW7eycOFCINLhau3atXi9XhYsWMANN9zAwoULcTgcjBkzhrlz\n53arTEKI/kPqvxADm6YS3Grx7LPPUl9fz2WXXYbD4YhNb95nor8pLCzksssuY8OGDYwaNaqviyPE\ngPHMM8+wdOnSvi5GtwyGbRCit3XnvJmwZeK1114D4J133olNa63PBEAoFMIwDDRNY9u2bRw8eJDz\nzjuv3wceQgghhOi6hMHExo0bO5TQ+++/z913382QIUP413/9V1atWsWFF17I6tWreeCBB1rtLCmE\nEEKIgS9hMFFfX88TTzzB+++/TygU4qKLLuKhhx5qcWvWr371K15++WUKCgr4wQ9+wHvvvcfIkSMp\nLi7me9/7ngQTQgghxCCVcNCqFStWEAgEeO6551i5ciWapvGTn/ykxXKmaTJhwgQuv/xyRowYERuo\nKicnB9M0k19yIYQQQvQLCVsm9uzZE/dcjccff5xrr722xXJOp5ONGzdy6aWXsm7dutj0t956SwaY\nEUIIIQaxhC0T4XA4rmXBNM1WB6J65JFH+NWvfhW37N/+9jeefvppHnnkkSQVVwghhBD9TcKWiZkz\nZ3L33XfHntL30ksvxQaGaWry5Mktngyan5/PNddck6SiCiGEEKI/ShhMPPjgg6xcuZKnnnoK0zS5\n+OKL+bd/+7cOJW61WrtdQCGEEEL0bwmDCYvFwl133cVdd93VG+URQgghxADTZjBx880389JLL3HO\nOefEPe5XKYWmafzjH//olQIKIYQQon9rM5h4+umnAVi7dm2nEty6dSsXXXRR7FUIIYQQg1ubd3Nk\nZ2cDkUeLl5aWMnLkSNavX8+zzz4b11LR3M9//vO4VyGEEEIMbglvDV22bBkbNmxg7969rF69mhEj\nRvDjH/84YcIJnh8mhBBCiEEiYTBx/Phx7r33XjZt2sTcuXNZunQpVVVVSS+IUopHHnmEhQsXsmTJ\nEo4fPx43f+/evSxatIhFixZxzz33EAwGk14GIUTfkPovxMCWMJiIVtr/+7//44ILLiAcDlNfX5/0\ngqxfv55AIMCaNWu49957WbFiRdz8hx9+mJ/+9Kf88Y9/ZObMmRQWFia9DEKIviH1X4iBLeGtodOm\nTeOaa67BMAymTZvGrbfeyoUXXtjm8l29vLFz505mzZoFwNSpU9m/f39s3tGjR0lPT+f555/ns88+\nY/bs2YwdO7ZL+Qgh+p9k13+/P0B1jYchaTKUvxC9IWEw8eMf/5hdu3YxadIkdF3n29/+Nl//+tfb\nXP6BBx6Ie+0oj8dDampqY8EsFkzTRNd1Kisr2b17N4888gijR4/mzjvv5Mwzz2x1JE4hxMCT7Pqv\nlKLKY2IYXtwpzt7YBCG+0hJe5ggEAlgsFlJTU3nhhRd49913KSoqanP5mTNncuTIEXJycjpVELfb\nTV1dXexz9EACkJ6ezpgxYxg7diwWi4VZs2bF/XIRQgxsPVH/rTY75dUBvD5/j5VbCBHRI3dzLFu2\njIceeqhTBZk2bRpbtmwBYPfu3UycODE2b/To0dTX18c6Ze3cuZPTTz+91XS2bdvGM888wwsvvNCp\n/IUQfScZ9T9a95955hlOnjoFgNXmoKS8XjpsCtHDEl7mOH78OL/61a94+umnmTt3Lt///veZN29e\nm8sfOnQIm82GruscPnyYSZMmdagg+fn5bN26lYULFwKwYsUK1q5di9frZcGCBSxfvpx77rkHgHPO\nOYdLLrmk1XRmzJjBjBkzKCwsZPXq1R3KWwjRt5JR/6N1H+DJJxvHubHaXZwsqWHMiMx2x8gRQnRd\nwmCi6d0cDzzwQMK7OV555RXmz5+PruusWbOmw48f1zSNxx57LG5a005WM2bM4E9/+lOH0hJCDCw9\nXf+tdjcniysZmZvZ5TSEEG1LeJnjnHPO4ZprrsHn8zFt2jT+5V/+pc27OXw+H5s3b+bqq6/mqquu\nYsuWLXi93qQXWgghOkvpTsoqavq6GEIMSglbJh5++OEWd3O0dYnhrbfeYtasWdjtdgAuvfRS1q5d\ny4IFC5JbaiGE6CRd16nza9hq60hLTenr4ggxqLQZTHz44YfMnDmTDRs2xD5HrVu3jiuuuKLFOmPG\njGHmzJmxz//6r//Kl19+mczyCiFEl1mtNiprfdisfhwOe18XR4hBo81g4q233mLmzJm8+OKLLeZp\nmtZqMHHeeefFfc7NzSU3NzcJxRRCiOSw2hwUldcxZrg1dvupEKJ72gwmHn/8cYBWgwkhhBjIbPYU\nThRXMSo3Q+7wECIJEvaZ+Oyzz3jxxReprq6Om/7000/3WKGEEKKn6RYXp0qqGJGT0ddFEWLASxhM\n3H333Vx88cUJx4soLS1l2LBhSSuYEEL0JE3TMDUHRSWV5GZLQCFEdyQMJhwOB8uWLUuY0B133MHr\nr78OwJ///Gfmz5/f/dIJIUQP0nWdYMhGeUUNWZlpfV0cIQashL2Pzj//fLZs2UI4HG53uaZPC/3j\nH//Y/ZIJIQaNUw3DW/dHhsVCXUCnusbT10URYsBK2DKRlZXFnXfeGeukpJRC0zQOHjwYt1zTTkxd\nfQy5EGJwevXVV/n+97/f18Vok8Vipcrjx2r14XI6+ro4Qgw4CYOJF198kVdeeYXRo0d3OFHpHS2E\naKo/BxJRVpud0op6RmQbWK3Wvi6OEANKwmAiMzOTs846K2FC5eXlPP/88y3eR912221dLKIQQvQO\nq93FqdJaRg+XW0aF6IyEwcQFF1zAXXfdxRVXXIHNZotNbz5o1UUXXcSnn37a4r0QQgwkhtUlDwUT\nopMSBhP79+8H4OWXX45Na20EzBUrViS5aEII0fs0TcPUnRSXVpEzLL2viyPEgNChPhNRSilCoVCb\n1xN37dpFZmYmeXl5rFq1io8//pizzjqL7373u1gs7WellOLRRx/l8OHD2Gw2li9f3mo/jYcffpj0\n9HTuueeeREUXQgwQ/a3+67pOIGSlsqqWjPTUHs1LiMEg4a2hH3/8MStXriQQCHDjjTdy7rnn8re/\n/a3Fcs899xx33XUXixcv5sEHH+T9999n9uzZ7N27l+XLlycsyPr16wkEAqxZs4Z777231ZaONWvW\nyOUTIQah/lj/DYuFmnpFvdfXa3kKMVAlDCaefPJJzj77bNavX8/QoUN56623+N3vftdiub/97W+8\n/fbbvPzyy7z11lv893//N7fccgvPPfdc3BNH27Jz505mzZoFwNSpU2OXV6J27drFvn37WLhwYUe3\nTQgxQPTX+h+5w8OL3x/o1XyFGGgSBhPhcJgLL7yQv//971x++eWMGjUK0zRbLGexWHC73QwfPpy8\nvDxSUyNNg1arFafTmbAgHo8ntk40vWg+paWlPPvsszz88MMyhoUQg1B/rv9Wu4uisjqCwWCv5y3E\nQJGwz4Rpmuzdu5fNmzfz3e9+l08//bTVStX0Ub5duUfb7XZTV1cXl280zXfeeYeqqiq+853vUFpa\nit/vZ9y4cdxwww0t0tm2bRvbt2+npqam02UQQvSNZNT/aN0HOJnkETetdhcnSzyMyHbLGBRCtCJh\nMPHd736Xe++9l/nz5zNq1CguvfRSHnrooRbL1dTUsG7dOpRS1NbW8t5778Xm1dbWJizItGnT2LRp\nE1dddRW7d+9m4sSJsXmLFy9m8eLFALz++uscPXq01UACYMaMGcyYMYPCwkJWr16dMF8hRN9LRv2P\n1n2AJ5/8edLLGB2DYmTOEAzDSHr6QgxkCYOJK664Iu420HXr1rVakUaMGBE7eQ8fPjzuLpDhw4cn\nLEh+fj5bt26NXRNdsWIFa9euxev1smDBgsRbIoQYsAZK/bfYUjhRXM2o3PS41lghvuoSBhOff/45\n//u//0tVVVXc9cpf//rXccs1DR66QtM0HnvssbhpY8eObbHc3Llzu5WPEKL/GUj132JL4WRxFSNz\nZZRMIaISBhP3338/06ZN47zzzpOKI4QQgGZxUVRSxfCcjL4uihD9QsJgIhgM8h//8R+9URYhhBgQ\nNE0jrDkoLq0kZ5gEFEIkvOiXl5dHSUlJb5RFCCEGDF3XCZp2Ssur+7ooQvS5Dt0aet111zFlyhTs\ndntsevM+E00VFxdTXl4e18diypQp3SyqEEL0L7ph4AsqKqpqyZRht8VXWMJgIj8/n/z8/A4n+Mtf\n/pLnn3+eoUOHxqZpmsaGDRu6VkIhhOjHDIsFjy+IJs/xEF9hCYOJ5r2nlVIUFBS0ufybb77Jxo0b\n44IJIYQYzCwWK7XeALruYUiau6+LI0SvSxhMrFmzhieeeAKv1xubNnToUD744INWl8/IyJBAQgjx\nlWOx2qjy+LEYXlJSEj9CQIjBJGEwsWrVKp5//nn++7//m7vvvptNmzbF9Z2IOnDgAABf+9rXePzx\nx/nGN74R99hx6TMhhBjsrDY7pVVedEPH6Wh5nBRisEoYTKSnpzN16lTOOOMMysvL+d73vseiRYu4\n/fbb45ZbunRp3OeNGzfG3kufCSFEbyrxj6agyMfoHDt6L4+PY7M7KSn3MjTdlBYK8ZWRMJiwWCxU\nV1eTl5fH3r17ueiiiygrK2uxXDR4KCoqIjc3N27eZ599lqTiCiFEYjXhDDbvqiTNZeFrY1MYP9KJ\nxei9oMJqd1JW7SdsmqSlpvRavkL0lYTjTNx0003ceeedzJ49m5dffpkbb7yR8ePHt1iuqqqKqqoq\n7rjjDqqrq6mqqqK6upqysjL+/d//vUcKL4QQrYmGDTX1IT46UM2rm0vY/VktPn+418pgtdmp8piU\nV8gTjMXgl7Bl4owzzuB3v/sdLpeLl19+mX379nHxxRe3WO7ee+9l69atALEn9wEYhtGpW0uFEKK7\nxjgOMnnMOXxeWE/IVPgCJns+97D/izrGjXByxmkuMlJ7/lHiFquN+kCIYEkludkyUqYYvBIGE/fd\ndx9vv/02ADk5OeTk5LS63P/+7/8CsGzZMlasWJHEIgohROdY9SAzpqRx9gQ3h76s51BBHb6ASdhU\nfFZYz2eF9eRm2pmc52J0th1d77lLIIbFQlgZnCiqYHi2PG1UDE4Jg4lJkybx5ptvMn36dFwuV2x6\nenp63HJHjhxh/PjxfOtb34rd2dGU3M0hhOhtdpvO1NPdnDk2haOnvHxytJ5KTxCAogo/RRV+UhwG\nk8a4OH2UE6fd6JFyaJoGlhQKi6rIzkzBIXd6iEGmzWAiEAhgs9nYsGED77zzTtw8TdM4ePBg3LSf\n/exnrFq1qsVdHdHlE93NoZTi0Ucf5fDhw9hsNpYvX87o0aNj89euXcvq1auxWCxMnDiRRx99tCPb\nJ4QYAHq6/huGxumjXIwf6aS4IsAnx+opLPGjUNT5wvzj01p2f+YhL9fBxNFOcjJtPfKUZIstheIK\nP2mugIyWKQaVNoOJb37zm7z++uvs27evQwk99dRTQPwtoZ2xfv16AoEAa9asYc+ePaxYsYKVK1cC\n4Pf7+a//+i/Wrl2LzWbj3nvvZdOmTcyZM6dLeQkh+pfeqv+appGbZSc3y46nPsShL+v5vNCLP2hi\nKsXRU16OnvKS6rIwcbSTcSOcuBzJba2w2ux4fCF8JZXkDkvvkaBFiN7WZjDR9CFdHXHBBRcwffp0\nZs+ezezZsxk7dmyn1t+5cyezZs0CYOrUqezfvz82z2azsWbNGmw2GwChUKjVgbOEEANTX9R/t8vC\nuZPTOGdCKseKfHz6ZT0lVQEAautD7Dxcyz8+9TBqmJ3TRzkZOcyOkaS+FYbFgqkMCk9VkjPUHds2\nIQaqNoMJv9/PJ5980mZQ0bwPxPvvv89HH33Ehx9+yB//+Ec0TeOSSy5h9uzZnH/++Qkri8fjITW1\nsdnPYrFgmia6rqNpGpmZmQC8+OKLeL1eLrzwwg5vpBCif+vL+m8YGuNHOhk/0kllbZDPjns5csJL\nIGSilOJ4iY/jJT4cVp2xI5yMG+Ega4i12y0KmqZh2FI4Weol1eUjKyMtSVskRO9rM5g4fvw4S5cu\nbTWYaK0PRGZmJtdccw3XXHMNACdOnODvf/87P//5zykoKGDXrl3tFsTtdlNXVxf7HD2QRCmleOKJ\nJygoKODZZ59tM51t27axfft2amrk3m4hBopk1P9o3Qc4eepUl8qRkWrl/K9ZmTYplS+LfXxe6OVU\nuR8AX9DkYEEdBwvqSEuxMG6Ek7HDHaSlJOzH3i6b3YE3aPLlyQqGZrhwOR3dSk+IvtBmLTj99NP5\ny1/+0ukECwsL2bBhA1u3buWTTz5hypQp3HTTTQnXmzZtGps2beKqq65i9+7dTJw4MW7+j3/8YxwO\nR+w6altmzJjBjBkzKCwsZPXq1Z0uvxCi9yWj/kfrPsCTT/68W+WxGBrjRkT6THjqQ3x+wssXJ33U\n1ocAqKkLsfuzWnZ/VktWmpXThjs4LdeB29W1wELXdXRbCqWVfuy1PoZlpWIYPXNniRA9oXshdRO/\n/OUv2bhxI3V1dcyaNYtbbrmFCy64AIejY1F2fn4+W7duZeHChQCsWLGCtWvX4vV6mTJlCq+99hrT\np09n8eIDhaclAAAgAElEQVTFaJrGkiVLuPzyy5NVfCFEH+rP9d/tsnD2hFSmnu6mtCrIFye9HDvl\nwx80ASivCVJeE2Tn4VqGDrGSl+tgTE7XWiysNjsmcPxUDW6XQWa6W8alEANCm3v7ueee26mEfvOb\n33DppZdyxx13cPbZZ3e6IJqm8dhjj8VNa9qJ85NPPul0mkKIgWEg1H9N08jOsJGdYeO8M9I4VRbg\n2CkvXxb7CYYjgUVZdZCy6khgke62MibHzuhse6f7WNgcLvxhxZenqnE7DYr8n/P3Y+sprStmWEoO\nF429nCm503pqU4XotDaDif/4j//oVELvvPMOmzZt4he/+AXHjh3joosuYvbs2Vx88cW43e5uF1QI\nIfoLQ9cYlW1nVLadcFhxoszPsVM+CksaA4sqT5AqT5C9Rzw47QajhkWWH55lw2pJ3NqgaRo2u4uD\nZbvZUrAmNr3Ec4rX970IIAGF6DeSdpnjtNNO47bbbuO2226jpqaG999/n3Xr1rFixQrGjRvH888/\nn6yshBCi3zAMjTE5kUsb4bDiVLmfL4v9HC/x4QtEAguvPxwbxlvXNHIybYwYamPEUDsZqZZ2Wy32\nFm9udfrWoxskmBD9RtKCiaZOnjxJRUUFgUAAq9UqHYmEEF8JhqExKtvBqGwHpkqjrCrIl8U+TpT6\nqfJEOm+aKhJwnCr3s/NwLU6bwfChNoZn2cjNsuN2xh8vq3wlreZV6imisqqW9CFuGfhK9LmkBROr\nV69m+/bt7Nixg/T0dGbNmsX8+fO54IILZIApIcRXjt6kj8W5kyN3gJws83Oi1E9ReYCQGbnt3hsI\n88VJL1+c9AKQ6rSQm2UjJzPyl+7IptJX1CJ9Qzco9nio9oRw2nXcKTZSXM5e3UYhopIWTHzwwQd8\n/etf50c/+hF5eXnJSlYIIbqk+VgVfS0txUJaioXJeSmEw4riygAny/ycKgtQURuMLVfrDVFbGOKz\nwnoAwtoVlIV2YHfWYHdVY7XXoWmKkBnkjcNPc+HouYx3nE15TYjSigrsNp0Up5VUt0taLESvSVow\n8dvf/jZZSSVVRVUtrpQaDEPHajGwWi0YhtGvDjJfJQeK/sHWo9IrXfQsu91GukvhD/gJmxAyFaap\nCIcVkQYBDU03sFq7P5JlVxiGxoihdkYMjbTa+vxhTpUHKKqI/NXUhRqXVRmkhM+nqrSYCtOHw2Lj\ntOwc0twhvvRsZsPnr3F82GEuHHM9NkcKCqjxmlRUV2G1atitOi6HFafTIcGF6DE90meiPwkpO76w\nHRVShOvDKNOHaYbRUOiGhq4R+dMj7w1dQ9fBYhhYLJE/Xddjw/qKrjtQ9I9YL3SQXumi52iahtvt\noq37yMLhMKFQCH8gRChsxgKNsIJQyMQ0QWk6hmHBYun5w6TDbjB2hJOxIyKXKep9YUoqI4FFSWUQ\nrTYDty09tnzAC2VesIavp6TuS4qOVfGPT95m+pizGZc9nIw0KykOJ5qmETChviaEWVGF1aJhs2o4\nbBZSUpzyo0okzaAPJqI0TWv3oGA2/IUa3pgBE9MMY4YDKBUZo1/XQIsGHrqGoUc+Gw2fNRSGYWAY\nOoauNbyXVpCoD754r9Xp7x95l8nZUzH0tjvqSouGSKZo3WyvP1coFCIQCBIIBgiFI8FGtIUjFFaA\njt4QbCT7h4bLYXDacCenDY8EF4GgSWlVkNKqAKWVQUqrgwRDJlbDzsjU06nwFVNVXczGfYfZ5agm\nw5GN3WaQkWolI9US+0t3WwgqHb/XpKy6BouhsFl0rBaNFJdd+reJLvvKBBOdFW2NIEEAAhA2Gz+o\noMI0I4GIMlsGIjQEH5qmRQITjRbvUSoSnGhg6DqGocdaSAzD6BctJIlO7kqZlNYVc7zqKIVVRyms\nOkalt6zVtMrrS/jZxvvJcA1lqCuHrJRhZKXkkOXKZmhKDkfKD0qLhuh1FkskUHC1MT8cDhMMBvH5\ng5hNLqVEAw9TgaYZGBZLt+9os1l1Rg6zM3JY5GSvlKK6LkxZVZDy6iDlNXZOVKRSVFNApa8Ib7CW\n7JQx+IMmRRX+uLTczkhQMcRtNLzqpKVYqPUGMM06bFYdm0XDatFJcTmwWq3dKrv4ahjUwUTq8Ol8\nfjJMpb8udnKOnrDjTuI6aEQub7R2ctf1xnV0DWjyXtM0NB10mrxvOOl3RDQgQTX8NZsRCUxMTDMU\nCUxME01TGLqOpqmGyy/QUKy4skOkR7nWEMhYLAYWozEgafpHQ/k7oq3LFaWeIgzdEgkeqo/hD/li\nyzgsTmyGnUDY3yI9h8VJVko25XUllNeVQGn8fE1rvWVH7rNvnVKRgFYpFfszTZOwqVCmwmwyHzRM\npYi0q0XWrav3JcpC0Ni60dYjA5RShMNh/IEgwWjrRkPAYSqNsKkIh0zQIi0cnTlpa5pGujsSFJw+\nKtJ6YZqZlFRn88GRzRSUFVLmLyJDm4Jdz4hb1+MN4fGGKGxWzxw2nTSXhdQUgzSXBbdLx2mtJdWh\ncDp0LIaGxYgEGU6HDZvN1i9+2Ij+YVAHE+ljZnGgIISjpPefIKrHTuqNJ3M9Fpg0C1hoCEKaLdu8\n5SIuDZoGOg1pNAmMYsvrWqxVRJmKSJSiGtZRaESWgYbWE0A34vM2GuYrwgSVl/8r3IQ/6EbTFGiR\nNNBMNh/+ALRIVJRmT2dc2tfIdo9kuHs0mc6hfFF5gI0FrwINrTQNvj7mWiZknYWpFN6gh2p/GZW+\nciq9pVT5yiio/qzV77jEc5LVO54lO2U4w9y5DEsZztCUHKxG4+Puo0+91TQtdmKNfo7ON01F2DQb\nvp9ITGeasTAP01QdOmiqVj5Ep0XTi5y/NVRDmmZDeZRq7BjY9Em9SjUsZ6pIWs0e4quUavjnNsSj\nqnHbNF2PBQnRPj+RPx1Na7vqB0ItnxQsOi96aTVRn4vGFo4AwZBJ2IRgyESp6P898j/TG4KXtvZF\nXdfIzUhj/vRv8HnFTj4s/CvB8E5GuM/l9LQrqKvXqaoNUeUJUV0XarhU08gXMPEFApRUtUzbZtFJ\ndRm4nQYpTgO71UuKXZHq0klzGbjsOhZLJOBw2K0yvtBX0KAOJvqSqSJH93Dzo38/o5QibAYJq1Dk\nzwwRNoOEmryPTA9ixtpRTm8zvZyU03BYUvDrFgqAgticamAUvsCtVHmLCZo+bBY7mc5c3i/N5AOt\nsknAlIGmZaBpE9A0KK/4B/6wBxoCoGjAomtQhELTPMBnoH2KBtgtTlwWN05rKi6bmxRrKg6LC13T\nGwKnSDp6w+WmaB+YxgAvGow1BmKapjVp/WnZetXYOqQ1WaZpYNfYGgZNP+tNlm2YHw3qtObpRYLO\nWDmID0DlV+LA1JEWjlAoRDAYIhgKRVo0wopQOBJ4hMORgBhNj92hMiHrXHLcp7H52Mscr/uYyuAX\nzD7tm5w5Li+WZp03THVdmCpPiNr6ENWeMLX1Iep84RZlCIRMymtMymuCLeZB5CmrLodBisPAYQW7\nVeGyRy6puJ06Jb7DHK76gOpAKZnOYcwYPZuzRp7bby7biu4b1MFE0d7VzP63fIYOG4qpFEqBqSK/\nNJVq/DVomo3vm05XDeuEG36xhs3GaZG0GtI0I7/bI+nEp2Uq1WJa5Ndw9BdwY3oK1ZBHYz7R9at9\nFVTUn8If9mPVnaQ7cuJ6dzfVPEAImUFMM5QgQGiboVmw6DYM3YqhW6gP1mCqlgccm+4gxTak3bTc\ntvQW5Y4+y6AtLn0MHm9Bi+nZKXm4rGkEwz78YS+BkI9A2EtNvY8qFQaqGv5AQ8dmOJr8ObEZDgx9\n8FQBjcZApGmAEb3k1TxAokkLVtNASPQvmqZhtVoTXgYxTZNgMIg/ECIYDJNhd/ON0/+Fj09tYl/p\n//HWp79mdOoZ1AbLqfKVku7IZmruHKaMnRqXTihkUusNU1MXCS483jC19WE83jB13nDseBi3TlhR\n4wlS7QlgKhOlTEwil9Lqg9VU+sqAM9D1CZywBPjk0GHGpAfIcefgsltwO224nTZcdgOnTcNp00lx\n6BhGJNC3Wi1YLY0d2iUA6X8Gz5G0FUFvBUNSdLKGDOwORF9U7mHzsZfIaTZ91JB/wmG4qQ/W4Q16\n8Da8Bs1gpDld6eiATUV6fqqGV5SOzXBhN1KwGynYjBTsDZ9thgub4cKqO7EbLiy6Aw091syuFJys\nOcKeoi0NaULDqYivDZtFtsvdLIhqHlQ1BEq0Pt1ssk40oBpGNkPqFcWeAnwhL3bDxVDXGNLsQxuW\ns6JUalwawXAAX8iLP/oXjgQa/nB93HdoaJZYYNE0yBiIByvVcK0j/nJI/24ZE8mj6zp2e8s7Mkbm\nzOPsyrP58+7f8WVt49NXK31FbD72EsW1R8lyDCcYDhImjNJMwipIyAwSsgYJGQEMVxC3GcQeDuAL\nKLw+DX9Ax+838AcM/AELoaCVcNBOOGwlckxoyTQtmAELwYCLT+sr+ZTKuPla9JIOOpqmYzHAalVY\nLQqLATYrWC1gt2jYbQZ2q47DauCwWXDaLDisFlx2G06bgcNqw26xY7PYseoWbBZ7XL+x7tzy35N3\nlw3UO9f6TTChlOLRRx/l8OHD2Gw2li9fzujRo2PzN27cyMqVK7FYLMybN48FCxYkPf+mr+3Na2/Z\nNhKPvGrtr6M1XPtu+lkBu09tbHX5gup9jcuiYTdcpDndOIwUHJYUnBZ346uRgtMa+Ww3Ik3+7Ra5\n4Zp+5Jp9/DX8sdnjGTnMy/7yrVT7SxliH8aZQy/itCHjopvZuD1N1jMbrv8rUzXe2tJkO9ujMR4Y\nH0047juKW67JZYTmHUtNFabKW0alr5jy+lNUeoup9JXgCVRDLOCCIDpuaxZDHDkMsQ9jiD2bNPsw\nXNYhgA5KUVB9iIMl26j1V+C2ZTEx63xGpk4CGlvAmgdOLT83Lhv9jhpbqVpvBYvuPo19LRrnN+03\n0Xz5poEazfJvXp7iBP+LntDX9X+wy8sYT4o9FW+ovsW8g+UfdSotXTOw2Kw47VZSh9iw6NYmfzY0\nrKiwAzPsJBS08XnpQUIhG2bIRjhkIxyOvJohG25bVqwVI9qi0fg5jC9o4g0mbj1tnULXw2h6OPZq\nMSJBiWHQ2KHU0LBYdGyGgdWiY7MY2CwGdosFm9WC3WLBbrFGXq1Wyr2F7CrZHDsMRTugV9dVMWHY\n12KXLpsOCaCav4u9xB/5Pis7wPpP/xr7PJDuXOs3wcT69esJBAKsWbOGPXv2sGLFClauXAlE7vf+\n6U9/ymuvvYbdbufmm2/msssuIzMzM2G6Nj2I3Qg0xskaoKId6hTx17BbOTE1uSsCtNg1986MHdGR\nOyaad7prqvpAafPFGzZF49sz7sVtT8NlS0kYICTTaaNmMZtZ3Uqj6d0GTae1t3xb71XcyVU1nGwb\nTrDKbDi5Kty2DEanZ6LUZKJ3MgRCfiq8xZR5i6nwllDuLaLCW0KR7wRFTW5ssOhWMhzZWDQrp+qO\ngQXsFghSyYGazxmWPo9xGf8UO6mbptkQkEUuZwGxgIqGFp1IcKHF7xea1vALLRJkNd9/unIHTmf8\nz4GkJ5lQT9V/0ai8vu3jyNVnzMdq2LDqNqyGFYthw6pbm7y3oSmNcEhhhiOXfKN3poTNhoG+FGia\n0aKj6OsH91PpO9oi33R7LtdNuKuh46eJP/oaNGPTAkGFLxDGHwzj9YfwB0OEzOZBR/znSH03UU3e\nm6aJCpuEgybBZss2P6FHmECg4a81syMd2PUwmmai6SYvfl6CphdFPjf9081IR/XYtMj7xmkKiE5X\noA1v0rk9ssxb//iI8Ol56JrCMCLBisXQIwMtWiLjGkWCIx1N1yJDCuhaQ4f9hnGRGtZr2ncrmfpN\nMLFz505mzYqcnKZOncr+/ftj844cOUJeXh5ud2Q8u+nTp7Njxw6uvPLKhOmmD0khKyO1Q2VQSnH4\n8GHKy8vJyspi0qRJLX5dHz58mJKSEj7++GO8Xi9Dhgzhjjvu4MiRIxw6dIiamhpSU1MJBoNUVVUR\nDodxOiMj0aWlpZGWlobFYmHbtm34/X6Uigx0ZZomVquVc845h/r6ek6cOEFtbS3BYBBII3avqKHA\nboKpoRsW/nr0LSoqKggGg1gsFgKByM5vs9k455xzOHz4MIFAAK/Xi1KRW0kzMzNRSkWur/r9hMPh\n2C2ose8tPR2v1xtLr8OtMD1M1/W4cjYXva5sGAY+X2MkEG3OtFqtpKSkEAqFItdz6+sJhUKx/7PL\n5WLy5MlUHQoTrquLrByJOwlrUGSpifwPcEdeAzqEIsHBhr9vYAMbcDgc2O120tPTyc7OZs+ePQSD\nwVjTajS/zMxMxo8fT11dXez/UVwcaRcYNmwYNpuN2tpa0tLSmDBhAmPG5LF69QvU1tZit9s566yp\nDB02jLFjx8VakUDDNE2OHv2CyspKMtIzyBs7Lva9Nd7FQlxLUfTOEk3TsPTBGGs9Vf97QnvHiei8\n0tJSPv30U+rr63G5XIwfP579+/fj9UYe5uXz+Vrsx06nk2uuuYbKysrYujk5OcybNy/uzgilFPv2\n7eMvf/lLi3rpcrn4wQ9+gK7rrFq1itLSpgFE9Dii4OxaDBeEa4D9Q3jz7+tiS40ePZq8vDxmzZrF\niy++SEVFBRkZGYwePZqysjJyc3OZM2cOAJs2baKwsJBjx461/6Xl2NHG+mn6G2xK1kxU2I/DouO0\n6RiGLeGPtMjttgp/SBEImgQbXgPR16DJwYMHKC4qot4bQGFgYqDQG14NlIq8t1gd5OWNJxLbq0jQ\n0SwIUTQEIg0tJUopKnynImVBQ5kNp9CGLmSptowmzRBtn7Cji5hKUVdaj5mg9aUExacfbQfAbrMz\ncvTYxn0udm2zSUtHBy5ztuhAroHP27LlqqP6TTDh8XhITW086VssltiDeprPS0lJoba2NullOHz4\nMAcORH6WFRVFntI3efLkFvM//fTT2Em2vLycZ599loyMDMrKyggGg0Tv64fGk3B07AmLxdLiQBIO\nR/bEQCDAtm3bsFgshEKNY/NHNOyYYQ3qIxXOBE7VnYotES1T07SaM02TsrLWB49qqqqqlfvD+oH2\nAgmgIfhqfG2+XjgcjgsyoqL/p7q6Onbu3NlsZpPXoAHBtipo5GTu8/nw+XzU1tZSUNDYcbT5/7S0\ntJTy8nIcDkdDb/3GMldXV8d+3ZWVlVFTU8P69etj+30oFGLnzo+ZOHEiLqcjbj89dOgQBce+AKCm\nuhK32xk3P5FUd+8/ebI/1P+Oau84EZ1XUFCAx+NB0zQqKys5efJkwn3X6/Xyxhtv4HK5qK2tRdd1\namtrefXVV7npppvi8n/jjTdaDfDr6+t59tlnGTJkSLNAIqrh8uLuVLiwFva3/KF1/Phxampq2LNn\nDx6PJ5buqVOnSE1NjUt33759VFdXt7tdAFqpC5crlcCICoam5HLR2MuYkjst9iMmFIoMax4KBRpa\n8bTGS3XEX4YDhdMCDiMS/IbNhnFSlMbH2z/Ce2ovzmAQZzt3pioFytRwe45x+TU3NYwBEulIaja8\nhpu8hhvGCAmFFTsKP8MT8KCUjjJ1UDpK6dj1NCZkjsdUilCYhlabSHphU8UuTTZ14stjCQOJhm8w\ndo3S7/dxsrCAkWNOa5jT2IKutRPAtPo90KSPlYKQ2fXbeftNMOF2u6mL/hIk/ol/brc7tlND5ICf\nlpbWajrbtm1j+/bt1NR0fmyJ8vLyDn1ufqKqq6vD5XLFDhaNAwLF969oHICq/Z2n9fnRnUTRXsQr\nekNDU0Wr0xt1pDWn6X7RfHnTNGOtVn6/P65+QON+2NH9tj9LRv2P1n2AwsLCHitre99v9L3fHxmc\nLfo/TVTno0KhUCzYja5bUlLSIr/20qurq4v9QGkp1qOJbPcITtH2L9Hm+1vT/TMaRHWYguHGaSy6\n/KG4ydFOkMl6/sn/1Vd3qN5Fx9nx1VUwbritSQfz6KXRxsuSEL/tw0eO5b0jf2qRZv64BUzISiM6\nNkx8v7NIi6GJFhsdNRRWvPjCLnyGn4ab1RtG/qHhtWEkIKXHfVZo2HBw8ZlnYJoNdyhG+181u0sw\neldidLwSZcb3t4rdtdjwuaam6y3Q/SaYmDZtGps2beKqq65i9+7dTJw4MTZv/PjxFBQUUFNTg8Ph\nYMeOHXz7299uNZ0ZM2YwY8YMCgsLWb16dafKkJWVFVdJsrKyWp1vtVrjWgFSUlKw2+3ouk44HG61\nKT4ynoEe+2vvYND+/PjOh/3l8oNoqSP/n6b7RfPrmNGTabSXfvNf5NFLOm3tp00/93fJqP/Rug/w\nzDPP9FhZ2/t+o/PsdjvBYDC2DySq81EWiwWHwxFbFyA7O7tF/u2ll5KSwpAhQ6ivbztQ0HULd8z8\nEY+991iby6SkpMQFcU33z9zcXIA2Wj9a0jQttk5Pys3N5fjx46207LYuKyurzfE92nL+kAtJcTnY\nenQDZXVFcS0tnTU8O62dwLfhh2Mrvx1H5Y5i1jnN7+/rvsLCQp7q4rr9JpjIz89n69atLFy4EIAV\nK1awdu1avF4vCxYsYNmyZdx+++0opViwYEGLCpYMkyZNAoi7Ftra/GHDhvVyn4kIwzCw2WxkZWXh\n9/vRtMjDxKTPRKNk9Zk4dOgQdXV1sZN9OBzGMAxSU1Ox2Wx4A/XUmR5MbxjCBk2H6+hKnwmIXJqK\n9pnIycmJ6zMxefJkxo8fz6pVq6iursbpdHLuueeSnZ3d5n7a1n7cH/WH+t9R7X2/0ffZ2dlJ7TPR\nPP/rr7++zT4T3//+99voM9Ho3nvvBeCHP/whv/zlL+PmdabPBNChPhPnn39+3Do9Zc6cOSilOHDg\nQMJLtVlZWdx6661dymdK7rSk3F1x66238vzzz3Py5MlW5jZGEbqux/7Xubm5XS53T9JUfzlLJFlh\nYSGXXXYZGzZsYNSoUX1dHCEGjGeeeYalS5f2dTG6ZTBsgxC9rTvnTXk2thBCCCG6RYIJIYQQQnSL\nBBNCCCGE6BYJJoQQQgjRLRJMCCGEEKJbJJgQQgghRLf0m3EmhBAiGUKhEDU1NT06CqYQg1F0MLa2\nR1Bt26APJl544YU2h95OlsLCwh4fy6I38pB8JB+IPG9hICsqKmL16tWdHgFXCBFRWlpKXl5ep9YZ\ntMFEbm4uS5Ys4YEHHkjauO9t6Y0BcnprEB7JR/LpyaGoe0O07vfUKIEvvPCCpC1pD8q0w+Ewv/71\nrznzzDM7ve6gDSYsFkts6GohxFdHtO73VAuOpC1pD+a0R4wY0ennlQAYjz766KPJL07/0VtDafdG\nPoNpWySf/p3PYBiCvie3QdKWtCXteIP22RxCCCGE6B1ya6gQQgghukWCCSGEEEJ0y6AJJjweD9/9\n7ndZvHgxCxcuZPfu3S2WeeWVV5g3bx4LFy5k8+bN3cpv3bp13Hvvva3OW758OfPmzWPJkiUsWbIE\nj8fTI/kkY3v8fj933XUXixYt4s4776SysrLFMl3dHqUUjzzyCAsXLmTJkiUcP348bv7GjRuZP38+\nCxcu5E9/+lOXyt+RfH7/+99z3XXXxcp/7NixLucFsGfPHhYvXtxierK2J1E+ydqeUCjE/fffz6JF\ni7jpppvYuHFj3Pxkb09P6sh+DFBRUcGVV15JIBBImGZP7r+J0gbwer3cfPPNHD16NKlpr127lptu\nuolbbrmFznaZS5T2u+++y/z587nppps6dWtuR74PgIcffpinnnoqqWXuTn1KlPbevXtZtGgRixYt\n4p577iEYDCYl7bKyMhYvXsySJUtYvHgx5513Hi+//HLSyr1u3TrmzZvHggULeOmllzqc6KDwX//1\nX+qFF15QSin1xRdfqLlz58bNLy0tVdddd50KBoOqtrZWXXfddSoQCHQpr8cff1xdffXV6p577ml1\n/s0336wqKyu7lHZH80nW9jz//PPqmWeeUUop9dZbb6nHH3+8xTJd3Z733ntPPfjgg0oppXbv3q2+\n973vxeYFg0GVn5+vamtrVSAQUPPmzVPl5eWdziNRPkopdd9996kDBw50Ke3mfvvb36rrrrtOffOb\n34ybnsztaS8fpZK3Pa+++qr6z//8T6WUUlVVVWr27Nmxecnenp7Wkf34gw8+UDfccIOaPn268vv9\nCdPsyf030T67b98+deONN6qLLrpIffHFFx1ON1HaPp9P5efnx7b/nnvuURs3bkxK2uFwWF1xxRXK\n4/GocDisrrzyyg4fNxJ9H0op9dJLL6lvfvOb6he/+EWHy9uRtLtTnxKlff3116svv/xSKaXUK6+8\n0qn/ZUe+E6WU2rVrl7r11luVaZpJS3vOnDmqpqZGBQIBlZ+fr2pqahKmOWhaJm677TYWLlwIRH5x\n2e32uPl79+5l+vTpWCwW3G43p512GocPH+5SXtOmTWszoldKUVBQwMMPP8zNN9/Mq6++2qU8EuWT\nrO3ZuXMnX//61wH4+te/zocffhg3vzvbs3PnTmbNmgXA1KlT2b9/f2zekSNHyMvLw+12Y7VamT59\nOjt27Oh0+RPlA3DgwAF+85vfcMstt7Bq1aou5RGVl5fHc88912J6MrenvXwgedtz9dVX84Mf/AAA\n0zTjbqNO9vb0tET7MYBhGPz+979nyJAhHU6zp/bfRPtsMBhk5cqVjBs3rsNpdiRtm83GmjVrsNls\nQOvHyq6mres6b7/9NikpKVRWVqKUwmq1djtdgF27drFv377YMb4zevL40F7aR48eJT09neeff57F\nixdTU1PD2LFjk1buqJ/85Cc89thjaJqWtLStVivV1dX4/X6ADqU9IAdh+POf/8wLL7wQN23FihWc\neeaZlJaWcv/99/PQQw/Fzfd4PKSmpsY+u1wuamtru5TP1Vdfzfbt21tdp76+nsWLF3PbbbcRCoVY\nsmQJ//RP/8TEiROTmk+ytmfo0KG43W4AUlJSWlzC6Mr2tFVGi8WCaZrout5iXkpKSsLydyUfgGuv\nvWwjdikAABocSURBVJZFixbhdrv593//d7Zs2cIll1zSpbzy8/M5ceJEwjJ0Z3vayweStz1OpxOI\nlP0HP/gBP/zhD2Pzkr09ydSV/Rhg5syZQCRA7oie3H8T7bPnnHNOp8ra0bQ1TSMzMxOAF198Ea/X\ny4UXXpi0cuu6zrp163jssceYM2cOLper2+mWlpby7LPPsnLlSv72t791uKwdLXN36lN7aVdWVrJ7\n924eeeQRRo8ezZ133smZZ57JjBkzklJuiFxqmzhxYqdHq0yU9u233868efNwuVzk5+fH6lZ7BmQw\nMX/+fObPn99i+uHDh7nvvvt44IEHOPfcc+Pmud3uuANMXV1dwmG228qnPU6nk8WLF2O327Hb7Vxw\nwQUcOnSo3ZNvV/JJ1vYsXbqUurq6WBpNd7Cubk/TMkbTBuJ21q6Uvyv5ANx6662xynDJJZfwySef\ndDmYaK8MydqeRJK5PadOneL73/8+3/rWt7jmmmti03tzezqrK/txUx39BdeT+2+ifbY7EqWtlOKJ\nJ56goKCAZ599NqlpQyQQzs/P54EHHuAvf/kLc+fO7Va677zzDlVVVXznO9+htLQUv9/PuHHjuOGG\nG5JS5u7Up/bSTk9PZ8yYMbHWiFmzZrF///4OBxMd+a7/+te/dmkkzPbSPnXqFH/4wx/YuHEjLpeL\n++67j3fffZcrr7yy3TQHzWWOzz//nLvvvpuf//znXHzxxS3mn3XWWezcuZNAIEBtbS1ffPEFEyZM\nSHo5jh49yi233IJSimAwyM6dO5kyZUrS80nW9kybNo0tW7YAsGXLlhZBWHe2p2nau3fvjgtAxo8f\nT0FBATU1NQQCAXbs2MHZZ5/d6fInysfj8fCNb3wDr9eLUoqPPvooKf+P5r8Yk7k97eWTzO0pKyvj\n29/+Nj/60Y9aHPB7ant6SqL9uKmO/trvyf23vbS7K1HaP/7xj2OXUaKXO5KRtsfj4Vvf+lasc6vT\n6exw4NZeuosXL+bVV19l9erV3HHHHVx33XUdDiQ6Uubu1Kf20h49ejT19fWxzo07d+7k9NNPT0ra\nUfv374+1YnVGe2n7/X4Mw8Bms8VasmpqahKmOSBbJlrz1FNPEQgEWL58OUop0tLSeO655/j9739P\nXl4ec+bMYfHixbET4z333NPpitSepvnMnTuXm266CavVyo033sj48eN7JJ9kbM/NN9/MAw88wC23\n3ILNZuMXv/hF0rYnPz+frVu3xq5zrlixgrVr1+L1elmwYAHLli3j9ttvRynFggULyM7O7nT5O5LP\nfffdF2tdmTlzZuzaendED5I9sT2J8knW9vzmN7+hpqaGlStX8txzz6FpGjfddFOPbk9P6ch+HNXR\nE1xP7r+J0u5sWTua9pQpU3jttdeYPn06ixcvRtM0lixZwuWXX56Ucl9//fV861vfwmq1MmnSJK6/\n/vqkpNsdPXl8SJT28uXLueeee4DIpavOtCAmSruioqLdFrjupH3DDTewcOFCHA4HY8aM6VDrkoyA\nKYQQQohuGTSXOYQQQgjRNySYEEIIIUS3SDAhhBBCiG6RYEIIIYQQ3SLBhBBCCCG6RYIJIYQQQnTL\noAsmTpw40aVBPDrq0ksv5eqrr2bu3Ln88z//M9/4xjd45ZVXOp3O7373O5YtW9bh5bdv387kyZN5\n8MEHW8xbvHhxj24z/7+9c41q6sri+J+EiAgWtDrT6oCO1gG0glbH8ggIimKASAAtPiBaig5L7NKo\n6KiIKCjDwyfiwrGd6uCL0orMMLYwoggFRqS4ULRUKQKCqCjIQ4KQZM8HVu4kJAGk1Xbs/X1Kzr33\nnL33edxz9j33bnTbdeLEifD29oZIJIJIJMLcuXMhFotRW1vb5/WJiYlMVMqDBw8iPT39hWXIzc3F\n/v37AajXg5eXF4RCIfbu3Qu5XP7C+Q4UIoJQKNSQx9vbGx4eHhAKhcjLyxtw/v2x08WLF7Fr164B\nl9HW1qb2Bb3Q0FBUVlYOOL/XDUtLSzx9+lQtLS0tDcHBwQPOU7UvaCM1NRUffPABPDw8MHfuXHz0\n0Ue4fv36gMvri8uXL+PgwYMAXrw9FRUVwcbGhmn3QqEQixYtQm5uLnNOWFiY1lgpqmzbtg23bt0a\nmAI/IRs3bkRFRYWaXj4+PhCJRFiwYAEuXbrUZx6pqak6I22eOXMGR48e7ff1/Tk/JiZGZ9iFV0q/\nw4z9n1BbW0tTp059afm7uLioRZirr68na2trevDgwQvl8+mnnzJR2/rDlStXiM/n04wZM6ijo4NJ\nr6urIz6f/1J1JtJt18jISJ3RU1Xx9/enzMzMAZff1tZGnp6eTLTDnvUglUpp1apVFBkZOeAyXpRv\nv/2WwsLCtMpDRPT111+Tg4PDK5NnINy7d0+tXmtqauiDDz74GSX6ZWFpaakR+fLs2bP0pz/9acB5\n9tYX9uzZQ0uWLKH6+nomrbCwkGxtbdXSfkoSEhIG3G+uXLlCnp6eamnfffcdOTg4UGlpab/zcXFx\nobKysgHJ8FNx/vx5xg669JoyZQo1Njb2ms+f//xn+tvf/jZgOV70emXU6P5Ewn2ZvDZfwOwPbW1t\n2LFjB8rLy6GnpwdHR0esX78eHA4Hly9fRnx8PPT19WFpaYmCggKcPn0ao0aN0siHVL7z9fTpUwwZ\nMoQJaFNcXIy4uDh0dHSAx+NhzZo1cHR0hEwmQ2RkJAoLCzF8+HCMGDECQ4cORX19PTw8PJCbm8t8\nH97NzQ0HDx6EhYWFWrnKb71fuHABHh4eAIBz585BKBTizJkzAACpVIqIiAhUV1fj6dOnMDIywp49\nezB27FhkZWUhKSkJHA4HXC4XoaGhmD59us70vujo6MCjR48wcuRIAEBVVRV27tyJ9vZ2PHr0CFZW\nVti3bx9SU1NRVlaG2NhYcDgcZGdn4w9/+AM+/PBDnfbqycmTJ+Ho6Kj2lU/Vehg8eDDCw8Ph6uoK\niUQCIyMjJCYm4vz589DX18fYsWMRHh6ON998k/HklJSU4P79+5g+fTpiY2MBABcuXEBiYiIUCgWM\njY2xadMmWFtb49ChQ7h27RoaGhpgaWmJ2NhYXLhwAbNnz9YqDwDU1tZi2LBhzP8vvvgCp06dYuoy\nLCwM48aNQ3t7OyIjI1FSUgIej4fZs2dDIpFg8+bNjJ0sLS0RFBSEgoICdHR0QCKRYM6cOUhLS0Nm\nZiaSkpLQ1taGXbt24fbt25DJZLCzs8PGjRvB4XBgbW2NlStXIj8/Hw0NDRCLxRCLxdiyZQs6Ojrg\n7e2Ns2fPwszMDG+88QYuXryIWbNm9dkGXnd61mlPehtTDh48iOzsbPB4PJiamiI6OhpZWVlqfUH1\ny5NPnjzB3//+d2RnZ+PNN99k0m1tbbF582a0t7cDAB4+fIjIyEjU19dDJpPBw8MDK1euRF1dHZYv\nX46ZM2eitLQULS0tWLNmDRNzJSkpCVlZWSAijB49GuHh4Xjw4AHOnDnDtPcxY8Yw7enx48fYvn07\nKisrweVy4efnh4CAgD5tZmlpiYCAABw/fhx79uxBQEAAAgICMHv2bOzcuRPXrl0Dj8eDmZkZdu/e\njSNHjuDRo0fYsGEDYmJiQESIi4tDV1cXGhoaYG9vj6ioqF71k8vliI2NRU5ODng8HqZMmYKIiAjo\n6+tr6L19+3ZmzFIlISGB8dDo0svQ0BD379+HQqFAeHg4njx5gsePH2PUqFHYv38/SkpKcPHiRRQU\nFMDAwACNjY24du0aHj9+DAsLC5ibm6OpqQnbtm3DqVOnkJKSgkGDBsHAwAA7duzA3bt3Na5Xnn/3\n7l1s374djY2N4HA4CA4Ohru7O4yNjTFt2jSkpKT0q35eGj/jROal0JtnYtOmTbRr1y4iIurs7KTA\nwED661//Sk1NTTRjxgz6/vvviYgoLS2NLC0tqa6uTiMPFxcXmjdvHolEIpo3bx5NnDiR9u/fT0RE\nTU1NZG9vT9evXyciojt37tD7779PtbW1dOzYMVq+fDnJZDKSSqXk6+vLeCZCQkLo1KlTRERUUFBA\nfn5+GuUqZ8pZWVkUFBTEpHt6etKtW7cYnb/++muKiopijoeHhzOzbVdXV2a1kJ+fT4mJib2m97Sr\nlZUViUQiEgqFZG9vTwKBgPbt20ft7e1ERBQTE0P/+Mc/iIioq6uLhEIhZWVlEVH3akz5Wznz7s1e\nPfHx8aGioiK1etC2krGzs6MbN27Ql19+SYsWLWK8OAkJCfTRRx8xsqxdu5aIuj0ejo6OdOXKFfrh\nhx/IwcGBKb+wsJAcHByora2NEhISSCAQkEKhYMoSiURqnhJlu5g5cyY5OTnR1q1b6d69e0REVFRU\nREuXLmXk+eabb8jd3Z2IiHbv3s14dzo7O8nf35+KiorUVigWFhZ0+PBhIiIqLy+n6dOnU2Njo9oq\nefPmzXTixAkiIpLL5RQaGkqffPIJc/3JkyeJiKisrIwmT55Mz58/19pfTpw48UJes9cZCwsLEgqF\nJBKJSCQSkZeXFzk7OzM237hxo9Yxpb6+nqZNm0adnZ1ERPTZZ5/RhQsXiEi3Z+Lf//43+fj49CmT\nWCymS5cuERHR8+fPSSwW01dffUW1tbVkYWFBOTk5RESUmZlJLi4uRNQ9pkkkEpLL5URElJKSQitW\nrCAidc+EansKCQmhuLg4Ivrf6rempkZNFm0reCKinJwcJl2p79WrV0kgEDDnxMfH07Vr14hI3bO3\nbt06pq8/e/aMbG1t6ebNm73qd/z4cfL392f6o0QiofT09F71VuX27ds0a9asXvXKzMwkBwcH6ujo\noOPHj9PRo0eZYytWrKDPPvuMiNQ9Cz3HDaWt5XI5vfvuu9TQ0EBEROnp6fT5559rvV5ZN97e3nT6\n9Gki6vaIz5kzh9ra2hh7+/v7a+j1KvlVeSZyc3OZFTyPx8PixYtx/Phx/P73v8eECROYYCcikQhR\nUVE689mzZw8mTpwIAMwqb8KECTAyMsKYMWMwefJkAMA777yDadOm4cqVKygsLISnpye4XC64XC68\nvLyYZ4RLlixBfHw8Fi9ejM8//xyLFy/WWbazszMzO62qqsK4cePUohW6ubnBzMwMJ06cQHV1NYqK\nipj9FB4eHli1ahWcnZ1hb2+PoKCgXtN7MnjwYKSlpQEAvvnmG2zcuBF8Pp8JZR0aGor8/Hx88skn\nqKqqQkNDg1pkOuqxyistLdWw13vvvYeioiKNb8FXVlbC3Nxcp11UMTQ0RG5uLnx8fGBgYAAAEIvF\nSEpKgkwmAwAmVoOyzpqbm1FRUQE7OzuMHj0aQPeKcMSIEbh58yYAwMbGhomV8MMPP8DMzEzNU6Js\nF3V1dfjwww8xfvx4/O53vwMA5OTkoKamBosWLWLs0NLSgubmZhQWFjL7Z3g8HpKTkwEAZ8+eVdPL\n398fAGBhYQELCwsUFxerHc/JycGNGzeQmpoKoDtgj2qUQaUXZdKkSejq6oJUKtVqP3Nzc/zzn//s\nj6l/FSQnJ8PExIT5r/QGAUBeXp7WMWXFihWwsrKCt7c3HB0d4eTkxIRA7w3VWBzPnj3D0qVLoaen\nh2fPnkEgECA4OBhXr15FS0sLs39IKpXiu+++w+TJk8Hj8Zj4DxMnTkRzczOA/7UNHx8fAN1RIp8/\nf96rLIWFhdi0aROA7iiTL9omBg8erPbfwsICXC4XCxcuBJ/Px5w5c2Btbc0cV/aLv/zlL7h8+TKO\nHDmCyspKdHR0oL29HSYmJjr1KywshJeXF9Mf9+7dCwBYu3Ztv/SurKzUCONdU1MDb29vEBHkcjne\neustHD58GAYGBhCLxSguLsaxY8dQVVWFiooK2NjYaLWD6rihhMPhQCAQwM/PD87OznBwcGD2X2mj\nubkZ5eXlTLTct956C1lZWcxxMzMz3L17V+f1r4Jf1WSi581MoVBAJpMxsdxV6S0csGo+I0eOxKxZ\ns3D16lXMnDlTowy5XA65XA4Oh6N2jMvlMr/t7e0hlUpRWFiI4uJixMTE6Cybx+PBzc0NGRkZqKio\nYDqJklOnTiE1NRX+/v4QCoUwMTFBXV0dgO6OtWDBAuTn5yMtLQ1Hjx5FWlqazvTe4PP5WL58OSQS\nCb766isYGxtDIpFAoVBAIBDAxcUF9fX1veZBRDrrpCdcLlejjnpSV1cHqVQKMzMzjXOV9aAsr+dA\np5RFW/0p5TEyMmLSs7Oz1R5xKPMAgNGjRyMmJgYBAQGYNm0arK2toVAo4OXlhfXr1zPn19fXw8TE\nBPr6+mqDzcOHD5lJkCqqbVLZpnrKeuDAAYwbNw5A92RF9RzVPLXpquSnDIf9OqDLTgA02plq+01O\nTkZZWRkKCgoQHR0NW1tbbNmyRWde1tbWqKysRHNzM0xMTGBkZIRz584BAA4dOoSnT58yG4yV7nEA\naGxshKGhIRobG8Hj8Zj89PT0GNkVCgVWrFjBBHbq6urS2FjaE3199dtDbW0tTE1NmcexvXHjxg2N\nKJdDhw5Feno6SkpK8J///AcSiQRisVgjhPaSJUtgZWUFJycnCAQClJaWMnro0q9nH2psbIRCodCq\nd1NTk4a8HA5HY/O2ubm5znEwLi4OZWVl8PX1ha2tLWQymc52ojpuqBIbG4uKigoUFBTg6NGj+OKL\nL3D48GGt53K5XOjp6anpWF1djbfffhuDBg2CQqFQu6f8HLyWI4auSuXz+Th58iQAoLOzEykpKXBw\ncMDUqVNRXV2N27dvAwAyMzPR2trar4h97e3tKCgogI2NDWxsbFBVVYUbN24AAO7cuYNvv/0WM2bM\nAJ/PR3p6Ojo7O9HZ2Ynz58+r5bN48WKEhYVBKBT2Gf3Ty8sLaWlpKC4u1thfkJ+fDx8fH/j6+mLs\n2LG4dOkSFAoF5HI5Zs2ahfb2dvj5+THPQru6unSm92XXwMBAvPHGG0hISGDKDgkJgUAgABGhtLSU\n6aD6+voak4Te7NWTsWPHMqF8tdHS0oKoqCj4+/tj0KBBcHR0xNmzZ5nVd3JyMv74xz+qDUY9sbW1\nRUFBAfN2SmFhIR4+fKi2elJy6dIltUiUPZk6dSq8vb0REREBAHBwcMC//vUvNDQ0AOi+GQQGBgIA\n7OzscO7cORAROjs78fHHH2t4HQAwN5abN2/i7t27Gnbi8/k4duwYk8/q1auZ9q4LbRPpe/fuMRMS\nlt7RNaaUl5fD09MT48ePx8qVK7F8+XKUl5cD0N4XAOA3v/kNxGIx1qxZozYRv3//PkpKSsDlcmFs\nbAwbGxt8+umnAIDW1lb4+/sjOzsbQO9jX2pqKtra2gB07w9QvhnG5XK19nd7e3vGO9ba2oply5ah\npqamT5tcv34dZ86c0Zgk5OTkYNmyZZg6dSpWr14NkUikYZOWlhbcunULoaGhcHV1xYMHD1BTU8OM\nI7r0s7OzQ0ZGBjo7O6FQKBAREYGMjAw4Ojpq6K3tLbq+xpee5OfnY9myZZg/fz6GDRuGgoICph/p\nsqcqTU1NcHZ2hqmpKcRiMdauXYvvv/9e5/XGxsaYNGkSM7mpr6+Hn58fWltbAfwy+uxr6Zno6OjA\ne++9B6C78enp6SElJQVbt25FZGQkhEIhurq64OTkhODgYOjr6yM+Pp7ZrPbuu++Cy+VqrF6B7tnw\nhg0bmGNSqRQeHh4QiUQAgAMHDiAyMhJSqRRcLhfR0dEYM2YMzMzMUFNTA09PTwwfPlzDZe/t7Y3Y\n2FhmBt0bU6ZMgVQqxezZszVWkIGBgQgPD8e5c+dgamoKV1dXXL58GVwuF1u3bsX69evB4/HA4XAQ\nHR0NHo+nM12b7qro6+tj27ZtCAoKwsKFC7Fu3TqEhIRgxIgRGDVqFObOncsMPi4uLoiJiUFnZydz\n/bBhw3Taqydubm7Izc1lbqCq9aBcobi5uTGv7C1YsAAPHjzAwoULQUQwNzdHXFycVj2U/8ePH4/t\n27dj9erVkMvlMDQ0RFJSksZKrKGhAQYGBmqPl7RNPNetWweBQICUlBT4+fkhKCgIgYGB4HA4MDY2\nxqFDhwAAq1evxq5duzB//nwQEdzd3eHq6srcIJSUlpbiyy+/hEKhwL59+zTCD4eFhWH37t2YP38+\nZDIZHBwcmEdWunQeOXIkrKys4O7ujtOnT8PExAR5eXk/70auXxB9LSjCwsJ0jikCgQA+Pj4YMmQI\nDA0NERYWBkC9LyjHDSVr165FRkYGNmzYAKlUiq6uLhgYGMDd3R1Lly4F0P04befOnRAKhZDJZPD0\n9ISnpyfq6up0yrtw4UI8evQIfn5+4HA4ePvttxEdHQ2g+0b88ccfg8fjYdKkScw127ZtQ0REBNMu\ng4ODmce7qty7d495LKmnp4ehQ4di7969jGdCKdPMmTORl5cHT09PDBkyBKampoiMjAQAZtNxVFQU\nVq5cCZFIhN/+9rd455134OTkhJqaGpiZmenUb9GiRbh//z58fX0BAO+//z7EYjH09PTw8OFDrXqr\nMmHCBAwePBiVlZX9uimHhIQgJiYGR44cwfDhw+Hm5obq6moAgJOTE6OXLoYNG4ZVq1Zh2bJlMDAw\nAI/HY17J1XV9fHw8duzYgeTkZHA4HOzevZvZqJuXl4d58+b1KfdL5WVvyvh/oLW1lWJjY5nNcTdv\n3iQ+n/9KZcjIyNC6MYilm9bWVvLw8FB7LfbXhIWFBT158uSll1NdXa11AzALy+tORkYGRURE/Nxi\nvDAtLS3sq6G/FIyNjcHj8eDr6wt9fX3weDwcOHDglZUfEBCAxsbGXl9L+rVjbGyMdevW4fDhw5BI\nJD+3OK8c1efDL5MDBw70uvmYheV1xcPDA9nZ2bhz5w4mTJjwc4vTbxITE7Fly5Y+H4+/bPToVYxQ\nLCwsLCwsLK8tr+UGTBYWFhYWFpZXBzuZYGFhYWFhYflRsJMJFhYWFhYWlh8FO5lgYWFhYWFh+VGw\nkwkWFhYWFhaWHwU7mWBhYWFhYWH5UfwXHHagiWMokngAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113c762d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (bmr_ax, hgd_ax) = plt.subplots(figsize=(7.5,3),nrows=1,ncols=2, facecolor='white')\n",
    "\n",
    "bmr_ax.plot(log_bmr_quantiles['x'], log_bmr_quantiles['y'], 'o-', color='#83BA70')\n",
    "hgd_ax.plot(hgd_quantiles['x'], hgd_quantiles['y'], 'o-', color='#83BA70')\n",
    "\n",
    "sns.regplot(data['log_bmr'], data['transmission'], ax=bmr_ax, \n",
    "                        line_kws = {'color':'#547BD3', 'alpha': 0.65, 'zorder': 10}, \n",
    "                        scatter_kws={'color': 'gray', 'alpha': 0.65, 'zorder': 10}, logistic=True) # Standard logistic regression plots, includes raw values\n",
    "sns.regplot(data['hgd'], data['transmission'], ax=hgd_ax,\n",
    "                        line_kws = {'color':'#547BD3', 'alpha': 0.65, 'zorder': 10}, \n",
    "                        scatter_kws={'color': 'gray', 'alpha': 0.65, 'zorder': 10}, logistic=True)\n",
    "\n",
    "hgd_ax.set_xlabel('Host Genetic Distance (Patristic)', size=12) # Label\n",
    "bmr_ax.set_ylabel('Transmission\\nWith BF >= 3.0', size=12)\n",
    "bmr_ax.set_xlabel('Log Body Mass Ratio (Donor/Recipient)', size=12)\n",
    "hgd_ax.set_ylabel(None, size=0)\n",
    "\n",
    "bmr_ax.tick_params(axis='both', which='major', labelsize=10, width=0.5, length=3, top=\"off\", right=\"off\")\n",
    "hgd_ax.tick_params(axis='both', which='major', labelsize=10, width=0.5, length=3, top=\"off\", right=\"off\")\n",
    "hgd_ax.spines['left'].set_linewidth(0.5)\n",
    "bmr_ax.spines['bottom'].set_linewidth(0.5)\n",
    "hgd_ax.spines['left'].set_linewidth(0.5)\n",
    "bmr_ax.spines['bottom'].set_linewidth(0.5)\n",
    "hgd_ax.spines['top'].set_linewidth(0.5)\n",
    "bmr_ax.spines['right'].set_linewidth(0.5)\n",
    "\n",
    "hgd_ax.set_ylim(-0.1, 1.1) # Format axes\n",
    "bmr_ax.set_ylim(-0.1, 1.1)\n",
    "\n",
    "sns.set_style('white')\n",
    "bmr_ax.set_title('A', fontsize=18, ha='left', x=-0.12, y=1.08, fontweight='bold')\n",
    "hgd_ax.set_title('B', fontsize=18, ha='left', x=0.0, y=1.08, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('../png/Fig4.png', bbox_inches='tight', dpi=600)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Figure 4: Cross-species transmission is driven by exposure and constrained by host genetic distance.  \n",
    "For each pair of host species, we **(A)** calculated the log ratio of their average body masses and **(B)** found the patristic genetic distance between them (from a maximum-likelihood tree of mtDNA). To investigate the association of these predictors with cross-species transmission, we treated transmission as a binary variable: 0 if the Bayes factor for the transmission (as inferred by the discrete traits model) was < 3.0, and 1 for a Bayes factor >= 3.0. Each plot shows raw predictor data in gray; the quantiles of the predictor data in green; and the logistic regression and 95% CI in blue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}