emp_nb_Pedicularis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook 10: \n",
    "\n",
    "This is an IPython notebook. Most of the code is composed of bash scripts, indicated by %%bash at the top of the cell, otherwise it is IPython code. This notebook includes code to download, assemble and analyze a published RADseq data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "### Notebook 10\n",
    "### Data set 10: Pedicularis\n",
    "### Authors: Eaton & Ree 2013\n",
    "### Data Location: SRP021469"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download the sequence data\n",
    "Sequence data for this study are archived on the NCBI sequence read archive (SRA). Below I read in SraRunTable.txt for this project which contains all of the information we need to download the data. \n",
    "\n",
    "+ SRA link: http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP021469"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## make a new directory for this analysis\n",
    "mkdir -p empirical_10/fastq/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For each ERS (individuals) get all of the ERR (sequence file accessions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    BioSample_s Experiment_s          Library_Name_s  MBases_l  MBytes_l  \\\n",
      "0  SAMN01998191    SRX271152           29154_superba        49        26   \n",
      "1  SAMN01999096    SRX840097            30556_thamno       102        55   \n",
      "2  SAMN01998194    SRX840089  41954_cyathophylloides       155        86   \n",
      "3  SAMN01998193    SRX840088  41478_cyathophylloides       155        82   \n",
      "4  SAMN01999094    SRX840095               39618_rex        58        32   \n",
      "\n",
      "                     Organism_s       Run_s SRA_Sample_s  \\\n",
      "0           Pedicularis superba  SRR1754715    SRS415733   \n",
      "1       Pedicularis thamnophila  SRR1754720    SRS816842   \n",
      "2  Pedicularis cyathophylloides  SRR1754721    SRS816834   \n",
      "3  Pedicularis cyathophylloides  SRR1754722    SRS816833   \n",
      "4               Pedicularis rex  SRR1754723    SRS816841   \n",
      "\n",
      "            Sample_Name_s  is_reference_s       ...       LibraryLayout_s  \\\n",
      "0           29154_superba  <not provided>       ...                SINGLE   \n",
      "1            30556_thamno            skip       ...                SINGLE   \n",
      "2  41954_cyathophylloides  <not provided>       ...                SINGLE   \n",
      "3  41478_cyathophylloides  <not provided>       ...                SINGLE   \n",
      "4               39618_rex            skip       ...                SINGLE   \n",
      "\n",
      "       LibrarySelection_s LibrarySource_s  LoadDate_s Platform_s  \\\n",
      "0  Reduced Representation         GENOMIC  2015-01-27   ILLUMINA   \n",
      "1  Reduced Representation         GENOMIC  2015-01-27   ILLUMINA   \n",
      "2  Reduced Representation         GENOMIC  2015-01-27   ILLUMINA   \n",
      "3  Reduced Representation         GENOMIC  2015-01-27   ILLUMINA   \n",
      "4  Reduced Representation         GENOMIC  2015-01-27   ILLUMINA   \n",
      "\n",
      "   ReleaseDate_s SRA_Study_s g1k_analysis_group_s  g1k_pop_code_s  \\\n",
      "0     2015-01-27   SRP021469       <not provided>  <not provided>   \n",
      "1     2015-01-27   SRP021469       <not provided>  <not provided>   \n",
      "2     2015-01-27   SRP021469       <not provided>  <not provided>   \n",
      "3     2015-01-27   SRP021469       <not provided>  <not provided>   \n",
      "4     2015-01-27   SRP021469       <not provided>  <not provided>   \n",
      "\n",
      "         source_s  \n",
      "0  <not provided>  \n",
      "1  <not provided>  \n",
      "2  <not provided>  \n",
      "3  <not provided>  \n",
      "4  <not provided>  \n",
      "\n",
      "[5 rows x 26 columns]\n"
     ]
    }
   ],
   "source": [
    "## IPython code\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import urllib2\n",
    "import os\n",
    "\n",
    "## open the SRA run table from github url\n",
    "url = \"https://raw.githubusercontent.com/\"+\\\n",
    "      \"dereneaton/RADmissing/master/empirical_10_SraRunTable.txt\"\n",
    "intable = urllib2.urlopen(url)\n",
    "indata = pd.read_table(intable, sep=\"\\t\")\n",
    "\n",
    "## print first few rows\n",
    "print indata.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def wget_download(SRR, outdir, outname):\n",
    "    \"\"\" Python function to get sra data from ncbi and write to\n",
    "    outdir with a new name using bash call wget \"\"\"\n",
    "    \n",
    "    ## get output name\n",
    "    output = os.path.join(outdir, outname+\".sra\")\n",
    "    \n",
    "    ## create a call string \n",
    "    call = \"wget -q -r -nH --cut-dirs=9 -O \"+output+\" \"+\\\n",
    "           \"ftp://ftp-trace.ncbi.nlm.nih.gov/\"+\\\n",
    "           \"sra/sra-instant/reads/ByRun/sra/SRR/\"+\\\n",
    "           \"{}/{}/{}.sra;\".format(SRR[:6], SRR, SRR)\n",
    "        \n",
    "    ## call bash script\n",
    "    ! $call "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we pass the SRR number and the sample name to the `wget_download` function so that the files are saved with their sample names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for ID, SRR in zip(indata.Sample_Name_s, indata.Run_s):\n",
    "    wget_download(SRR, \"empirical_10/fastq/\", ID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read 696994 spots for empirical_10/fastq/29154_superba.sra\n",
      "Written 696994 spots for empirical_10/fastq/29154_superba.sra\n",
      "Read 1452316 spots for empirical_10/fastq/30556_thamno.sra\n",
      "Written 1452316 spots for empirical_10/fastq/30556_thamno.sra\n",
      "Read 1253109 spots for empirical_10/fastq/30686_cyathophylla.sra\n",
      "Written 1253109 spots for empirical_10/fastq/30686_cyathophylla.sra\n",
      "Read 964244 spots for empirical_10/fastq/32082_przewalskii.sra\n",
      "Written 964244 spots for empirical_10/fastq/32082_przewalskii.sra\n",
      "Read 636625 spots for empirical_10/fastq/33413_thamno.sra\n",
      "Written 636625 spots for empirical_10/fastq/33413_thamno.sra\n",
      "Read 1002923 spots for empirical_10/fastq/33588_przewalskii.sra\n",
      "Written 1002923 spots for empirical_10/fastq/33588_przewalskii.sra\n",
      "Read 1803858 spots for empirical_10/fastq/35236_rex.sra\n",
      "Written 1803858 spots for empirical_10/fastq/35236_rex.sra\n",
      "Read 1409843 spots for empirical_10/fastq/35855_rex.sra\n",
      "Written 1409843 spots for empirical_10/fastq/35855_rex.sra\n",
      "Read 1391175 spots for empirical_10/fastq/38362_rex.sra\n",
      "Written 1391175 spots for empirical_10/fastq/38362_rex.sra\n",
      "Read 822263 spots for empirical_10/fastq/39618_rex.sra\n",
      "Written 822263 spots for empirical_10/fastq/39618_rex.sra\n",
      "Read 1707942 spots for empirical_10/fastq/40578_rex.sra\n",
      "Written 1707942 spots for empirical_10/fastq/40578_rex.sra\n",
      "Read 2199740 spots for empirical_10/fastq/41478_cyathophylloides.sra\n",
      "Written 2199740 spots for empirical_10/fastq/41478_cyathophylloides.sra\n",
      "Read 2199613 spots for empirical_10/fastq/41954_cyathophylloides.sra\n",
      "Written 2199613 spots for empirical_10/fastq/41954_cyathophylloides.sra\n",
      "Read 17540645 spots total\n",
      "Written 17540645 spots total\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "## convert sra files to fastq using fastq-dump tool\n",
    "## output as gzipped into the fastq directory\n",
    "fastq-dump --gzip -O empirical_10/fastq/ empirical_10/fastq/*.sra\n",
    "\n",
    "## remove .sra files\n",
    "rm empirical_10/fastq/*.sra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 1.1G\n",
      "-rw-rw-r-- 1 deren deren  42M Dec 17 20:49 29154_superba.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  86M Dec 17 20:49 30556_thamno.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  80M Dec 17 20:49 30686_cyathophylla.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  58M Dec 17 20:49 32082_przewalskii.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  40M Dec 17 20:49 33413_thamno.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  61M Dec 17 20:49 33588_przewalskii.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 108M Dec 17 20:49 35236_rex.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  84M Dec 17 20:49 35855_rex.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  82M Dec 17 20:49 38362_rex.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren  50M Dec 17 20:49 39618_rex.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 100M Dec 17 20:49 40578_rex.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 128M Dec 17 20:49 41478_cyathophylloides.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 136M Dec 17 20:49 41954_cyathophylloides.fastq.gz\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "ls -lh empirical_10/fastq/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a params file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pyRAD 3.0.63\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "pyrad --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\tnew params.txt file created\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "## remove old params file if it exists\n",
    "rm params.txt \n",
    "\n",
    "## create a new default params file\n",
    "pyrad -n "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note: \n",
    "The data here are from Illumina Casava <1.8, so the phred scores are offset by 64 instead of 33, so we use that in the params file below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## substitute new parameters into file\n",
    "sed -i '/## 1. /c\\empirical_10/          ## 1. working directory ' params.txt\n",
    "sed -i '/## 6. /c\\TGCAG                  ## 6. cutters ' params.txt\n",
    "sed -i '/## 7. /c\\20                     ## 7. N processors      ' params.txt\n",
    "sed -i '/## 9. /c\\6                      ## 9. NQual             ' params.txt\n",
    "sed -i '/## 10./c\\.85                    ## 10. clust threshold  ' params.txt\n",
    "sed -i '/## 12./c\\4                      ## 12. MinCov           ' params.txt\n",
    "sed -i '/## 13./c\\10                     ## 13. maxSH            ' params.txt\n",
    "sed -i '/## 14./c\\empirical_10_m4         ## 14. output name      ' params.txt\n",
    "sed -i '/## 18./c\\empirical_10/fastq/*.gz ## 18. data location    ' params.txt\n",
    "sed -i '/## 29./c\\2,2                    ## 29. trim overhang    ' params.txt\n",
    "sed -i '/## 30./c\\p,n,s                  ## 30. output formats   ' params.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==** parameter inputs for pyRAD version 3.0.63  **======================== affected step ==\r\n",
      "empirical_10/          ## 1. working directory \r\n",
      "./*.fastq.gz              ## 2. Loc. of non-demultiplexed files (if not line 18)  (s1)\r\n",
      "./*.barcodes              ## 3. Loc. of barcode file (if not line 18)             (s1)\r\n",
      "vsearch                   ## 4. command (or path) to call vsearch (or usearch)    (s3,s6)\r\n",
      "muscle                    ## 5. command (or path) to call muscle                  (s3,s7)\r\n",
      "TGCAG                  ## 6. cutters \r\n",
      "20                     ## 7. N processors      \r\n",
      "6                         ## 8. Mindepth: min coverage for a cluster              (s4,s5)\r\n",
      "6                      ## 9. NQual             \r\n",
      ".85                    ## 10. clust threshold  \r\n",
      "rad                       ## 11. Datatype: rad,gbs,pairgbs,pairddrad,(others:see docs)(all)\r\n",
      "4                      ## 12. MinCov           \r\n",
      "10                     ## 13. maxSH            \r\n",
      "empirical_10_m4         ## 14. output name      \r\n",
      "==== optional params below this line ===================================  affected step ==\r\n",
      "                       ## 15.opt.: select subset (prefix* only selector)            (s2-s7)\r\n",
      "                       ## 16.opt.: add-on (outgroup) taxa (list or prefix*)         (s6,s7)\r\n",
      "                       ## 17.opt.: exclude taxa (list or prefix*)                   (s7)\r\n",
      "empirical_10/fastq/*.gz ## 18. data location    \r\n",
      "                       ## 19.opt.: maxM: N mismatches in barcodes (def= 1)          (s1)\r\n",
      "                       ## 20.opt.: phred Qscore offset (def= 33)                    (s2)\r\n",
      "                       ## 21.opt.: filter: def=0=NQual 1=NQual+adapters. 2=strict   (s2)\r\n",
      "                       ## 22.opt.: a priori E,H (def= 0.001,0.01, if not estimated) (s5)\r\n",
      "                       ## 23.opt.: maxN: max Ns in a cons seq (def=5)               (s5)\r\n",
      "                       ## 24.opt.: maxH: max heterozyg. sites in cons seq (def=5)   (s5)\r\n",
      "                       ## 25.opt.: ploidy: max alleles in cons seq (def=2;see docs) (s4,s5)\r\n",
      "                       ## 26.opt.: maxSNPs: (def=100). Paired (def=100,100)         (s7)\r\n",
      "                       ## 27.opt.: maxIndels: within-clust,across-clust (def. 3,99) (s3,s7)\r\n",
      "                       ## 28.opt.: random number seed (def. 112233)              (s3,s6,s7)\r\n",
      "2,2                    ## 29. trim overhang    \r\n",
      "p,n,s                  ## 30. output formats   \r\n",
      "                       ## 31.opt.: maj. base call at depth>x<mindepth (def.x=mindepth) (s5)\r\n",
      "                       ## 32.opt.: keep trimmed reads (def=0). Enter min length.    (s2)\r\n",
      "                       ## 33.opt.: max stack size (int), def= max(500,mean+2*SD)    (s3)\r\n",
      "                       ## 34.opt.: minDerep: exclude dereps with <= N copies, def=1 (s3)\r\n",
      "                       ## 35.opt.: use hierarchical clustering (def.=0, 1=yes)      (s6)\r\n",
      "                       ## 36.opt.: repeat masking (def.=1='dust' method, 0=no)      (s3,s6)\r\n",
      "                       ## 37.opt.: vsearch max threads per job (def.=6; see docs)   (s3,s6)\r\n",
      "==== optional: list group/clade assignments below this line (see docs) ==================\r\n"
     ]
    }
   ],
   "source": [
    "cat params.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Assemble in _pyrad_\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "pyrad -p params.txt -s 234567 >> log.txt 2>&1  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "sed -i '/## 12./c\\2                    ## 12. MinCov           ' params.txt\n",
    "sed -i '/## 14./c\\empirical_10_m2       ## 14. output name      ' params.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "pyrad -p params.txt -s 7 >> log.txt 2>&1  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "We are interested in the relationship between the amount of input (raw) data between any two samples, the average coverage they recover when clustered together, and the phylogenetic distances separating samples.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Raw data amounts\n",
    "The average number of raw reads per sample is 1.36M."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count         13.000000\n",
      "mean     1263500.153846\n",
      "std       494927.550159\n",
      "min       580917.000000\n",
      "25%       901009.000000\n",
      "50%      1312128.000000\n",
      "75%      1614305.000000\n",
      "max      2083670.000000\n",
      "Name: passed.total, dtype: float64\n",
      "\n",
      "most raw data in sample:\n",
      "11    41478_cyathophylloides\n",
      "Name: sample , dtype: object\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "## read in the data\n",
    "sdat = pd.read_table(\"empirical_10/stats/s2.rawedit.txt\", header=0, nrows=14)\n",
    "\n",
    "## print summary stats\n",
    "print sdat[\"passed.total\"].describe()\n",
    "\n",
    "## find which sample has the most raw data\n",
    "maxraw = sdat[\"passed.total\"].max()\n",
    "print \"\\nmost raw data in sample:\"\n",
    "print sdat['sample '][sdat['passed.total']==maxraw]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at distributions of coverage\n",
    "pyrad v.3.0.63 outputs depth information for each sample which I read in here and plot. First let's ask which sample has the highest depth of coverage. The std of coverages is pretty low in this data set compared to several others. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "summary of means\n",
      "==================\n",
      "count    13.000000\n",
      "mean      6.665923\n",
      "std       2.428387\n",
      "min       3.623000\n",
      "25%       5.075000\n",
      "50%       5.935000\n",
      "75%       7.577000\n",
      "max      12.438000\n",
      "Name: dpt.me, dtype: float64\n",
      "\n",
      "summary of std\n",
      "==================\n",
      "count     13.000000\n",
      "mean      60.777000\n",
      "std       32.695511\n",
      "min       23.321000\n",
      "25%       37.463000\n",
      "50%       57.140000\n",
      "75%       77.513000\n",
      "max      116.124000\n",
      "Name: dpt.sd, dtype: float64\n",
      "\n",
      "summary of proportion lowdepth\n",
      "==================\n",
      "count    13.000000\n",
      "mean      0.728338\n",
      "std       0.069037\n",
      "min       0.590968\n",
      "25%       0.698748\n",
      "50%       0.735235\n",
      "75%       0.743905\n",
      "max       0.866382\n",
      "dtype: float64\n",
      "\n",
      "highest coverage in sample:\n",
      "8    38362_rex\n",
      "Name: taxa, dtype: object\n"
     ]
    }
   ],
   "source": [
    "## read in the s3 results\n",
    "sdat = pd.read_table(\"empirical_10/stats/s3.clusters.txt\", header=0, nrows=14)\n",
    "\n",
    "## print summary stats\n",
    "print \"summary of means\\n==================\"\n",
    "print sdat['dpt.me'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of std\\n==================\"\n",
    "print sdat['dpt.sd'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of proportion lowdepth\\n==================\"\n",
    "print pd.Series(1-sdat['d>5.tot']/sdat[\"total\"]).describe()\n",
    "\n",
    "## find which sample has the greatest depth of retained loci\n",
    "max_hiprop = (sdat[\"d>5.tot\"]/sdat[\"total\"]).max()\n",
    "print \"\\nhighest coverage in sample:\"\n",
    "print sdat['taxa'][sdat['d>5.tot']/sdat[\"total\"]==max_hiprop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "highest prop coverage in sample:\n",
      "8    38362_rex\n",
      "Name: taxa, dtype: object\n"
     ]
    }
   ],
   "source": [
    "maxprop =(sdat['d>5.tot']/sdat['total']).max()\n",
    "print \"\\nhighest prop coverage in sample:\"\n",
    "print sdat['taxa'][sdat['d>5.tot']/sdat['total']==maxprop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Means for sample 38362_rex\n",
      "10.2151397661 101.259826132\n",
      "22.714617861 157.486621215\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "## print mean and std of coverage for the highest coverage sample\n",
    "with open(\"empirical_10/clust.85/38362_rex.depths\", 'rb') as indat:\n",
    "    depths = np.array(indat.read().strip().split(\",\"), dtype=int)\n",
    "    \n",
    "print \"Means for sample 38362_rex\"\n",
    "print depths.mean(), depths.std()\n",
    "print depths[depths>5].mean(), depths[depths>5].std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the coverage for the sample with highest mean coverage\n",
    "Green shows the loci that were discarded and orange the loci that were retained. The majority of data were discarded for being too low of coverage. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<toyplot.mark.BarMagnitudes at 0x7ff250987e90>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<div align=\"center\" class=\"toyplot\" id=\"t38d8a6346130409684866126fc30dbdb\"><svg height=\"300.0px\" id=\"tb754be2a5d9f4ff89270b1000c9c287d\" preserveAspectRatio=\"xMidyMid meet\" style=\"background-color:transparent;fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:helvetica;font-size:12px;opacity:1.0;stroke:rgb(16.1%,15.3%,14.1%);stroke-opacity:1.0;stroke-width:1.0\" viewBox=\"0 0 350.0 300.0\" width=\"350.0px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:toyplot=\"http://www.sandia.gov/toyplot\"><g class=\"toyplot-axes-Cartesian\" id=\"te2ac3316e4d94b2c97a8b25cec1766a5\"><clipPath id=\"ta21867796812428283b171c97966fc5e\"><rect height=\"200.0\" width=\"250.0\" x=\"50.0\" y=\"50.0\"></rect></clipPath><g class=\"toyplot-coordinate-events\" clip-path=\"url(#ta21867796812428283b171c97966fc5e)\" style=\"cursor:crosshair\"><rect height=\"200.0\" style=\"pointer-events:all;visibility:hidden\" width=\"250.0\" x=\"50.0\" y=\"50.0\"></rect><g class=\"toyplot-mark-BarMagnitudes\" id=\"t4d18450398204b2e876bb687211e2ec3\" style=\"stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\"><g class=\"toyplot-Series\"><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"60.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"177.08760000000001\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"67.666666666666671\" y=\"62.912399999999998\"></rect><rect class=\"toyplot-Datum\" height=\"48.031199999999984\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"75.333333333333329\" y=\"191.96880000000002\"></rect><rect class=\"toyplot-Datum\" height=\"18.331200000000024\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"83.0\" y=\"221.66879999999998\"></rect><rect class=\"toyplot-Datum\" height=\"10.288800000000009\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"90.666666666666657\" y=\"229.71119999999999\"></rect><rect class=\"toyplot-Datum\" height=\"7.9992000000000019\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"98.333333333333329\" y=\"232.0008\"></rect><rect class=\"toyplot-Datum\" height=\"7.0020000000000095\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"106.0\" y=\"232.99799999999999\"></rect><rect class=\"toyplot-Datum\" height=\"6.6887999999999863\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"113.66666666666667\" y=\"233.31120000000001\"></rect><rect class=\"toyplot-Datum\" height=\"7.466399999999993\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"121.33333333333333\" y=\"232.53360000000001\"></rect><rect class=\"toyplot-Datum\" height=\"8.2727999999999895\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"129.0\" y=\"231.72720000000001\"></rect><rect class=\"toyplot-Datum\" height=\"8.4527999999999963\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"136.66666666666666\" y=\"231.5472\"></rect><rect class=\"toyplot-Datum\" height=\"9.3348000000000013\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"144.33333333333331\" y=\"230.6652\"></rect><rect class=\"toyplot-Datum\" height=\"10.018800000000027\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"152.0\" y=\"229.98119999999997\"></rect><rect class=\"toyplot-Datum\" height=\"10.515600000000006\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"159.66666666666669\" y=\"229.48439999999999\"></rect><rect class=\"toyplot-Datum\" height=\"10.900799999999975\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"167.33333333333334\" y=\"229.09920000000002\"></rect><rect class=\"toyplot-Datum\" height=\"10.443600000000004\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"175.0\" y=\"229.5564\"></rect><rect class=\"toyplot-Datum\" height=\"10.400399999999991\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"182.66666666666666\" y=\"229.59960000000001\"></rect><rect class=\"toyplot-Datum\" height=\"9.529200000000003\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"190.33333333333334\" y=\"230.4708\"></rect><rect class=\"toyplot-Datum\" height=\"9.0\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"198.0\" y=\"231.0\"></rect><rect class=\"toyplot-Datum\" height=\"8.0315999999999974\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"205.66666666666666\" y=\"231.9684\"></rect><rect class=\"toyplot-Datum\" height=\"7.3692000000000064\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"213.33333333333331\" y=\"232.63079999999999\"></rect><rect class=\"toyplot-Datum\" height=\"6.3972000000000264\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"221.0\" y=\"233.60279999999997\"></rect><rect class=\"toyplot-Datum\" height=\"5.4971999999999923\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"228.66666666666666\" y=\"234.50280000000001\"></rect><rect class=\"toyplot-Datum\" height=\"4.68719999999999\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"236.33333333333334\" y=\"235.31280000000001\"></rect><rect class=\"toyplot-Datum\" height=\"3.927599999999984\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"244.0\" y=\"236.07240000000002\"></rect><rect class=\"toyplot-Datum\" height=\"3.0671999999999855\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"251.66666666666669\" y=\"236.93280000000001\"></rect><rect class=\"toyplot-Datum\" height=\"2.5776000000000181\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"259.33333333333331\" y=\"237.42239999999998\"></rect><rect class=\"toyplot-Datum\" height=\"2.0735999999999706\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"267.0\" y=\"237.92640000000003\"></rect><rect class=\"toyplot-Datum\" height=\"3.1968000000000245\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"274.66666666666669\" y=\"236.80319999999998\"></rect></g></g><g class=\"toyplot-mark-BarMagnitudes\" id=\"t74a3d7a9e1a0423c966b4bdf6b72e910\" style=\"stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\"><g class=\"toyplot-Series\"><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"60.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"67.666666666666671\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"75.333333333333329\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"83.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"90.666666666666657\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"98.333333333333329\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"7.0020000000000095\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"106.0\" y=\"232.99799999999999\"></rect><rect class=\"toyplot-Datum\" height=\"6.6887999999999863\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"113.66666666666667\" y=\"233.31120000000001\"></rect><rect class=\"toyplot-Datum\" height=\"7.466399999999993\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"121.33333333333333\" y=\"232.53360000000001\"></rect><rect class=\"toyplot-Datum\" height=\"8.2727999999999895\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"129.0\" y=\"231.72720000000001\"></rect><rect class=\"toyplot-Datum\" height=\"8.4527999999999963\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"136.66666666666666\" y=\"231.5472\"></rect><rect class=\"toyplot-Datum\" height=\"9.3348000000000013\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"144.33333333333331\" y=\"230.6652\"></rect><rect class=\"toyplot-Datum\" height=\"10.018800000000027\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"152.0\" y=\"229.98119999999997\"></rect><rect class=\"toyplot-Datum\" height=\"10.515600000000006\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"159.66666666666669\" y=\"229.48439999999999\"></rect><rect class=\"toyplot-Datum\" height=\"10.900799999999975\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"167.33333333333334\" y=\"229.09920000000002\"></rect><rect class=\"toyplot-Datum\" height=\"10.443600000000004\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"175.0\" y=\"229.5564\"></rect><rect class=\"toyplot-Datum\" height=\"10.400399999999991\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"182.66666666666666\" y=\"229.59960000000001\"></rect><rect class=\"toyplot-Datum\" height=\"9.529200000000003\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"190.33333333333334\" y=\"230.4708\"></rect><rect class=\"toyplot-Datum\" height=\"9.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"198.0\" y=\"231.0\"></rect><rect class=\"toyplot-Datum\" height=\"8.0315999999999974\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"205.66666666666666\" y=\"231.9684\"></rect><rect class=\"toyplot-Datum\" height=\"7.3692000000000064\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"213.33333333333331\" y=\"232.63079999999999\"></rect><rect class=\"toyplot-Datum\" height=\"6.3972000000000264\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"221.0\" y=\"233.60279999999997\"></rect><rect class=\"toyplot-Datum\" height=\"5.4971999999999923\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"228.66666666666666\" y=\"234.50280000000001\"></rect><rect class=\"toyplot-Datum\" height=\"4.68719999999999\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"236.33333333333334\" y=\"235.31280000000001\"></rect><rect class=\"toyplot-Datum\" height=\"3.927599999999984\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"244.0\" y=\"236.07240000000002\"></rect><rect class=\"toyplot-Datum\" height=\"3.0671999999999855\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"251.66666666666669\" y=\"236.93280000000001\"></rect><rect class=\"toyplot-Datum\" height=\"2.5776000000000181\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"259.33333333333331\" y=\"237.42239999999998\"></rect><rect class=\"toyplot-Datum\" height=\"2.0735999999999706\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"267.0\" y=\"237.92640000000003\"></rect><rect class=\"toyplot-Datum\" height=\"3.1968000000000245\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"274.66666666666669\" y=\"236.80319999999998\"></rect></g></g></g><g class=\"toyplot-coordinates\" style=\"visibility:hidden\"><rect height=\"14.0\" style=\"fill:rgb(100%,100%,100%);fill-opacity:1.0;opacity:0.75;stroke:none\" width=\"90.0\" x=\"200.0\" y=\"60.0\"></rect><text style=\"alignment-baseline:middle;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"245.0\" y=\"67.0\"></text></g><line style=\"\" x1=\"60.0\" x2=\"282.3333333333333\" y1=\"250.0\" y2=\"250.0\"></line><g><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"60.0\" y=\"250.0\">0</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"136.66666666666666\" y=\"250.0\">10</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"213.33333333333331\" y=\"250.0\">20</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"290.0\" y=\"250.0\">30</text></g><text style=\"alignment-baseline:middle;baseline-shift:-200%;font-weight:bold;stroke:none;text-anchor:middle\" x=\"175.0\" y=\"250.0\">Depth of coverage (N reads)</text><line style=\"\" x1=\"50.0\" x2=\"50.0\" y1=\"62.9124\" y2=\"240.0\"></line><g><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 240.0)\" x=\"50.0\" y=\"240.0\">0</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 204.0)\" x=\"50.0\" y=\"204.0\">10000</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 168.0)\" x=\"50.0\" y=\"168.0\">20000</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 132.0)\" x=\"50.0\" y=\"132.0\">30000</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 95.99999999999999)\" x=\"50.0\" y=\"95.99999999999999\">40000</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 60.0)\" x=\"50.0\" y=\"60.0\">50000</text></g><text style=\"alignment-baseline:middle;baseline-shift:200%;font-weight:bold;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 150.0)\" x=\"50.0\" y=\"150.0\">N loci</text><text style=\"alignment-baseline:middle;baseline-shift:100%;font-size:14px;font-weight:bold;stroke:none;text-anchor:middle\" x=\"175.0\" y=\"50.0\">dataset10/sample=38362_rex</text></g></svg><div class=\"toyplot-controls\"><ul class=\"toyplot-mark-popup\" onmouseleave=\"this.style.visibility='hidden'\" style=\"background:rgba(0%,0%,0%,0.75);border:0;border-radius:6px;color:white;cursor:default;list-style:none;margin:0;padding:5px;position:fixed;visibility:hidden\"><li class=\"toyplot-mark-popup-title\" style=\"color:lightgray;cursor:default;padding:5px;list-style:none;margin:0;\"></li><li class=\"toyplot-mark-popup-save-csv\" onmouseout=\"this.style.color='white';this.style.background='steelblue'\" onmouseover=\"this.style.color='steelblue';this.style.background='white'\" style=\"border-radius:3px;padding:5px;list-style:none;margin:0;\">Save as .csv</li></ul><script>\n",
       "(function()\n",
       "{\n",
       "  if(window.CSS !== undefined && window.CSS.supports !== undefined)\n",
       "  {\n",
       "    if(!window.CSS.supports(\"alignment-baseline\", \"middle\"))\n",
       "    {\n",
       "      var re = /\\s*alignment-baseline\\s*:\\s*([^;\\s]*)\\s*/;\n",
       "      var text = document.querySelectorAll(\"#t38d8a6346130409684866126fc30dbdb text\");\n",
       "      for(var i = 0; i != text.length; ++i)\n",
       "      {\n",
       "        var match = re.exec(text[i].attributes.style.value);\n",
       "        if(match)\n",
       "        {\n",
       "          if(match[1] == \"middle\")\n",
       "          {\n",
       "            var style = getComputedStyle(text[i]);\n",
       "            var font_size = style.fontSize.substr(0, style.fontSize.length - 2);\n",
       "            var dy = text[i].dy.baseVal.length ? text[i].dy.baseVal[0].value : 0;\n",
       "            dy += 0.4 * font_size;\n",
       "            text[i].setAttribute(\"dy\", dy);\n",
       "          }\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "    if(!window.CSS.supports(\"baseline-shift\", \"0\"))\n",
       "    {\n",
       "      var re = /\\s*baseline-shift\\s*:\\s*([^;\\s]*)\\s*/;\n",
       "      var text = document.querySelectorAll(\"#t38d8a6346130409684866126fc30dbdb text\");\n",
       "      for(var i = 0; i != text.length; ++i)\n",
       "      {\n",
       "        var match = re.exec(text[i].attributes.style.value);\n",
       "        if(match)\n",
       "        {\n",
       "          var style = getComputedStyle(text[i]);\n",
       "          var font_size = style.fontSize.substr(0, style.fontSize.length - 2);\n",
       "          var percent = 0.01 * match[1].substr(0, match[1].length-1);\n",
       "          var dy = text[i].dy.baseVal.length ? text[i].dy.baseVal[0].value : 0;\n",
       "          dy -= percent * font_size\n",
       "          text[i].setAttribute(\"dy\", dy);\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "  }\n",
       "})();\n",
       "</script><script>\n",
       "(function()\n",
       "{\n",
       "  var data_tables = [{\"data\": [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 49191.0, 13342.0, 5092.0, 2858.0, 2222.0, 1945.0, 1858.0, 2074.0, 2298.0, 2348.0, 2593.0, 2783.0, 2921.0, 3028.0, 2901.0, 2889.0, 2647.0, 2500.0, 2231.0, 2047.0, 1777.0, 1527.0, 1302.0, 1091.0, 852.0, 716.0, 576.0, 888.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"t4d18450398204b2e876bb687211e2ec3\", \"filename\": \"toyplot\"}, {\"data\": [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1945.0, 1858.0, 2074.0, 2298.0, 2348.0, 2593.0, 2783.0, 2921.0, 3028.0, 2901.0, 2889.0, 2647.0, 2500.0, 2231.0, 2047.0, 1777.0, 1527.0, 1302.0, 1091.0, 852.0, 716.0, 576.0, 888.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"t74a3d7a9e1a0423c966b4bdf6b72e910\", \"filename\": \"toyplot\"}];\n",
       "\n",
       "  function save_csv(data_table)\n",
       "  {\n",
       "    var uri = \"data:text/csv;charset=utf-8,\";\n",
       "    uri += data_table.names.join(\",\") + \"\\n\";\n",
       "    for(var i = 0; i != data_table.data[0].length; ++i)\n",
       "    {\n",
       "      for(var j = 0; j != data_table.data.length; ++j)\n",
       "      {\n",
       "        if(j)\n",
       "          uri += \",\";\n",
       "        uri += data_table.data[j][i];\n",
       "      }\n",
       "      uri += \"\\n\";\n",
       "    }\n",
       "    uri = encodeURI(uri);\n",
       "\n",
       "    var link = document.createElement(\"a\");\n",
       "    if(typeof link.download != \"undefined\")\n",
       "    {\n",
       "      link.href = uri;\n",
       "      link.style = \"visibility:hidden\";\n",
       "      link.download = data_table.filename + \".csv\";\n",
       "\n",
       "      document.body.appendChild(link);\n",
       "      link.click();\n",
       "      document.body.removeChild(link);\n",
       "    }\n",
       "    else\n",
       "    {\n",
       "      window.open(uri);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function open_popup(data_table)\n",
       "  {\n",
       "    return function(e)\n",
       "    {\n",
       "      var popup = document.querySelector(\"#t38d8a6346130409684866126fc30dbdb .toyplot-mark-popup\");\n",
       "      popup.querySelector(\".toyplot-mark-popup-title\").innerHTML = data_table.title;\n",
       "      popup.querySelector(\".toyplot-mark-popup-save-csv\").onclick = function() { popup.style.visibility = \"hidden\"; save_csv(data_table); }\n",
       "      popup.style.left = (e.clientX - 50) + \"px\";\n",
       "      popup.style.top = (e.clientY - 20) + \"px\";\n",
       "      popup.style.visibility = \"visible\";\n",
       "      e.stopPropagation();\n",
       "      e.preventDefault();\n",
       "    }\n",
       "\n",
       "  }\n",
       "\n",
       "  for(var i = 0; i != data_tables.length; ++i)\n",
       "  {\n",
       "    var data_table = data_tables[i];\n",
       "    var event_target = document.querySelector(\"#\" + data_table.id);\n",
       "    event_target.oncontextmenu = open_popup(data_table);\n",
       "  }\n",
       "})();\n",
       "</script><script>\n",
       "(function()\n",
       "{\n",
       "  var axes = {\"te2ac3316e4d94b2c97a8b25cec1766a5\": {\"x\": [{\"domain\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 30.0, \"min\": 0}, \"range\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 290.0, \"min\": 60.0}, \"scale\": \"linear\"}], \"y\": [{\"domain\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 50000.0, \"min\": 0.0}, \"range\": {\"bounds\": {\"max\": -Infinity, \"min\": Infinity}, \"max\": 60.0, \"min\": 240.0}, \"scale\": \"linear\"}]}};\n",
       "\n",
       "  function sign(x)\n",
       "  {\n",
       "    return x < 0 ? -1 : x > 0 ? 1 : 0;\n",
       "  }\n",
       "\n",
       "  function _mix(a, b, amount)\n",
       "  {\n",
       "    return ((1.0 - amount) * a) + (amount * b);\n",
       "  }\n",
       "\n",
       "  function _log(x, base)\n",
       "  {\n",
       "    return Math.log(Math.abs(x)) / Math.log(base);\n",
       "  }\n",
       "\n",
       "  function _in_range(a, x, b)\n",
       "  {\n",
       "    var left = Math.min(a, b);\n",
       "    var right = Math.max(a, b);\n",
       "    return left <= x && x <= right;\n",
       "  }\n",
       "\n",
       "  function to_domain(projection, range)\n",
       "  {\n",
       "    for(var i = 0; i != projection.length; ++i)\n",
       "    {\n",
       "      var segment = projection[i];\n",
       "      if(_in_range(segment.range.bounds.min, range, segment.range.bounds.max))\n",
       "      {\n",
       "        if(segment.scale == \"linear\")\n",
       "        {\n",
       "          var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "          return _mix(segment.domain.min, segment.domain.max, amount)\n",
       "        }\n",
       "        else if(segment.scale[0] == \"log\")\n",
       "        {\n",
       "          var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "          var base = segment.scale[1];\n",
       "          return sign(segment.domain.min) * Math.pow(base, _mix(_log(segment.domain.min, base), _log(segment.domain.max, base), amount));\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "  }\n",
       "\n",
       "  // Compute mouse coordinates relative to a DOM object, with thanks to d3js.org, where this code originated.\n",
       "  function d3_mousePoint(container, e)\n",
       "  {\n",
       "    if (e.changedTouches) e = e.changedTouches[0];\n",
       "    var svg = container.ownerSVGElement || container;\n",
       "    if (svg.createSVGPoint) {\n",
       "      var point = svg.createSVGPoint();\n",
       "      point.x = e.clientX, point.y = e.clientY;\n",
       "      point = point.matrixTransform(container.getScreenCTM().inverse());\n",
       "      return [point.x, point.y];\n",
       "    }\n",
       "    var rect = container.getBoundingClientRect();\n",
       "    return [e.clientX - rect.left - container.clientLeft, e.clientY - rect.top - container.clientTop];\n",
       "  };\n",
       "\n",
       "  function display_coordinates(e)\n",
       "  {\n",
       "    var dom_axes = e.currentTarget.parentElement;\n",
       "    var data = axes[dom_axes.id];\n",
       "\n",
       "    point = d3_mousePoint(e.target, e);\n",
       "    var x = Number(to_domain(data[\"x\"], point[0])).toFixed(2);\n",
       "    var y = Number(to_domain(data[\"y\"], point[1])).toFixed(2);\n",
       "\n",
       "    var coordinates = dom_axes.querySelectorAll(\".toyplot-coordinates\");\n",
       "    for(var i = 0; i != coordinates.length; ++i)\n",
       "    {\n",
       "      coordinates[i].style.visibility = \"visible\";\n",
       "      coordinates[i].querySelector(\"text\").textContent = \"x=\" + x + \" y=\" + y;\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function clear_coordinates(e)\n",
       "  {\n",
       "    var dom_axes = e.currentTarget.parentElement;\n",
       "    var coordinates = dom_axes.querySelectorAll(\".toyplot-coordinates\");\n",
       "    for(var i = 0; i != coordinates.length; ++i)\n",
       "      coordinates[i].style.visibility = \"hidden\";\n",
       "  }\n",
       "\n",
       "  for(var axes_id in axes)\n",
       "  {\n",
       "    var event_target = document.querySelector(\"#\" + axes_id + \" .toyplot-coordinate-events\");\n",
       "    event_target.onmousemove = display_coordinates;\n",
       "    event_target.onmouseout = clear_coordinates;\n",
       "  }\n",
       "})();\n",
       "</script></div></div>"
      ],
      "text/plain": [
       "<toyplot.canvas.Canvas at 0x7ff2508b2ed0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import toyplot\n",
    "import toyplot.svg\n",
    "import numpy as np\n",
    "\n",
    "## read in the depth information for this sample\n",
    "with open(\"empirical_10/clust.85/38362_rex.depths\", 'rb') as indat:\n",
    "    depths = np.array(indat.read().strip().split(\",\"), dtype=int)\n",
    "    \n",
    "## make a barplot in Toyplot\n",
    "canvas = toyplot.Canvas(width=350, height=300)\n",
    "axes = canvas.axes(xlabel=\"Depth of coverage (N reads)\", \n",
    "                   ylabel=\"N loci\", \n",
    "                   label=\"dataset10/sample=38362_rex\")\n",
    "\n",
    "## select the loci with depth > 5 (kept)\n",
    "keeps = depths[depths>5]\n",
    "\n",
    "## plot kept and discarded loci\n",
    "edat = np.histogram(depths, range(30)) # density=True)\n",
    "kdat = np.histogram(keeps, range(30)) #, density=True)\n",
    "axes.bars(edat)\n",
    "axes.bars(kdat)\n",
    "\n",
    "#toyplot.svg.render(canvas, \"empirical_10_depthplot.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print final stats table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "\r\n",
      "44945       ## loci with > minsp containing data\r\n",
      "44875       ## loci with > minsp containing data & paralogs removed\r\n",
      "44875       ## loci with > minsp containing data & paralogs removed & final filtering\r\n",
      "\r\n",
      "## number of loci recovered in final data set for each taxon.\r\n",
      "taxon\tnloci\r\n",
      "29154_superba         \t22283\r\n",
      "30556_thamno          \t34175\r\n",
      "30686_cyathophylla    \t28448\r\n",
      "32082_przewalskii     \t15422\r\n",
      "33413_thamno          \t19521\r\n",
      "33588_przewalskii     \t18022\r\n",
      "35236_rex             \t35638\r\n",
      "35855_rex             \t35508\r\n",
      "38362_rex             \t36300\r\n",
      "39618_rex             \t29769\r\n",
      "40578_rex             \t36407\r\n",
      "41478_cyathophylloides\t33687\r\n",
      "41954_cyathophylloides\t30314\r\n",
      "\r\n",
      "\r\n",
      "## nloci = number of loci with data for exactly ntaxa\r\n",
      "## ntotal = number of loci for which at least ntaxa have data\r\n",
      "ntaxa\tnloci\tsaved\tntotal\r\n",
      "1\t-\r\n",
      "2\t-\t\t-\r\n",
      "3\t-\t\t-\r\n",
      "4\t5740\t*\t44875\r\n",
      "5\t4506\t*\t39135\r\n",
      "6\t3988\t*\t34629\r\n",
      "7\t3585\t*\t30641\r\n",
      "8\t3779\t*\t27056\r\n",
      "9\t4729\t*\t23277\r\n",
      "10\t5527\t*\t18548\r\n",
      "11\t5863\t*\t13021\r\n",
      "12\t4629\t*\t7158\r\n",
      "13\t2529\t*\t2529\r\n",
      "\r\n",
      "\r\n",
      "## nvar = number of loci containing n variable sites (pis+autapomorphies).\r\n",
      "## sumvar = sum of variable sites (SNPs).\r\n",
      "## pis = number of loci containing n parsimony informative sites.\r\n",
      "## sumpis = sum of parsimony informative sites.\r\n",
      "\tnvar\tsumvar\tPIS\tsumPIS\r\n",
      "0\t1756\t0\t10432\t0\r\n",
      "1\t3512\t3512\t10233\t10233\r\n",
      "2\t4749\t13010\t7935\t26103\r\n",
      "3\t5195\t28595\t5574\t42825\r\n",
      "4\t5153\t49207\t3946\t58609\r\n",
      "5\t4704\t72727\t2609\t71654\r\n",
      "6\t4150\t97627\t1631\t81440\r\n",
      "7\t3631\t123044\t1072\t88944\r\n",
      "8\t3016\t147172\t602\t93760\r\n",
      "9\t2385\t168637\t390\t97270\r\n",
      "10\t1915\t187787\t220\t99470\r\n",
      "11\t1452\t203759\t126\t100856\r\n",
      "12\t1061\t216491\t47\t101420\r\n",
      "13\t747\t226202\t30\t101810\r\n",
      "14\t515\t233412\t16\t102034\r\n",
      "15\t356\t238752\t7\t102139\r\n",
      "16\t229\t242416\t2\t102171\r\n",
      "17\t142\t244830\t2\t102205\r\n",
      "18\t92\t246486\t0\t102205\r\n",
      "19\t54\t247512\t1\t102224\r\n",
      "20\t26\t248032\t0\t102224\r\n",
      "21\t19\t248431\t0\t102224\r\n",
      "22\t8\t248607\t0\t102224\r\n",
      "23\t5\t248722\t0\t102224\r\n",
      "24\t2\t248770\t0\t102224\r\n",
      "25\t0\t248770\t0\t102224\r\n",
      "26\t1\t248796\t0\t102224\r\n",
      "total var= 248796\r\n",
      "total pis= 102224\r\n",
      "sampled unlinked SNPs= 43119\r\n",
      "sampled unlinked bi-allelic SNPs= 37740\r\n"
     ]
    }
   ],
   "source": [
    "cat empirical_10/stats/empirical_10_m4.stats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "80709       ## loci with > minsp containing data\n",
      "80639       ## loci with > minsp containing data & paralogs removed\n",
      "80639       ## loci with > minsp containing data & paralogs removed & final filtering\n",
      "\n",
      "## number of loci recovered in final data set for each taxon.\n",
      "taxon\tnloci\n",
      "29154_superba         \t25218\n",
      "30556_thamno          \t38412\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head -n 10 empirical_10/stats/empirical_10_m2.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Infer ML phylogeny in _raxml_ as an unrooted tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## raxml argumement w/ ...\n",
    "raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \\\n",
    "                      -w /home/deren/Documents/RADmissing/empirical_10/ \\\n",
    "                      -n empirical_10_m4 -s empirical_10/outfiles/empirical_10_m4.phy\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## raxml argumement w/ ...\n",
    "raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \\\n",
    "                      -w /home/deren/Documents/RADmissing/empirical_10/ \\\n",
    "                      -n empirical_10_m2 -s empirical_10/outfiles/empirical_10_m2.phy\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "This is RAxML version 8.0.16 released by Alexandros Stamatakis on March 21 2014.\n",
      "\n",
      "With greatly appreciated code contributions by:\n",
      "Andre Aberer      (HITS)\n",
      "Simon Berger      (HITS)\n",
      "Alexey Kozlov     (HITS)\n",
      "Nick Pattengale   (Sandia)\n",
      "Wayne Pfeiffer    (SDSC)\n",
      "Akifumi S. Tanabe (NRIFS)\n",
      "David Dao         (KIT)\n",
      "Charlie Taylor    (UF)\n",
      "\n",
      "\n",
      "Alignment has 148988 distinct alignment patterns\n",
      "\n",
      "Proportion of gaps and completely undetermined characters in this alignment: 36.33%\n",
      "\n",
      "RAxML rapid bootstrapping and subsequent ML search\n"
     ]
    }
   ],
   "source": [
    "%%bash \n",
    "head -n 20 empirical_10/RAxML_info.empirical_10_m4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "This is RAxML version 8.0.16 released by Alexandros Stamatakis on March 21 2014.\n",
      "\n",
      "With greatly appreciated code contributions by:\n",
      "Andre Aberer      (HITS)\n",
      "Simon Berger      (HITS)\n",
      "Alexey Kozlov     (HITS)\n",
      "Nick Pattengale   (Sandia)\n",
      "Wayne Pfeiffer    (SDSC)\n",
      "Akifumi S. Tanabe (NRIFS)\n",
      "David Dao         (KIT)\n",
      "Charlie Taylor    (UF)\n",
      "\n",
      "\n",
      "Alignment has 114159 distinct alignment patterns\n",
      "\n",
      "Proportion of gaps and completely undetermined characters in this alignment: 56.87%\n",
      "\n",
      "RAxML rapid bootstrapping and subsequent ML search\n"
     ]
    }
   ],
   "source": [
    "%%bash \n",
    "head -n 20 empirical_10/RAxML_info.empirical_10_m2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the tree in R using `ape`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAMgCAIAAABUEpE/AAAgAElEQVR4nOzdeVxTV94/8JOwSEQr\nkLAIyiKLWFYJsqit4IJorWLdWnXcUGyro1SFp+KIWm0da61YQamIVltAoYWiMlJZNYKBJBSQUazF\nR8GlyiYiO8n5/XGf5scoUKdeSMHP+9U/br4559wvy7zm470nFw6llAAAAAAAe7iqbgAAAACgv0HA\nAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAs\nQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACAZQhYAAAA\nACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgA\nAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUI\nWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACA\nZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAA\nAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAEL\nAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAM\nAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAA\nsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsgJ517do1Gxub\nnp5eXV1tYGBAKVUe/OkzAgDAy0PAgr7kyZMnK1euNDQ0NDY2/uSTT5T1/Pz80aNH6+vrnzx5Ull8\n/PjxggULhg4d6uHhcezYMUJIXV3d+vXrzczMdHR0tm3b1s3Il9fS0jJz5kxKqVQqtbe3/9PrvOB0\nmUzm6urK4XCUB3/6jAAA8PIQsKAvWbx4sVwuLykpuXz58q5duxQKBSGkvr5+7ty5x48fj4mJ2bRp\nEzNSoVDMnDlz6tSp9+7d++KLL86dO0cIWbRokYGBgUwmE4vFu3fvZqZ3OvLlFRUVVVZWcjgcqVRq\nZ2f3p9d5welSqVQoFHY8AAAAFULAgj6jvb19xYoVERERcrk8JSVlzJgxXC6XEBIfH//22287OzuP\nHz++pqaGiU2ZmZkDBgxYsWIFl8sdP358YmIiISQnJ8fOzk5NTS0hIcHV1ZWZ3ulIJYlEMmXKFCMj\nI1NT09jYWIVCoaWlVVNTQwh5/PixmZnZvXv3rl69OnHiRIFAYG5uvmXLFkJIYWHhrFmzysvLQ0ND\nZTJZXV3dxIkTjYyMtm7dyiybn58/adIkfX39kSNHXrhwgRDS0NAgEAhCQkJ8fHyMjY1jYmKYkc9P\n77QH5sIV+f1SFiHk+a4AAKD3UIA+RSaTMb+6qampTMXHx+fMmTOU0rKyMgsLC6a4ePHilStX2tnZ\n8fn8Y8eOMcUNGzYof/MzMjK6GcmQSCRGRkapqant7e3nz58fNGgQpdTGxqaoqIhSGhwcvH37dkqp\nl5dXVFRUW1tbUVGRpqamXC6nlC5cuPDUqVNtbW08Hu/vf/97a2traWnpgAED5HK5VCo1NTXNzs5u\na2s7cuSIoaGhQqEQiUQ8Hi83N5dSev78eV1dXYVC0en0TnswNTWtqKjoeNBpVwAA0DvUVZTretDc\nuXOZKxPQn3C53H379pmYmIwePbq2tvb48eMrV66sqKgghNy8edPBwYEQUlxcbG1tzYwvKChwcnLK\nyspKS0v74IMPli5dGhkZWVhYWFpayufzv/vuuzVr1ly/fr3Tkcrfn5CQkM8//3zq1KmEEB8fH2a8\nubn53bt3hwwZ8v333xcXFxNCkpOTCwsLY2NjMzMzR48ezUzPycnZvXt3aWkpj8f74osvNDQ0NDQ0\nLCwsuFxuaGjonj17JkyYQAiZN2/eRx99RCmVSCT+/v6enp6EEE9Pz9bWVkJIp9Of76GysrKlpcXE\nxER50FVXAADQO/phwPrtt9867nSG/kFdXT0oKCg6OprH4+no6MyaNevgwYOEEIVCUVFRYWhoSAhJ\nTU2dP38+U7x161ZWVpa+vr6Pj8+AAQM4HE54ePjZs2ctLS0JIX5+fmFhYV2NVJ5UJpOdOHGCOa6p\nqeHxeIQQc3PzioqKuLi40NBQbW3t6OjoyMjIGTNmuLm5DR06lM/nE0Lu37/f3t4+fPjwEydOeHt7\na2pqEkKuXLnC5Ke8vLzo6GhmWZFIJBQKuVyuRCJhkhwhJCcnx8XFhdm/9fz053tgFmF2uDMHnXYF\nAAC9ph8GLA0NjREjRqi6C2BfWlpaXFzc0qVLf/vttzVr1qxbt44QwuFwXnvttZycHD09vXPnzu3d\nu5cp8ni8vLy8qVOnhoaGLl++nMPhVFdXi0QiCwuL8vLygICAoKCgrkYqzzho0KCkpKSAgIDi4uK3\n3norPj7+jTfeMDMzO3PmzMOHD5nsFRwcnJmZ6eDgIJVKT58+vWfPHkLIlStXlJ/pGzt2LLNabm4u\nk5C0tLQSExNXr14tlUrXr1/P/HtAIpG0t7fPnTu3vLw8MDAwIiKCENLpdELIMz08vwGr064AAKD3\nqPoeJfu8vLxU3QL0iISEBHNzc11dXScnp+joaIVCwdSPHz/O5/OdnJzy8vKUg2NjY42MjIYNGxYU\nFNTc3EwpTU5OtrKy0tPTEwqFcXFxyunPj1RKT08fOXLkoEGDbG1tv/32W6bIbD/PyspiXm7YsGHY\nsGEjR44MCQmZNm3a9OnTKaVisVggEISHh3t4eDDbqiilzs7OV69epZSmpqYynbzxxhuXLl2ilNbU\n1Nja2gYGBgoEAmdn56SkJGZKp9Of78HPzy85ObnjQaddAQBAr+HQfvdAQm9v76ysLFV3Af1WTEzM\n999/n5SUxOKaaWlpUVFR8fHxKuwBAABY1A9vEQL0nLt37+7cufPs2bPsLiuRSBwdHVXbAwAAsAgf\nLAJ4UceOHRs3btynn36q/KwiW6RSKfNBSBX2AAAALMItQgAAAACW4QoWAAAAAMsQsAAAAABYhoAF\nAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiG\ngAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAA\nWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAA\nAABYhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCw\nAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAqExIS8umnn6q6CwAAYB8CFvQx\nV65cGTt2rJ6e3ogRI86cOUMI+eWXX9566y0DAwNLS0umotTa2jp27NjS0lLmJaXU3t6e87vly5d3\nHJydnW1sbCyXy3vta5FKpXZ2dr12OgAA6DXqqm4A4L/w4MGDGTNmxMTE+Pj4nDhxYvXq1ePHj/f1\n9T106NCUKVPi4+OXLVtWVVXF5XLb29u3bt2ampp6584dGxsbZnpOTg6fz6+vr2deamhoKFe+dOnS\nggULXFxc1NTUeudroZQiYAEA9Fe4ggV9iUgkWrZsma+vL5fLdXZ2VlNTS0xMnDBhgq+vr5qamre3\nd0tLCzOyra1t3LhxkydPdnd353L/7/f866+/bmhoMDExcXBwEIlEAwYMYOo5OTl79+718/MbN27c\nM2fMzs52dHQ0MDDYsmULIeT27dvW1tbMWxKJxMPDgxDS0NAgEAhCQkJ8fHyMjY1jYmKU3U6YMEFf\nX9/CwuLkyZNMce3atdu2bXNxcVm3bl1DQ0NsbKyxsbGZmVlqaioz4OrVqxMnThQIBObm5sxJAQCg\n76H9jpeXl6pbgB538+ZNT0/P2NjYiIgIZ2fnioqKsrKyadOmCYXCjsOCgoJ27NjBHDc0NHh4eCQm\nJlZXV4eGhhobGzP1K1euzJo16+nTp25ubpmZmR2nKxQKIyOjixcvVlVVubu7t7W1xcXFLVmyhHk3\nLCxs/fr1lFKRSMTj8XJzcyml58+f19XVVSgUubm57u7uJSUlcrk8OztbS0tLLpdTSoVC4ZtvvimV\nSk+ePKmmpnby5MnW1tbDhw+bmJgwy3p5eUVFRbW1tRUVFWlqajKzAACgb+FQSlWd8Vg2fvx45dUC\n6E/Mzc25XG5NTU1oaGhBQUFkZKSjo2N9ff38+fMvXbrk5ORkbm6uqan5zTffKKe88cYbW7du9fHx\neWap2tpaExOTp0+fSqXSPXv2fPvttxwOh8/nV1ZWamtrK4dRSt3c3HR0dObPnz9v3jwdHZ3AwMBR\no0atXr2aELJgwYLZs2e/++67+/fvv3Xr1sGDBwkhdXV1JiYm9fX1U6dOTUtLUy7F5/OrqqoaGxv1\n9PTKyspMTEyCg4NramqOHj1KCHn8+PHQoUMbGhq4XO6TJ08KCwtv376dmZlZWloqFot78psKAAA9\noh/uwXrw4MGePXtU3QWwLyQkpKSkZM+ePR999NHBgwc5HA4hpK6uLiUlhdl0ZW9vHxYWphzf2toq\nk8nc3NyYl0lJSU5OTiNGjCCE/Prrry4uLlwud8uWLenp6YmJicyYQYMG1dbW6ujoMC85HI5YLM7M\nzIyJidm1a9ft27fFYrG/vz8h5PHjx2lpacxvmkQimTp1KjMlJyfHxcWFw+Hk5eU9fPjQwMCAENLU\n1MR0W1BQ4OzsbGJiQgiRSqVLlixhZuXl5bm6unK53Ojo6MjIyBkzZri5uQ0dOpTP5/f0dxUAAHqE\nqi+hsQ+3CPurzMxMV1fXysrK+vr6+vr6p0+fUkqHDBny448/Njc3b9y40dfXt+MNtby8vNdff135\ncv78+YsWLXr69OmdO3c8PDyeuRsYFhYWEBDwzBl1dHRSU1PlcnlcXJylpSWldODAgdeuXWtra1u1\napWhoaFCoaCUWllZzZs37+nTp9euXbO2tr5w4QKl1NHRcdeuXc3Nzfn5+SYmJjk5OZTSvXv3fvTR\nR5RSuVyuo6OzePHipqamoqIiGxub9PR0Sqmenl5hYaFcLs/Ly7OwsIiPj++Z7yUAAPQsbHKHPuPz\nzz+XSqX6+vqDBw8ePHiwr68vISQsLGzVqlUjRoxobGw8ffq0cj87ISQ3N9fT01P5cufOnb/++quJ\nick777yzfft2b2/vjovn5+crr3Up7dix43/+53/4fP7+/fvj4uIIIWvWrJkwYYKXl9eQIUM8PDw4\nHE5tba26urqJiYm5ufnChQs///zzKVOmEEK++eabxMREfX39VatWhYeHjx07lhAiFouZln799Veh\nUDh48ODhw4f7+/vv379/0qRJhJBly5bNmDHj9ddfT05OtrW17Xi7EwAA+pB+uAfL29s7KytL1V3A\nqyItLS0qKio+Pl7VjQAAwF8IrmABvBSJROLo6KjqLgAA4K8FAQvgpUilUgcHB1V3AQAAfy398FOE\nAL1J+QlEAAAAJVzBAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAEL\nAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAM\nAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAA\nsAwBCwAAAIBl6qpugE3vv/9+TU1NWVmZqhsBAACAV1q/uoKVJZZ6B2wYKDBUdSMAAADwSutXAUtj\nwADD4WbqGpqqbgQAAABeaf0qYAEAAAD8FSBgAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIW\nAAAAAMsQsAAAAABYhoAFAAAAwDIELAAAAACW9auApZC3P31Sp1DIVd0IAAAAvNI4lFJV98Ca2bNn\nP336tKSk5MGDB6ruBQAAAF5d/SpgMby9vbOyslTdBQAAALy6+tUtQgAAAIC/AgQsgJ517do1Gxub\nnp5eXV1tYGBAKVUevPgpkpOT58+f/6c7BACA5yFgQZ+UnZ1tbGwsl///DzS0traOHTu2tLS0m5EH\nDx7k/KfS0tK6urr169ebmZnp6Ohs27aNrQ5bWlpmzpxJKZVKpfb29n96nRecLpPJXF1dORyO8qDT\nYdnZ2Twer729vWOxpKTE1dX1T3cIAADPQ8CCvufSpUsLFixwcnJSU1MjhLS3t2/evNnd3b20tPSZ\niz3PjAwICKj/3Z49e95///2RI0cuWrTIwMBAJpOJxeLdu3crFApWmiwqKqqsrORwOFKp1M7O7k+v\n84LTpVKpUCjseNCpUaNGiUQidXX1jsWSkpJupgAAwJ+AgAV9TE5Ozt69e/38/MaNG8dU2traxo0b\nN3nyZHd3dy6X283IAQMGDBo0aNCgQbm5uVlZWV999RWHw8nJybGzs1NTU0tISHB1de24AiFEIpFM\nmTLFyMjI1NQ0NjZWoVBoaWnV1NQQQh4/fmxmZnbv3r2rV69OnDhRIBCYm5tv2bKFEFJYWDhr1qzy\n8vLQ0FCZTFZXVzdx4kQjI6OtW7cyy+bn50+aNElfX3/kyJEXLlwghDQ0NAgEgpCQEB8fH2Nj45iY\nGGbk89M77YG5cEV+v5TV6SkIITt37pRKpYSQM2fOuLi4CASCwMDAkpISFxcXQohIJJowYYK+vr6F\nhcXJkyeZKdnZ2Y6OjgYGBsyXBgAAL4T2O15eXqpuAXrKlStXZs2a9fTpUzc3t8zMzI5vBQUF7dix\n40VGXrt2zd3dvaamhnm5YcMG5f8cMjIyOo6USCRGRkapqant7e3nz58fNGgQpdTGxqaoqIhSGhwc\nvH37dkqpl5dXVFRUW1tbUVGRpqamXC6nlC5cuPDUqVNtbW08Hu/vf/97a2traWnpgAED5HK5VCo1\nNTXNzs5ua2s7cuSIoaGhQqEQiUQ8Hi83N5dSev78eV1dXYVC0en0TnswNTWtqKhQHnR6CkqpUCgs\nKChIT08fMWJEcXFxfX393LlzraysKKW5ubnu7u4lJSVyuTw7O1tLS0sulysUCiMjo4sXL1ZVVbm7\nu7e1tbH7AwUA6K8QsEDFgoKC/Pz85v2RBQsWnDt37p133mloaGhsbOTxeE+fPu24zvjx43/66Sfm\nOC8vr6uRNTU19vb2N27cYF5GRERMnDixtLS0srJy//79tra2HdecMmXKyZMnmWO5XM4kGB8fn5SU\nlNu3b48YMYJZua6u7uLFiydOnFi6dKm7uzsz3szM7M6dO1evXtXT02tpaaGUlpWVMetPnz49Li6O\nGVZbW6utrS2Xy7/88su1a9cyxcePH2traysUik6nP9/Do0ePmAilPOj0FA0NDa+99lpra6uXl9eP\nP/7IvHv8+PEFCxYwX2zHf3rx+XxKqUKhcHV1nTx58pEjR2pra1/mBw0A8EpRf+6SFkCvqqqqWrt2\nrYWFRffD1NXVV6xYkZGRkZiYyFQGDRpUW1uro6NDCGltbZXJZG5ubsxbW7ZsSU9P73Tkjh07/va3\nvym3aoWHh589e9bS0pIQ4ufnFxYW1vGkMpnsxIkTzHFNTQ2PxyOEmJubV1RUxMXFhYaGamtrR0dH\nR0ZGzpgxw83NbejQoXw+nxBy//799vb24cOHnzhxwtvbW1NTkxBy5coVT09PQkheXl50dDSzrEgk\nEgqFXC5XIpFMnTqVKebk5Li4uDD7t56f/nwPzCLMDnfmoNNTFBQUODk5aWhoFBYWenh4MO+KxWJm\nA1ZeXt7Dhw8NDAwIIU1NTcw2eQ6HIxaLMzMzY2Jidu3adfv27a62zwMAQEfYgwUqpqGhMXz48BF/\nxNTUND09nflnQVhYWEBAAKWUyUyEkMLCQgsLC+XLtLS0TkeWlpb+9NNPH330kfLs1dXVIpFIoVDc\nvn07ICAgKCioY2+DBg1KSkpqb28vKChwcHC4du0aIcTMzOzMmTPXr1//29/+RggJDg4+evTo1q1b\n+Xz+6dOnmeBy5coV5Wf6xo4dy6yWm5vLJCQtLa3ExES5XJ6Xl7d+/fpPP/2UECKRSFJSUhoaGq5f\nvx4YGMhst+p0+vM9PL8Bq9NTiMVid3d3QgiPx2O+rvj4+OPHjzNTzM3No6KiWlpaJBKJtbV1QUEB\nIURXVzc9PX3SpEm+vr4aGhpIVwAAL0pVl856Dm4R9i0BAQHKG3YvaOHChUePHu1Y2b9/v7+//x+O\nfOuttw4fPtxxQHJyspWVlZ6enlAojIuLYzYqKaWnp48cOXLQoEG2trbffvstU2S2n2dlZTEvN2zY\nMGzYsJEjR4aEhEybNm369OmUUrFYLBAIwsPDPTw8mG1VlFJnZ+erV69SSlNTU5mTvvHGG5cuXaKU\n1tTU2NraBgYGCgQCZ2fnpKQkZkqn05/vwc/PLzk5uePB86eglM6ZMychIYFSmpCQYGhoaGFhsWvX\nLn19/cePH1NKCwoKXFxcBg8e7OTkpGzgwIEDTk5OOjo6bm5u+fn5f/STAQCA/4M/lQMqtnr16o0b\nN77Mozh7WUxMzPfff5+UlMTimmlpaVFRUfHx8SrsAQAAWIQ9WAD/hbt37+7cufPs2bPsLiuRSBwd\nHVXbAwAAsAh7sABe1LFjx8aNG/fpp59aW1uzu7JUKnVwcFBtDwAAwCLcIgQV63O3CAEAAP4QrmAB\nAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYh\nYAEAAACwDAELAAAAgGUIWAAAAAAsU1d1A9Bv3b9//6uvvvrDYfn5+V988YWenl5XA8aOHTtz5kxW\nWwMAAOhZCFjQU65fv15458GkeQu7H7bQzbubd5/UVJ8+fRoBCwAA+hYELOhBugaGlnaOL7NCzaOH\nt7JS2OoHAACgd2APFgAAAADLELCgl7S1tn4wyf1x1SNl5XF15fbl85e6jfrnh8sa6590UwQAAOhb\nELCgN6ScPLp18exH9yo6Fk/8c4eZzaijl4sG8AZ+f/hAN0UAAIC+BQELeoPZyFFzPwzsWKGU5mf8\n5OHzloampvuUaZKsC10VAQAA+hxscv8Dc+bMqaysHDBggKob6Xtqampsxk9iju3dxz3zbnNjQ3Nj\ng8kIK0LIMEvr2kcPuyoCAAD0OQhYf4DP5//zn/+0trZWdSN9T0ZGxtEzP3U/hsPhEEIopQq5vPsi\nAABAH4JbhKAaWgO1tQZq3//fMkLIg9u3+EZDuyoCAAD0OQhYoBocDmfMRB/ZxQxK6c+irDGTfLsq\nAgAA9DkIWKAyyz7eduNnaYCX8ElN9bwP1ndTBAAA6FuwBwt6zw+l9zu+1BEY7DiR8MyYZ4pNjY29\n0RkAAACrELCgp/D5/LsFVw6v/bn7YdevXzczMxs4cGCn7yoUijfffLMHugMAAOhBCFjQU5ydnUUi\n0R8OW7169caNG21sbHqhJQAAgN6BPVgAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAAAFiGgAUAAADA\nMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAAAAAAWIaABdCz\nrl279jJ/yvoFp1dXVxsYGFBKlQcvfork5OT58+f/6Q4BAOB5CFjQJ2VnZxsbG8vlcmWltbV17Nix\npaWlzMtffvnlrbfeMjAwsLS0PHPmDFOklNrb23N+t3z58u7XfBktLS0zZ86klEqlUnt7+z+9zgtO\nl8lkrq6uHA5HedDpsOzsbB6P197e3rFYUlLi6ur6pzsEAIDnqau6AYD/2qVLlxYsWODi4qKmpkYI\naW9v37p1a2pq6p07d5iLPTU1Nb6+vocOHZoyZUp8fPyyZcuqqqq4XG5OTg6fz6+vr2fW0dDQ6GrN\nl1dUVFRZWcnhcKRSqZ2d3Z9e5wWnS6VSoVDY8aBTo0aNEolE6ur/8T/8kpKSlStX/ukOAQDgebiC\nBX1MTk7O3r17/fz8xo0bx1Ta2trGjRs3efJkd3d3LpdLCElMTJwwYYKvr6+ampq3t3dLSwsz8uuv\nv25oaDAxMXFwcBCJRAMGDOhqTSWJRDJlyhQjIyNTU9PY2FiFQqGlpVVTU0MIefz4sZmZ2b17965e\nvTpx4kSBQGBubr5lyxZCSGFh4axZs8rLy0NDQ2UyWV1d3cSJE42MjLZu3cosm5+fP2nSJH19/ZEj\nR164cIEQ0tDQIBAIQkJCfHx8jI2NY2JimJHPT++0B+bCFfn9UlanpyCE7Ny5UyqVEkLOnDnj4uIi\nEAgCAwNLSkpcXFwIISKRaMKECfr6+hYWFidPnmSmZGdnOzo6GhgYMF8aAAC8ENrveHl5sbjaqlWr\nfvnlFxYXhGcEBATcuHHjBQdfuXJl1qxZT58+dXNzy8zM7PhWUFDQjh07mOOIiAhnZ+eKioqysrJp\n06YJhUJKaUNDg4eHR2JiYnV1dWhoqLGx8R+uKZFIjIyMUlNT29vbz58/P2jQIEqpjY1NUVERpTQ4\nOHj79u2UUi8vr6ioqLa2tqKiIk1NTblcTilduHDhqVOn2traeDze3//+99bW1tLS0gEDBsjlcqlU\nampqmp2d3dbWduTIEUNDQ4VCIRKJeDxebm4upfT8+fO6uroKhaLT6Z32YGpqWlFRoTzo9BSUUqFQ\nWFBQkJ6ePmLEiOLi4vr6+rlz51pZWVFKc/OB+iMAACAASURBVHNz3d3dS0pK5HJ5dna2lpaWXC5X\nKBRGRkYXL16sqqpyd3dva2v7Uz9kAIBXDgLWH0DA6mnLly9PT08vewFXrlx55513GhoaGhsbeTze\n06dPO64zfvz4n376iTl+8uSJr6/vwIEDPT0933vvvaVLlz5z0pqaGh6PJ5fL8/LyullzypQpJ0+e\nZI7lcjmTYHx8fFJSUm7fvj1ixAhmfF1d3cWLF0+cOLF06VJ3d3dmvJmZ2Z07d65evaqnp9fS0kIp\nLSsrs7W1pZROnz49Li6OGVZbW6utrS2Xy7/88su1a9cyxcePH2traysUik6nP9/Do0ePmAilPOj0\nFA0NDa+99lpra6uXl9ePP/7IvHv8+PEFCxYwX2zHf3rx+XxKqUKhcHV1nTx58pEjR2pra//0TxkA\n4FWDPVigYjo6Ot9++63ybl03rl27dvny5cTERObloEGDamtrdXR0CCGtra0ymczNzY15q66uLiUl\nhcvltre329vbh4WFEUKSkpKcnJxGjBhBCPn1119dXFy4XO6WLVvS09M7XZMQIpPJTpw4wRwzmYwQ\nYm5uXlFRERcXFxoaqq2tHR0dHRkZOWPGDDc3t6FDh/L5fELI/fv329vbhw8ffuLECW9vb01NTULI\nlStXPD09CSF5eXnR0dHMsiKRSCgUcrlciUQydepUppiTk+Pi4sLs33p++vM9MIswO9yZg05PUVBQ\n4OTkpKGhUVhY6OHhwbwrFouZPVt5eXkPHz40MDAghDQ1NTHb5DkcjlgszszMjImJ2bVr1+3bt7va\nPg8AAP9B1QmPfbiC1e+FhYUFBAR0rOTl5b3++uvKl0OGDPnxxx+bm5s3btzo6+vL3FabP3/+okWL\nnj59eufOHQ8Pj2fuBj6/JqXU1NQ0IiKira1NJpMZGRldunSJUvrpp59Onz5dKBQyy+rp6RUWFjIX\nwywsLOLj4yml33///axZsyila9eu3bdvH7Pahx9+eOTIEUqpiYlJREREe3u7WCy2sLAQiUSUUisr\nq3nz5j19+vTatWvW1tYXLlzoavrzPezatSs0NLTjQaen2Lt376ZNmyilQ4cOPXz4cFtb2+nTpzU1\nNZlvhaOj465du5qbm/Pz801MTHJyciilOjo6qampcrk8Li7O0tLyJX9wAACvDgSsP4CA9Re0cOHC\no0ePdqzs37/f399f+fL48eP6+vrGxsYffPBBXV0dU7xx44a7u/uQIUOEQmFqauofrkkpTU9PHzly\n5KBBg2xtbb/99lumyGw/z8rKYl5u2LBh2LBhI0eODAkJmTZt2vTp0ymlYrFYIBCEh4d7eHgw26oo\npc7OzlevXqWUpqamWllZ6enpvfHGG0xoq6mpsbW1DQwMFAgEzs7OSUlJzJROpz/fg5+fX3JycseD\n509BKZ0zZ05CQgKlNCEhwdDQ0MLCYteuXfr6+o8fP6aUFhQUuLi4DB482MnJSdnAgQMHnJycdHR0\n3Nzc8vPzX+znAwAAlEP/mwcS9gne3t5ZWVlsrRYQEBAUFGRtbc3WgtDXxcTEfP/990lJSSyumZaW\nFhUVFR8fr8IeAACARdiDBfBfuHv37s6dO8+ePcvushKJxNHRUbU9AAAAi/AcLIAXdezYsXHjxn36\n6aesX9GUSqUODg6q7QEAAFiEW4R/ALcIAQAA4L+FK1gAAAAALEPAAgAAAGAZAhYAAAAAyxCwAAAA\nAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAAAADLELAA\nAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQsAAACAZQhYAAAAACxDwAIAAABgGQIWAAAAAMsQ\nsAAAAABYhoAFAAAAwDIELAAAAACWIWABAAAAsAwBCwAAAIBlCFgAAAAALEPAAgAAAGAZAhYAAAAA\nyxCwAAAAAFiGgAUAAADAMgQsAAAAAJYhYAEAAACwDAELAAAAgGUIWAAAAAAsQ8ACAAAAYBkCFgAA\nAADLELAAAAAAWIaABQAAAMAyBCwAAAAAliFgAQAAALAMAQvgr6i6utrAwIBS+vxbxcXFo0eP7v2W\nAADgxSFgQV/y5MmTlStXGhoaGhsbf/LJJ90UKaX29vac3y1fvpwQUldXt379ejMzMx0dnW3btjEj\nOy2qnEwmc3V15XA4z78llUqfCVjZ2dk8Hq+9vV150FttAgBA59RV3QDAf2Hx4sV8Pr+kpKS+vt7W\n1vYf//gHl8vttJiTk8Pn8+vr65mJGhoahJBFixZ5enrKZLKqqipHR8dt27ZxudxOiyr9KgkhRCqV\nCoXCrt5ycXHpWBk1apRIJFJXV1ce9EqPAADQJdX/HwnAC2pvb1+xYkVERIRcLk9JSRkzZgyXy+20\nSAj5+uuvGxoaTExMHBwcRCLRgAEDCCE5OTl2dnZqamoJCQmurq7MyE6LjIaGBoFAEBIS4uPjY2xs\nHBMTw9TXrl27bds2FxeX/fv3L1myxMrKysrKSiAQvP7664QQkUg0YcIEfX19CwuLkydPMlNsbGwi\nIyMJIf7+/sx1qaamJlNT07KysqtXr06cOFEgEJibm2/ZsoUZz1zBIoRkZ2c7OjoaGBgo32KuYLW0\ntLz//vvLli1ramrauXOnVColhCgPAABAxWi/4+XlxeJqq1at+uWXX1hcEF6STCZjfnVTU1O7KjY0\nNHh4eCQmJlZXV4eGhhobGzPDNmzYoPzNz8jI6KbIEIlEPB4vNzeXUnr+/HldXV2FQkEpFQqFb775\nplQqbW5uZkaWlJS4ubmVlJTk5ua6u7uXlJTI5fLs7GwtLS25XE4p9fT0PHz4cHl5ORPCKKVhYWEB\nAQGUUi8vr6ioqLa2tqKiIk1NTWa8qalpRUWFQqEwMjK6ePFiVVWVu7t7W1tbS0sLj8e7fv36+PHj\nw8PDlf0UFBR0PAAAANXi0M520fZp3t7eWVlZbK0WEBAQFBRkbW3N1oLwjI8//vjhw4eampp/OHLr\n1q3Dhg2jlNbV1R0/fvzLL7+sqKhg3uq0yKitrTUxMXn69GlkZOQPP/xw6NAhPp//3Xffff3119ev\nXz906NDzReXc/fv337p16+DBg4SQuro6ExOT+vr6pqYmPT29srIyExMTZlheXt66detOnTplYWHh\n4+OTlpamXIHP51dVVRFCZs6cOX369OvXrzs5Ofn7+zc1Nb3++usXL14cPnz4kydPCgsLb9++nZmZ\nWVpaKhaLKysrHRwcHjx4QAhxc3PT0dGZP3/+vHnzdHR0CgoK3njjjREjRrS0tNy4cYPD4TQ2Ng4d\nOrSqqqqtrY05YO6HAgCACmGvBqjYw4cPFy9ebGFh8Ycjg4KCjh07xuPxdHR0Zs2axeSeRYsWHT16\n9JliUlKSk5PTiBEjCCG//vqri4sLl8sNDw8/e/aspaUlIcTPzy8sLIwQ0mlRSSKRTJ06lTnOyclx\ncXHhcDgFBQXOzs7KdJWWlrZ9+/bk5GQjIyNCSF5e3sOHDw0MDAghTU1Nyl3qAoGgsrLywoULe/fu\n9ff3P3r06MyZM4cPHx4dHR0ZGTljxgw3N7ehQ4fy+XxCiEwmEwqFzFyxWJyZmRkTE7Nr167bt29L\npVI7O7t//etf3t7eGRkZkydPLigocHJy0tDQyMvLYw5Y+KkAAMDLeXUD1tWrVz/88EMtLa3uh12/\nfv3f//73wIEDO31XoVBMmjQpJCSkBxp8VWhqag4fPpwJQ91LT0+Pi4tbunTpb7/9tmbNmnXr1hFC\nLly48Hzx1KlTP/zww9dff11dXb1u3brPPvuMEFJdXS0SiSwsLMrLy5kLk10VlSQSSXt7+9y5c8vL\nywMDAyMiIgghYrF47NixzICEhITDhw+fO3dOV1eXqZibm0dFRW3atKm4uHj27Nnx8fHMYIFAEBkZ\nuXnzZuZaXXh4+MWLFwkhwcHBmZmZDg4OUqn09OnTe/bsIR02YOnq6p46dWrKlCnV1dWXL1/mcDhS\nqXT58uUCgWD16tURERGTJ08Wi8Xu7u5MY8wBAAConopvUfaAF9yDlZaW9t764B9K77/Mf1EXCxYt\nWtTTX1H/FhAQcOPGjRcZmZCQYG5urqur6+TkFB0dzWw/6rR448YNd3f3IUOGCIVC5Vat5ORkKysr\nPT09oVAYFxfHjOy0yKipqbG1tQ0MDBQIBM7OzklJSUx9zpw58fHxlNJ79+5xOBwTExNLS0tLS8vx\n48dTSgsKClxcXAYPHuzk5KScQinds2ePgYFBY2MjpZQQsnnzZqa+YcOGYcOGjRw5MiQkZNq0adOn\nT6eU+vn5JScnU0oPHDjg5OSko6Pj5uaWn59PKR09evTPP/9MKa2trR08ePCdO3fmzJmTkJDANMYc\nAACAyr26e7DS09OPnUub+0Hgy5yr5uFv2Uf2fffddy+zyCtu9erVGzdutLGxUXUjz0pLS4uKioqP\nj1d1IwAA0PfgMQ0AnZNIJI6OjqruAgAA+iQErP/T1tr6wST3x1WPlJXH1ZXbl89f6jbqnx8ua6x/\n0k0R+iWpVOrg4KDqLgAAoE9CwCKEkJSTR7cunv3o3n98vP/EP3eY2Yw6erloAG/g94cPdFOEfikx\nMXHWrFmq7gIAAPqkfvgpwpaWFuVjJ7tx8+ZN5bHZyFGGpma731+qrFBK8zN++kdUjIampvuUaXEH\nPl8SvPX54oylq3rkawAAAIC+rB8GLHNz8yNHjvzhsHv37r1mZccc27uPe+bd5saG5sYGkxFWhJBh\nlta1jx52VQQAAAB4Rj8MWLGxsS8yjPkUYfdjmCc9UkoVcnn3RQAAAAAl7MHqnNZAba2B2vf/t4wQ\n8uD2Lb7R0K6KAAAA8NdUXV1tYGDQ6ROpiouLR48e3XOnRsDqHIfDGTPRR3Yxg1L6syhrzCTfrooA\nAAB9y5MnT1auXGloaGhsbPzJJ590U6SU2tvbc363fPlyQkhdXd369evNzMx0dHS2bdvGjOy0qHLM\nH8ZQ/tWyjqRSKQKWaiz7eNuNn6UBXsInNdXzPljfTREAAKAPWbx4sVwuLykpuXz58q5duxQKRVfF\nnJwcPp9f/7vIyEhCyKJFiwwMDGQymVgs3r17NzOy06LKSaVSoVDY1VsuLi49d+pXOmDVPHpY9u9i\n5X+f/5Ba/fA35cvqh78tCd768aFv5n740f07/9tp8faNa6r+IgAAAP4L7e3tK1asiIiIkMvlKSkp\nY8aM4XK5nRYJIV9//XVDQ4OJiYmDg4NIJBowYAAhJCcnx87OTk1NLSEhwdXVlRnZaZHR0NAgEAhC\nQkJ8fHyMjY1jYmKY+tq1a7dt2+bi4rJ///4lS5ZYWVlZWVkJBILXX3+dECISiSZMmKCvr29hYXHy\n5Elmio2NDRPy/P39metSTU1NpqamZWVlV69enThxokAgMDc337JlCzNe+adds7OzHR0dDQwMlG8x\nV7BaWlp66PvcDze5vyBbW1tnU6Pq/D/4ozqpqanu7u7KP+X7vLlz57LdGgAAQE9RV1f38/MrKChg\nLu2kpqZ2VWxsbPz111+3bNkyYcKEAwcOrFy58t69e4SQFStWzJ49m1ktIyODOei0yPj5558bGxvf\nfvvtzz77LDU1deHChQsXLuRwOGKxWFtbOyoqyt7enolu//73v1esWHHs2LErV64EBQVFR0ePGjVK\nJBL5+vouXryYy+UKBAJCSEVFhVgsZhY/cuTItGnTLC0tvb29Fy1adOHChWvXro0ZM2bnzp1cLlcq\nlR44cIBS+t57750+fdrOzu6tt97asWOHQqEoKSnh8/mTJ08WiUQ98o1W1R9B7CtWrVr1yy+/qLqL\n/uzF/9gzAACwRaFQ1NbWfvnll8OGDeu+yKipqeHxeHK5PCIiYuLEiaWlpZWVlfv377e1taWUdlpU\n+vLLL9euXcscP378WFtbW6FQNDQ0DBgw4O7du8phYrHYzc3t1q1blNIpU6Z0zCp8Pp8Z8/bbbx8+\nfHjdunXR0dGEkKamJgsLi/LyckppXV3dxYsXT5w4sXTpUnd3d0rpo0ePDA0NFQqFQqFwdXWdPHny\nkSNHamtrKaUymWzgwIH29vbW1tZsf2v/zyt9ixAAAOBVs2jRoqamJg6Ho6OjM2vWLA0Nja6KSUlJ\nt27dYmb9+uuvLi4uXC43PDz8yJEjI0eOFAgEfn5+TU1NhJBOi0oSiYS5T0cIycnJcXFx4XA4BQUF\nzs7OJiYmTD0tLW3Dhg3JyckWFhaEkLy8vIcPHzJJpbGx8e7du8wwgUBQWVl54cKFxYsXE0KOHj06\nc+bM4cOHR0dHT5o0KSsrS19ff+jQoePGjSOEyGQyoVDIbM8Xi8XBwcE5OTlOTk6UUqlUamdnl5WV\nxVw56wkIWAAAAK+QCxcuxMXFyeXye/furVmzZt26dV0VT506FRoa2tDQUF5evm7dup07dxJCqqur\nRSKRQqG4fft2QEBAUFBQV0UliUSSkpLS0NBw/fr1wMDArVu3EkLEYvHYsWOZAQkJCbt37z537pyR\nkRFTMTc3j4qKamlpkUgk1tbWBQUFTF0gEERGRq5Zs0ZTU5MQEh4evnnzZkJIcHDw0aNHt27dyufz\nT58+7eHhQTpswNLV1U1PT580aZKvr6+GhgaHw5FKpcuXLxcIBKtXr+6h7/OruwcLetr9+/e/+uqr\nPxyWn5//xRdf6OnpdTVg7NixM2fOZLU1AIBX1+HDh4OCgjZt2mRqarpu3TrmyQudFnfu3LlkyRIT\nExMrK6tPP/3U29ubEBIVFbVx48aNGzdaWFhs2rRpwYIFXRUZtbW16urqJiYm5ubmw4YN+/zzz5nb\nf2KxmBl2//79BQsWGBsbjxkzhhAydOhQkUj0zTffrFy5cs+ePSNGjAgPD1dGMYFA0N7e7u/vz7x8\n5513DA0NCSHLli2bMWOGtrb2nDlzbG1tv/nmm3nz5jEpihCyY8eO//mf/3n33XdtbGzi4uIIIVKp\n9MMPPySEMFfCegKHdvb0LVBikri1tbWqG+l7MjIy9h49OWnewpdZ5ElN9a2sFOVHTgAAoJc9efJk\nw4YNZ8+eVVNTe//990NDQ5n6lStXNm7cWFpaqqOjExYWxvxLOD8/f/Xq1Xfv3t23b9+SJUsIIT/+\n+GNwcHBLSwvzoKwdO3Z0M70/wRUs6EG6BoaWdo4vs0LNo4e3slLY6gcAAP5bixcv5vP5JSUl9fX1\ntra2//jHP7hc7oMHD2bMmBETE+Pj43PixInVq1fPnDmzvr5+7ty5Z86cefTo0eLFi5mA9fHHH5uZ\nmcXFxVVVVTk6Om7btq2r6ar+QlmGPVgAAADQua6ejyUSiZYtW+br68vlcp2dndXU1Agh8fHxb7/9\ntrOz8/jx42tqaphnjd66dWvMmDHPPB+r0+mM3n9oVk9973ro04n9Bh7T8Kelp6e/uy7oh9L7zH+n\nim8bmAyPvlyorETnFDl4jh/02pAxE32+lZR2Woy69PPChQtV/aUAALwqmpub29raninKZDImM6Sm\npj7z1s2bNz09PWNjYymlPj4+Z86coZSWlZVZWFgwAzZs2KCMHBkZGd1MZ4hEIh6Pl5ubSyk9f/68\nrq6uQqGglAqFwjfffFMqlTY3NzMjS0pK3NzcSkpKcnNz3d3dS0pK5HJ5dna2lpaWXC6nlHp6eh4+\nfLi8vJwJYZTSsLCwgIAASqmXl1dUVFRbW1tRUREr37fn4QoW9IaUk0e3Lp796F5Fx+KJf+4wsxl1\n9HLRAN7A7w8f6KYIAAC947333ispKXmmOHr0aOb5WCtXrlQWa2pq1q5du2TJksjIyPfee48QcvPm\nTQcHB0JIcXExs3f50KFDhYWFyudjrVmzppvpDIlE4u/v7+npSQjx9PRsbW0lhDQ2NpaUlMTGxgqF\nQubBCnl5eStWrDh16pSdnd22bdvy8vLs7e3V1NS8vLy0tbWZ62TMU0m/+OKLjRs3EkKam5sPHDjw\nj3/8gxCSnJxsY2MTGxv75Zdf9sw3ErcIoVeYjRw198PAjhVKaX7GTx4+b2loarpPmSbJutBVEQAA\neseRI0esra2dnZ2VlU6fj0UISUlJ8fPzmzRpUk5OjqOjIyFEoVBUVFQwn+lLTU2dP38+6fr5WM9P\nV+r9h2b10DcTm9x71b59+2JiYvh8vqob6Q01NTU24ycxx/bu4555t7mxobmxwWSEFSFkmKV17aOH\nXRUBAKAX3Lx585tvvsnK+o+/IMc8H2vp0qW//fab8vlYWVlZ27dvP3/+vJaWVkNDA4fD0dbW5nA4\nr732Wk5Ojp6e3rlz5/bu3Ut+fz4W87B15fOxOp2uPKNEImlvb587d255eXlgYGBERAR57qFZhw8f\nPnfunPKv2DEPzdq0aVNxcfHs2bPj4+OZwcxDszZv3qx8aNbFixcJIcHBwZmZmQ4ODlKp9PTp07t3\n7+6J7ycCVq/S1dUNDAxkPljR72VkZBw981P3Y5hdh5RShVzefREAAHqOXC5fvXp1VFRUxyeb379/\nf+PGjTt37nzm+Viff/65VCrV19dnho0fP14kEnE4nH379r377rvDhg1LTEwcPHgw6eL5WJ1OZ45V\n8tCsHvqWImCBamgN1NYaqH3/f8tGjnZ9cPsW32hoV0UAAOhpe/funTp1qp2dHfOypqZmz549+fn5\nn3322ccff/zM4PPnz3e6yLJly5YtW9axMnPmzOefv9DVdEKIVCp1cHDYv3///v37O9a///575sDY\n2Jj5cGJHo0ePVm7D7yg4ODg4OJg5ph2e+rlv3759+/Z11QNbsAcLVIPD4YyZ6CO7mEEp/VmUNWaS\nb1dFAADoUVKpND09nbl/V1tb+/HHH8+ePXvKlClZWVnMZvNeI5FIntmS1XchYIHKLPt4242fpQFe\nwic11fM+WN9NEQAAekhzc/PatWujoqKePn368ccfv/322+PGjcvOzp48eXLvN8Ncwer98/YE3CKE\n3vND6f2OL3UEBjtOJDwz5pliU2Njb3QGAPCqCg0NXbJkSXx8fFJS0saNGz/77DPmGQcqkZiYqKpT\nsw4BC3oKn8+/W3Dl8Nqfux92/fp1MzOzgQMHdvquQqF48803e6A7AAAgGRkZFy5c0NLSWrFihUgk\nUj6FAV4eAhb0FGdnZ+UHQ7qxevXqjRs32tjY9EJLAACgVFFR4efnt3Llyu3btw8ZMkTV7fQ3CFgA\nAACvHKlU+uGHH4rFYuUnB4FdCFgAAACvlpSUlN27d//444/Gxsaq7qXfQsACAAB4hXz99ddnzpw5\nf/488yxQ6CEIWAAAAK8ESunmzZsrKiqSkpKYvx4DPQfPwQIAAOj/2traVqxYoaWl9d133yFd9QIE\nLAAAgH7uyZMnb7/9touLy/bt25m/9wo9DbcIAQAA+rOHDx/Omzdvw4YNfn5+qu7lFYKABQAA0G/d\nuHFjyZIlX331lbu7u6p7ebXgFiEAAED/dPny5SVLlsTGxiJd9T5cwQIAAOiHEhISwsLCkpOTjYyM\nVN3LqwgBCwAAoL85dOhQamrqhQsXtLW1Vd3LKwoBCwAAoP9QKBQfffRRc3NzYmKiujr+X15lsAcL\n4K+ourrawMCAUvr8W8XFxaNHj+79lgDgr6+1tXXZsmW6urqRkZFIV6qFgAV9yZMnT1auXGloaGhs\nbPzJJ58o61euXBk7dqyent6IESPOnDnDFPPz80ePHq2vr3/y5EmmUldXt379ejMzMx0dnW3btnU/\nXbVkMpmrq2unj6uRSqXPBKzs7Gwej9fe3q486K02AeAvpK6ububMmd7e3njY1V8BAhb0JYsXL5bL\n5SUlJZcvX961a5dCoSCEPHjwYMaMGaGhoVVVVVu3bl29ejUhpL6+fu7cucePH4+Jidm0aRMzfdGi\nRQYGBjKZTCwW7969u5vpKieVSoVCYVdvubi4dKyMGjVKJBKpq6srD3qlRwD4C7l79+60adM2bNiw\nfPlyVfcChCBgQR/S3t6+YsWKiIgIuVyekpIyZswYLpdLCBGJRMuWLfP19eVyuc7OzmpqaoSQ+Pj4\nt99+29nZefz48TU1NUyWysnJsbOzU1NTS0hIcHV17WY6o6GhQSAQhISE+Pj4GBsbx8TEMPW1a9du\n27bNxcVl//79S5YssbKysrKyEggEr7/+OrPghAkT9PX1LSwslBfPbGxsIiMjCSH+/v7MvyybmppM\nTU3LysquXr06ceJEgUBgbm6+ZcsWZjxzBYsQkp2d7ejoaGBgoHyLuYLV0tLy/vvvL1u2rKmpaefO\nnVKplBCiPACAV0pxcfHMmTMPHDjg4+Oj6l7gdxS6tWrVql9++YWt1aKjo0+cOMHWav1DQEDAjRs3\nXny8TCZjfnVTU1OfeevmzZuenp6xsbGUUh8fnzNnzlBKy8rKLCwsmAEbNmxQ/uZnZGR0M50hEol4\nPF5ubi6l9Pz587q6ugqFglIqFArffPNNqVTa3NzMjCwpKXFzcyspKcnNzXV3dy8pKZHL5dnZ2Vpa\nWnK5nFLq6el5+PDh8vJyJoRRSsPCwgICAiilXl5eUVFRbW1tRUVFmpqazHhTU9OKigqFQmFkZHTx\n4sWqqip3d/e2traWlhYej3f9+vXx48eHh4cr+ykoKOh4AACvjosXL3p4eNy6dUvVjcB/wK0EULH2\n9vaKioo/vKvF5XLNzc0JIaNHj66trT1+/PjKlSsrKiqYd2tqakJDQwsKCiIjIx0dHQkhN2/edHBw\nIIQUFxdbW1sTQg4dOlRYWFhaWsrn87/77rs1a9Zcv369q+kMiUTi7+/v6elJCPH09GxtbSWENDY2\nlpSUlJWVmZiYMMPy8vLWrVt36tQpCwsLHx+fvLw8e3t75i0+n89cJxMIBISQL774YuPGjf7+/s3N\nzQcOHLh48SIhJDk5ubCwMDY2NjMzc/To0Vwut7KysqWlhVl/2LBhO3funD9/fmpqqrq6ekFBAYfD\nmTdvXktLy4cffsjhcBobG2/evGlvRWQBrAAAIABJREFUb688YOPHAgB9w6lTp6KiolJSUvT09FTd\nC/wHBCxQsUGDBh08ePAP/7S7urp6U1NTbGwsj8fT0dGZNWvWwYMHmbdSUlL27Nnz0UcfHTx4kLn7\nplAoKioqDA0NCSGpqanz588nhISHh589e9bS0pIQ4ufnFxYW1tV0JYlEMnXqVOY4JyfHxcWFw+EU\nFBQ4Ozsr01VaWtr27duVj/LLy8t7+PChgYEBIaSpqUm5oEAgqKysvHDhwt69e/39/Y8ePTpz5szh\nw4dHR0dHRkbOmDHDzc1t6NChfD6fECKTyYRCITNXLBZnZmbGxMTs2rXr9u3bUqnUzs7uX//6l7e3\nd0ZGxuTJkwsKCpycnDQ0NPLy8piDl/2RAEAf8fnnn8tkspSUFC0tLVX3As9CwAIVO3DgwAuO1NfX\nj4uLW7p06W+//bZmzZp169YRQrKysrZv337+/HktLa2GhgYOh6Otrc3hcF577bWcnBw9Pb1z587t\n3buXEFJdXS0SiSwsLMrLywMCAoKCgrqarjyjRCJpb2+fO3dueXl5YGBgREQEIUQsFo8dO5YZkJCQ\ncPjw4XPnzunq6jIVc3PzqKioTZs2FRcXz549Oz4+nhksEAgiIyM3b97MRMnw8HDm8lVwcHBmZqaD\ng4NUKj19+vSePXtIhw1Yurq6p06dmjJlSnV19eXLlzkcjlQqXb58uUAgWL16dURExOTJk8ViMfNH\nMJQHANDvMQ+7am9vj42N7bhzFP5CVH2P8q8Oe7D+OhISEszNzXV1dZ2cnKKjo5ntR76+vh1/n8eP\nH88MPn78OJ/Pd3JyysvLYyrJyclWVlZ6enpCoTAuLq776ZTSmpoaW1vbwMBAgUDg7OyclJTE1OfM\nmRMfH08pvXfvHofDMTExsbS0tLS0ZOYWFBS4uLgMHjzYyclJOYVSumfPHgMDg8bGRkopIWTz5s1M\nfcOGDcOGDRs5cmRISMi0adOmT59OKfXz80tOTqaUHjhwwMnJSUdHx83NLT8/n1I6evTon3/+mVJa\nW1s7ePDgO3fuzJkzJyEhgWmMOQCA/q2pqWnu3LlhYWGqbgS6w6GdPckQlJhLHcwmnpd37NgxdXX1\nJUuWsLIa9Ki0tLSoqKj4+HhVNwIA8P89fvx4wYIFixcv/tvf/qbqXqA7uEUI0DmJRNJxwzsAgMqV\nl5e/++67O3f+P/buPB7K7X8A+DNja1Bh7CODslRkK2MrQnKLIiFrckurqNteqNRtcaMSlSXpZomi\nUCmFSMYacVG3Fd2K7GasM8/vj+f3nZcrtE2W7uf9un88znO2ee68+HTOec7xMzY2Hu2+gM+AfbAA\nGFxRURH2HiIAAIwFpaWly5cvP3/+PERX4wKMYAEwuMTExNHuAgAA/L87d+74+fklJSWxXmEGYxwE\nWAAAAMCYFh0dffHixZSUFNYLy2DsgwALAAAAGLv279//7Nmz1NRUHh6e0e4L+AoQYAEAAABjEYPB\n8PLy4uXljY6OHrANMhj7YJE7AAAAMOZ0dnba2tpOmzbt2LFjEF2NRxBgAQAAAGNLQ0PDokWLHB0d\nPT09R7sv4BvBFCEAAAAwhjx//tzJySkwMBA7aR6MUxBgAQAAAGNFcXHxhg0bIiMjZ8yYMdp9Ad8F\npggBAACMD5WVlQoKCj+6eGNjo6ioKIqirIsvb+LGjRu2trbf3MNbt255eXldv34doqufAARYAADw\nn5aVlSUpKclgMFgpPT09urq61dXV2I/Pnj1bvHixqKjo1KlTk5OTsTAFRVFlZWXc/6xatWr4Olm+\nIcoRERFZsmQJiqJFRUXKysqfLTtUlPOFxYuLi2fPno3D4VgXg2bLysoiEAh9fX39EysqKmbPnv3Z\nJgYVFhZ29uzZ27dvS0hIfFsNYEyBAAsAAEbfZ6OcQXMGBQXh/q26urq1tdXT05NMJgsICPj6+g7f\nbnZ2tp2dnaqqKgcHB4IgfX19u3fvplAo1dXVWBjU1NRkZmbm4eHx+vVrAQEBV1fXwsJCZWXl3Nxc\nIpHY/j/nzp0bqs4BviHKUVBQaGhowOFwRUVFM2fOHPA0vjzK+bT4oIqKijQ1NftfDGr69Ok5OTmc\nnP9aaVNRUTFMkWHs378/Ozv72rVr/Pz831AcjEEQYAEAwCj7bJQzVE53d3dWiHPs2LF169YpKio6\nOjqKiooWFxdTqdQjR44wmcyh2s3NzfX397e0tNTT08NSent79fT0TExMKBQKHo9HECQxMdHAwMDM\nzKyiogKHw3V3d2Nhyvnz52k0GolEUlFRycnJYe2B+WmdLIWFhQsWLFi/fv39+/djYmKYTOaECROa\nmpoQBGlpaSGTyW/fvi0vLzcyMhIWFrazs2tra0MQJCUl5cmTJzU1NT4+PsXFxa2trUZGRuLi4t7e\n3giCTJ8+/dy5cwsXLhQREVFUVLx79y6CIKWlpYcPH96zZ4+pqamkpGR0dDTWgU+LD9oHLLZD/hfk\nIQhSUFBgbGzcvwkEQfz8/IqKihAESU5O1tDQEBYW9vLyqqio0NDQQBAkJyfHwMBAREREVlb20qVL\nWJGsrKxZs2aJioru3buX9Vj6+vrc3NwQBLl06RI3N/fXfHHA2IaCYa1Zs+bZs2fsqi0iIiIqKopd\ntQEAfgIPHz40Nzd3d3f38/PDUuh0ekpKyrZt28zMzIbPyXLnzh0zM7Oenh4URQUEBJKSkpqamg4e\nPKijozMgZ0FBgYmJiZiYmJiYmKamZltbGw6Hu3HjBoqizc3N0tLSdXV1rq6usrKyRCKRTCYvXLhQ\nTU0tLS1NRESEh4dHQkJCV1d3/fr1kyZNEhAQ+O2333x8fCQlJfPz842MjAQEBPj5+W/cuKGlpXXr\n1i0ikbh79+4FCxZISEj4+fmJi4unpaXp6OgsWbKEg4NDTExMSEiorKyMwWBwcHDs3LkTRVF9fX0h\nIaHXr18bGxtzcnIyGIxly5bNmzcvLi7u0aNHeDyeQCDIy8tHRETw8PAwGAxbW1shIaGsrKzExERp\naWk8Hu/p6SknJzdhwoRHjx5lZ2dj85gyMjIXLlwgEAgeHh7p6ekKCgo4HG737t0oiiooKJSVlaEo\numPHjv3796MoKi0tXVtby7ooKiqSlpbOysrq7e0NDQ0VExNjMpkoimpqapaUlNy7d09OTu7Jkyft\n7e3Lly+fNm0aiqKPHj2iUCgVFRUMBiMrK2vChAkMBoPJZIqLiz948ODjx48UCqW3txdF0ba2tkWL\nFoWHh7P7awVGHwRYnwEBFgDgx8nLy1u6dGlHR4eWllZGRkb/W9u3bz9w4MCX5KysrKRQKE1NTdiP\nW7duZf0T+v79+/1zFhYWYlFOXl6enp4ePz8/nU7H4XBUKhXtF2FMnjzZy8urt7e3rKyMm5t74cKF\nvLy8wsLCenp6Li4uWJjS09NTXV3Nw8Pz8eNHHh4eaWnps2fPWllZnTlzRlRUlEAgpKenEwiER48e\noSh6+/ZtLi6uqKio3t5eAoGwadOmV69eVVdX4/H41NTU169fc3FxYX3w8vJatWpVVFQUHx+furo6\niqLS0tJSUlKpqakSEhKTJk2i0WihoaHCwsJKSkpYV48cOYJFObm5uby8vNbW1sLCwps2bcKiHCqV\nysfHl5mZycPDIyQk1NXVJS4uHhsbKy8vj0U5pqamN2/efP36tZycXEdHR319PRZCsS4WLVoUGxuL\nPcDm5mY+Pj4Gg0Gj0SZNmtTT02NoaHj9+nXsbmRkpJ2dHYqiCxYs6D+QQSQSURRlMpmzZ882MTEJ\nDQ1tbm5GUfTDhw8GBgas4uAnAwHWZ0CABQD4QfLz85ctW0aj0eh0OoFA6Ojo6H9XX1//zp07n83Z\n1NSkrKz89OlT7Mfg4GAjI6Pq6uqGhobAwEAsCmFZsGDBpUuXUBQ1MTEZMJvx5MkTLMLo7u4mEAg3\nb96MiopauXKluro6g8FAUZRMJsvJyZ09e1ZISCg+Pv7FixcvXrxQUlIqKCgQFBSMjY39tM41a9Zg\n7ba0tOBwOGwGUEhI6O3btx8/fnzx4oWAgMC5c+ecnJyUlZVv3rx57NgxHh6ePXv2xMbG8vHxbd26\ntb6+XkREhEQiLVq0aN26ddbW1iiKNjc38/DwuLq60mg0HA5XU1ODRTnJycnz5s2LjIyUlpa+ePHi\ngCiHn5/f2toai3KUlZX19PSwKMfd3R3rw8WLF1EUvX379qJFi/pfEInEd+/eYR8EawJF0ZycnLlz\n56IoKiAg8P79e+zu2rVrjx8/jqLopEmTPnz4gCXS6fTOzk7suq+v7+7duytXrpSWlq6ursbiv+/+\nHoExCtZgAQDA6Ni7d29iYiIfHx8vL29nZyc/P39LSwt2q6enp7i4WEtL67M5Dxw44OzszFqqdebM\nmdDQUEVFRWFhYUtLy87Ozv4tFhcXY2FQenp6Q0PDx48fT548iS1jOn78uI+PDx8f38GDBxEEKSws\nFBERkZCQqKysTElJef36dWNjo4KCAg8Pz/z5869everj44OtKNq8eTODwTA0NExPT0dRNDk5eerU\nqe7u7vb29qxlWLm5udzc3NevXy8oKFBXV9fU1KysrMzLy8NeS6yqqtLR0amtrfX29t63b5+fnx8W\nnWhraxcXF0+ZMmX27Nn5+fnYqn8EQXJycohEoq6ubklJCRcXV0pKSmlpKTc3t6en5+HDh6lUKo1G\nu3nzJpVKffjwoby8/N27d+l0urOzs66uLjZcJysry2QyVVVVURQlk8lYH5ydnZF+665YFxMmTEhM\nTGQwGPn5+VgTCIJQqVQKhYIgCIFASEpK6uvri4+Pj4yMxIrIyMiEhYV1d3cXFhbKy8uXlJQgCCIo\nKHjv3j1jY2MzMzMGg+Hq6hodHY1VAn5OoxzgjXkwggUA+NFOnjzp7u7ePyU/P3/GjBmfzVlVVaWk\npIQtvcKIiopGRkYyGIxXr14tWLDgzJkz/YtLS0sHBwf39vYWFxeLi4tnZ2c7ODhYWVktWrRIU1MT\nG6ni5eW1srLC4glZWdkNGzaIiIgICQnJysq2trZu2rTpxIkTT58+pVAo3NzcZDI5LS2NRCIFBwf3\n9fVhscuCBQvCw8OnTZtmY2PT0dFRWVkpLy9/7NgxRUVFLi4uUVHRP//8E0XRDRs2rF69GkGQzMzM\nw4cPL1q0SExMTEpKSlFR0cDAYNq0aYsWLTp06NCvv/4qLCwsICAgKyv78OFDrAl5efny8nJ/f39s\n2RMej582bVpmZuaVK1e4ubmlpaW9vLw4OTklJCTi4+MLCgpIJNLMmTMfPXokICCQlpampqZ2/Pjx\nqVOnoiiKLYHPzMzEHpGlpSW2Io11kZaWNm3aNCEhoblz52ZnZ2PZrK2tExISUBRNSEgQExOTlZU9\ndOiQiIhIS0sLiqIlJSUaGhoTJ05UVVVNSkrCipw6dUpVVVVAQEBBQUFdXZ017gV+VhBgfQYEWACA\nH83BwWHAMufAwMBff/31szkXL1589uzZ/hlu3LiBRQOampqxsbHYcmyWe/fuKSoq8vPzKykpYVEO\n+kmEsXXrVizK2bNnzy+//IJNk1GpVGFh4TNnzmhra2PLqlAUVVNTKy8vRwcLQZqampSUlLy8vISF\nhdXU1FhBxoDiR48etbS0/LQP3xPlTJ48Gatz+ChHS0uroKAARdHLly9j+UfG2bNnlyxZMmA6GPyU\ncOjX7FH7H+Tu7r59+3Z5eXm21HbhwgVOTk4XFxe21AYAAN8vOjr66tWrSUlJbKwzPT09LCwsPj5+\nmDx1dXUmJiYpKSny8vJs7MPvv//OZDL37dv3JZn79+H7mx4eiqK7d+9uaWkJDg4edIcw8JOBNVgA\nAPDfVVdX5+fnd/z4cfZWW1hYOGvWrGEyXLhwQU9P7/Dhw/Ly8uztQ1FRkYqKypfk7N8HtjQ9jN7e\n3lWrVgkKCp47dw6iq/8ICLAAAOA/6sdFGJ+Nctzc3N68eWNtbc32PiQmJi5duvRLcrL6wJZ2h9HW\n1rZ06VJDQ8OdO3f+6LbA2AFThJ8BU4QAADDCUBQ9ePDgZ8/5GRfevn1rZ2fn6+s7YNsI8NODESwA\nAABjS09Pz+nTp0e7F2xQXV1tbW198uRJiK7+gyDAAgAAMLbw8PB8yanMY1xOTo6bm1tsbOygJ0+D\nnx7n57MAAAAAI4tIJDY1NQkJCY12R75RfHx8aGhoamrq+P0I4DtBgAUAAGDMIZPJb968GafRyZkz\nZ3Jzc2/evMnDwzPafQGjBqYIAQAAjDlYgDXavfhqTCbT09Ozurr68uXLEF39x0GABQAAYMwZjwFW\nV1eXvb29nJzcmTNnYLMrAAEWAACAMWfcBVitra2WlpYWFhaenp6j3RcwJsAaLLapq6vr6ekZPk9D\nQwMnJ+fLly+HysDPzy8qKsrurgEAwDgjIyMzjgKs2tpae3v7gwcPGhkZjXZfwFgBARZ7dHZ2qmtp\nq8+d/yWZ7xSUDnXr499/lZSUsK9fAADw7fbs2cPHx7d3796Rb5pIJDY2No58u9+grKxs7dq14eHh\nysrKo90XMIZAgMUeTCaTrDh93cHvPUsraJ0DW/oDABgLkpOTlyxZ8uzZsy1bthQWFk6cODEwMHDJ\nkiWsDD09PYaGhhcuXFBSUkIQBEVRFRWVv/76C7vr6uoaGRnJypyVleXg4FBbWzti63uKioo2bNgw\nMm19alwcNPLgwQNfX9/ExERJScnR7gsYWyDAAgCAH2Xt2rX6+vpmZmYhISELFiyIj493dXX9+PEj\nHo/v6+vz9vZOS0t78+aNgoIClj83N5dIJLa3t2M/cnFxsarKzs62s7PT0NAYsegKRdGioqJR3PCT\nj4+vvb194sSJo9WBz4qNjb148eL169cFBARGuy9gzIFF7j9Eb0/PemNKy8d6VkpLY8P+VbYrtaYf\n3eBKb28bJhEA8NPg4OBITEw0MDAwMzPj4OCYP39+d3c3dqu3t1dPT8/ExIRCoeDx//+r+Pz58zQa\njUQiqaio5OTksN7zz83N9ff3t7S01NPTG9BEVlbWrFmzREVFsYm8169fs85OLSws1NbWRhCERqMJ\nCwvv2bPH1NRUUlIyOjoay5CTk2NgYCAiIiIrK3vp0iUscdOmTb6+vhoaGps3b6bRaDExMZKSkmQy\nOS0tDctQXl5uZGQkLCwsIyPzQ2cPyWRybW3tj6v/O+3fvz81NTU5ORmiKzAoCLDY7+alcG8nq/q3\n//q9EHX0AFlhevjDMh4C79Wzp4ZJBAD8NPz9/Xt6ekpLS+vq6l6+fOnm5jZ9+nQsnCIQCObm5jgc\nTkdHB8tMp9OfP3++d+/eV69eubi4rF69GkunUqn+/v5xcXGlpaUDAiwURe3t7c+cOVNVVXX//v2+\nvj4qlaqrq4vdffToERZgPX78mE6nW1hY3L1798KFCx4eHiiK5uXlbd++PSQk5MOHDxcvXly7di2T\nycSay8rKCgsL09LSYjAYcnJyb9682b17N6s/mzdvdnBweP/+fXJy8h9//IGV+hHG7IuE2GZXXV1d\nsNkVGAYEWOxHVpy+fINX/xQURQvu39E2XczFzU1Z8Eth5t2hEgEAPxN7e3tnZ2dxcXFFRUUnJycB\nAYEB66Dz8vKwGAhBEF5e3ry8PCsrKyEhIS8vr+bmZiaTWVBQ4O/vHxMTg8fjy8vLtbS0BjQhJSXl\n5+eXmJiYlpbGyck5aIBVWFj466+/YpGcjo4O9r6zr69vfn6+srIyBweHoaEhHx8fHo+n0+kVFRUx\nMTGamprl5eWurq7Ozs5cXFwrVqxobGzEYqkbN24oKCjExMQEBASoq6uzht/YbmwGWJ2dnba2tnJy\nckePHsXhcKPdHTB2wRqsz3j37l1oaKiwsPDw2fpv0KBMGTiG30WnddFpJLlpCIJITZVvrv8wVCIA\n4CfT2tp68+ZNbNGVsrLyyZMnWbd6enqKi4tZMVNSUpKqqqqcnByCIM+fP9fQ0MDj8Xv37r13715i\nYiKWh5+fv7m5mTUnhcPhqFRqRkZGdHT0oUOHXr9+TaVSf/31VwRBWlpa0tPTjx07hiBIYWHhwoUL\nsSK5ubkaGho4HC4/P//Dhw/YvjCdnZ1YrFBSUqKmpkYikRAEKSoqcnFxwUrl5+fPnj0bj8dHRESc\nO3fO3NxcS0tLQkKCSCT+uEdHJpNTUlJ+XP3f4OPHj7a2th4eHlZWVqPdFzDWQYD1GV5eXr29vZyc\nn3lQXV1dNzKyh8+D/f5CUZTJYAyfCAD4aSgrK0dFRZmZme3du1dWVtbU1JR1q7S0VFZWlhUtxcXF\nXbt27fz5842NjZs3b/79998RBElPT8funjp1qrKy8vz58/0rFxQUjIuLW7BgQWNj48OHD3E4XHl5\nOScnZ19f344dO7i5uclkMoIghYWFfX19y5cvr6mp8fLyCg4ORhBERkYmLCxs27ZtT548sbKyio+P\n19XVZQ2AMZnMx48fk0ikFStWPHv2bPPmzSEhIQiC7NixIyMjQ0VFpaio6MqVK1gA94OMtRGsN2/e\nODs7Hzt2jDWrC8AwIMD6DGNj4y/JRqPR9p8YchHVBF6+Cbx8/7x6oag++93rl0RxiaESAQA/mZMn\nT65Zs4aLi2vp0qVXrlzpP6H26NGj/n+q/fz8XFxcSCTStGnTDh8+PH/+v/bVKygo+HQTywMHDuzc\nuXPFihUKCgqxsbEIgmzcuNHAwEBBQUFHR0dbWxuHwzU3N3NycpJIJBkZGSkpqePHjy9YsABBkIsX\nL65evfrYsWPY0S5YXEWlUu3s7BAEef78uaam5sSJE6dMmSIjIxMYGIj9MnR1dTU3N+fj47O2tlZS\nUrp48aKNjU3/LtXX17Nrt2QJCYn379+zparv9/jx4w0bNly4cGH69Omj3RcwPuDGxUYjYx+NRjMw\nt9xz7hIrxVpJMuJhqYDw//+iOblto6iUtL3njnM+2/kmCbhs3zdoYtA6h8zMzNH5DACAn1F6enpY\nWFh8fPzINOfh4SEsLOzr68uW2gwNDbOysthS1fdIS0s7cuRIXFychAT8Sxh8KVjkPkJcd/k+fVzk\nbqjZ1tRos95zmEQAAGCjwsLCWbNmjVhzQUFBEyZMsLGx6erq+v7auLi4PnsE2Y92+fLls2fP3rx5\nE6Ir8FUgwPpRrlX/wxq+QhBEQFj0QFRC2IOSncGRBP6JwySCsamyspK1GyQAI+k7v3uZmZlnz579\nbLbGxkZRUVEURVkX39zizp07Fy1atHTp0tbW1m+uBDNlypS6urrvrOR77N+/Pzs7+9q1a/z8/KPY\nDTAewRostnn78vnlE4e/s5KWlha2dOZn1dbWtnXr1pSUFA4OjnXr1vn4+AyVyPLp6SIDDicZvnh3\nd7eNjc2NGzeKiorgoDEwktj13XN2dr5+/fpnsxUXF8+ePRuHw7EuvrlFBEFWrVpFIpF++eWXK1eu\nTJky5Zvrwda5Y29WjrC+vr5169ZJSUmdP38etmMA3wACLPbg4+O7fSOJ8bk3AZOTkzk4OBYvXjxU\nBv7fPNjdtZ+Kk5MTkUisqKhob29XUlLat28fHo8fNBHLP+B0kUEPJxmmOIIgZWVlDQ0NOBxudM8M\nAf9B7PrufWHxoqIiTU3N/hffydTUVFBQcOnSpZGRkaqqqt9WyWi9SEin011cXBYtWuTm5jbyrYOf\nBApGUERERFRU1Gj3Yrzq7e1NSkqi0Wjv3r07ffq0rq7uUImYhw8fmpubu7u7+/n5YSl0Oj0lJWXb\ntm1mZmbD1Mny+PFjcXFxSUlJb29vXV1dDw+P+fPni4mJ7du3D8vw5MmT+fPnE4lEMpm8Z88eFEX7\n+vqEhYU3b968ePFiWVnZ0NDQZcuWUSgUd3d3rMjixYvXrFnj5OSkqqq6Y8cOLDE/Px87e0RBQeHO\nnTs/8CGCEVdQUGBiYiImJjZlypTo6GgGg8HDw9PY2IiiaHNzs7S0dF1d3adfpM9+9z79znR0dBCJ\nxN27dy9YsEBCQuLy5ctYzk+LD9qHZcuWXb9+HUVR1sWnvfoGL1++pFAoWVlZ31Y8IyNj//7931b2\nm9XX1xsaGt64cWOE2wU/GXiLcERduHCBk5OTtXffd6qpqVm7dq20tDRbahv7vL29paSkSkpKsH9e\np6WlsfZO/DSRSqUePXo0OjrayMjo6NGj/d9437FjBz8/f/+pwEHrxDg6Oi5ZssTa2nrSpEmrV68+\nceLEy5cvVVVV6XQ6Ho+fP3++o6Ojq6trZWXlnDlzOjs7KysrNTQ0Hj9+PHPmTFtb29bW1tTU1JaW\nFikpKWwvR1FRUR8fHw8Pj3fv3pHJ5K6ursePHy9btuzSpUt6enqRkZHe3t7v3r2DKYmfQ1FRkYWF\nxcWLF01MTNLT021sbNrb2xUVFRMSEmbNmrVz505eXl5fX99Pv0h4PH6Y796g35nc3FxTU9P79+/r\n6OikpaU5ODg0NjYyGIxBv7qf9oFMJufm5kpJSbEuBu3VNzyExsZGGxubNWvW2Nvbf23Zly9fHjp0\n6MKFC9/Q7rf5+++/XVxcTp069emm+QB8FZgiHMfa29uJROLOnTtHuyMjRExMDEEQdXX15ubmyMjI\n1atXsw6CHZDIOl0E23dxwC/KvLw8b2/v/imD1onJzc09cuRIdXU1gUD4448/uLi4uLi4ZGVlsb80\nN27cKC0tjYmJycjIwM4MoVKp9vb22IxMbW1tQEAAFxdXbW2tqqoqHo9/8eIFDw/Pxo0bEQRpamrC\njqXz8fE5duyYgYEBgiA2NjZbtmxBURQCrJ/Dnj17jh8/jkXtpqamVVVVCILIyMjU1dVNnjz56tWr\nT548QQb7IiHDfvcG/c4Meh7OUF/dAX1oaGjo7u4mkUisi6F69Q2IRGJqaqq9vX1NTc3X/r4a4UXu\nRUVFW7dujY6OHpVVX+AnAwHWOIbH4ydNmvTf+UXg6OgYHh5OIBAEBASWLl0aFBQ0VOIwp4sMOJxk\n0OIs//zzT19f35QpU6KiouY2NDkiAAAgAElEQVTPn8/NzY0gSF5eHvY3bNAzQ6hUqr6+PoIg3d3d\nlZWV2MBYfn4+hULBLgwNDbE/VKx68vPzIyIisBZzcnI0NTV/3OFuYIQVFxdHRUVh101NTQQCAUEQ\nGRmZ2tra2NhYHx8fPj6+Qb9Iw3/3Bv3ODHoeTlFR0afFP+0DVgm2wh27YO+ROLy8vImJiR4eHp6e\nnoGBgV/+Defi4urt7f2epr/czZs3T506lZSU9EPP/wH/HfB7HIwbd+/ejY2NZTAYb9++3bhx4+bN\nm4dKTE9Px6bAT548iS1+Yp1GMuBwkkGLs+Tl5bHeq+p/gC72V2rHjh3h4eHe3t5EIvHKlSvYqbpU\nKhW7ePz4sbKyMvaHjZXIuuh/PWHChMTERAaDkZ+f7+npefjw976LCsYOfn7+pKSkvr6+kpISFRWV\nyspKBEHIZHJycnJVVZWzszMyxBdp+O/eoN+ZwsLCmzdv0mi0qqoqLy8vbJh20OKf9gF7c7D/xaC9\n+h4cHBwhISFycnK2tradnZ1fXhCHw2GHTP9Q58+fDwsLg+gKsNPoLf/6L2LvIvfKysr169ezq7ax\nLyEhQUZGRlBQUFVVNSIigslkDpXI4uDgEB4e3j8lMDDw119/Hb5OFiqVKiwsfObMGW1t7UePHmGJ\nampq5eXlKIpu3bpVSkpKUVFxz549v/zyy6JFi1paWiZPntzX14c19Ntvv2FFFBQU/v77bxRFtbS0\nCgsLscSZM2dWVVWhKJqWljZt2jQhIaG5c+dmZ2ez85GB0Xbv3j1FRUV+fn4lJaU///wTS4yOjkYQ\nJDMzE/vx0y8S+rnv3qffmaamJiUlJS8vL2FhYTU1taSkJKzIoMU/7YOlpSW2ppt1MWivBsCG2b72\nmURFRc2bN6+hoeEL8zs4ONTV1X1tK19l586d69at+4bPAsAwYJH7iGLvIveqqqqgoCDsBFYAwHgR\nHR199erVpKQkNtb5tefhsKUPR48eTU1NFRISmjdv3rx58zQ0NDg5v2jZyb1797y9va9cufIl7+js\n2bPH3NycNQ7HXr29vWvXrlVUVPzvLGYFIwamCAEAYOTU1dX5+fkdP36cvdV+1Xk47OrDrl27Hj58\nGBsbq6amduvWLQsLC11d3V27dt27d2/4SUATE5PQ0FBra+vHjx9/tpUftxVWe3u7paWlgYEBRFfg\nR4AACwAARsiFCxf09PQOHz4sLy/P3pqLiopUVFRGpQ98fHwmJib79++/fft2enq6iYnJw4cPly1b\npqent2vXrpSUlEEPzFFRUUlISHB3d09LSxu+/h8UYH348MHS0nLLli0rV65ke+UAIAgCU4QjCqYI\nAQD/EXQ6vaSkJDc3Nysrq729XV9fX09Pb+7cuaxXTBAEaW5utrGxWbVqlaOj41D1VFVVnT59+kuO\nU/xyT58+Xb16dXBw8Egegw3+a2CbBgAAAOzHy8urr6+vr6+/c+fO3t7eJ0+e3Lt3Lzg4mEajaWho\n6Ovrm5iYCAoK3rx5c+XKlX///ff+/fsHrYdMJtfU1LCxYw8fPty5c2d0dLSMjAwbqwVgAAiwAAAA\n/FhcXFyampqampo7d+7s6+srKyu7d+/eqlWr6uvrVVRUrKys0tPTPTw8Tp48yTqUnYWXl5dGo7Gr\nJwkJCWFhYSkpKUJCQuyqE4BBQYAFAABg5HBycmLBFoIgdDr90aNHDx48ePXq1atXr27dupWTkyMp\nKTmgCLsONggJCXn48GFKSgoPDw9bKgRgGBBgAQAAGB28vLwmJiYmJiYIgnR2dp4/f97R0ZFCoezc\nuVNQUJCVjUgkNjY2fs8WoEwmc8uWLUwm888///x0kAyAHwHeIgQAADD6CASCl5dXZmamiYnJ0qVL\nd+3axXr98DtfJOzp6XFzc5OTkwsKCoLoCowYCLAAAACMISYmJg8ePNDU1DQzMzt27FhnZ+f3BFit\nra1WVlbm5uaenp7s7ScAw4MACwAAwNiCw+FsbGyys7MFBQUNDQ2fPn366tWrb6inrq7O3Nx827Zt\ny5cvZ3snARgeBFgAAADGIi4uLnd39/T09O7u7uDg4K8tXlVV5eDgcO7cufnz5/+I7gEwPAiwAAAA\njF2TJk0KDw93cHD4qhgrOzt7w4YNcXFxM2fO/HF9A2AYEGABAAAY63x9fa9evfr06dMvyRwXF/f7\n778nJSV9uuMDACMGAiwAAABjHScnZ0hIyOrVq3t7e4fPefz48dTU1OTk5P5n8gAw8iDAAgAAMA5M\nnz7d2tr62LFjQ2VgMpleXl50Ov3PP//k5uYeyb4B8CkIsAAAAIwPnp6e+fn5hYWFn97q6upasWLF\n1KlT9+/fz66d3wH4HhBgAQAAGB9wONy5c+e2bNnS2dnZP72lpWXZsmUODg4eHh6j1TcABoAACwAA\nwLhBIpHWrVvn4+PDSqmpqbG0tPTx8bG0tBzFjgEwAARYAAAAxhMnJ6e3b99mZWUhCFJaWurk5BQe\nHq6trT3a/QLgX+CwZwAAAONMSEiIhYXF1q1bT506FR8fLy4uPto9AmAgCLAAAACMMwICAjY2Nhs3\nbqysrITtGMDYBFOEAAAAxh8nJycdHR2IrsCYBQEWAACA8UdISKipqWm0ewHAkCDAAgAAAABgMwiw\nAAAAjEt8fHw0Gm20ewHA4CDAAgAAMC5JS0vX1NSMdi8AGBwEWAAAAMYlMpn85s2b0e4FAIODAAsA\nAMC4BAEWGMsgwAIAADAuQYAFxjIIsAAAAIxLEGCBsQwCLAAAAOOSuLj4u3fvRrsXAAwOAiwAAADj\nEh6PZzKZo90LAAYHARYAAIDxipubu6enZ7R7AcAgIMACAAAwXk2ZMqWurm60e8EGlZWVCgoKo90L\nwE4QYAEAABiv2LvOva2tbfXq1WJiYpKSkgcPHhwmkSUrK0tSUpLBYLBSenp6dHV1q6urv6R4d3f3\nkiVLUBQtKipSVlZm1wcBYwEEWAAAAMYr9gZYTk5ODAajoqLi4cOHhw4dwhZ4DZqIyc7OtrOzU1VV\n5eDgQBCkr69v9+7dFAqlurqaNRw1THEEQcrKyhoaGnA4XFFR0cyZM9n1QcBYAAEWAACA8YqNAVZf\nX5+bm1twcDCDwbh58+acOXPwePygiVj+3Nxcf39/S0tLPT09LKW3t1dPT8/ExIRCoWDZhimOIEhp\naenSpUtramp8fHyKi4tbW1uNjIzExcW9vb2xDOXl5UZGRsLCwjIyMnv37kUQhMFgiIiIeHp6mpub\ny8nJhYWFWVtba2trr127Fitibm7u7u7u7Oyspqa2c+dOLLGgoMDY2FhERERRUfHu3btseVzg81Aw\ngiIiIqKiothVW2Vl5fr169lVGwAAjDsvX75ctWoVGyssLi7G/jimpaUNk5iXl7d06dKOjg4tLa2M\njIz+NWzfvv3AgQOfrRPj4OAQFxfX29tLIBA8PDx6enqqq6t5eHgYDAaKooaGhmFhYb29vWVlZdzc\n3AwGo7y8nIuLq6KiAkVRGxsbU1PTnp6e+vp67C6TyRQWFj59+jSKov/88w8XFxeDwSgqKpKWls7K\nyurt7Q0NDRUTE2MymWx8YmAoMIIFAABgvJKSkmLvInd1dfXm5uaAgIDVq1cPlVhQUODv7x8TE4PH\n48vLy7W0tPrXkJeXp62t/dk6Mbm5uTo6OtXV1QQC4Y8//uDi4uLi4pKVlcUGum7cuKGgoBATExMQ\nEKCuro7H46lUqr29PTaZWFtbu3//fi4urtraWlVVVTwe//LlSx4eno0bNyII0tTUNH36dDwe7+Pj\nc+zYMQMDA05OThsbm46ODhRF2fjEwFAgwAIAADBecXFx9fb2sqUqR0fHzs5OHA4nICCwdOlSLi6u\noRL37t2bmJjIx8fHy8vb2dnJz8/f0tKCVdLT01NcXMwKuQYtzvLPP//09fVNmTKlqKho/vz53Nzc\nCILk5eXp6OggCBIREWFsbJyZmSkiIiIhIYFNRFKpVH19fQRBuru7KysrNTU1EQTJz8+nUCjYhaGh\nIRacserBErEWc3JyNDU1+09Tgh8HnjIAAIBxDIfDsWW70bt378bGxjIYjLdv327cuHHz5s1DJaan\np2NzQCdPnnR3d0dRVEBAAKuktLRUVlaW9eOgxVny8vJmz56Nw+GKi4t1dXWxxEePHmGB0Y4dO8LD\nw729vYlE4pUrV7BRMSqVil08fvxYWVkZi8lYiayL/tcTJkxITExkMBj5+fmenp6HDx/+/mcFvgQE\nWAAAAMYxCQmJ9+/ff389Z8+e9fPzExERWbx4sY2Njaen51CJLAUFBQPmB1nh0TB1skhJSeXm5gYH\nBxcVFbFKsWpwdXU1NzefMWPGjRs3lJSULl682NraWldXN2PGDARBqFQqqwiVSmWNYH0aYEVERAQG\nBoqKim7fvj0qKgobAAMjAAdzsSPpwoULnJycLi4ubKmtqqoqKCgoJCSELbWB4VVWVlpaWj579my0\nOwIA+Jc5c+YwmUwhIaHvrCckJEReXp4tXQIAQRDO0e4AAF+hra1t69atKSkpHBwc69at8/HxwdIL\nCgrWrl1bV1d34sQJLH5FUVRFReWvv/7CMri6ukZGRg6aiCBIS0vL2rVrs7OzyWSyu7u7m5sbq8Xu\n7m4bG5sbN27ANoAAjE1ycnJBQUGioqKj3REA/gUCLDCeODk5EYnEioqK9vZ2JSWlffv24fH49vb2\n5cuXJycn19fXOzk5YQFWbm4ukUhsb2/HCmJrSwdNZDKZS5YscXV1jY2NffToUUBAQP8AC7YBBGCM\na2hoEBERGe1eADAQBFj/CbW1tex60Wa04PF4KSkpNzc3U1PTtra2/lv2xcfHW1hYqKmp0en0pqYm\nJpOJx+PPnz9Po9FIJJKQkFBISMgvv/yCIMigiRkZGTw8PFhQpa+v33+BArYNIPaec3FxsaamppGR\nUWVl5Zo1a/z8/BAEKS8v9/T0fPLkCT8/v6Oj4+HDhxkMhri4uIODw4sXLyorK3fv3p2Wlvb27VtV\nVdXz588jCGJubi4pKdnZ2VleXr5w4cJjx44hCFJQULB79+4nT54ICQkFBQWZmpqOxjMGYFxCURSH\nw412LwAYCAKsnx+NRtPX18eWQI5fnJycx44ds7S0LCkpwd5MTktLw27Fx8dv2rQJQZD3799LS0vj\n8Xg6nf78+fO9e/caGBicOnVq9erVb9++HTQRQZCoqCgZGRllZeX379/7+/uvWrWK1aiampqRkdGS\nJUusra3/+OMPTU3NO3fuvHz5UlVV9cCBA3g8fvPmzY6Ojnfv3q2srJwzZ46fn19VVVVra6u7u/vM\nmTNtbW2vXr2ampra0tIiJSV19uxZHA6Xn5/v4+Pj4eHx7t07Mpl85MiRx48f29jYXLp0SU9PLzIy\n0sXF5d27d/AHA4Av0djY+Pjx48zMzPnz5492XwD4F1jkPqK+fJH7zp07S0pKhs9Do9Fqa2uVlJSG\nykAgEOLi4phMpr29fUpKyld3d0xCUbS1tTUyMjIgIKC2thZBEDk5uYyMDBkZmevXr589e/bOnTv9\n8zc3N5NIpI6Ojv5bv/RPnDlzpqqq6qlTp9LT09evX9/c3Nw/p4yMTHZ2dltbm4GBwbt377i5uV++\nfLl48eKqqioEQdra2kpLS1+/fp2RkVFdXU2lUsPDw3NycqKiohAE0dHRCQgI0NHRKSkpWbduXUFB\nwYsXLwwMDGpqavB4/F9//eXg4FBWVrZ48WJnZ+cVK1YgCIKFYm1tbbBRDQBfIigoCEGQt2/fFhQU\n7Nu3z8jIaLR7BMD/gxGsMaqgoMDjXMx3VhKwdV17ezsfHx9bujTqHB0dw8PDCQQCtmUf9ouVyWTW\n1taKiYkhCJKWlmZra4sgSFJSkqqqqpycHIIgz58/19DQwOPxgyYymcyXL19iW/mZmpry8PD0Hzpi\nbQMYFRU16DaA586dMzc319LSkpCQIBKJyLduAxgREYG1CNsAAvBVrl69euvWLT4+vrq6Oj8/P39/\nfz8/v9mzZ492vwCAfbDA+DHoln04HG7SpEm5ubklJSWpqalYgBUXF+fj40Oj0WpqajZv3oytlxo0\nEYfDEQiE/Pz8np4eHx+fVatW9Q+wYBtAAMay4ODg9+/fY/+GlJKSOn/+/NmzZ8PCwszMzMrKyka7\nd+C/DgIsMG4MumUfDoc7ceLEihUr3NzcEhMTJ06ciCCIn5/f8+fPSSTSsmXL9u/fjy3OGDQRh8MF\nBwe7u7tPnTqVn5//4MGD/VuEbQABGLP279+fl5cnIiJCp9NZiTIyMufPnz9x4sSRI0fs7OyeP38+\nij38KpWVlQoKCqPdC8BWo3LE9H9WREREVFTUl+Q0NDS8Vv0P67+4J69FSVMiHpayUiJyy1R09Pkn\nTZ5jZPpnYfWgiXqLlrx//769vd3c3PxHfzQAABgZb968UVRUtLGx6evrww5dHjQblUr95Zdf3N3d\n3759+4U1t7a2/vrrr6KiohISEgcOHGCl5+fnq6mpCQsLs36BM5nM/vu2uLq6DpWIomhzc7Otra24\nuDiFQomIiOjfYldXl4WFBZPJjIqKsrKy+tpHAcYyGMEaB25eCvd2sqp/W9s/MeroAbLC9PCHZTwE\n3qtnTw2TCAAAP4e3b9+uXbt2165dPDw8cXFxHBwczs7Oly9fHjQzhUK5deuWs7Ozm5ubp6fnhw8f\nPlu/k5MTg8GoqKh4+PDhoUOHsCMOsZ32IiMjo6Ojt23bhuVkbaqHOXfu3FCJ2E57CxcufPv27R9/\n/JGamtq/Rdhp7ycGAdZIo9FozV+gr6+PVYSsOH35Bq/+laAoWnD/jrbpYi5ubsqCXwoz7w6VCAAA\nP4EPHz6sXbvW1dXVxcXF2dnZwsICexdETEyMi4urpqZmqIL6+vppaWkWFhY2Nja7du1qbW0dKmdf\nX5+bm1twcDCDwRh0pz19fX1spz2k36Z6KioqOTk5PDw8QyWydtrD4/H6+vqJiYmsFrGd9mpqarCd\n9lpbW42MjMTFxb29vbEM5eXlRkZGwsLCMjIye/fuRRCEwWCIiIh4enqam5vLycmFhYVZW1tra2uv\nXbsWK2Jubu7u7u7s7KymprZz504ssaCgwNjYWERERFFR8e5d+NMwQuAtwpF29erVzMzMz2Z7+vQp\n61qZojfgbhed1kWnkeSmIQgiNVW+uf7DoImy0+HfQwCA8a21tfXIkSO5ubm+vr4mJiYIgtjZ2R09\nepSVYeXKlZcuXdq3b98wlZiYmBgZGV27ds3MzMzS0nLz5s0EAmFAHk5OTthpD7DTaM9RgsENWIN1\nrfofBEFYa7AuFz9DEOQi9a9r1f8EJN/nmUAYNBHWYAEAxq/Ozs6jR49qa2vHx8czGAwssaGhwdjY\nuH+2np4eLS0tJpP5JXX29PRgr5KcPHmyq6vr0wxMJrO5uTkgIEBKSgpLkZWVffXqFYqiSUlJpqam\nA/I3NTURCARW9z5NnDFjhr29fX19fXR09KRJkwbkJJPJb968KS8vFxIS6u7uRlH0xYsXSkpK2N3W\n1tYHDx5ERUWtXLmSQqGgKBoWFubi4oLd1dbWfvToEYqixcXFc+bMQVEUe48Ha6KiomLWrFkoii5a\ntCg2NhYr0tzczMfHN6AP4AeBKcJxaQIv3wRevn9evUAQ5N3rl0RxiaESAQA/t7a2ttWrV4uJiUlK\nSrJeg3327NnixYtFRUWnTp2anJzMytzS0mJnZychIaGtrX3hwgVWek9Pj66ubnV1NfZja2urp6cn\nmUwWEBDw9fUdyY+D6evrCw0NNTAwEBQUzM7OtrGxYW0Od+XKFQcHh/6Zubi4dHV1c3JyvqRmLi4u\nFxeX+/fvEwgEAwOD0NBQBoOB3XJ0dOzs7MThcNhOe6yzSgfdae/ly5dYqf477X2aiO20d/LkyeF3\n2isqKhp0pz1jY2Nslz4JCQk9PT3kW3faMzQ0xFqEnfZGEjzlcQmHw80xMi1+cB9F0cc5mXOMzYZK\nBAD83D5dl93U1GRmZoZNEh06dMjV1RVbNjToauu+vr7du3dTKJTq6mrWNgGOjo6ioqLFxcVUKvXI\nkSNY8ZGBomhCQoKRkVFnZ2dmZqa7uzsW6LBcvXoVi3L6c3Z2xpYofSFubm53d/f79+83NzfPmzcv\nISEBRVHYaQ+w2egOoIGhDD9FeK36n4iHpcoUPSExcS3jhZeLng6aCFOEAPzcent7k5KS3r9/7+Dg\nMHHiRG5u7gMHDoSFhbm6ura0tGzevJlEIuFwOG9vbxRF09PTDQ0NB2xDQKfTU1JSPDw8xMXFWfsI\nCAgIJCUl3blzZ8qUKZycnLKysjdu3GA1unHjRh8fH3V19YCAgOzs7Hnz5mGrsFlbGBgbG//5558o\nisbHx8+bN49Go33hx4mPj9fV1T169GhbW9ugGUpLS1euXDnoLUFBQSqV+qUPrp/37997eHgsWLDg\nxIkTMjIygoKCqqqqERERrDnHyMhIIpGoqqqan5+PpTx9+pRCoUyePFlTUzMtLW2YRBRFY2JixMXF\npaSktm/fPmBSkkqlCgsLnzlzhjXZh6KomppaeXk5iqJbt26VkpJSVFTcs2fPL7/8smjRopaWlsmT\nJ/f19aEoGhgY+Ntvv2FFFBQU/v77bxRFtbS0CgsLscSZM2dWVVWhKJqWljZt2jQhIaG5c+dmZ2d/\nwyMC3wbOIhyj5s2bZ+9z/DsriTjsfTPuMh8f3890FiEAYAADA4Ps7GwEQSIjI93d3QMDA8PDw7GY\noKys7P3798+ePevq6lq5ciV2sCYOh9u1a9e+ffu6urqwaSwymaypqZmYmPjo0aOAgABZWdmAgACs\n8nv37tXU1OzZs+fdu3dYyuzZs/n4+AICAtrb23ft2hURETF9+vScnBwzMzMajYbH4+/fv+/l5XX6\n9Ondu3enpaUJCAh89iM8ePDg4MGDc+bM2bZtm7Cw8FDZtm7dunTpUgMDg09v6erq4nC4tLQ0bLfh\nr/XmzZuGhgY4YwewEbxFOEYZGhpmhgYMn6etrQ2b6R8qw4wpEhMnThzJ4X0AwAjr6+vz8vKKi4sL\nCwvbsmXLnDlzXFxcUlNT79y5U1dXp6ys3NjYOHv2bDweX1xcTCaTL1++nJqa+ttvv7G2IcjIyOjo\n6NiwYQO2j8CTJ0+uXbsWGBj47NkzBQWFTZs2xcTEcHBwYM3R6fSKigp7e/vk5OTAwMD29nZlZWXs\nFh6Pv3z5souLi5GRUU1NjbOzc3FxcXp6+pkzZ27fvs3Lyzto/x8/frx7924ymRwZGSktLT3MJ+3t\n7aVSqSdOnBj0rri4uIODw7Zt286fP/8Nj5FMJpPJ5G8oCMCQRnsIDXy7ysrK9evXfzYbTBEC8LNy\ncHCg0+koihYXF2O/0tPS0mpraxkMxtatW1m/5+/fv89gMCZMmPDhwwdWztu3b7Mq4eTkVFJSIhKJ\nFy5cmD59+vPnz7Fbr169IpFIOjo6rN3Sc3JyKBSKpqbmvHnz+Pj4kpKSKBRKRUVFR0fH3bt3J0yY\ngC0IExcXJxAI9+/fp1Aozc3Ng3b+2bNntra2dnZ2lZWVX/Jhr1+/3n939QG8vLyKi4vXrVuXkpLy\n5Q8QgB8HAqxxDAIsAP7jhIWFIyIi+vr6amtrjY2NLS0tpaSkJk+evHbtWkNDw1WrVhkZGQUEBCgp\nKTGZTEFBweTk5K6uLjc3NyMjo/7bEEyePJm1j4CoqGhkZCSDwXj8+PGUKVNkZWXLyspYLfr7+3t4\nePDw8NTV1c2aNWvatGn9/8VOJBJfv349c+bMJ0+eEAgEGRmZ9+/fD+hzX19fenr6kiVLjI2NHzx4\n8OUf1srK6vXr10PdDQwMTExM7Ojo0NPT+7RRAEYevEUIAADj1dmzZzdv3iwsLGxubu7g4PDHH39w\ncXGdPHkyPDy8srJywoQJSUlJVlZW2O4DM2bMWLNmzbRp04hEYnBwcP9tCBYvXszaRyA0NPTw4cMC\nAgJz5861sbF5/vz5rFmzWC1SqVQRERE1NTUSiXTx4sVXr17x8/OrqqomJSXR6fTHjx9bWFgEBwfj\n8Xg+Pr6mpiZRUVGsYGVlZWBg4JIlS3R0dFxdXRUUFO7duzdv3rwv/KTv37+n0+nDzOJhG0rx8fH5\n+/tv3LjxOx4qAOwBARYAAIxXy5cvJxAIJ06cKC4uXrhw4aZNmzZv3uzq6kokEo8dO3bmzJmmpiZ3\nd/ft27cjCPL06dPff//99evXnp6eW7ZsYW1DMHHixBUrVrD2EVi6dGloaKiiouKrV68OHDhAp9Np\nNBqrxatXrxIIBGxPAXV19ZkzZ+7atSs/P59EIsnLy9fW1j558kRGRsbOzi4zM5NMJmdkZGAFc3Nz\nSSSSp6cnHo+/e/euv7//V33SmJgYZ2fnYTJgARaCIDo6OoqKipcuXfrKZwkAu432EBr4djBFCABI\nSEj4dHOBGzduYG/ma2pqxsbGYomD5kQH20fAzOxfu+jp6+v3b9Ha2jo+Ph67Likp0dDQmDhxIjaI\nhaLohw8fVFRUsrKyUBT9888/+++6npqaqq+v/88//3zDx9TT08NWmw3l48ePVlZW2HVPT4+xsfEw\n84kAjADYpmEcq6qqCgoKCgkJGT5bR0cHbNMAABh17u7uR44cIRKJX1uwqKgoLCzss68HGhoaZmVl\nYdeVlZVeXl5paWmwazkYLfDNAwAAMBJCQ0O/IbpCEOTixYuurq6fzdZ/S5oZM2YsXLjw1KlT39Ac\nAGwBARYAAICxq6urq7S0FDs9ZngTJ05sb29n/aihobF79+6ysjLsx4KCAnV1dREREdbyLBRFlZWV\ncf+zatUqBEGCgoJw/1ZdXT3qhzOC8Qg2GgUAADB2JScnW1hYfElOMplcU1Mzc+ZMBEGys7NXrFhB\noVD8/PyuXr3a3t6+fPny5OTk+vp6JycnFxcXBEFyc3OJRCIrJsNeq3R3d8ciLQRBQkJCXr16paio\naGFhoaOjU1xc/PHjx07XSSgAACAASURBVFmzZvn6+sLMI/gs+IoAAAAYuy5fvrxy5covycl6kTA3\nN9ff39/S0nLBggVXr15FECQ+Pt7CwkJNTU1fX7+pqQmbTDx//jyNRiORSCoqKjk5OTw8PAiC8PDw\n8PPz8/PzP3r0KDMz8/Tp0zgcLjc3d+bMmRwcHAkJCdi2+P3b3bRpk6+vr4aGRmBgYE5OjoGBgYiI\niKysLGuozMTE5PLlywiCJCQkGBgY0Ol0tj4hMFaN9ip78O3gLUIAwM+ktbV1wFnUNTU1FhYW2KAU\nxtXVFUVR7ChraWnpyZMn+/j4YMVjYmIkJCRYOYWFhTMyMphMZv/i1tbWsrKyKIrSaDRtbe3ExMTG\nxkYfHx9JScn+PamsrKRQKE1NTdiPA7bFH9BtbF/7oqKizMxMbF97BoORlZWF7WuPoui9e/eUlZUz\nMjKG2dce/HxgihAAAMCY4OTkRCQSKyoq2tvblZSU9u3bFx0dPWfOnNbW1gETeY6Ojp/O2WH7dS1Z\nsiQsLAyHw5HJZC0tLWwekEwm37p1S1pa+vbt21hVvLy8eXl5WJ1eXl7+/v5MJhMbmmpubra1tb12\n7ZqgoCCCICEhIaWlpdXV1UQi8fLlyxs3bqyqqmL1GTuc8cWLFyQSydTUND8/n3U4I5FIxCo0MjIi\nEAgeHh7379//kqOvwc8BpggBAAB8Xltb2+rVq8XExCQlJQ8ePIglooOtE0cQJC8vT1dXV0hISE5O\nLjk5eajiCIK0tLTY2dlJSEhQKBQymRwcHMxgMG7evDlnzhwcDpecnFxdXf3pRN6gc3YPHjyg0WjJ\nycliYmKioqKdnZ38/PxBQUE0Gu3NmzeLFi3Kycm5f/++ra0tgiBJSUkvX77EZvd0dHTExcVzc3Ox\n2T0ZGZkZM2YoKCggCGJiYnL48OHQ0NAnT55YW1svXLiws7Oz/2MpKSnB9rVHECQ/P//Dhw/Y6AWd\nTq+rq8PyVFZWNjQ0EAgE1r724D9hVMfPwHeBKUIAwIixsLBwdXWtr69/8eIFFxcXNvmVk5Mzb968\n9v/B9in9559/hISEbt++zWAwLly4IC4uPlRxBoMxd+7ciIgIBoORk5NjZWXV/9Tq1tbWEydODDqR\n9+mc3YApv8mTJ7u7u2OJ2FiUk5OTiIgIiURqa2tDUdTW1tbR0VFdXZ1CoaioqAQEBGCze3/99Ze0\ntHT/2T1OTs4dO3ZQKJSysrIFCxacOXOm/2Px9/ffsmULdj1r1qxDhw51dXUVFBSQSKTc3FwURbHD\nGcvLy1VUVO7duzdS/7vA6IMAaxyDAAsAMDJ6e3uTkpJoNNq7d+9Onz6tq6uLpTs5OWlqak6aNElG\nRubWrVtY4pUrV7Zu3Ypdl5SUkEikoYqnp6ebmJj0b4jJZDY3NwcEBEhJSW3cuNHHx0ddXT0gICA7\nO1tXVxeHw8nIyLi4uBgZGVVXV8+bN8/R0VFJSSk+Pn7evHk0Gg2rpKmpiYODIywsjFVtZGSkkJAQ\nHo/Py8vDUp4+fTpnzhwEQWbNmpWWlrZgwYL+Qw9EIpHVH3l5eW5ubgEBgf7b4rN887724KcHAdY4\nBgEWAGAk9R9eQj+3ThxF0b///ltHRycmJmbQ4iiKOjk5rV69eubMmUQiUVtbm3UYzosXL2RlZTU1\nNWfMmHHjxg1s8XhcXJyenl5WVhYOh3v27BmKovfu3VNQUBATE6NQKH/++eeLFy+w4gUFBXp6eiiK\nJiYmfprIkpOTQ6FQsOtJkyb1n93r7OzErisqKmRkZGbPnj0grgLgsyDAGscgwAIAjKT+w0sDbjU1\nNREIBGxaDUXRxsbGjRs36ujolJWVDVN8xowZ9vb29fX10dHROBwuPDy8r6+vrq7OzMzs2LFjPDw8\nFhYWjo6ORkZG/YeX8Hh8ZGQkg8F4+fLlpEmTJCQk3r9/j035dXR0vHnzRltbOyMjA/3fPOCARBaY\n3QM/FARY4xgEWACAkeHg4DBgeAkdenwoNTV17ty5iYmJrFGfQYszGIwJEyZg40YNDQ3YPCPrLOrs\n7GwKhfL06VMKhYLD4VRVVbFxLzqdnpCQgB1lPWPGDBEREWx4Ccs5efJkTU1N1gjZoIksMLsHfigI\nsMYxCLAAACNDWFg4IiKCNbwUGBiIDjE+lJGRMXv27IaGBmzZe0dHx1DFmUymoKBgcnJyd3f3+vXr\nd+zY0b9FGF4C4x0EWOMYBFgAgJGRkJDQf3gJG5oadHzIzMys/3Sevr7+UMVRFI2JiREXF5eSktq+\nfTv2BiILDC+B8Q6Hoii7d34AI6SqqiooKCgkJGT4bB0dHfb29ikpKSPTKwAAAADARqMAAAAAAGwG\nARYAAAAAAJtBgAUAAAAAwGYQYAEAAAAAsBkEWAAAAAAAbAYBFgAAAAAAm0GABQAAAADAZhBgAQAA\nAACwGQRYAAAAAABsBgEWAAAAAACbQYAFAAAAAMBmEGABAAAAALAZBFhgXMrKypKUlGQwGNiPBQUF\n6urqIiIily5dwlJQFFVWVsb9z6pVqxAECQoKwv1bdXV1a2urp6cnmUwWEBDw9fUdtY8EAADgJ8I5\n2h0A4KtlZ2fb2dlpaGhwcHAgCNLe3r58+fLk5OT6+nonJycXFxcEQXJzc4lEYnt7O1aEi4sLQRB3\nd3cs0kIQJCQk5NWrV4qKihYWFjo6OsXFxR8/fpw1a5avry8eD//wAAAA8F3gDwkYZ3Jzc/39/S0t\nLfX09LCU+Ph4CwsLNTU1fX39pqYmJpOJIMj58+dpNBqJRFJRUcnJyeHh4UEQhIeHh5+fn5+f/9Gj\nR5mZmadPn8bhcLm5uTNnzuTg4EhISJg9e/aA6GrTpk2+vr4aGhqBgYE5OTkGBgYiIiKysrKsoTIT\nE5PLly8jCJKQkGBgYECn00f0cQAAABiTYAQLjLJdu3Z9+PCBm5v7szm9vb3r6ur8/f3j4uKMjIxW\nrFiBpcfHx2/atAlBkPfv30tLS+PxeDqd/vz587179xoYGJw6dWr16tVv375l1VNVVeXj43P79m1s\nWMvNzc3Kygq7df/+/QGNUqlUPj6+sLCw9vb27du3R0RETJ8+PScnx8zMzMnJCY/H796928vLi0Qi\nnThxIi0tjZeXly2PBQAAwLgGARYYZR8+fHBycpKVlf1szpqamhMnTsTExOBwuPLyci0tLSz977//\nVlFRQRDkyZMn8vLyCILw8vLm5eVhd728vPz9/ZlMJjY01dzcbGtre+3aNUFBQQRBQkJCSktLq6ur\niUTi5cuXN27cWFVVxWqRTqdXVFS8ePGCRCKZmprm5+crKytjt4hEIlahkZERgUDw8PC4f/++gIAA\nG58MAACA8QsCLDDKuLm5p0yZIicn99mca9euvXfvXmJiIvYjPz9/c3PzpEmTamtrxcTEEARJS0uz\ntbVFECQpKUlVVRWr8/nz5xoaGqyJvwMHDjg7OysoKGA/njlzJiUlZerUqQiCWFpanjx5sn+LJSUl\nampqJBIJQZD8/PwPHz6IiooiCNLZ2YnD4bA8lZWVDQ0NwsLC2C0AAAAAgTVYYBxJT09HURRF0ZMn\nT7q7u6MoKiAggMPhJk2alJubW1JSkpqaigVYcXFxPj4+NBqtpqZm8+bNfn5+WA3V1dV37tzZsmUL\nq87GxsacnBwmk/n69Wt3d/ft27f3b5FKperq6mLXMjIyYWFh3d3dhYWF8vLyJSUlCIK8efPGzs4u\nJSWlu7s7IyNjhB4EAACAMQ8CLDD+FBQUsOYHcTjciRMnVqxY4ebmlpiYOHHiRARB/Pz8nj9/TiKR\nli1btn///vnz52OZt23b5unpiS29woSFhR0+fFhERGT58uVubm4bNmzo3xCVStXR0cGuL168mJiY\nKCIismbNmjNnzujq6tbX11tYWAQHBysrK+/YsePIkSMj8eEBAACMBzgURUe7D+AbVVVVBQUFhYSE\nDJ+to6PD3t4+JSVlZHr1tdauXfvbb7+x5uwAAACAnwCMYAEAAAAAsBkEWAAAAAAAbAYBFgAAAAAA\nm8E2DeNbd3f3y5cvh89Dp9PpdPow2XA4nIyMDGvfAQAAAAB8JwiwxreMjIzXrXQOjs/8f+QQlly3\ne8hjjJ+VFSfFRqurq7O7dwAAAMB/FARY4xuKoqv3HeIhfNfxLJeO+zEYDHZ1CQAAAACwBgsAAAAA\ngM0gwPqp9Pb0rDemtHysZ6W0NDbsX2W7Umv60Q2u9Pa2YRIBAAAAwC4QYP08bl4K93ayqn9b2z8x\n6ugBssL08IdlPATeq2dPDZMIAAAAAHaBAOvnQVacvnyDV/8UFEUL7t/RNl3Mxc1NWfBL4f+1d7dR\nUV0Ho8f3DCqCUZERUEfLIGpKQccF0jFoikZFkhijtcpKpI0QKstA1cTyPHqtkBpb69WIJtFIEY1K\nJJUKKyJPiQQwGb3BGEa7Sn1JwJcYVBBFQEZAZs79MDeEq4CmbsDB/+/TmT17ztlz5MN/nTnOFBxq\naxAAAEjETe52rL6+/ubNm80P/QwT7p5grqs312mHjxBCDPUeWVVR3tYgAACQiCtYdszX19fb2/u+\n02xfcKUoirXFfxVsdRAAAEhBYNmxXr16DR48uJ0JvZ379Hbuc/l8qRDiyoVzmkGD2xoEAAAS8RFh\nd6ZSqQKfCSn6LG/U2IATxoLAKaGtDypKRxy9oqJi06ZNVqu1/Wlffvnlhg0bXF1d25owadKk0NBQ\n2asDAKADEVjd3ILlCYnLYhZOChjhp1/w3wmtDqZv3dQRh/7Xv/5VeObcL174ZfvTXnzSv51na29W\n7dmzh8ACANgXAqu72X/mcsuHLgPd/7gr/a45rQ52BI9hPxkT9PTD7OFGRfnVrz6XtR4AADoHgWXf\nVCpV2jvre/R4qH/HU18VqlSvyloSAAAgsOzb1q1bKyoq2p9jNptXrlyZmJjY5ozoCL1eL3ll97jT\n2Lj42afX/i3LZaC7beTm9Wubfh9z/t//8hlnWLzuHee+/e4d7OhVAQDQEQgs+zZs2LBhw4a1P+fW\nrVv9+/cPCAjonCW1Knv3duPBzFa/ZX5lUup7K5b+/f3Nv/mvVfcOzliwsIuWDADAf46vaUBn4Fvm\nAQCPFa5g4Qe+vr5OTk6y9lZbW+sf+qJtm2+ZBwA8Vggs/MDd3b2goEDW3vLy8rYf+KT9OXzLPACg\nW+IjQnQNvmUeANCNEVjoGs1fKK8oyr3fMt9yEAAAu0NgocssWJ5w9sRXCycF1Ny4PnfRknYGAQCw\nL9yDhc7zH3zL/G2zuTNWBgCAVAQWOopGozn6PweO593nPvf2WRVr+JzZspYEAEDnILDQUcaOHfvt\nN2fuOy06OnrZsmWjRo3qhCUBANA5uAcLAABAMgILAABAMgILAABAMgILAABAMgILAABAMgILAABA\nMgILAABAMgILAABAMgILAABAMgILAABAMgILAABAMgILAABAMgILAABAMgILAABAMgIL9qSmpiYq\nKsrDw2PIkCGrV6+2DSqK4ufnp/peRESEEKK6unrJkiWenp4uLi4JCQntzGx1EACAh9GjqxcA/Ajh\n4eEajaa4uLi2tvanP/3pH/7wB7VaffToUY1GU1tba5vTs2dPIcT8+fOfeuqpoqKiysrKMWPGJCQk\ntDWz1UEAAB4GV7BgN5qamiIjI7ds2WKxWLKzswMDA9VqtRAiKSmprq5Oq9WOHj3aaDQ6OjoKIY4e\nPerr6+vg4JCenj5u3Lh2ZrY62Cw2NjYhIcHf3z8xMdFoNAYHB7u5uXl5ee3evds2YerUqampqUKI\n9PT04OBgs9ncyacFAPAI4goWupjVaq2pqamqqrrvzAEDBsyaNctkMgUEBAghcnJyhBBms7mkpGTl\nypXBwcGbN2+OiooqKysTQkRGRs6ePdv2wry8vLZmtvXyZoWFhX369ElOTq6trY2Li0tJSfHx8TEa\njaGhoeHh4Wq1esWKFUuXLtVqtW+//XZOTo6zs7P0UwQAsDsqRVG6eg3oWLdu3XrppZeysrLuO3Py\n5MkFBQWdsKSWoqOjr169etd1o1atW7fOy8tLUZTq6uqdO3du3Ljx0qVLLSdUVVVptdpbt25t27Zt\n//79W7du1Wg0qampSUlJp0+fbnWm7cpWW4Nms9nV1bW0tFSr1YaEhOTm5jZP1mg0lZWVQghFUQwG\ng9lszsvL8/DweJhTAQDoNriChS6WlJT0gDPnz5+/fft2JycnFxeXF1988d133xVCZGZm6vX64cOH\nCyFKSkr8/f3VavV7772XlZXl7e0thJg1a9amTZvamtnqYPMRTSbT2LFjtVqtEOLYsWPl5eXu7u5C\niNu3b6tUKtucU6dOXbt2beDAgbanAAAQ3IMFO3Lo0KG0tDSLxVJWVhYTE7N48WIhxEcffRQfH19X\nV/ftt98uXrz4rbfeEkJcv37daDRardYLFy4sXLgwLi6urZmtDjYrLCwMCgqybet0uuTk5IaGhuPH\nj48cOdJkMgkhLl68GBYWlpWV1dDQkJ+f38knBADw6FLQ3dXW1s6YMeNBZk6aNKmjF/Mw0tPTdTrd\ngAED9Hp9SkqK1WpVFOXs2bMGg6F///4BAQE5OTm2mR9//PGIESNcXV0DAgLS0tLamdnqYLM5c+bs\n27fPtm0ymfz9/fv27avX6zMzMxVFKS8vHz169OHDhxVF2bNnz5QpUzrlNAAA7AD3YHV/j/g9WAAA\ndD98RAgAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZ\ngQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUAACAZgQUA\nACAZgQUAACBZj65eADpVfX39mjVrmpqaWn323Llzy5cvv2vQ29v7t7/9bccvDQCA7oPAerxUVVVl\nf3b05dfvriib2J9Pvndwz6Y/EVgAAPwoBNZj54n+Lt6+Yx58voODQ8ctBgCAbol7sAAAACQjsB5T\nN69f2/T7mIig0a9OHLv1D8tu192yjf/z6Gevv/DMfP+Rcb+cfrroy65dJAAAdorAekxt+V9vPNHf\nJTErf33GJ2q1et97G4UQitW68Y1Fc19buuvYqSm/emnT72O6epkAANglAusxdep44fO/iXLRuLm6\ne7ywILow93+EEEKlcn6ib11tbf1tc73Z3Kdfv65eJgAAdomb3Ls/BweHU6dOTZs2TQjR0NAg+mmE\nEN6+Yw7s2DYr6rU7jY2ZyVuqKq4KIVQqVcyfNya8MndbfJxDj56bsvK7eOkAANgnrmB1f05OTqWl\npbm5ubm5uX/7299sg79bt7ny6uU3Xpy6dtEr2uHeT/QfIIS4VVOd+MZrb2zctt144plfhu34c3yX\nLhwAAHvFFazHlKOT83+9s72no6MQ4otPDg71HimE+PpkkcdPPCc8N1MIMf2lV1aFz+7iVQIAYJ+4\ngvWY2pu4dlNcbF1NTeWVy+lbEkNffkUI8ZORT5adK/nn//m84bb5swN/93zSp6uXCQCAXeIK1mMq\nfNnKLSvfiH4mcMBAt+d/EzU+5HkhxMDB2ti1m3b+OeHa5e+8fHxj/rSxq5cJAIBdIrAeU0/0d/nv\n93bcOx74TEjgMyGdvx4AALoTAuvx4uzsbK2qeD/2N60+W1RUFBAQcNegTqfr8GUBANC9qBRF6eo1\n4FExefLkgoKCrl4FAAB2j5vcAQAAJCOwAAAAJCOwAAAAJCOwAAAAJCOwAAAAJCOwAAAAJCOwAAAA\nJCOwAAAAJCOwYE9qamqioqI8PDyGDBmyevVq26CiKH5+fqrvRURE2Ma/+OKLoKAgV1fX4cOHHzhw\noK2XCyFu3rwZFhY2ePDg8ePH79jRyi8IAQDwoxBYsCfh4eEWi6W4uPjIkSNr1qyxWq1CiKNHj2o0\nmtrvbdu2TQhx5cqVGTNmxMfHV1ZWrlq1Kjo6uq2XW63WmTNnTp8+vaysbMOGDQcPHuza9wgA6AYI\nLNiNpqamyMjILVu2WCyW7OzswMBAtVothEhKSqqrq9NqtaNHjzYajY6OjkIIo9G4YMGC0NBQtVo9\nduxYBweHtl6en5/v6OgYGRmpVqsnTpyYkZHR8qCxsbEJCQn+/v6JiYlGozE4ONjNzc3Ly2v37t22\nCVOnTk1NTRVCpKenBwcHm83mzj4vAIBHDz/2jI6ydOlSs9msUqmk7E2tVickJMyaNctkMtl+kTon\nJ0cIYTabS0pKVq5cGRwcvHnz5qioqLKyMiHEvHnz5s2bJ4QoKSmJiYlZv359jx497n25EGLXrl06\nnc7Pz+/q1avr169v/oTRprCwsE+fPsnJybW1tXFxcSkpKT4+PkajMTQ0NDw8XK1Wr1ixYunSpVqt\n9u23387JyXF2dpbyfgEAdo0fe8YP5P7Y89SpU995553evXvL2uGwYcN69uypKEp1dfXOnTs3btx4\n6dKllhOqqqq0Wu2tW7dsl6Zu3LgRHx9vMpm2bds2ZswY25x7X+7r66vX6zdv3pybm7to0aKqqirb\ny4UQZrPZ1dW1tLRUq9WGhITk5uY2H0uj0VRWVtp2aDAYzGZzXl6eh4eHrDcLALBrXMFCR+nVq5eX\nl5eTk5OsHc6fP3/79u1OTk4uLi4vvvjiu+++K4TIzMzU6/XDhw8XQpSUlPj7+9vyKDs7e926da+/\n/vq7775ru4rW6sutVuu5c+cKCgrc3NxCQkIcHR1bXnIzmUxjx47VarVCiGPHjpWXl7u7uwshbt++\n3Tzt1KlT165dGzhwoO0pAAAE92DBjhw6dCgtLc1isZSVlcXExCxevFgI8dFHH8XHx9fV1X377beL\nFy9+6623hBAFBQVvvvlmRkbGtGnT6urq6urq2nq5SqVycnI6duxYY2NjfHx8REREy8AqLCwMCgqy\nbet0uuTk5IaGhuPHj48cOdJkMgkhLl68GBYWlpWV1dDQkJ+f3/nnBADwiFKA702aNEni3p599lmz\n2Sxxh+np6TqdbsCAAXq9PiUlxWq1Kopy9uxZg8HQv3//gICAnJwc28zQ0NCWf+QTJ05s6+WKouzd\nu3fQoEFDhw6Ni4urr69vecQ5c+bs27fPtm0ymfz9/fv27avX6zMzMxVFKS8vHz169OHDhxVF2bNn\nz5QpUyS+WQCAXeMeLPxA7j1Yzz333P79+yV+RAgAgL3gI0IAAADJCCwAAADJCCwAAADJCCwAAADJ\nCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwA\nAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJ\nCCwAAADJCCwAAADJCCwAAADJenT1AmB/Tp48+f777993WnFxcWxsbI8ebf6NzZ49OzQ0VOrSAAB4\nJBBY+NFOnjxZ3999fMjz7U8bOzeinWcvnDl1+PBhAgsA0C0RWPhP9HPVeAzzfJg91Ny4XvXd17LW\nAwDAI4V7sAAAACQjsPCw7jQ2LppiuFlZ0Txy8/q1NyPmvfJzn7+8tsBcW9POIAAA3RKBhYeSvXv7\nqvDZFWWXWg7u+ssfPUf5bD/yT0cn57+/v7mdQQAAuiUCCw/F80mfX722tOWIoihf5n0yPuT5nr16\nGaY9e7zgUFuDAAB0V9zkjh9UVFRMmzbtvtMuX778zK9/a9v2M0y469l6c129uU47fIQQYqj3yKqK\n8rYGAQDorggs/ODf//73g0z74IMPiq7ebH+OSqUSQiiKYrVY2h8EAKD74SNCSNbbuU9v5z6Xz5cK\nIa5cOKcZNLitQQAAuisCC5KpVKrAZ0KKPstTFOWEsSBwSmhbgwAAdFcEFuRbsDzh7ImvFk4KqLlx\nfe6iJe0MAgDQLXEPFiTYf+Zyy4cuA93/uCv9rjmtDgIA0C0RWPjRdDrdhx/+6UxedvvTTpw4odfr\n1eo2r5L++te/lr00AAAeCSpFUbp6Deiennvuuf379zs5OXX1QgAA6GzcgwUAACAZgQV7UlNTExUV\n5eHhMWTIkNWrV9sGq6urlyxZ4unp6eLikpCQ0M5MIcTNmzfDwsIGDx48fvz4HTt22Aa/+OKLoKAg\nV1fX4cOHHzhwoJPfFACg+yGwYE/Cw8MtFktxcfGRI0fWrFljtVqFEPPnz3d3dy8qKiosLFy7dq1t\nsNWZVqt15syZ06dPLysr27Bhw8GDB4UQV65cmTFjRnx8fGVl5apVq6Kjo7v2PQIAugECC3ajqakp\nMjJyy5YtFoslOzs7MDDQdgf90aNHfX19HRwc0tPTx40bp1ar25qZn5/v6OgYGRmpVqsnTpyYkZEh\nhDAajQsWLAgNDVWr1WPHjnVwcGh50NjY2ISEBH9//8TERKPRGBwc7Obm5uXltXv3btuEqVOnpqam\nCiHS09ODg4PNZnNnnxcAwCNIATrGs88+azab7zutrq7uxoOxzS8qKrL96ebk5NhG3njjjea/57y8\nvOY93zszPDw8KirK19dXo9Hs2LHjrpV88803Tz311N69e1sOBgQE/OIXv/jqq68KCgoMBkNxcbHF\nYjl8+HDv3r0tFouiKJ9++qmfn19+fr7BYKiqqvpPzxYAoFvhfxGiozzg/yIMCQnp2bNnnz597rvD\ndevWeXl5KYpSXV29c+fOjRs3Xrp0aevWrfv379+6datGo0lNTU1KSjp9+rRt/l0zhRC+vr56vX7z\n5s25ubmLFi2qqqqyXdm6ceNGfHy8yWTatm3bmDFjmo9oNptdXV1LS0u1Wm1ISEhubm7zUxqNprKy\n0nYUg8FgNpvz8vI8PDx+/HkCAHRHXdt36MYe8ArWwoULz549+yA7fPnll5t3WFpaaostHx+fkpIS\n2+D58+c9PT3bmmmxWHr37l1eXq4oyrVr19zc3KxWq6IoBw8efPrppzMyMmwPWzIajQaDwbbdr18/\n22sVRTGbzbdvuP9oqwAACixJREFU37ZtFxcX63S6cePG3ftyAMBji3uwYDcOHTqUlpZmsVjKyspi\nYmIWL14shLh+/brRaLRarRcuXFi4cGFcXFxbM1UqlZOT07FjxxobG+Pj4yMiIlQqVUFBwZtvvpmR\nkTFt2rS6urq6urqWRywsLAwKCrJt63S65OTkhoaG48ePjxw50mQyCSEuXrwYFhaWlZXV0NCQn5/f\n2WcEAPDI6urCQ7cl/QpWenq6TqcbMGCAXq9PSUmxXTH6+OOPR4wY4erqGhAQkJaWZhtsdaaiKHv3\n7h00aNDQoUPj4uLq6+sVRQkN/f9+dnrixIktjzhnzpx9+/bZtk0mk7+/f9++ffV6fWZmpqIo5eXl\no0ePPnz4sKIoe/bsmTJlyoOfHABA98Y9WOgoD3gPVnR09LJly0aNGtU5qwIAoBPwESEAAIBkBBYA\nAIBkBBYAAIBkBBYAAIBkBBYAAIBkBBYAAIBkBBYAAIBkBBYAAIBkBBYAAIBkBBYAAIBkBBYAAIBk\nBBYAAIBkBBYAAIBkBBYAAIBkKkVRunoN6J4mTJhw7rvL7kOHPcxO7jQ2Pu2vT0pKkrUqAAA6QY+u\nXgC6raampkmz5s597fWH2cmNivKCbetlLQkAgM7BR4QAAACSEVgAAACSEVjoJHcaGxdNMdysrGge\nuXn92psR8175uc9fXltgrq1pZxAAAPtCYKEzZO/evip8dkXZpZaDu/7yR89RPtuP/NPRyfnv729u\nZxAAAPtCYKGj1NXVNW97Punzq9eWtnxWUZQv8z4ZH/J8z169DNOePV5wqK1BAADsDv+LEB1l4sSJ\n1d9v+xkm3PVsvbmu3lynHT5CCDHUe2RVRXlbgwAA2B2uYKGjzJ07975zVCqVEEJRFKvF0v4gAAB2\nhMBC1+jt3Ke3c5/L50uFEFcunNMMGtzWIAAAdofAQtdQqVSBz4QUfZanKMoJY0HglNC2BgEAsDsE\nFrrMguUJZ098tXBSQM2N63MXLWlnEAAA+8JvEaKj5OXlrfjz/w6YNPVhdlJXU2O5fO7DDz+UtSoA\nADoBgYWOUlVV9emnn9532l//+tcXXnhh8OA2b7fy9vb29/eXujQAADoWgYUuFh0dvWzZslGjRnX1\nQgAAkIZ7sAAAACQjsAAAACQjsAAAACQjsAAAACQjsAAAACQjsAAAACQjsAAAACQjsAAAACQjsAAA\nACQjsAAAACQjsAAAACQjsAAAACQjsAAAACQjsAAAACQjsGBPampqoqKiPDw8hgwZsnr1atvg119/\n/fzzz7u7u3t7ex84cKB58s2bN8PCwgYPHjx+/PgdO3Y0jzc2NgYFBZ05c8b2sLq6esmSJZ6eni4u\nLgkJCZ35dgAA3RWBBXsSHh5usViKi4uPHDmyZs0aq9V648aN0NDQ3/3ud1euXFmzZs2CBQusVqsQ\nwmq1zpw5c/r06WVlZRs2bDh48KAQoqmpacWKFQaD4cyZM6NGjbLtc/78+e7u7kVFRYWFhWvXrrW9\nHACAh0FgwW40NTVFRkZu2bLFYrFkZ2cHBgaq1eqMjIzg4ODQ0FAHB4fJkyc3NDTYJufn5zs6OkZG\nRqrV6okTJ2ZkZAgh7ty5M2HChKlTpxoMBrX6//3xHz161NfX18HBIT09fdy4cc3jNrGxsQkJCf7+\n/omJiUajMTg42M3NzcvLa/fu3bYJU6dOTU1NFUKkp6cHBwebzebOOyMAgEdVj65eAB53ly5dio+P\n79+/f/vT1Gp1QkLCrFmzTCZTQECAECInJ0cI0djYePLkye+++66xsTE2NtbHx8dWSLt27dLpdH5+\nflevXl2/fn1ERIQQwsnJacaMGZ9//vlTTz3VvOfIyMjZs2fbtvPy8u46bmFhYZ8+fZKTk2tra+Pi\n4lJSUnx8fIxGY2hoaHh4uFqtXrFixdKlS7Va7dtvv52Tk+Ps7Czv3AAA7JVKUZSuXgMea9988011\ndbVKpWp/Wo8ePX72s5/17NlTUZTq6uqdO3du3Ljx0qVLtbW18+bN+/zzz/V6vU6n69Wr1wcffCCE\n8PX11ev1mzdvzs3NXbRoUVVVVfOlqaeffnrVqlUhISFCiK1bt+7fv3/r1q0ajSY1NTUpKen06dPN\nBzWbza6urqWlpVqtNiQkJDc3t/kpjUZTWVkphFAUxWAwmM3mvLw8Dw8P2acHAGCfFMBOvPzyy2az\n2bZdWlrq5eWlKMqlS5csFouiKHfu3HnyySf/8Y9/KIpisVh69+5dXl6uKMq1a9fc3NysVqvthQ0N\nDU5OTlVVVbaHPj4+JSUltu3z5897enq2PKLRaDQYDLbtfv362XaoKIrZbL59+7Ztu7i4WKfTjRs3\nrvkQAABwDxbsxqFDh9LS0iwWS1lZWUxMzOLFi4UQfn5+WVlZDQ0Ny5cv9/Lysl2XUqlUTk5Ox44d\na2xsjI+Pj4iIaL5CdvLkSS8vLxcXF9vD69evG41Gq9V64cKFhQsXxsXFtTxiYWFhUFCQbVun0yUn\nJzc0NBw/fnzkyJEmk0kIcfHixbCwMNsC8vPzO+1UAAAedV1deMCDSk9P1+l0AwYM0Ov1KSkptitG\nO3fudHNzGzJkyKJFi6qrq5sn7927d9CgQUOHDo2Li6uvr28eT0xMfPXVV5sffvzxxyNGjHB1dQ0I\nCEhLS7vrKtScOXP27dtn2zaZTP7+/n379tXr9ZmZmYqilJeXjx49+vDhw4qi7NmzZ8qUKR357gEA\n9oR7sAAAACTjI0IAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJ\nCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwA\nAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJ\nCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwA\nAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJ\nCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwA\nAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJ\nCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwA\nAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJ\nCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJCCwAAADJ/i+A\nyQzMQhQ4mgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -h 800 -w 800\n",
    "library(ape)\n",
    "tre <- read.tree(\"empirical_10/RAxML_bipartitions.empirical_10_m4\")\n",
    "ltre <- ladderize(tre)\n",
    "\n",
    "par(mfrow=c(1,2))\n",
    "plot(ltre, use.edge.length=F)\n",
    "nodelabels(ltre$node.label)\n",
    "\n",
    "plot(ltre, type='u')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get phylo distances (GTRgamma dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 0.02792806\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "mean(cophenetic.phylo(ltre))\n"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}