plot-H5N1-HA-beast-tree.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot H5N1 beast tree using baltic for fellowship applications"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, subprocess, glob, os, shutil, re, importlib\n",
    "from subprocess import call\n",
    "import imp\n",
    "bt = imp.load_source('baltic', '/Users/lmoncla/src/baltic/baltic.py')\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.patheffects as path_effects\n",
    "import matplotlib.lines as mlines\n",
    "from matplotlib.font_manager import FontProperties\n",
    "\n",
    "import numpy as np\n",
    "from scipy.special import binom\n",
    "\n",
    "import datetime as dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rcParams['font.serif'] = \"Helvetica\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create empty baltic tree object\n",
    "myTree = bt.tree()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<baltic.tree instance at 0x117ad4d40>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "treeString='((A:1.0,B:2.0):1.0,C:3.0);'\n",
    "bt.make_tree(treeString, myTree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load trees\n",
    "tree = bt.loadNexus(\"/Volumes/gradschool-and-postdoc-backups/post-doc/stored_files_too_big_for_laptop/H5N1_Cambodia_outbreak_study/beast/dta-2018-08-23/HA_2018-08-23.mcc.tree\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeface='Helvetica'\n",
    "mpl.rcParams['font.weight']=300\n",
    "mpl.rcParams['axes.labelweight']=300\n",
    "mpl.rcParams['font.family']=typeface\n",
    "mpl.rcParams['font.size']=14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make color dictionary to set colors based on geo or on host species  #other cambodia option, lighter: #ecbc91\n",
    "#colors = {'wild':'#fdb771', 'domestic':'#13257a', 'unknown':'#47c5b8', 'human':'#e7e17e'}\n",
    "colors = {'wild':'#CEB540', 'domestic':'#4E83AE','human':'#DE4428','not yet annotated':'#CBCACA'}\n",
    "color_key = {'wild':'#CEB540', 'domestic':'#4E83AE','human':'#DE4428','not yet annotated':'#CBCACA'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in colors from annotations file \n",
    "annotations = \"/Users/lmoncla/Documents/Fellowships/data_query/H5N1_tree_annotations.txt\"\n",
    "colors = {}\n",
    "alphas = {}\n",
    "sizes = {}\n",
    "\n",
    "with open(annotations, \"r\") as infile: \n",
    "    for line in infile: \n",
    "        sequence = line.split(\"\\t\")[0]\n",
    "        host_group = line.split(\"\\t\")[3].strip()\n",
    "        \n",
    "        if host_group == 'wild':\n",
    "            color = color_key['wild']\n",
    "            alpha = 1\n",
    "            size = 3\n",
    "        elif host_group == 'domestic':\n",
    "            color =  color_key['domestic']\n",
    "            alpha = 1\n",
    "            size = 3\n",
    "        elif host_group == 'human':\n",
    "            color =  color_key['human']\n",
    "            alpha = 1\n",
    "            size = 3\n",
    "        elif host_group == 'unknown' or host_group == \"\":\n",
    "            color =  '#CBCACA'\n",
    "            alpha = 0.6\n",
    "            size = 3\n",
    "        else:\n",
    "            color = '#CBCACA'\n",
    "            alpha = 0.6\n",
    "            size = 3\n",
    "\n",
    "        colors[sequence] = color\n",
    "        alphas[sequence] = alpha\n",
    "        sizes[sequence] = size\n",
    "        \n",
    "#print(colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def convertDate(x,start,end):\n",
    "    \"\"\" Converts calendar dates between given formats \"\"\"\n",
    "    return dt.datetime.strftime(dt.datetime.strptime(x,start),end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for k in tree.Objects: ## iterate over objects in tree\n",
    "    #print(k)\n",
    "    #if isinstance(k, bt.leaf) == True:\n",
    "        #print k.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to remove the top group and reduce length of tree:\n",
    "ones_to_remove = [\"a/northernshoveler/mississippi/11os5805/2011|a|h5n1|amnongsgd|ha|4|cy166736|usa|2011-12-08\",\"a/mallard/illinois/12os5080/2012|a|h5n1|amnongsgd|ha|4|cy186667|usa|2012-10-20\",\"a/mallard/ohio/11os2216/2011|a|h5n1|amnongsgd|ha|4|cy132445|usa|2011-08-19\",\"a/mallard/ohio/12os2218/2012|a|h5n1|amnongsgd|ha|4|cy187005|usa|2012-11-16\",\"a/ruddyturnstone/delawarebay/317/2016|a|h5n1|amnongsgd|ha|4|mf613724|usa|2016-05-24\",\"a/canadagoose/delawarebay/601/2016|a|h5n1|amnongsgd|ha|4|mf046544|usa|2016-05-25\",\"a/mallard/wisconsin/2576/2009|a|h5n1|amnongsgd|ha|4|cy097350|usa|2009-10-06\",\"a/mallard/california/2535p/2011|a|h5n1|amnongsgd|ha|4|cy133981|usa|2011-08-12\",\"a/mallard/california/2536v/2011|a|h5n1|amnongsgd|ha|4|cy134005|usa|2011-08-12\",\"a/mallard/california/2542v/2011|a|h5n1|amnongsgd|ha|4|cy134021|usa|2011-08-12\",\"a/blackduck/tennessee/17os0306/2017|a|h5n1|amnongsgd|ha|4|mg279989|usa|2017-01-13\",\"a/mallard/minnesota/ai114429/2011|a|h5n1|amnongsgd|ha|4|cy206115|usa|2011-09-15\",\"a/mallard/california/2531v/2011|a|h5n1|amnongsgd|ha|4|cy133957|usa|2011-08-12\",\"a/gadwallduck/tennessee/17os0341/2017|a|h5n1|amnongsgd|ha|4|mg280333|usa|2017-01-14\",\"a/mallard/wisconsin/10os2672/2010|a|h5n1|amnongsgd|ha|4|cy133077|usa|2010-10-02\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n"
     ]
    }
   ],
   "source": [
    "# remove from tree \n",
    "for k in tree.Objects: ## iterate over objects in tree\n",
    "    #print(k)\n",
    "    if isinstance(k, bt.leaf) == True:\n",
    "        if k.name in ones_to_remove:\n",
    "            print(\"yes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/lmoncla/miniconda3/envs/mumps/lib/python2.7/site-packages/matplotlib/font_manager.py:1320: UserWarning: findfont: Font family [u'Helvetica'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAEsCAYAAAB679d0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsvXlgW9d57fvDDJAAeDiDAziC1EBR1mhTgyVZsmzZcRwPsRMndZypSdqmt+lt09zXl6bvJm1u26Rp07Rpk/a6SdM0o+3aceJBlizJki1b1EhKpChwHsEBIzFP7w8KCEmRFEkckCcy1l82AS6u/Qkfzt7f/vbasng8HieDDN7FkK+2gAwyWG1kkiCDdz0ySZDBux6ZJMjgXY9MEmTwrkcmCTJ41yOTBBm865FJggze9VCmStDf3w+A2WwGwOfzkZWVNed7F3ptNs9ytfj9furr65f896e/loqWhXikHJuhoSGGhobQaDQ3vOb3+9HpdABcu3YNgLq6uhnviQU6kKlKkSn082qYzjMyMgKAyWSa8Zrb7YaJfyA3N5f1d/7LEke5PKScBEajcVm/53Q6sVqtyaB7PJ7kz5cLj8dDIBAgGAzO+fr0f4SFXktFy0I8i/37s7ESsTl79ixarZaGhoYFeUZHR7HZbMTjcRQKBdFoFNdEJ5vLXyKCgNb85UXpyc7Onvc1lXySSDC2KB4xkHISBAIBAHJycpb0e1arld7e3uQ3UygUSlUKoVBo3n/kpfKIASnxpBobt9sNQFVVFX6/n4GBAXbu2MHwyAiDA0NsKgO5LEyO4i0q1n9qTo7pTzubzQZAcXHxDa/Fos8tW+dykHISRKPRZf2eRqOhvr6exsZGYOpxDFBaWrpsLUNDQ/j9fmpra+d8fbHTkVS0LMSz3OnQSsQmPz8/+XfmwpUrV8jJyaGqshK5XE4wGOSa1Yp3chKtoYxIzj5kCjXDw+dxRVtItKTJZLIkx/Sn3ejoKJOTk/j9flQqFaGgj0rDS5jXPI65/qFlj3M5WLXpUDp4jEYjSmXKQ1rVMc2eJsJUggBMTEwsW4vP58Pv9ye5ZmP2dCzxzW8w6JGFe7DU1jLYfYSwe4Cmpl1kZWVx7tw5HE4ntbW1eENGJux2WlvDZGW9xB23305WVhYarXZGIiSg1+uJxWL4/X4AVDIPOsUwk2Ovw29aEjgcDmBqUKvN43A48Pv9FBUVrbiWxIc3FovddPE3FxKvnTt3DpvNNuMbO/Hhz8/PX/pgpnEEg8F5nyZ2ux2Hw5H8Eunq6sLlcrFprYz9m/vp7b+dUqEXue8CF7v15OfnJz/cvb299Pb24vf7mZyc5OrVq1RXV+OVZRP1XWJTJajz3jPj7yWSrKGhgX8+3o+ldAv37dyPVl+x7DEuFyknQUFBgRg6ROEpKCiY95su3VoSa5zpFRxBEJbMY7FYsFgsM54iiXWXVqvF7XbT39+PSqVaEq/BYCASidDZ2Tnn09LlcuF2u1EqlbhVBSgL64ma8rHaLxK6oEanK8DpyycaCRGJOZicnASYwaXT6RgZGUGpVPLOO+9wzXgHSnkNm02/BGYmwfTYVOXrqCzQky2sXdKYxELKSTA+Pg6kVr4Ti2d8fBy/35/SN+ZytSTWOLW1tSmXSGdPiRJVIZPJRF9fHzKZbN5SZwKzk8Xj8STLm3ONS6fTkZOTg9/vxxouJyZTI0NJSd1ezkyE+PMDtchlshuebnP9XZiaCh67aqekuIjNux9DJpPNWz5ubESUL6/lIuUkmO8fdTV4pKQlFZ7ZlTOtVrtkjtnJ4nQ6KSgooLi4OPmUmZ4owWCQcDhMW1sbMWUYNFrixDFE3eyqLWJjYyMymSw5nVrMF831msecECvGYiD1VWQGomN25Wz6msBsNmO1WpfFazAYZkyzrFYr3d3d1NTUJKc1Wq2WvOA4RkWQrNgkVaYG1pTpefbIabQqBYR8GLWKOad6kWgcpUKWTK6RiAG9wcjH37sLuVy6zQkpJ0HiMZbqFEQMnkQFJFVIaUyzeQRBYNu2bTf9ndnJEggEkmuLBCwWCz6fj5GREcrLy9Hr9TQ2NpL4ArdanVgvvcO5K/lExq4xWLgHgyLEngL3DWNyh+CFThkbhQChuIxqfYQOdwSlxztndUis2IiBzMI4TVpWm2d2siQ+8LM/dA0NDTdUtRIIhUJcvHiRrrPvAKDI34FKFcdgMMz6a3F0Chl5OLnoNBJHwR1b1/MPd9QRhzmTQKzYiIHMwjhNWqTGM19sEsky16LdbDZz9epVcnNzMQoCA3IVvmgEm22IUCiEUqlEr3GxtfIknY49rNXqOevTEFZo8QUCKBXzT4HEio0YSDkJllMGTBePIAio1WpJaFlJnrk22GbD6/Uuqq9qdlXJbDbT0NCAXC7nbKeNHI0Mg96QXKcEI1pGPSX0DExy/lI7d+69i7N+A8++3c2H9t+27DGtJFJOgkS9+MZH5MrzTE5O4vP5kv0oq6klnTyzP/Stra3IZLJk89tcewkOh4NQKJSMjdfrTdb0YWqfIBqNMjExwXBBEwZNiM1Z43SGcjlgyqWiKIf6+nqutrdTVVWFSq2etth9AGOFG1NZJ+vWrsXkCKERFt6wFCs2YiDlJJhrvrdaPDKZTDQeMZAuntkl1NmYay9hfHycQCCA0Wikv78fr9dLf39/svza29ubTJzcgwfxyqC2Nodz7VGeueLnr++4g/ICA9nZ2YRCIdasWXODrl27dhEOh9mmVN507GLFRgys2j5BMBikt7c3+f+Jb4axsbFla5mcnCQQCOD1eud8fbGtzMvR0tHRQWVl5YyfpWu/YXYJdTFl08SHPZEgGzduJDc3N/m0qK+v5/Lly3R2drLNnIVWKae4yMBdvhFOj8rwev1QYEguaOf7EC92J/uW2idwuVzA0pvFLBbLjP9PfPBS6R1KJEGqlQcxtMDyY7NUntmVoLmSIhEbtVpNZWXlDWVWn8/H97//fa5du4bHaUdp0OJ0xagsNHCbJYe3r3SyprJoxca0kli1EulcJTxg3m/qxWC+MmACi21bEEMLrF6JdK69hERsuru76e3tvWERnTgjEAgEaBBClJYWoFSpCF45T1fcQOGWHcvSMh8yJdI08WRKpAtz+P3+G57A0+Hz+QiHwxQWFhIIBHjz9de5O2xna/s7jK2t5NQpF50DNo73w+++ZzNb11WnpAdukRJpqlMGMXn0er0o2/NSGpNYPInYzLfj7PP5MJvN2Gw2XG43d+3bx3sffJDmM2e43H8PjsFhotEoblUO4wE5Ry/0JpNgrgM0KzEmsZByEkQiETF0iMITiURE4xEDUuJZTGxKS0upr69HJpMxNj7O2Pg4vX19FBUVsWf/AV555RU0kTfJVt1JlizA0aNH+dnpbtRCMR+8vZwd2zYtSY9UkHISiHGmVyyeYDAoynlcKY1JLJ7FxEalUlFdXT1jyhQKhejs7EQQBMrLy2m/3M0+43miDhPXnDLycgoIdrbx1xMhfrZ5I0r54srUYsVGDKScBHl5eWLoEIUnLy9PlN4hKY1JLJ7Fxkav199QaWptbQWmKjmNm/Zz5coVLly4wN0HDmBpPUOpeowfeOW0Hj+C8ptfouwzf0zefY/eVI9UkHISJLbPU637isEzMTGB3+9PufIgpTGJxbOY2Mze54CpStOWLVtmHIbJysqivb2dH/7Xf7GpYT0HHvk4weYRfnBuhAO//zVychTc7CMuVmzEQMpJIFYrrBg8+fn5ojwJpDQmsXgWE5uFeo+m8+Tn5/O5z32Ozs5Ompub+fc3ehjwKsiKx/nemXH+6ZO7FsUjFaScBHa7HUg9o8XgsdvtojwJpDQmsXiWGpuRkRE6OzuJouA/z4yxtb6MT9+/JamlvLwck8mEx+Ph6JtjoMomJ1uPPSJf1N8QKzZiIOUkEKNrUywetVotStVhOVoSbSDT2y8Sj/yEe8Vi2zZmt2CsRmxeeeUVzpw5Q83mvRQbjQyMOufU0tTUxJD/LJNyA//jfbcTj8cXtTAW63MjBjJJIJKWREUlFvu1feBSHSHE1DMXx81iM7ufq7i4mKpcBXvMeeivd3vO1iIIAh9/+EDy/xe7V3BLJUHCJzPV/nAxeDwejyjHK5ejZa7DKWIZ8q5UbKaXRtVqNS6Xa8bvNDc3J59qJpOJsbExftUyhiOmY+f6Sj5x78zzA+0DE+RmaynOvdF3VKzPjRjIHK9Mkxap8SwmNtN3k81mM06nk4GBAWqqq7lw4QKTk5MMTsroC2TTVBlG5h3C1j+GV7DwwpkutphAOW3H/ssvdnJbdRFf/ODOtIxJLKTcYzA+Pp7sA1ltnvHx8ZSsCsXUIjWepcZGEASeeOIJDhw4gEKpJB6PMzAwQGtHF66gjJeueSkoLOT3fC08lO9lR5UepVxOJBpNtlE8taOUJ/fNbaglVmzEQMaLNE1apMaznNhMn+LV19dTXFxM8PBhyp1vUF+3hrdbXBhNFazbsg6zz8eLR97g4K7t9LiiPHx3EwvYDkmihTqBlJ8Ec1l5rBaPlLRIjSdVDkEQ2L9/Pzt37kSxZgND2UZ0Oh295jXI5HLO9TgYD2fx4mU7nuDNncrFio0YWDVr9nTwRKPRGdWZ1dQiNR6xYmOxWFAqlVy6dImenh6ys7OprKzk9IicaMzI/85xUFGauyg9UkFmOpQmLVLjESs2ZrMZs9lMYWEhPT09+F55lqxv/r+s3XiIgqCd7H8+jLtmLTn/8uxN9UgFKU+HHA5Hsmy22jwOhyOlK43E1CI1nnTEpqWlherBa7xacZAW9QYG7QG+sevz9D78Gc42N9PV2ZlcJC/Es9rIlEjTpEVqPGLHxmAwsHfvXo6GAthcbop9HXQY64grBF5qG6Fpg45YLIZKrZ7z9FimRJomnkyJdGEOMWOTWCjf/cCDGI1aVHKI59VSFB1BEQvz9NvjOHxhJubRfUuVSFM9jC4mj06nm/fxu9JapMaTrtg0NDTwn//5n3i8I1TU52CIjFJSsYcaWRhBpyQSiRCLxW449ipWbMRAyk8CMQ2vUuWRkhap8aRLiyAINDU10diwjvfv3Ugou4RX2l30emQ0brwNS139nOe+xdIjBlJOAp/PJ8pcUwweMa3ZpTImsXjSGZvE9VKvnzjHtUgRj/a/Sq6zny/8++t89dlzi+ZZLWQWxmnSIjWedMamoaGBS5cuUXLyJ2wUSlnj6eKwKkJQv5mwXA7cBTBjOialhXHGdyhNWqTGk87YJPqM/q/TjkGj5C3jvTQfO8Zj44eJylW0tOwiEIryf17upmlNCX/y/qZby3dISjbm71Zr9sVypDM2giCw9+C9XLlyhYnhYbRaLaf0DUxqDAQHnNSZDMjlMlTX7yyQQgt1Ahlr9jRpkRrPSsVmfHyc7du309/fT3t+PTJklJhK2LbZwn+sqU+emchYs6eJ591izb5cjnTHxmKxJI2Ay8vL6Q8BMgj88//HxXwDr218lKbGOnZtqJRMZQgyV7i+a3hWQosgCFgsFo4cOcLk5CQVA6/ysdGznNNXErFPMBo8y/cnIuzaUCmJA/YJrJo1ezp4EscB57toeiW1pMIz+6zv7Au057oLYTFaViI2FouFK1euEI/HGSnI5Zmccu4aOMU7JVvZpxih7uGPL4pnJZEpkaZJSyo8s52jxVhErlRsBEHgwQcf5Otf/zp39bdxsewO8oIOXIVV7P38h5NPgEyJNE08t0qJdLZz9OwD+8vVstjY9AzaGOjrJUc/1drg9/vpccaorTJTpIvdVIsgCHz2s5/lxHe+Sf61c7xVvI12eT5vXunn7m1rknpSHZNYyFizp0nLSvAs5tbKBNxuN36/f4YRbuKCv1AoRDAYTPqD/vhKmMfankFjgL57PkgwGORX/VryrHY+cFvO1N+emACFApnHherVn1P6od9G6feiu+12YGrqds8ffIHTp09zqrkLnVJOad6vXScy1uxp4nm3WbPf7AK/2Rzznebq6+tjZGSEmpoaAHSOSQLxCGHd1Afe6/WyhmGi9jD/9m9XWT/eyyOBYRxCEY6NO6g99SqDbgealncw/+g4yoKpmysFQeDQoUMcOjT1d6ZPxzLW7GniebdZs8++wG8hDAwM4Pf7qaurSz5BjEYjDQ0NmM1mLl68yOTkJEqlknhOOaXvuxuFRk6tLgu3dQiL3kk4HMJtH2XT+FWQycjxe1CuacB778PIIiFM73simQDLHdNqIGPNniYt6eZJTGUWun4Jfj1lisViycPt0+89TvCMjIzg8Xgwm82UKN0YtQo8Hg9Hznfx9oSOJ5uq2FCczdq1a/Ha72Ms5KF11E73mQFqT5zlY//wDWRLmIpmrNnTxPNusmbv6+tjeHj4pkmQmDIZjUYCgQC5uTMPwSeudN27dy9arTZZsnS73Zw8dYoJp5vqstu41mylT63g7GQ+k3Ij71uXDdjpzq7AJ8+CaATki2/LuKWs2aX0rfluehIkrmK92aI4MWWqqanB5/NRWFg444pXtVpNbm4uarV6Rs0+MVX6+ZtW/AEDQ6g4WCojFlIhj8IDdzVRlHcv90140KiUyFRL60u6pZ4EicPb2dk3+k2uNI/T6cTv91NYWLjqWlaCJxQK0dHRseDvJjbWpsdmegk2GAzy8ssv43A4MBqNGAwGhsJ6BmMFNGjHcHZfRJNbT1CmJZBr5i/3FBGNxSnKm1o0l+Uvr/dHrNiIgZSTQCznZTF4VCoV4XBYElpWgqeiomJR5dEEx1yxsVgs7Ny5k1/+8pdcvnyZ/Px8hkx3gQKqqmswGVR0eb3k5OZhdFv53i/aiRtK+ErN0nasFzum1UDGmj1NWlaCx2g0LqoylOBIxMbpdHL+/HlsNlvyw1hbWwtAT08PdWordoOJjTUWZLWlbPT7CYfD9Pd7sHp1eF1uItEYSsXy92Qy1uxp4llNa3ap80yPjdVq5Z133kGlUiV7kBLNb+Pj49w1cpoJVwE/eLuUD+xpZENtbnIBu67NitMXSikBxBqTWMj0DqVJi1R4Es14FRUVBAKB5P9XVlZiNpvRaDRoNBr8fj+dnZ2UlJTwf+PbMeYV4LD5+ccX3uRT9T+kdN1vk118F/VVZas+JrGR8R1Kkxap8FgsluTCeLYDXV9f34xu1VAoRGFhIetVIyisx4gH3OSFO9CpXJw5/ys+/U+v0D9kW/UxiY2MF2matEiFJ1EJcrlc+Hw+SkpKCAaDnDhxInmDfWNjIz6fj8nJSZqbm/F4PDg05aAxMBwu44fXHqdnsoCSvGyKC25utpvuMYmNjDX7u4RnOofFYqG0tPSG92zZsoXt27fjk+mIl29BJpMhCLk8vmMrT+6o4H8+uEVSYxILGWv2NGmRGs/02AiCQENDww17DImzAMfbRwnKZECcSkFFc3Mzer2e77QNMuyFapUDr6GS//XYTgpzlrfjm7FmTxNPZjq0MMdiYiMIAh/Yv5kfvvQmFblqGIZz1mtT+wzGCoIyAxf1peD14PQGl50EUpoOpfyJSdhrp9ofLgaPw+HA7/dTVLS4TsZ0apEaz+zYTK8Szcbt27ailMvQaDRYrVbi8RjtQ25w2QiW1AAytlTkULeIyzgW0gPSOFeQKZGmSYvUeGbHZrozxGwkHKdhKlmOvv46ztr3IgPi0Qgy4MmDm1PWIxVkSqRp0iI1ntmxSawLbtZ2sWXLFrZs3sx6WQ+RwYuEA14gzr8fbiEYXv7uvJRKpCkngU6nE806PFUenU6HVquVhBap8Sw3NoIg8Nhjj/Hoob3kertRXO8WlQXsPP7Xv+Cff9m8bD1SsWfPmG+lSYvUeOaKTeJADYDNNrUJplKpGHT4UcriqLVZvD4g57d2VnBg23q+8pWvcOz0ebIMOfzk/ASgQKteXiPcLWW+lZhnpurwIAaPmPbjqWqRGs9csbFarbS2tqJSqbDb7TgcDpRKJSc9pahiPupUVlTKKuIeBwMP34F6224OfPpPcblcHPKFCcXlbDPJOX76HA5fmLKyMravKV+xMYmFzMI4TVqkxjNXbO69997k3QJ+vz85PamxudEo5Yz2KjDaenFNhAjJ5Pj8QfKVSl555RXa29uprKzk37uK6YwUoFdE+ew7X2Bg3RrKv/qdFRmTWMj4DqVJi9R45oqNyWRKOtL5fL5kp2iiO9vpXJc8n3z+I/8LmUzGxLPPcvnyZYaGhhjv64Vdn0Imj/PBLfkUHhsnOpazYmMSCxlr9jRpkRrPcmKT6Dvy+XzU19djtVrx+XxJS8e96jB+RztFRQLvvfcRYrtPI1vk4lsKLdQJZKzZ06RFajypxiaREHq9nrNnzwJgGrvKnV0vXP/2/yzy7MVvfEnJmj1zcd+7hEcsLYIgsGbNGqLRKNqtu1H4PMTXLX3jTEoX92Ws2d8lPGJqycvL45577kEe9BKzvoOnpIafHm/h8b2LO+opph4xkPKTwOVyJW22V5vH5XLhdrsloUVqPGLGRhAE1q5dy7gvQCdaDtsU/PR4G5e//Ef0dHcv6r5ksWIjBjIl0jRpkRqPmLERBIHDhw9PLZJ3PMCbAyPsdJ1CVhTH3tlBbl4eBoNhwSmPlEqkmd6hNGmRGo+YsREEAZlMhlKp5J6DBzlYZ8CYp0HbcgZXVyff+M+XaD5zhsuXL8/43WGHdwaPVHqHMtbsadKSTp6F2qAX4hArNnq9ns985jNcuHCBM83NeH0+1uy9h5/kl3CxzY4/v4QDwRDDDi/vdB1ne3UeJzrsvNpm53fv38S9W6pR93agrFzYQnKlkLFmT5OWdPIs1Aa9EIeYsTGZTBw6dGjGHQkmkwnjxUucbzvFc89FGS28A71GRU54HJkPslVyfvWzH3D4aQ/dRdv5Y9/z7Pirf0hZU6pIOQnEOicq1rlVMSy/pTSmuXgSbdBzYb6LO0ZGRggEAjMWx4n/nt024XG7GRoeRiaTEZGp6Alk87FD2yjKy5lTy/RbdRoaGlh37hwajQZdth4ZMpQKGfVAVX4/Tz99ia2BSQQ0ZLuWlsjpQmZhnCYtq8Uz38UdgiDc8AG2Wq3YbDZqa2sJBoOEw2G6u7sZFrZgUKsoDXQyTAEdvih7RpwU5eVQUFCA0+mkubkZjUZDPB4l7j2DpfERsvX5Mw7kzEZtbS1btmyh+8JZSlrfRvOev116MNKATO9QmrSsFs/0izumPxWcTieBQGBGu4LFYpnxxOjv78dms6FQ2RB0ApWVlQQdoAnJwGenpcXNyMgIHR0daDSaqfWBZowyzbO0nw+jFu5cUFviaeOKx/hvtZ6BX7Tww3XrlhkV8ZCxZk+TFinwTH8qGI3GG3qHEvbrMPUBTSy0px+CrwcO+P3oNMrka1lZWQwNDaHX66m07CMeNKNQV93w98fGxjhx4gQulwuZUkWwfCd3VBoI2drYWniZmrAdeHxZYxMTGWv2NGlZSZ7p3/jT7zie/lQYHBzE7/fPe6nH9C7ShV4bHBzEYDDQ3PzrE2UyTe2cv9fZ2UlfXx8ajYaoVoc6NkDPsJ1HDzzOqDXCpqrdiwtAmpGxZk+TlpXkWcwFfqnEJhCKkMgPlUpFQUEBDz744E0rVOvXrycSieDxeCguLqbFE+HsxAYekitp2PbJZWlJBzLW7GnSspI8i7nAb67YTH+C+K/br7e1taFWq5NPnSvDfrr82Ty4o46Hd29IapldFZoPDQ0Nyb8xdGGYXo+XmtJyIiFpuM9Bxpo9bVqkxON0Ojl58iTxeHzGrvG5c+eS1SG73U5nZyft7e189J67KDn/OqdL1pLjCkL+Zkq08WVpmZ4slZWVRKMx9Fla+ifGUhqTmMiUSNOkZbV4ErvJMHVVU25uLufPnycejyfvKU7AYrEkj1fabDYGBwcRBAF5fye5g530jvmpRsaf9pzAnv8R2NaY0piUSiVGY9aSx5RuZEqkadKyWjyzF74DAwMMDw+j1+vp6urC5/Ph8XhusGQMhUKoVCp0Oh2eikZcD32YrF/8AldeOS19bZjWbVq1MaUbGS/SNGlZLZ7Zc/WOjg52795NOBzm7NmzNDc3YzKZZnz4/H4/w8PD6HQ69Ho9p93ZeC/3UVZWhr7jApsvHydytgxu3yi52IgB0domcnIWd8A6nTxi2oanqkVKPEajkaKiInJzcyksLEStVs9ou2htbaWsrIytW7cCcLHlMiPD/VivdVCogM133M3bdgWBlhYmJiYYcIVR55bx+J71qzYmMZFyW2E0GhXNOjxVHjGt2aUyJrF4otEoer2eBx988KYXgOs0Kvbt3cMf/uEfYosq+KJ2P+eN9fi8XqLRKJdtQZ45bSUSXX6sxYqNGMhMh9KkRWo8idjMVdo0m820trbe8Dsej4eCfAG3d4Ct69chl8vp7u6mY1QgX6dI6fK+W2o6JGX78dXUIgWe6ecOFoqNIAhs2bIluSuc+D273Y5SoeB9Gwpou3CCy6e8+Hw+7j7wGMpIgPPnzlFWXr6seGes2dPEkymRzsT0cwdLiU1iunTp0iVqqqups1ioqqqit7cXv9+P2RhFp9UTjUbxer03YZsbmRJpmngyJdKZSJw7cLvdHDlyhHg8TlVV1Zzv9Scv7O7H5rezttrC7/zO7/D60aMAaNRqzGYzza+f4FuDAh9ZD488fOeybVNuqRKpWPbaYlmYL8bpYCW0SImnr6+PQCBAWdnMO4gTrtShUIje3l6Ghobo7++n6Y6t9J28SNZHBO7av5+2tjYKCgrwBwIIOgX39rzJlvs/mJJvkFRs2UGEJBDj3KpYPHK5XDQeMSAVnlAoRCQSueFD29fXh0wmS/YV5efnE4lEGNPVEyvN4egbvfzZwCVqXv4JZd/8L5SNjay73v+fapu4WLERAykr8Xq9y54Xis3j9XpFWRNIaUxi8IyMjDAwMDDvHkowGCQSiVBWVkZ+fj77N5qJ16wHpRxnIEQwywgKpShaEhCLRwxkFsZp0iIVnmAwyPj4OBs2bGDnzp0z1ktmsxmr1cqRI0eYmJhAJpMxNDREiaClyKBmb60em1NL9KP/i/qc3JS1TIeUFsYZ36E0aZEKj8ViobCwkMnJyRtik9gzsFgsrFu3jrKyMkpLS7nYNcKwK0hvdycn+mN84+gA4UhUMmMSGxlr9jRBec+zAAAgAElEQVRpkQqPIAjU1tbi9/vnbVGorKyku7ubS5cuAXCXqQCnL4LFYqGzzUmOVoni+saYFMYkNlJ+EkxOTiZttlebR0papMYTCAQW5LDb7dhsNq5du8abnW5aOkcxGo1sqy0kEI3TOewUTYuYPGIgc3FfmrRIjWeh2FgsFqqrqwHo7u5mT2GE/JFBZO4iVAod1UUGcrI0omkRk0cMZKzZ3wU8lZWV5Ofnz8sjCALvec97+PM//3Pa29sxmooIFeaj0Gk5fXGCjuFJ8o06UbQkICVr9pSTIGGvnWpDlBg8iWuEEvdwraaWdPPMPkG2kC+pRqPB4XDgdrvnjY0gCDz66KNTfyseQxbwgU7HrpoctteVoJDL0j6m1UKmRJomLenmuVk79GwIgnDTdpKmpiYA4q89T/6x5xn89JfQFpaQp4VYLIZcLpdcbMRApncoTVrSzbNYtweYao+4cOECIyMj8/YOTb/Me6wsC9XDDzPpnGS91onhx/9IX+B3qLpzv+RiIwZSToJUjaXE5MnOzhZlwSWlMYnB09fXx9DQEDDVQjEdk5OTdHd309PTQzQaxWQyYSvYgEHt5zN33oGrvQX5cB8jndfwCIXJaUz30BjBcAxTzq9tHBN9WzKZjHjMR8z1EjWNH8OYe+NUTazYiIGUk0Cs00FSOYEllhYp8VRUVNDa2kosFrvBoMvtdtPd3Y1CocBisZCVlUWOXo1CnUs8HidnbSPH3vNJHKMOis+do7Ozk0gkQpfSQlhtoDhbxs5SGTolvNwjI0ev46md5RAZJz55goBn15xJIJVTZZCxZk+bFinxNDQ0cP78eUpKSuY06CorK+P5559Pxm5XXQlvDwT5wVuDfGRnOTU1NWiv30/s8XiIRCKoA25iOXo63GpKzeUUGtTYZA7KSoppbFwDNBIJ70GpmvuKVrFiIwYyC+M0aZESjyAIbNq0CavVSktLyw2vt7a2MjQ0hNvtpry8nJ+2uFHJ5WRpFLjdbowGA5a6OnJzc9myZQswVZ2yWq3IFCrGPCHcciO/+8jMBJsvAcQYk5jILIzTpEVqPIIgLHhxtsvlQiaToVKpqM/2IYuECEwM8crLXh544AGuXr1KNBrl9OnTrKmvZ3xigtHRUbpUFuxBOZHYIE3rzCgW2SJ9Sy2MV9t+fDZHxpp9blRXV1NcXDynY7fZbKazsxOTycSmTZumTSlrOXHixJTbdyBAS0sLoVAIk8mEIAjU1dVRMSlHrRfYVFuy6AQQa0xiIWPNniYtUuNZKDaCILBv3z60Wm2y7Or2BfGHIjQ0NNDa2kpJaSkFBQWcPXsWAI1Wy/bt21PSA9KoEmWs2dOkRWo8N4tNPB7nzJkzXLx4EYPPzYsTBUyq9Dz9e3fROR7k7vpyGhtzkyVo3fWFcip6pIKMNXuatEiN52axqa2t5ciRI7hcLu7MlnPb2CAdFbfxl1/7Bh3aTbzUcYIf/8n7WL9+PaFQKOWnrVixEQMZa/Y0aZEaz81iIwgCBw8epL+/n2A8jmZkhMKJq1xqaaHu/iY6xwI888wz+N3XqDe7maz+KDW1818Kshg9ib+72siUSNOkRWo8i4lNQ0MDGo2Ga9euMe5yEYlE2LZ7P/dvLeTV85309g5Rlt2FMthJ1LYbUkgCKZVIM8cr06RFajyLiU2iH+m+++5jz549VFZWElNMzf3zs7REIhHau4tQ/psbX3s3g4ODKem5ZY5XStFvUwpapMazmNg4nU7OnTuHwv0DKio2sOWJT/D9/z5CJBbHmCuQo6zgcjb8b+NW/mRDeUq2KVJooU4gY82eJi1S45krNrPvLGtra6OtrY3H7+zBM9JLQc0n2FRbjNsfJkefRXHpOoqddnZF4mTptDO+zSPROMc67LxvdyM1pTffrJSSNXvKSSBGNUYsnkgkIkpjlpTGJBbPXLGZfeulxWKhuLgYt6cVFLk8//zzdHd3k2syc3Q8j3XGICVxGxtLXqZvrJQO2wbWr1+PNktPXKnleIeD6vLxRSWBWLERAxm3iTRpkRrPXLGZfuvl9LuKnc4mrFYrlmAQn8/H+NgQFmWQar0BtcxIv03LmHOS0ckRmnbs5LlLdv7mU7u4rbEBIVsz159Py5jEQspJYLfbgdR3/sTgsdvtouwYS2lMYvEsJTbTD+wkdoz7+/uRyWS43X5iWQ+QpQxRlhPmrw4PgEyWbMy74nbj8XimdoTlKoJ+L1mKHqorTcRUteiyi4jH49jtdkpLSxe8e3mlkCmRpkmL1HiWG5vE3QUbNmzAarVy+fJlurq6WLduHQ0NDZgGnfjDMdxuN9aefnr7B8k3aOnsHaDTuB2TRs29lcXENRsgNjUdi8fj6LOzUUtk1zjTRZomLSvJM/3Q/XwYGRkhEAjMOF55swP605F4Omi1Wmw2G+3t7VNGv8psIjF47WI7t91+JwXqIm6vMtC4cRP/8tY42xrX8MihzcTjcfx+f3LKZbVakSsUybPLq4mMNXuatKwkz2IO3Ws0mkXFxul0Yr12DYVCQV19/Q03yZhMJt7znvfwox/9iBdffJH6pnoaiocYoI5xXxStSoZKrSZXq+X/uW+qDJqYKvn9/uRYHA4HMpmMUZsNs9mMYRVLphlr9jRpWUmexRy6t9vt+Hw+ysvLZ/zc7XbT3NyM0+mkp6eHU6dOcc+B3dSt2cjRt85zrMvLWnMBn7p/W1JLSUkJn/3sZ3nxxRdRc4GaXA8NG+4hEleiWuQ9ZokxyRUKtKt8V0HKSZCw1061P1wMHq/XK0rvkJTGJBbPfLHp6+tDrVajUCjo6elhYmKCurwXiLm7eKv3TvQaJYLafYMWs9nMJz/5ySVpmF6B6u/vT/IkXlstZBbGadIiNZ75YqNWqzGZTFy8eBGlUsmaNWuQxZqR+d4mZyKHxqoWZKGtjPTpMFXskNSYxEJmYZwmLVLjWSg2fX19qFQqysvLOX36NP/6fAFGYz7Doy1YZXfhjuWRY8rGVDG3lsTOs0qlori4eFEOgLfU8Uqpbwitlhap8cwVm0RVqbS0lNtuuy1ZWLDb7cnF9hvXHLjjoNYX3qAl8eG/fPkyAwMDHLz7bvr7+vC43Vjq6hb0gLqlNssS9toLHeJeKZ7JyUl8Ph/FxcWrrkVqPHPFZvoVrwnk5eXR1NREY2MjTqeTcPhljEY1MvcgLS3jjI2NAVBYWMi5c+fo7u7GYDDg8XgYHh6mtLQUl8uF2+1esC9IrNiIgZTLF2LaoafKIyUtUuOZiyNxxet8Rx0TvUWnT59OJst0nkSv0djYGMZYGNO548giEVRq9U0/3GLFRgyknARZWVmi2aqnyiMlLVLjmY8jMSWay7RMo9Gwe/duamtrKSoqYmxsjO7ubgJDfeS88mOabmvkiSeewGKxkD/cTXZODqhU5OXl3bSsK1ZsxEDGmj1NWqTGM19s5poSTYfRaKSmpobXXnuNoaEhDAYDG594gtGaOkYvXqKs1sLuA/czVF3NZDhEtkZDUVHRioxJLGRKpGnSIjWe+WKTmBIttL9SWzvlP5Sfn487LOeN02fYurEBgKvXuvn74yP86+/fS2HO4r/ZpVQizRyvTJMWqfGkEhtBEFi7di0ajYZxw1recgi88ItfcOT1Y4SicQ6uMdJnbbtp/9JsPbfM8Uqp2I8nODLW7PNzpBIbk8nEsWPHqMrTUlFeyVtHhpGVFGF2hKjJU2NzTOKLyBbdkCcF060EUk6CWCwmhg5ReGKxmGg8YkBKPKnGpqGhgfz8fEIhDxfePIrP5yPYfZFjxWvI0ymY8EUpLYhx5x2L1yMVpJwEYvTqiMXj9/tFOWMspTGJxbNQbBIO0xqNZl6XOkEQ0Ol0jIyMsGbNGuRyOWNjY2i0bnY11HBpxM+nHty5JD1SQWZhnCYtUuNZKDYWi4VAIIDNZluQo7KykpKSEmpqajh8+DBjY2Ps3yAjT69mb62KHms7G2+7bdF6pILMwjhNWqTGs1BsEqfH6uvrF/QIValUBAIBBEFg7969FBYWJkuc0WgUpVLJ5OTkos4t3FILYynZj2es2RfmECM22dnZVFdXs2HDlNPEc889R44gIJfJcLvdtLe1Ya6ouGnripSs2VN+EjidzqTN9mrzOJ3O5CbMamuRGo9YsfH5fEktOTk5lJWV8Y1vfIN+d4whVxiNRrMoLyGxYiMGMtbsadIiNR4xYlNZWYlMJuO1117D7/ezdcsWTCYTn/nMZ3i23YVSqeSRQ3sWrUcqyFizp0mL1HjEiI1Go8Fut3P8+HFcLhc1azbwX+fsRJHzhffvoKJw8S0QUrJmT3k65PF4kjbbq83j8XiSLbqrrUVqPEuJTeLccUdHxw2NdbFYjJKSEvLz8/mvH/4Ahy+C2x+md9SNMWtxxlsJPWLERgxkSqRp0iI1nqXEJnHueC4YDAaqqqoIhUL4lLlolDKiMdhZm7tkPVJB5nhlmrRIjWcpsVGr1XM6wwWDQdrb2/H7/Wg0Gh44tA+9TktMJsdgMBCPxxfdmnFLHa8U62SQGDwGgwGFQiEJLVLjESM2FouFiooKvF4vdrsdx5gNY1UVbUOTBJsv8U/H+nlq/3ree8fNrRWlcKIsgZSTQIwb5MXiCQaDovGIASnxLCY2088czwVBELBYLPh8PjQaDaFQiPZBF0VGHUdffYn8nI382+HLNFQWUWNa+AyxWLERAxlr9jRpkRrPYmJzswM2AGqVCo0g0LRjBwBtbW2cPHkSJVE+uLWAbI2SbNnNS7EZa/Y08WTcJhbmuFlsEgdsOjo65n1PIBhkfGyM7Ov2jKdPn2Z0dJT7Dh0iGo0il8nw+2++AJeS20TKJVK73Z60Dl9tHrvdjsPhkIQWqfGIFZvCwkJ0084G22w2rly5wtmzZ3EE4rSP+Dh8qW9ResSIjRjIlEjTpEVqPGLFpra2ltra2qSxbmFhIV/96lcpKSnB67CxtsKM3X/zKamUSqSZLtI0aZEaT7piYzKZWLd2LSaTicryMsLRGJvW1iyZZzWRsWZPkxap8aQzNg8/8gjf+ta3aGi8jZ91qqi7OMEXn9i9oJOEWLERAxlr9jRpkRpPOmNjMpkoLy/nzIl/Rm3+NNmabDquXkWuUGAwGKirq1sUz2ohZSVerzdp173aPF6vV5R5r5TGJBZPumNTXV3NiKeQbQWXGHV7aO/qo2PQQefI3O3bYsVGDGQWxmnSIjWe5cbG7XZz7ty55BmBtp4RLgwH+cyD2VSYfn0wpqmpia7OLq6FKynPCvHza3Hy9T6qTFoOzaNHKsj0DqVJi9R4lhObYDDIq6++CkBj40ae65QRCMeIIOdsxwCjA12cH/TT54Yv/9Zu7rn3HjYOD6PWaMkWfNy/owFL2dwf9luqd+g3bUNopbRIjWc5sbFYLKjVamKxGGvW1LMhOEZ+TjbVxTkoPEMcOXKGAVUB9pCanz87hltm5MJwiL/+ZBO3377woRkpbZZlrNnTpEVqPMuJjSAI3HvvvQQCARobG2lshOHhYQBefPEN2tsu8fuPDRGOqvHo/phJrwLkv27Si8VixOPxORv3pGTNnnISZHDro6+vL3kDZeJ+gitXrmAbGQRAoZSz/8AB7HY7jT09vPnWaXr6B6kzm9DpdHR1d6NSqQgGgwwENFSWFHHS6mTfhjIeLSlZtXElkLFhTJMWqfEsNzYJT6IEEvV9pVLJlm27uTqehVJbgubyZVpbW3E6nbw1psMrz0YZ7UWr1XHq1Cm0Wi0OWQ5jxvW0TEygVSmIy6XxHZyxZk+TFqnxLDc2CU+i2bdOllz/Br/vfR+f0Xna29uLNqigseE2zth1TPjjFF5/zZ5dC0COTkmRkM0dNUs7jZYuZEqkadIiNR6xY5M4czD9DmWz2UxXVxdabRizPsbwyChMTJ0jVmcZiMrVKMNehp0R+p1hKovHeUoC1aFM71CatEiNR+zYqFSqG2xTBEFg3759yZ8/etdWdlYbyMnJIeh1E7P3EFFmE0GJhjCukT5ePPYOb16+eddpOpFZE6RJi9R4Vio2DQ0N9Pf3Y7FY2LZtG0JBMV999hy58n6UV88jMxYRVek5VJ+FZ7iTnw6FCTDC0589kLK25SJjzZ4mLVLjWanYJI5g9vb2otFoeLv5HPKYkk0bN7DepMOvlGMsyiUy2kFXVxeevCJiCgUn3jmPyagjLy8PtUYz7xHPdCBjzZ4mLVLjETs2C13GkbgDGWBdvYV19VOLen9FbrK6dGIwjN/vRzt+GUXxGt45cYm79u4lHA4jl8spKipCqVyZ6lFmYZwmLVLjWWxsEoft5/uQJ7Qs5CM6fbE8HT6fL1llKigooL+/n/iEjTq1EZvNzT99+9t8/o//GIC+3l5qamtvqlcMZBbGadIiNZ7FxsZisSw4FRFrTJFIhOLiYsrKypIO1RqNhkAggEwmw7GCZr0Za/Y0aZEaz2JjkzhsfzMtqTpK5+Xl0djYyLFjx3jxxRcJFW1g/Y73c35gkp11OrJX8I7jlJMgEYxUKxhi8DidTvx+P4WFhTd/c5q1SI1H7NhMRzwep3fUTVXxzS3Zp/OUlpaye/duAoEAfdl1fHB7KeHrVixr1q5NSedSkLFmT5MWqfGIHRu3201/fz+Dg4Mcu3CZa7EaPrDFwCvtAR68o4737ViDddiBVq2kPP/GJjmVSkVBQQHr16+nqamJU2++idPpTCZ6LBYTxU1wMchYs6dJi9R4xI5NX18fw8PDDA4MYNaHefT2MsZGrejQ4p6w0dIS4muv9lBWYOT/PLV3Xh6YmoLdtW8fWVlZRKNRIpHIiiUAiJAECXvtVPvDxeDxeDyilBOlNCaxeMSOTXZ2Nn19fZg37+f2dSZc9pMUFe3g05X65Hs/tN1EXW31gjyzx6RQKFY0ASBTIk2bFqnxiB2bwcFB3njjDe6rMdM+4OKdXi2/fdeU39DAqIOXrjg4uFbAOzFEi30YbzDKgE/Jw7vXo5DLJXW8MlMiTZMWqfGIHZuNGzdOtUZk/RKj/1c05Bp48/IYR6+Mcq2rl357AL1KRui68W7r0CQ/OHGNrusH728p36FbzX5cLC1S4xE7NhqNhsHBQdq7olhqYaS/k3DpVq7ZPfzlozvY0N+HVqOisqqKwsJC1q2PcccGG3WluTN4pICUkyAUComhQxSeUCgkGo8YkBJPOmJTVlZGVB5FKNpGeTyH2loD92l0aDVq9u/fP+P3lAo5taacOXlWGylPh8LhsCilNzF4wuGwKBUQKY1JLJ50xCYvL4+NFgdZnKGstBS3awytSk4wePMeJbFiIwYybhNp0iI1HrFjk52dzcaNG3numAuLxUK2vpWmtR2MuJ6iQPADC99Wc0u5TSTstVPdFRWDx263i7IrKqUxicUjdmzMZjPbt2/H5XJxxpWDPgLqYT0vD47zIZWFqhudF+fkEevMRSrIlEjTpEVqPOmIjcVi4cc//jE+9yRao56sgJX7zXECY91A46J5VhuZEmmatEiNJx2xEQSBu+++m9LoMJ6ON+n21bAp7zIh5Y0t0BMeP+0Djjl5VhsZa/Y0aZEaT7pi09TUxNmzZ/F6vTTbt3H+wj4OWpRJn6IEnj1v40K/h5984UHUSkXGmj1dPBlr9oU50hGbxOH64eFh3LEICpmMTkec2d1C+9fmcWBTNWqlYk6e1UTKSZCw1061510MHq/XK0p/jJTGJBZPOmPT0NBAeXk58lgPLiXEY0ac/jBGtRxzRUXyvdPXJGLFRgxkFsZp0iI1nqXGxu1209zcjEajwe/3J2v6CqUKd1RNaWlpcgdaEAQ++MEP8tVvfwtfjo7isI+vv9rLH91VgmpkZM4PupQWxhlr9jRpkRrPUmPT19dHOBzG6XTS19fHQzvOEo1pOe78BG/bFBQWFrK5riz5/qqqKr78Pz9PZ2cnY5Nhes/ZiWgELJa57y+7pazZE5c3SIEnJydHlAMoUhqTWDxLjU1iY62trY2R4T6y7gIIkO9ycm9dBWH7IC0tduLxGDHnM1SuexyTaUPS5vGhe2+uRyoQbU2Qqt+mGDxiXkmUqhap8SwnNmq1msceewyPx8OF7l+SrYswMjjIjiYzQ8NDuFwuHKMX2NtwgTPHu/jXjkewlObx5d+6E9X1BbDT6cRqtSanVYmq0MDwKBqlDDfZ3LbOQnHO6lWLpGELnIFkUVhYSFNTE6HQXqxWK9qCQQLBILW1tbjdbg53d3HkTDa2SSNKIcS1IQfnLraQpZ5KgtbWVmQyGQ0NDbjdbi5dusTguJtz4UoOmqNYKrX85I0O/scDt63aGFNOAr1ef/M3rRCPXq8XpfQmpTGJxbPU2ExOTmKz2ZL+QwkvoTVr1gC/boVuaGjg3LlzaDQajEYjsXgc+Tx2jx0dU65zI+RTVaPj3ESYDRY5H92/LsXRpQbR3CZS7Q8XgyfhqJDqTTVSGpNYPEuNzcjIyJw7urO1CIJwQ9v0dJjN5hnW7S6XC1V2Fr+1rYC+CR9aBdhtAxgqKwkEAmi12qUMSxRkSqRp0iI1nqXEJhgMMj4+PtVsF7HTe/6L1G79U9TagiVrme5Gl5+fz9q1a7l8+TIul4uQYwxlQR3dvb28fnmI26sEbrttI2qNZsnjSwWZ3qE0aZEaz1JiY7FYkt2m7qEfEXCd561f/R5nX/0YfT0tjI+PE/KPEo9Fl6Qh8dR48sknGZ+YQKPVEgqHcbvc9I17cUXVqERy6FgKMtbsadIiNZ6lxEYQBGpra7FarVzqNLK5SsbIhIt67SjO7r/HITMyHu+goPrDmNd9YslaBEHgnnvuAaZ2jHfu3Mn7pvmUrjRSfhKIafmdKo+UtEiNZzkcJpOJLK0cY3acvLw8QhEdPYPgdNiAOG0d3XzzP/6Ojv6xFdGTLmSs2dOkRWo8y4mNXq/n3nu/wDM/+w90in6yVXDyvIIKcy3dIw6scS2P1L3CcLeRevPSnghixUYMZBbGadIiNZ7lxsZkMvHwo7/F0NAQR44c4X17T5Kt8/Daxa1kxa30eMrRy5e+Xsn0DqWJJ9M7tDDHcmMjCAKlpaUEg0FOvt5LXq4Rr9eLpWSc+hw/J0dLiEUDyBWLL2/eUr1Dv4n24yuhRWo8S43NXJd0NDQ0cOXKZvr7+yktLaC1340zy0CrvYrAZC9ZOWuWpEcqSHlh7HQ6U/aqF4vH6XQm7/xdbS1S41lqbDQaDZpZ9XpBENi1axf5+fmMjo4S947Q4aohGNcuKQESesSIjRhIOQmUSqUod0uJwaNUKkVxWZPSmMTiESs2+fn5HDp0iGAwyN33fYzN5kJuy2snHluap5FYsREDKauY/W2xmjwajYZodGkbOOnSIjUeMWNTXFyM2Wxm3OFmKJLPpdFNhMOTqDWL9xISKzZiIOUngcfjSdpsrzaPx+NhcnJSElqkxiNGbDxuD06nE4/Hw+OPP87RCz3EfBfYWXxxSQmQ0CNGbMRApkSaJi1S41lObNxuN0ePHqW9vR3b6CjVW/aBTM4H7t5OMBjkPU1reM4aolrwJN0l4vE4BQUFlJSULLhDnSmRpsjzuc99jr/927/lK1/5Crt37+buu+9Ocvzwhz/EZDLx0Y9+dEW0/KbwLKdE2tfXlzTOHbXZ2KWXM2Z38Z3vfIerV68yOTnJxvoc7q1s49ypK0z4KiktLcVSW8v5s8cIhDQUFRaSIwg0Ns4047qlSqTztfcGfUN0Nf8pAW8/2mwzNdu+Csz/yFxKm/Df//3fz8shxlxTSpbqYvEs1Zo9GAwyPDzMnj17WL9+PaFQiBdeeIH+/n7C4TBKpZI1a9agz1bTNxEGdS3uETudnZ3U5PyKgvgwp3rvo8JsJuD309LSQiAQSJ4sc7vdwFSi9ff302JXUltVwcfvvz3lsS4VKSXBwJV/xDXeBsCYeur8ajQaQ6GQ43W2E49NXdAQmOyl7cTH0RrqUShmLkN0Rgvl6z87p1X35z//eZ566imGhob4whe+wPnz53nqqac4c+YMly5dSr4vFAoljwHq9XoeeeSRVIYlKUt1sXiWas1usVi4cuUKVquVWCyGTqcjOzub7du3o1aryc/PT367O51OXnjhBcrLsygvL2d49DVkyipqamo4f/48k9ePh+qzsyksKWPCE0QenCrXDg0NMTExgdu8j/DABNFo9DfvuqZ4fO4mqEQCzPf/szGXTffOnTs5deoUw8PDlJSU4PF4sNlsFBUVzXjff//3f3P77bfz0Y9+lD/6oz9a4ggWp+U3nWep1uyCIPDggw8mkyCBSCRywwk1q9XK8PAwlZWVNDQ00NqazdWrV9nYGCUYCtFxzcpIwQ7Wq1S0D4DHa+cj684S0rw3eVKtxxmhoLBoxRMAUkyC8vWfJfd6lifafX3XW2LbTnyUwGQfEAdkaPUVVGz+q3nbZeey6t61axef//znicVifPjDH+b555+nuLiY3t7eGe/r6upi8+bNMw5wpAIpWaqLxbMca/ZEPBP/psFgkKtXr96gSaPRsHv37uRUVK1Ws2fPHgBi8TgXL11C5h1luGcSh7CBHUX9nB/Vs616FIOxgng8Tr0qQFlZaoYEy0XKJVK73Z602Z6Omm1fRauvAJkcrb7i+ppgaTxFRUUMDw+jUqnYtWsXX//619m5c+cNv1tdXc3Fixex2+00NzenNqB5tPym89jtdhwOx83fuAAsFgs5OTlJ94vpMBqNbNu2jcbGRiorK3E4HBQWFuL1einIz2d3GZRqQ+SPvY21X8uRoZ309L3Gtcs/5dq1a/gDAZ451c6Pj19OSeNykLYSqSarlHV7vjfjZwuV6ObjKSkpobGxkaqqKsbGxti5cyc/+tGPZrznoYce4v3vfz9Hjx4V5VC6lEqbYvGIUT4WBIEDBw7Q2tq6YAHCYrHM+G+fz8eVK1emLGO8E1xrOcPme/0U6n1MxmzU1dfj8/mweYOsM6689YosnqJVcX9/P/DrUpdvgRNCC702m2e5WnI6axYAABc7SURBVPx+P/X1c9+SslhtqWhZiOdWiU17ezsmkym5ME7sEcwugybgdDp5+eWXKSsrIxwO881v/SOxxscwZGu4o7aAP3j/ngX/frqRsWZPkxap8YgZm6WWoQVBwGKx0NHRgUaj4QOPvZ9nW5x4wyZCkRi+4OreXZaxZk+TFqnxiBmbcDhMR0dH8mcdHR1ztl5Ph8ViIRAI0N/fz5UrV1DFiwgDJ60O+v79OH/9kRvXeiuFlKPi9XrnXCitBo+YNoxSGZNYPGLGpri4+KYf+tkQBIEtW7Zw3333AWB753ni0SjE4zgnV/eoZaZ3KE1apMYjZmwKCgqWPUUTBIHdu3fT2dmJN+wDrZ5QOPVz4angN7J3aCGOzPHK+TnSGZuxsTHePHUKlUqF9nqCRKNRqqurb3Cgbmpqwul0Mux0c2xcj065umuClKdDOTk5olmHL5enp6eHo0ePkpOTsyzn5qefflo0LfPx9PT0cOzYsXnf63Q6efbZZ2/Kk+BaaoPgcmNzMy0JHD9+PNmW8cILL3D48GHCoRCd0ywYExAEgfvvv59PfOiD/MH9a/nzD92dsq5UkHIS+Hy+OR+zwUCA1pYWms+cobWlheBN7D7m41kMEkmwXI7ZSZCKlvl4enp6OH78+Lzvdblc8yaBGHrSMabpMBgMDA0N4fP5uHjxIufPnwem/IXmqkolePZvWU91WdENr68kUpoO9fX14by+C5k4LxqLRpErFHi9XuLXe04CgQCtly+j02qRz+oNycrKoqKiYs5AHTt2jL/7u78jFotht9t55ZVX0Gq1PPnkkwwODlJWVsYPfvADvvvd73Lq1CmOHz/O9773PWCqC/KBBx7g8OHDAOzfv5/nnnuOF198kb/5m78hEonwpS99iVAoREtLC4cOHeLP/uzPOHjwIPF4nJGREZ544gnC4TAbN27k29/+9px6mpubZ/zsueeeQ61W8+STT9Ld3Y3JZOLnP/853/3udzl58iTNzc28/PLLHDp0iHA4TGFhIT/96U95+umnOXz4MPv27eNnP/sZzz//PN///vcB+OIXv8jGjRv50pe+xOuvv87mzZuX/G8Vj8dFKZHOxREMBunt7UWpVJKt11NaWnq9jBoH7xvEYltu6AkSQ4tYSPlJoFAo5mx6is9yF5v9/7Oh1+vn3e39xS9+wf3338+RI0d47rnnWL9+PSdOnKChoYFnnnmGT33qUzz55JP86le/SnJoNBoqKyuxWq1cvXoVi8WCQqHg61//OkePHuXYsWN87Wtf48EHH6SxsZGXX36ZgwcPJrVUVFRw+PBhTp48idvt5tq1a3Pqmf2zY8eOJTUePnyYxsbGpMYnnniCI0eOoFQqefHFFzlx4gTr1q3j6NGjfPzjH+fgwYMcO3YMmUzGCy+8wIkTJ3j++ef51re+xeTkJO+88w5vvPEGTU1Ni/vHWWR8U+WxWCyYzWa8Xi8tLS1Eo1FsIyPI/G8h8x/hpTfO0mezE4/Hk814YukRAyk9CSoqKubdFW293j+egFarpaq6et5dwfnsxzds2ABAWVnZ1GJqeJgtW7YAsG3bNs6ePZu0G59tP/7hD3+YH//4x0SjUZ544gnGx8dpa2tLHsIZHR2d8xvJ6XRis9n4i7/4C5xOJz09PQwNDc2pJycn54af2e12tmzZgtPppL6+HqvVOsMS3ev18qlPfYrBwUFsNht1dXWUl5cnX+/q6uLixYvcddddwNQ37ZUrV9i4cSMAW7du5dVXX53332W++KbLtj7hXdrb28sdd9xBZ2cnTqeT4MTzyOUy/q25H93JazyxOY7ZZGHz1ibRbOvFQMpPgkTJbDbq6uqSXvNarZa6urpl8Uw/ohePx6mpqeHs2bMANDc3U1tbi0qlIhqNUlBQMKP6sXfvXt544w3eeOMN9u7dS0FBAY2NjRw5coRjx45x8eJFZDLZDccACwoKeO2113jooYc4duwYu3btSibLbD0LaSwoKMBqtc7QCPDKK69QX1/P8ePHefTRR4nH4zNer67+/9s716goy36N/2ZgYBQGBhnkpCGGgJxfFQ0KOXhaeWhbHsggK5eyqm1pfnCt3toul6W1ltUX65V2B82Utdur0sJdWiCkeD6LpaIpBxlEgYFRBgZmePYHmomzA3NgeJnrkzzKNdfzf7zneZ77vu7rH0x8fDyFhYUUFhZy8OBBoqOjjfYEw/N2f9C1NgNFb9fJAA8PD5566ikSExP59eLj/Os7X/Q3i2gURvD1RVeqG+4at2AqFArK7qqprB3cvcZWi2Z3lUqJio5mSnw8UdHRuD6k+YKp8eNPP/00v//+O9OnT6e4uJhFixYRFRXF0aNHWbp0aaf4cbFYTExMDNHR0cYV03Xr1jFjxgxSU1NZu3YtAFOnTiU9PZ0jR44YtcTGxvLhhx+ycOHCfi9UGTSmpKRw7tw5o8YTJ06Qnp7OtGnT+PHHH5k/fz6lpaUA+Pr6UldXx+LFi3FycmLevHlMnz6d1NRUNm3ahEQiYfLkySQlJVFUVNQvPYZzsmZsfVBQ+9ZK+NuCPf7RELTaVsT15YhaGmlpc6L2wUjOnjlDUVERRUVFbNpzmI/2nTZblzkw20BnsPh2bNg8EJNYV56BatFoNJ0eLdavX8+SJUuIj483WZs5WvrisbfamPr5ptSmuLiYkpISQkNDjUa6o0ePkpubS1VVFa0iFxpD5pIU7EZKqBfN2vZNVvfbXAkeF8QjCveha6CzVLy2NSLMN2zYQFlZGfHx8TbXYm88loyJ7wmGGaKOdorIyEh+/vln9Ho94VETmTrNl8tVGqbExxvvJoZHK0tM3w4U/9bR7Js2bRo0LfbGY+3Y+pCQkE69yeBvi0RlZSWh4wLxdnNmUqCUgqIT+Hh54OkhQ6vte9utLeDwDllJi73xWLs2crnc2Kb1zJkzSCQSRo8ezbRp0wBQqVTcu3ePC6V1xE4cT1lVDdI7VZRcv8H2E/XEjfPm7eemm61vIHB4h6ykxd54bFWbP/74A5FIRGxMDJW3bxMyYQKpqam4urqiVCrxGV1Ji+CMpqEG9f0HPGhspk1wRdc2eItnjmh2K2mxNx5b1aa0tBRXFxdioqNpaWnhekkJYrGYyVOmGPscAMYZNzc3NzIFYVA71wyZaPadO3caLRF9cQw0mt1gUQBYvXq13USqW4rHVrH15eXlnL9wgV/z8nB1cUHUdBRn/dU+eSzRbNEcWC2a/Y6qkdXZeTy9eR+rs/O4o+p7rn2w48e//vpr45/ff/99u4lUtxSPLWLrtVoter0eiUSCTCajtvos4pZiWlQ/odZouVRWYxKPrWGWis9/ucT1yvYFGGen9u12+jY9TmInritVaHXtK6AVNfd5/dN8xvt54CTufCGC/TxZOTumx07mhmQ5rVaLl5cXc+bM4fXXX+fChQt4eHiwZ88eVCoVK1asQKFQcPXqVdasWcM333yDRqPh4MGDuLm5sWnTJg4dOoQgCHz11VfodDqWL1+Oq6srs2bNQqFQGBe3tm3bRlZWFgcOHODq1au8/PLLQPsC2Jo1a/pVH0t1Z7cEj1QqtcgUaV9aQkJCCA0Npbq6GqVSSeulEibN03FHmMvh3OMcvq7iy9d98PYYMSid63uD2UPRkGrm7NSZyjAAevu5KwzZlB296oZkubfeeousrCzu3btHY2Mjhw8fZvfu3WRnZ5Oeno5KpSIvL49PPvmEnJwc8vPz2bJlCwcPHmTChAlUVlZSWFjIuXPneO+994iPjycrK4sXX3wRQRAQiUTs3LnT6PfX6XSo1WrefPNNsrOzCQ8PH9B/oJ7OaSCwBI9araapqQl/f3+raZH/FbwbHBzMnVt/srD8Knd8v+RRZzHBYhcmhQYioRWdTmKx2lgCZg2ClbNjjC80hu12htXF1dl53K69jyCASARjvGVsWDq111XBnqbeDMly0G4aU6lUncxzBn9+REQEYrGY0NBQYmNjAQgICEClUnHlyhUKCwtJSUmhra2NwMBAli5dysaNG8nIyCAzM9O479UAiUSCQqGgpqaG8PBwYGCb3YfTFGlHuLu7MypwLLX/+Q66Zg0tIhHXbl2jVSJDpqtnxMiRjB8/3mwtloLVvENvpycwxluGWCRijLeMt9MT+s1jSJaDdtOYWCzuZp6Dv1+sGhoaOs0yCIJAWFgYs2fPprCwkAMHDrBr1y4kEgkfffQRO3bsYMOGDZ04oD23s6amBh8fH2OqwkDuBKb6oWzBY23vkAFBQUEEBQWxYMECWhV+lFZUsHv3bqpbpdS3jaRV34ZcLrdYbSwBq0Wz+3m58fHLnbfN9fVN1BOPIVluzpw5eHl5MXr0aMrKykhKSkImk5GTk9NppmLkyJHd8jZjY2Px8/MjJSUFQRDIyMjA09OTjz/+GI1GQ2ZmJtBug160aBGbN2/GyckJmUzGli1bWLVqFSKRaEDvBEM5mn2gWrRaLRUVFWi1Wu42aAiNjGPUqFFcrhVz+nYjL82ZRKCfwm6a9oEFDHR3794FMCZFD9Qk1pVnoFo0Gg3jxo3r9+d3/DtztPTFMxxq8/nnn3P27FkCAgKY/Hga6mY9E0aPxD8ggPLKaqZOjkUsFverNtaG2XeCoRw/bk0t9sZjq9pIJO19KpRKJWFaHVKJE7UPWpgcEIBcLje+W1mqd4MlYPYgGOrx49bSYm88tqpNZGQkgljCpT+uoblXwZlmb2J8JXyd8z8cqXJmevhoYh/xormpCemIEb1au20JsweBwV9u6E8wmDx1dXU0NTXh4+Mz6FrsjceWtdl/uRbBM45nwwKYKAjo9Hr+fADNOhEtuvanb01T06A9/nSFw0VqJS32xmOr2oSEhOBHLmN9fBGJRAiCgEwm4x8xkYwOuI9W6k1UVCi376pwl0rM1mMJOFykVtJibzy2qo1cLmdq5KOUlV2jYZwXI0aM4Nq1azQ1NXGkVMNdsS9tdWXsPP+AAE8p21+ba5YeS8DsQWBta0BhYSF5eXm8++67JnFY2xowVHlsWZuIiAgiIiJwdXXl5MmT3L59G5lMRmODFrx82XXhAdK2Zka7ijl59gJiEQi6lm75pnqdDpmHh3E9yFowexD0NvfcWlVB9X+9QmtFKZKx4/B9Zzt49v4tZIk57N4ykAbCYwnYE4+ta9O1r9yNGzdoErfHvfhLdUyU1nP7z4uU+HmCSESYb/f3A0EQaLaBxdqsQVD7ry2o/2hf0dX+1bhB39ZGg1iM9loxgrZ9O19r2Z9UrnoK5wmRNHSxH7g8Go73q/80+st78qsXFxezYMEC6urqWLZsGVKplJUrV7Jx40ZSUlKA9t7GGo2G2tpa1q1bx5dffolUKmX//v1UVVWRkZGBVqslLi6uxzS5vXv3Gl/U+tLSH9gTT2Njo0U8+6ZoMbRrKisrw93dnVmzZrVnJ125ilL1JzGBY2htlVGi0SCTOtPT/UmtVvPgwQObNPs22zYhkUiMc8MdYRgAvf3cFQ/LszEkvBmyebp9niDwww8/8OSTT3Lq1Cny8/MJDAzk/PnzKBQKfv31V/Ly8npNk+sYlvswLabCnnhslTsEf0euBAUFoVQqcXV1JS0tjWeeXkj4I+2LY83Nzfj6+vK/52v48VItV65coaSkhJKSEs6fP09ubi4eMhkaC/R3eBjMuhN4v/pPYwKdf5cEutsr59NafguENhCJkTwSzKjN/93rtFhfL1wdE95EIpFxwabjYndUVBQ1NTW4u7sbEw8MJrra2lpeeeUV6urqKC8v7zVNzhQt/YE98QzGpEHXzfcdeyOXl5fj7e2NU00L1VVKypseMC4oiLCwMBoaGrh586bN8krNfifozQrr+872bu8Efa019mWp7WhuU6lU3L/fnlhWXFxsjCoUiUR4enoyYsSIbolwOTk5LFy4kKVLl7Jq1ao+0+QepqU/sCceT0/PHu/Y1tRi2HzfEdrmZvx8fY0NBNPa2jh37hyCIHBfrUYsFuPt7c3atWu5cf06Y2zwOGT2IDDMPXfNvpf4j2XM5/s7HWvtY566N56ucHJy4sCBA5w+fbrbziSNRoNWq+22kJOWlsby5cv57rvvTHqxM1XLUOKxZDR7f7QYclQrKyu5c+cOE8PDEYvFnPsrSrJW8EChUDDjHyFUVFRQr1Lh5+/f/kjVz7yogcJsA11VVRWAcbPGQE1iXXkGqkWj0fQ6pWaqNnO09MUzHGtTX1/Pt99+y8WLFxGPGofeO4RFsaP4v9wfUKvVpGe8CCJImBI3dA10lorXtlRsuCU6NNrTOVmKZ7BqI5fLWbx4MS4uLpy6WYeqsYXiy7/j4uKCQqEg51wtUomYhClmSxswzB4ElorYtgSPNePHhzrPYNbG8ELsc/w4paWlVFdXIxaLUSqVzEj0QCwe3LQJh3fISlrsjWewayOXy0lOTiYhIYEbN26wd+9elEolR27eR68X+I8Us6UNGA7vkJW02BuPvdTGsIYgEolQKpXEBY7EeajfCcy1CVuSx83NzSJBTvZ0TpbisbfaBAcHExoaSotOoIXB7V9m9puSJSO/zeWxJy32xmNPWgw8c+bMoe2vNJLBxL91NPtgarE3HnusjY+PDxljxw56DKPZ6wQOODDUYf7EsQMODHE4BoEDwx6OQeDAsIdjEDgw7GGTQaBUKpk0aRJSqRSdTodOp+PZZ58lNTWV9evXA+1TZpmZmSQnJzNz5kzjoszWrVt54oknyMjIsFiYlSl6AF577TVSUlJYsWKFcTPPnj17SExMZP78+cZkZXNw8uRJEhMTSUpK4o033gB6PmdTj9lCT2trKwkJCbi7u3faL2BpPaZouXXrFklJSUyfPp3nnntuQNfJJoNg1KhR5Ofn89hjjwGwd+9eYmNjKSgooKmpiYsXL3LhwgVcXFz47bffeOmll9izZw/37t2joKCAoqIiYmJi2Ldvn830nD59mpaWFgoLC4mMjGT//v20traSnZ3N4cOHef755/n000/N1hIUFMShQ4c4cuQId+/e5ciRI93Ouac6WKs2puhxdnZm3759LF682Ph71tBjiha5XE5ubi6HDx8mODiYn376qd/XySaDQCqV4uXlZfz55s2bxMTEABAXF8fx48eNu8ag3aTl7e3NqVOnjHuIZ86cyYkTJ2ymp6djJSUlREdH4+zsbDE9fn5+xgQHZ2dnLl261O2ce6qDtWpjih6RSNTNiGcNPaZo8fLyMqbiGbrx9Pc6Dco7QVhYmLG3QEFBASqVCoVCgVarZeLEiWzfvp1nnnmG+vp64+YNT09PVCqVzfR0PHbo0CFUKpVV9Vy6dImamhrkcnm3z+jpc61dm7709ITBqo0BSqWSvLw8Zs+e3W8tgzIIFixYQFNTEzNmzMDV1RVfX19++eUXPD09uXLlChs3buSDDz5ALpcbn+fUarXFsj1N0RMXF0dUVBSpqamo1Wp8fX2tpqeuro7Vq1fzxRdf9PgZph6zFB6mpycMVm2gfffaCy+8wGeffYazs3O/tQzKIHBycmLbtm3k5+fj5OTE7NmzEQTBGOOhUChoaGggPj7e+G2cl5dnfIa3hR6ADRs2UFBQgLe3N/PmzSM0NJTLly+j1+stpken05GZmcnWrVvx8/Pr8ZxNPWYJmKKnJ1hDj6lasrKyePXVV4mIiADo/3USbICWlhZhxowZglwuF9LS0oTCwkIhOTlZSE1NFXbs2CEIgiC0trYKixYtEpKTk4WkpCThxo0bgiAIwvvvvy88/vjjwrJlywStVmszPXq9XkhOThbS0tKEzZs3G393165dQkJCgjB37lyhvr7ebC05OTmCQqEQkpOTheTkZOHYsWM9nrOpx2ylZ8mSJYK/v7+QmJgo7Nu3zyp6TNFy7Ngxwd3d3fhvvv/+e0EQ+nedHN4hB4Y9HItlDgx7OAaBA8MejkHgwLCHYxA4MOzhGAQODHs4BoEDwx6OQeDAsMf/A+OMYNytyq8GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118b12590>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots(figsize=(3,5),facecolor='w')\n",
    "\n",
    "# this sets the grey vertical shading \n",
    "every=12\n",
    "xDates=['%04d-01-01'%(y) for y in range(1975,2018)]\n",
    "\n",
    "[ax.axvspan(bt.decimalDate(xDates[x]),bt.decimalDate(xDates[x])+1/float(12),facecolor='k',edgecolor='none',alpha=0.04) for x in range(0,len(xDates))]\n",
    "#ax.set_xticks([bt.decimalDate(x)+1/24.0 for x in xDates if (int(x.split('-')[1])-1)%every==0])\n",
    "\n",
    "#ax.set_xticklabels([convertDate(x,'%Y-%m-%d','%Y') if x.split('-')[1]=='01' else convertDate(x,'%Y-%m-%d','%b') for x in xDates if (int(x.split('-')[1])-1)%every==0])\n",
    "ax.tick_params(axis='x',labelsize=8,size=0)\n",
    "\n",
    "# this sets the vertical dashed lines on the tree\n",
    "[ax.axvline(i,ls='--',lw=0.45,color='grey',zorder=0) for i in range(1975,2018,5)]\n",
    "\n",
    "\n",
    "branchWidth=0.35 ## default branch width\n",
    "\n",
    "cmap=mpl.cm.viridis\n",
    "\n",
    "# k objects are tips, nodes, branches\n",
    "for k in tree.Objects: ## iterate over objects in tree\n",
    "    #   x=k.x ## or from x position determined earlier\n",
    "    x=k.absoluteTime ## or use absolute time instead\n",
    "    y=k.y ## get y position from .drawTree that was run earlier, but could be anything else\n",
    "    \n",
    "#     xp=k.parent.x ## get x position of current object's parent\n",
    "    xp=k.parent.absoluteTime ## get x position of current object's parent\n",
    "    if x==None: ## matplotlib won't plot Nones, like root\n",
    "        x=0.0\n",
    "    if xp==None:\n",
    "        xp=x\n",
    "        #c='indianred' ## colour can be fixed\n",
    "#     c=cmap(k.height/ll.treeHeight) ## or be a function of something else\n",
    "#     c=[cmap(k.traits['posterior']) if k.traits.has_key('posterior') else cmap(1.0)][0]\n",
    "    #c=colors[k.traits['geo']]\n",
    "    #c=['#63C78A' if k.traits['geo']=='cambodia' else '#6A34CA'][0] ## can be discrete too\n",
    "    \n",
    "    if isinstance(k,bt.leaf) or k.branchType=='leaf': ## if leaf...\n",
    "        #x=decimalDate(k.name.split('_')[-1],variable=True) ## get x position from name\n",
    "        s=sizes[k.name]\n",
    "        #s=sizes[k.name] ## tip size specified in a dictionary\n",
    "        #s=50-40*k.height/tree.treeHeight\n",
    "        \n",
    "        \n",
    "        #ax.scatter(x,y,s=s+0.8*s,facecolor='k',edgecolor='none',zorder=10) ## plot black circle underneath\n",
    "    \n",
    "    if isinstance(k, bt.leaf):\n",
    "        if k.name in colors:\n",
    "            c = colors[k.name]\n",
    "        else: \n",
    "            c = '#CBCACA'\n",
    "            print(k.name)\n",
    "        a = alphas[k.name]\n",
    "        ax.scatter(x,y,s=s,facecolor=c,edgecolor='none',zorder=11) ## plot circle for every tip\n",
    "    \n",
    "    if k.branchType =='leaf':\n",
    "        c=\"#000000\"\n",
    "        #c = colors[k.name]\n",
    "        \n",
    "    elif isinstance(k,bt.node) or k.branchType=='node': ## if node...\n",
    "        c=\"#000000\"\n",
    "        ax.plot([x,x],[k.children[-1].y,k.children[0].y],lw=branchWidth,color=c,ls='-',zorder=9)\n",
    "        \n",
    "        \n",
    "    ax.plot([xp,x],[y,y],lw=branchWidth,color=c,ls='-',zorder=9)\n",
    "        \n",
    "    # add in a legend\n",
    "    han_list = []\n",
    "\n",
    "    for key in color_key:\n",
    "        marker = mlines.Line2D(range(1), range(1), color = color_key[key], marker='o', markerfacecolor = color_key[key], label = key, markersize = 4)\n",
    "        han_list.append(marker)\n",
    "\n",
    "    ax.legend(handles = han_list, markerfirst = True, frameon=False, loc=3, fontsize=8)\n",
    "\n",
    "        \n",
    "# set axis limits, remove border lines         \n",
    "ax.spines['left'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['bottom'].set_visible(False)\n",
    "\n",
    "ax.set_ylim(-5,tree.ySpan+5)\n",
    "ax.tick_params(axis='y',labelsize=0,size=0)\n",
    "#ax.set_xticklabels([])\n",
    "ax.set_yticklabels([])\n",
    "\n",
    "#ax.set_xlim(bt.decimalDate('2008-01-01'),bt.decimalDate('2016-12-30'))\n",
    "ax.set_ylim(-5,tree.ySpan+5)\n",
    "\n",
    "\n",
    "plt.savefig('/Users/lmoncla/Documents/Conferences/LSRF-annual-meeting-2019/poster/HA-annotated-tree-2019-10-01.pdf', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot tree highlighting the Cambodia 2013 samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a dictionary for color and size\n",
    "color = {}\n",
    "size = {}\n",
    "\n",
    "for k in tree.Objects: ## iterate over objects in tree\n",
    "    #print(k)\n",
    "    if isinstance(k, bt.leaf) == True:\n",
    "        if \"cambodia\" in k.name and (\"2012\" in k.name or \"2013\" in k.name or \"2014\" in k.name):\n",
    "            size[k.name] = 2\n",
    "            color[k.name] = colors[k.name]\n",
    "        else:\n",
    "            color[k.name] = '#CBCACA'\n",
    "            size[k.name] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to change the colors \n",
    "for c in color: \n",
    "    if color[c] == '#13257a':\n",
    "        color[c] = '#368FE3'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# color dictionary for key \n",
    "color_key2 = {\"human\" : '#E3363B', \"poultry\": '#368FE3'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIIAAACJCAYAAADgxOA+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHR1JREFUeJztfWuMHGe55lNd1VXV1V19n5ked3rGzoyxExPHIRexBMd2EkEUKZCFg2RjBy0Ru8CG3YWDxArBQdYKrVjFOYjD8Un+JCCBfY5WWBBgIUCc2E5E7hDf1gp2PB6P7enp+6Wq69JVXefH+Psylx67pqdt99j1/LFn6plv3p5666v3/jGO4zjwcMPDd60F8NAb8BTBAwBPETxchKcIHgB4iuDhIjxF8ADAUwQPF+EpggcAANfJDxUKBQBAMplse73RaECSJBiGgVqthnq9DgCIxWJt+ZqmIRAI0K/L5TLlz73WDm755NpM/qWwWH47eQguXLiAVquFcDhMv6frOiqVCprNJprNJliWRaPRQCgUQiQSmcUlqNVqMAwDkUgEDMPA7/dDkiTEYjFwXEe3E0CHiuDzudtI6vU6Go2Ga/5i118M3zRNVKtVNJtNNBoNAADLsgvym80mqtUqRFGEqqoAphXDDdrxa7UaHMeBqqowTRO2bcO2bZimCY7jwDAMdF2HruswTROlUgk+nw/ZVh9WJCSsHQyCZVkwDAOfz4dGowHDMJBMJhEOh5ekBECHikA+aDwevySPZVmEw2EoigJg4SdKEARIkkS/nsmfe60d3PALhQJ8Ph/C4TBVhHZPHIFpmtB1Ha1WC81mEwBgWRYMwwDDMPP4tm3Tz1ytVinXtm00m02YpgnDMOhaoVAI4XAYkUgEmqaBZVmEQiH6OxuNBmzbxp73UhgsMvinNRLdMRiGQSqVgqIoCIfDl1Rot+hIERZ6JfQyn2VZyLKMWCwGURQB4LKvHNM00Wq16Po+nw+KorR9+iqVCgzDoEpomiZ4noeqqmi1WrN2LY7jwLIsms0mHMeBJE3fZI7jEAqF0Gw2IQgCgsEg/nEIWJnuQ0yellUURTQaDSQSCQwMDLj7A7nAkmyETCZzXfPJU0p2nEwmQ//4xP4hSiGKIlUSXdchiiKSySRM04RlWdA0DY7jwDAM+P1+VCoVFItF8DyPer0O27apEjWbTSiKgmg0ilWrVkHiP1CiQqEATdOQSCRcfVa36EgRotHoDcsnClCtVqFpGsLhMEzTpIYbAASDQTSbTUxNTaHRaNAtndzwvr4+iKKIRCJBubZt053HcRwUCgX09/cjk8lAEIRZsvA8vyj53aAjRSBPiCzLl+TZtg1VVWEYBoDpd2w7aJpGOQCQz+cpf+61dnDDr9VqdHt2Kz/BTD4xgEOhEGRZBs/zKBaL0DSNKoKmaVAUBTMz/MRWMAwD+XwePM/jtttuQyKRwNTUFAKBAIaHhymfeF7tZGk0Gl19LQAdKoJbyLJMt8PrBcQAnmn4SpI06zVhGAYEQUAikaBPbzabxdmzZ1EsFiGKIhiGgWEYmJycRDabxdDQ0DX5PAQdKUIwGHTFEwQBfX19dGtbyEqfa+kTKzgcDrvyGtzyiTvnVn6Cy/HJ52wnDwHDMKjX6xAEAZlMBoqiQNd1MAwDy7IwPj6O48ePQxRFrFmzBjzPg+f5eYZpMBhs67UsFR0pQrVaBXBp98vjz+bzPI9Vq1YBAPr6+mYZm36/H2fOnEG9XoeiKGi1WrjpppuQSCTmKQKxTVKplCtZ3OKGcR9t20a9Xke5XKZ/XBIBXAj1eh08z3dNHoZhUKlUUCqVkM/nMTw8jMHBQYyOjiKVSmFychK5XA5+vx/BYLBtfCCZTNI4SDdxw7iPsixD13UA0z4/ANdPVbfkkWUZxWIRY2NjUBQFoVAI5XIZhUIB4XAYwWAQfX19SKVSVHHJzxHl7Sn3sdvv2KvBFwQB8XgckiRRa95NroFl2a7JIwgCVq5cCY7jUKvVkEgkMD4+TsPP5XIZjUYDuVwODMNg7dq189bqKRuh1Wp5/A75giBgZGSEfp1IJDA2NoZsNotqtQrDMMCWy/jQb34Hfvs2xD7zH+etvVh53KAjRXCbfPH4l4ckSQiHw9A0DZZlgWEYMKoK6CbOTuWQP358Orp40RPSNI2+4rqJG8ZY7HV+OBzGwMAAjh49iqos4/B/2o6YFERtagqxWIwqwpUyFjsqTCkUCtQg8vjd4ddqNWQyGQwNDSEWi2H4zDhG/+lf8NLvJvGrp1/CVL6In7xRQT5fQLFYdL22W3S0I1wu/ezxO+fffPPNYBgG5zQd1WAQoXeKyNf8eONvCp55i8U9fxdGf0xc1Ppu0JEiEPfLrTW9XPm2baNWq9F3/sDAAGq12mV/binySJKE0dFR2LaNcY7F8/o63P/+Afh/tAu+j/1PFEpVCJI5K5LZDXSkCCS5cr3zZVmmZWSdrk8iiCThZNs2HMehcQFZlufJIwgChoaGoCgK/u53/w93jL+DE31rcGcKGI6xgN39FFFHKy42Dbpc+TzPIxqN0htF0tGmaSKfzy9YHqZpGprNJg0ZG4YBy7JooQvJhzSbTRpT8Pl8qNVqNKooCAI2bNiAPtPC/gNh/Ku8Ec/cl0BUMnonskgiXm7z+tcbv16vg2GYebaAaZpQFAWVSgWKoiAQCNDKI0VRYFkWWJaFoij0ppumCVVVIQgCWJadlV6WJAmrN9+H1Zvvw1cufm9iorBod9YNPPexA76maW0LZlVVBcMwWLFiBS03q9frNOtI6h9ZlqUV0pIkged5ZDIZcByHm266CYIgLPjUe+5jD/EDgUDbekdSq0ACQ5Ik0WIcTdNoXsE0TUSjUcRiMcTjcQiCAI7jMDIyMqsaaSFZesZ9dFvZc73y/X7/vOIU4IOKLILTp09DVVVayUTsgnK5DNu2EYvFEAwGwXGc69eULMtdqVqei44UwTRNj98GpCKL1COEQiEEAgEIgkDL+k+dOkWTTOTJNk0TpmkiEolctj+BcLuNjhRhse7UjcInlUrkRgmCgKmpKZr1TKVSEEURuVyONrOMj4+j0WggmUziQx/6kCtZFqr9XAq8KuYryOc4ju4OwAeZx5tuugm5XA6nT59GNpuFYRiuq596qoq5VCoBcB85u1H5mUwGoVAIU1NTcBwHrVYLDMPQpzoejyMcDmNsbAyVSgXHjh2D3+9HIBCgu0i7tTVN643IYq+7d73El2UZiqK09f0ty4IsyxgaGsLY2BhKpRKGh4cv2YHllaotU74kSRgZGWnbp0DCz2NjY7AsC5ZlwbZtDA0NLWgQ9lSp2uV6Bj2+Oz4xLlesWEG9jWq1ijfeeAPDw8NQVXVeu3sgEMCVGI15Rdvir0e+bdtQFOWS9sJi108kEtSuaLVaqNVqOH78OIDpAttQKETti3w+33X7ALjCbfHXI1+W5VmeANneW60WVFUFy7L01VCv1+k8BsMwZoWWBUGgJWiqqiIYDCKZTKLZbEIURdrGR5Sq0WigXC6jVqtB1/V5PZFLhWcsLpJPqqFJBpIkpIDpMDLpiZyLRqOBVqtFt3XbtmkrPclOkqQUy7JIJpMQRRGGYdAKJjJxpWfcx142zq4Wnzz9pNydZVlaW3jy5EmadzAMA7IsIxKJ0I7oueXolmWB4zhIkkTtARI0UlUVpVIJlmWB53mwLItoNNrV3QDoUBEikcgNzyc3emZ+we/3I5fLoVwu017HRqMBy7Lo2BzDMCCKIniex/DwMGKxWNteSeJlpFIpDAwMUKUDFu4qXwqWZCO4jYZdz3ySX2g0GnQWwsxBF2TCiWmaYFkWfr+fzks6d+4cfD4fSqUSLWufC0EQsGLFCvr15ORk78QRPHwA4gIWCgVMTU1hYGAAPp8P9XodqVSKzk+oVqtIJpOoVCqoVCqo1Wo4f/48jRRmMhmkUqklD8XqFB391lAo5PHngHgT5MknsxFIeZpt2ygWi9B1HeFwmIaZHceB4ziu3/mhUGjR7qkbLKmK2W1e/0bgE28in89jamoKkiTRHYCM2tF1ndYhiKKI0dFRBAIBjI+Pt52FsJAsmqb1xsSUXnbvriXftm2Uy2UEg0EMDg5CFEUIggDLspDP5+nuoOs6QqEQEokEZFl2basQWbxStR7ny7KMQCAA27ZhWRZisRhSqRTS6TR4nofP54PP50M8HkcoFKJewGLk6alStV5oc7+W/IXCzIIgYHBwcF4TTCwWQyQSgaqq1C6YGUtYjDxeW3wP8eeGmedySTk7iToahoFWq4VAIADLsuj4m1AohGg0Sm2QmTd45nS4mYMyvLb4HuITw7AdLMtCsViEz+ejN69UKqFWq9Hpr+RGlkol2LaNiYkJ2LaNgYEBGqVcyDPw2uKXCT+ZTGL16tX0pgOgTS0sy0IURTqIG5ie46TrOnieh6IoEEWRTodr56V4xuIy4ZMeBYZhEIvFEIvFsGrVKlqhzHEcRFGkc5d9Ph8d1tnf349oNIpoNLpgyXpPGYu93Fbei3xBENDf34+pqSnYtg1d12EYBi1hD4fDWLVq1ayk1kIVSvF4vHd2BBIm9fiX55ORuxzHUWMxGAzSlvt6vU7tCjfJpEqlsqChuhR0tCMsNh5+I/NJBpLneTiOA8uykE6nwXEcNE2jPZFTU1NQVRXDw8O0+2mhtXum02mxufAbmU/S1YIgUFeR5BeKxSL1EAzDoK1wZPzuzB2EuJC6rl+RxJTXFt8hf26f40xks1kA0+P2Sfs7STjV63UMDQ1hYGAAQ0NDUFUVJ0+epCfDOI5DPYyZMQMyXKNQKIBhGFoJ1S2l8NzHDvmkDqEd5ipYrVbDxMQEKpUKbNuGpmngOA7r169HX18feJ5HLpeDoiiQZZmGotsVxPA83ztxhOVQSnal+XMnss/ExMQEgOnQsm3byGaz8Pl81Ej0+XxQVRW5XA6O49AdhFQzk3T23XffPU+pJiYm0Gw2XZ845xZeW/wV5pNEFMdxtCyNRBPL5TIdu0vcwnw+j0ajMavKae7aPWMs9kob+nLgk0QUMQxJypl4CuVyGaqqolarIRAI0Mhhq9Vq+2rw2uKXMZ9EGmeC5Bw0TYOqqlBVFY1GA8FgEOl0mtYqtFu7Z4pXe7UNvVf5xMOYOUaPVDEXi0VqdJZKJUxNTV2yytpri1/GfFmW27a4k9Az2UECgQBUVcWJEyegqio2bNgwL4bhtcUvYz4ZnUNOhFUUhT715DBSclRwMBiE4zgLpqG97ON1wldVdVYcwDAMVKtVFAoFqKqKZDKJhx56COvWrWsb0eyp7GOvtJUvJz4pbwsEArPqFSVJQn9/PwRBoHOfX3/9dfj9fmzYsGHeK8Vri1/mfFmW27p9sVgM2WwWyWQSsiwjl8uhXq/D5/NhaGhoniKQqGO30dGKxN3x+O75PM8jHo/PCwaRSiSe5zE0NESjj47j0BPo567dMy1vvW6cLTc+6W0MBAIIBoO0TuHs2bOwbRvpdJraC94MpeuIb9s28vk8zSQWCgW6K5DJ7bquL2gs9swMpeXQtt7L/GAwSCubRVGEz+ejFUrkCOXR0VHIsgzDMDAxMTErGLXYvgw36EgRyNbktlXL488Gz/O4+eab6ZHAhUIBExMT8Pl84HkekUiEjslpNptgGIYamoFAAOl0Gul02pUsbtGRIizWffH48zEzja3rOqLRKILBIOLxOCKRCD0Q1HEcml+wbZuO3Dl58iTS6fRlD1B3C68t/iryF2qVC4VCkCQJoihS5SADNkqlEiqVCgKBAHRdh2ma0HUdp06dwrlz52BZFlqtFu66664l2Q1eW/xV5C/UKlepVGBZ1qynmxwjlEgk6CsEmC6Dm5ycpAUt5IAQRVGuviL0mvu1XPgLtcolk0lan3D06FEMDAyA4zjE43Fa0EIUQRRFqlBnzpyBpml0WNdS4LmPPcJ///33aU2jz+ejhqZpmujv7wfP85BlmbqP69atw+joKC14Wezray68tvge4ZNIIsuydCgX2RFmhpRntsXPHbS1FHht8T3Cj0aj0HUdiqLQ9367KiWvLf4655PzHsmxgIqitE1SeW3x1zk/HA4jGAwin8/T+c1kzgKBLMteruF64ZMDPzVNQzabpbtANBrFyMgIJEnC22+/jWazScf8hsNh6jX0VK6h19vQe5mvqirq9ToKhQKy2Sy19nmex9/+9jcUCgWk02k6c9nv98O2bXAch2KxCEVRwHEcyuVyV1vertu2+DNnzmDHjh0AgOeee+6ay0NA6hFI2TpxDcm8JMMwoGkaotEowuEwbXqp1Wo4e/Ysjh07hmPHjiGXy7WtV+gUV6wt/nzVwjefL2C8bCEtA/+wiYPbNEm329afe+45PP7447P4i7G8uy0PGZUjSRI9P9rv92PlypUYHh6mySYAtNmVtNVPTEzANE3UajWMj49jeHi4K5PaO1IEUnn7jwcqOJmfb9narRbey1vQL/ZhTNSA//Y7C+tSuXnc1X08vnLP7Dp9sj4AHDp0CLt374ZlWTAMA7/4xS+wc+dOvPvuuwiHw9izZw/eeustHDx4EE8++SR+9rOfQRAEbN68GQDw61//GkePHsXmzZvxzW9+E0899RRuv/12ZLNZ6LqOn/70p4jFYvj617+Oxx57DHfeeeeCn3exf592kCQJ0WiU3kjiGUiShPHxcdx5553o6+trewZULBaD3+/HhQsXut772NGroVaroVarXZKjz2nGMWz3Gby56+u6jt///vf48pe/jN27d0NVVRw6dAhbt27FM888Q09IaYdPfepTuO2223DgwAE88MADAIAtW7Zg165d+NznPod9+/ah1Wrh8OHDbZXA7ed1y2cYBpZlIZVKYXh4mHY3V6tV5PP5y05DIa1wo6OjGB0d7dq5DUtyH/8+0766t9Fo4PFf1HCmbMFxAIYBhqIsnv5c/4L8dusT3HHHHQCADRs24Dvf+Q6+8Y1vAADuuusuHDx4EHfccQetHHaTAt6yZQt4nsejjz6K7du3Y/Xq1bjvvvsu+3nd4lJ8kicIBAL0XMiTJ09C13UwDHPZLqZl19fw1KeTWBnj4GOAdMjBt+5xb9jMXf/w4cP038ceewzvvPMOAODtt9/GyMgIWq0WxsfHAYAeijUTc/sOS6USCoUCQqEQwuEwfvSjH2Hbtm2u5Vms/DNBEk+kXI1kGX0+nyvjr6f6GtykY9MRDv/2hRQALMribre+3+/HQw89BF3XsW/fPnzve9/Dxo0bIcsy9u7di1arhV27duHhhx9um4W755578Oijj+KJJ56g6xODbuvWrfjud7+LtWvXupZnsfLPhWmauHDhAh3QTY7o8fl8OHHiBN5//31YlgVRFBGJROD3++l4PnIuVLlcnvdZloJl0Ra/YcMGfP/736df7969e9b1XC6HvXv3or+/f5aR9fOf/xwA8NRTTwGYfgUdOHAAuVyO/g6GYfD5z3/+iso/F2QecyAQoMO2fD4fMpkMVFVFuVyGZVkQBAGFQoG2wPn9fgiCAL/fT5NTo6Oj104Rer0N3S1/3759+OEPf4jnn3/+qsqTTCbp+QyiKOLChQuIx+O49dZbMTU1RUfokZPfyGxnMmSDKAS55vash0uBcTronyLNGwulW+e6Pkvht3OjLiXPpfjk2uXk6ZTfTp65a82FYRg4deoUfD4fUqkUnezOcRyazSY4jkO9XsfY2BgURaHZyXXr1iEcDiMSiSAajS5ZEby2+GvMFwQBoijOy1gSIxKYjh8kk0kUi0UcP34c+Xwe586dw7p16xbt0SwEL/vYA/yZnU6maaJUKuH8+fPQdR0+nw9r166lMxZEUcRf/vIXaJp27UPMyynbt9z4qqqiUqmgVCpBURTYtj2rQIUcB5xIJBY9yeVSWFKI+WrwDx06hFdffXWW13At5bnSfEmS0Gw2wfM8QqEQ6vU6Tp06hb8c+f8oBkawMTM9tZW4nN1CRwGlmaeSLgTr/Hlkt+7Auf9wH8z/+t/B5PJdXf965cdiMciyjMHBQWzZsgW33HILVFXFeNnG/z0ZQEU1eme8HrGKfT/9GcyTJ+ddb9ktWO+9B+jT8f/W2QloX/kach++dR6XX70a/Ff+c9v1SV7/6NGjeOSRR1AqlbBt2zaIoogvfelL2LlzJzZv3oxcLodnn32WnrL+xBNP4LnnnoMoivjtb3+LyclJbN++HYZhYMOGDfj2t7+N1157DXv27EGr1UKpVMIf/vCHBSuBr+Zp9GQcX7lcxsTEBGKxGBKJBFqtPHaOTkCvN1G2dPj9/llnSS8VHe0IyWTy8gaR/kESiAHALCIo02793/zmN3j44YfbGkgk+vb888/jk5/8JN58803s378f6XQaf/3rX5FMJvGnP/0JL774Imq1GiqVCo1AknX379+/tM/bRT7HcTh//jyOHTuGI0eOwO/3Y9WqVUgkEgiHw121Dejv7OSHqHHz9/+j7fVGo4Ha4/8F1plxwHHgMAyQXoH+p/95QX7b9S8aWx/+8IcBAOl0GgzD0IANCYFUq1WsWrUKAJBKpagVvmLFCpTLZRSLRXz1q19FqVTC2bNn8elPfxqO48xa91Jh8KttXJLWOFLNRMbx9vf3Y3JyEpqmYWRkxNXabtHRjhCJRC7bWZN86v+AWzkM+HxghzKQ//f/6nj9mUmjcrmMyclJANOvDGDaX59pcM3kO46DvXv34tFHH8ULL7yAe++9F5IkzTs271JxNTeft5t8WZbR19eHeDxOFf/48eN47733EAgEFnVgqFtcsbZ4Lp1G6t+mY/3ZbBaLKcC+1Posy+KFF17AW2+9Rd+PhmFcMqx7//334wtf+AL27dsHlmUXrF3oRJ4rwRcEgaabyRSVVquFer0OWZavyOTVjkLM5IkcHBxse31uOHUpfDchZrd8cu1y8nTKbyfP3LXc8g8fPoyxsTE0Gg1kMhkkk0lEIhEoigJd13Hrrbd29QAPry2+R/nk/EeWZdFsNpHL5ZDNZhGPxxcdp3ADry2+R/nJZBLDw8M4d+4cPWuaRBoZhsHk5CQ9iByYTn0zDINMJtNR+ZqXa+hRPs/zWLNmDf1/q9VCNpsFx3EwDIMO+SYFLeSwj1gs1pEieCN4e5hPhnaTs55kWYbf74ff70cqlUIymQTLsuA4DoODg4jH4x1Pbvfa4nucT44DOn36NFRVBc/zqNVqKJVKyGQyWLlyJYDFn0Q3F15bfI/zBUFAOp1GLBZDsVjE5OQkSqUSDMNAvV537clcDl5b/DLhS5JE3U+O4+j8pG6hoziCh+sP3R/z7WFZwlMEDwA8RfBwEV1RhAsXLuAjH/kIRFGktfhbt27Fli1b8K1vfQvAtGW8Y8cObNq0CQ8++CD1m5988kl8/OMfx/bt27tq/LjFG2+8gY997GPYuHEj7alsJ5Pb7y1XdEUR4vE49u/fj49+9KMAgF/+8pe4/fbb8fLLL0PTNBw+fBjvvvsueJ7HwYMH8cUvfhF79uxBPp/Hyy+/jFdffRXr16/Hr371q26IsygMDw/jpZdewiuvvIJcLodXXnllnkzt5OwF2buJriiCKIqz+vVPnz6N9evXA5huV3vttddoUQkwHVtPJBJ488036RyDBx98EK+//no3xFkUUqkUjddzHIcjR47Mk6mdnL0gezdxRWyENWvW4ODBgwCAl19+GeVyGclkEoZh4JZbbsHTTz+Nz3zmM6hUKjQHT0bTXyscOXIEhUKBjqyZKVM7OXtJ9m7giijCI488Ak3T8MADD0AQBAwMDOCPf/wjIpEITpw4gZ07d2LXrl10cggwPVziStTiuUGpVMLXvvY1PPvss21lcvu95Ywroggsy+LHP/4x9u/fD5Zl8YlPfAKO49Cq3mQyiWq1irvvvpvuHC+++CK1Ma4mLMvCjh078OSTTyKVSrWVye33ljWcLsA0TeeBBx5wotGoc//99zsHDhxwNm3a5GzZssX5yU9+4jiO4zSbTeezn/2ss2nTJmfjxo3OqVOnHMdxnB/84AfOvffe62zbts0xDKMb4iwKe/fudZLJpLNp0yZn06ZNzp///Oe2Mrn93nKFF2L2AMALKHm4CE8RPADwFMHDRXiK4AGApwgeLsJTBA8APEXwcBGeIngAAPw7Mfksim3FtosAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1175a4fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots(figsize=(2,2),facecolor='w')\n",
    "\n",
    "# this sets the grey vertical shading \n",
    "every=12\n",
    "xDates=['%04d-01-01'%(y) for y in range(1975,2018)]\n",
    "\n",
    "[ax.axvspan(bt.decimalDate(xDates[x]),bt.decimalDate(xDates[x])+1/float(12),facecolor='k',edgecolor='none',alpha=0.04) for x in range(0,len(xDates))]\n",
    "#ax.set_xticks([bt.decimalDate(x)+1/24.0 for x in xDates if (int(x.split('-')[1])-1)%every==0])\n",
    "\n",
    "#ax.set_xticklabels([convertDate(x,'%Y-%m-%d','%Y') if x.split('-')[1]=='01' else convertDate(x,'%Y-%m-%d','%b') for x in xDates if (int(x.split('-')[1])-1)%every==0])\n",
    "ax.tick_params(axis='x',labelsize=8,size=0)\n",
    "\n",
    "# this sets the vertical dashed lines on the tree\n",
    "[ax.axvline(i,ls='--',lw=0.45,color='grey',zorder=0) for i in range(1975,2018,5)]\n",
    "\n",
    "\n",
    "branchWidth=0.35 ## default branch width\n",
    "\n",
    "cmap=mpl.cm.viridis\n",
    "\n",
    "# k objects are tips, nodes, branches\n",
    "for k in tree.Objects: ## iterate over objects in tree\n",
    "    #   x=k.x ## or from x position determined earlier\n",
    "    x=k.absoluteTime ## or use absolute time instead\n",
    "    y=k.y ## get y position from .drawTree that was run earlier, but could be anything else\n",
    "    \n",
    "#     xp=k.parent.x ## get x position of current object's parent\n",
    "    xp=k.parent.absoluteTime ## get x position of current object's parent\n",
    "    if x==None: ## matplotlib won't plot Nones, like root\n",
    "        x=0.0\n",
    "    if xp==None:\n",
    "        xp=x\n",
    "        #c='indianred' ## colour can be fixed\n",
    "#     c=cmap(k.height/ll.treeHeight) ## or be a function of something else\n",
    "#     c=[cmap(k.traits['posterior']) if k.traits.has_key('posterior') else cmap(1.0)][0]\n",
    "    #c=colors[k.traits['geo']]\n",
    "    #c=['#63C78A' if k.traits['geo']=='cambodia' else '#6A34CA'][0] ## can be discrete too\n",
    "    \n",
    "    if isinstance(k,bt.leaf) or k.branchType=='leaf': ## if leaf...\n",
    "        #x=decimalDate(k.name.split('_')[-1],variable=True) ## get x position from name\n",
    "        s=size[k.name]\n",
    "        #s=sizes[k.name] ## tip size specified in a dictionary\n",
    "        #s=50-40*k.height/tree.treeHeight\n",
    "        \n",
    "        \n",
    "        #ax.scatter(x,y,s=s+0.8*s,facecolor='k',edgecolor='none',zorder=10) ## plot black circle underneath\n",
    "    \n",
    "    if isinstance(k, bt.leaf):\n",
    "        if k.name in color:\n",
    "            c = color[k.name]\n",
    "        else: \n",
    "            c = '#CBCACA'\n",
    "            print(k.name)\n",
    "        ax.scatter(x,y,s=s,facecolor=c,edgecolor='none',zorder=11) ## plot circle for every tip\n",
    "    \n",
    "    if k.branchType =='leaf':\n",
    "        c=\"#9F9F9F\"\n",
    "        \n",
    "    elif isinstance(k,bt.node) or k.branchType=='node': ## if node...\n",
    "        c=\"#9F9F9F\"\n",
    "        ax.plot([x,x],[k.children[-1].y,k.children[0].y],lw=branchWidth,color=c,ls='-',zorder=9)\n",
    "        \n",
    "        \n",
    "    ax.plot([xp,x],[y,y],lw=branchWidth,color=c,ls='-',zorder=9)\n",
    "        \n",
    "    # add in a legend\n",
    "    han_list = []\n",
    "\n",
    "    for key in color_key2:\n",
    "        marker = mlines.Line2D(range(1), range(1), color = color_key2[key], marker='o', markerfacecolor = color_key2[key], label = key, markersize = 4)\n",
    "        han_list.append(marker)\n",
    "\n",
    "    ax.legend(handles = han_list, markerfirst = True, frameon=False, loc=3, fontsize=8)\n",
    "\n",
    "        \n",
    "# set axis limits, remove border lines         \n",
    "ax.spines['left'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['bottom'].set_visible(False)\n",
    "\n",
    "ax.set_ylim(-5,tree.ySpan+5)\n",
    "ax.tick_params(axis='y',labelsize=0,size=0)\n",
    "#ax.set_xticklabels([])\n",
    "ax.set_yticklabels([])\n",
    "\n",
    "#ax.set_xlim(bt.decimalDate('2008-01-01'),bt.decimalDate('2016-12-30'))\n",
    "ax.set_ylim(-5,tree.ySpan+5)\n",
    "\n",
    "\n",
    "plt.savefig('/Users/lmoncla/Documents/Fellowships/LSRF/figures/H5N1-Cambodia-only-long-2018-10-01.pdf', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes on TMRCA for human Cambodian samples \n",
    "\n",
    "I used this website to convert dates: http://sopac.ucsd.edu/convertDate.shtml\n",
    "\n",
    "date that all Cambodia 2013 human samples share a common ancestor: 5.2992 \n",
    "date of youngest tip: 2017.6356 (August 20, 2017)\n",
    "\n",
    "2017.6356 - 5.2992 = 2012.3364 (May 3, 2012) \n",
    "\n",
    "95% posterior density: 5.1193, 5.4974\n",
    "\n",
    "2017.6356 - 5.1193 = 2012.5163 (July 7, 2012)\n",
    "\n",
    "2017.6356 - 5.4974 = 2012.1382 (February 20, 2012)\n",
    "\n",
    "Probability is very high that this node was in Cambodia at that time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mumps",
   "language": "python",
   "name": "mumps"
  },
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}