emp_nb_Finches.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook 6: \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 6\n",
    "### Data set 6 (Finches)\n",
    "### Authors: DaCosta & Sorenson (2016)\n",
    "### Data Location: SRP059199"
   ]
  },
  {
   "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",
    "+ Project SRA: SRP059199\n",
    "+ BioProject ID: PRJNA285779\n",
    "+ Biosample numbers: SAMN03753600 - SAMN03753623\n",
    "+ Runs: SRR2053224 -- SRR2053247\n",
    "+ SRA link: http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP059199"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## make a new directory for this analysis\n",
    "mkdir -p empirical_6/fastq/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For each ERS (individuals) get all of the ERR (sequence file accessions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    BioSample_s Library_Name_s  MBases_l  MBytes_l  \\\n",
      "0  SAMN03753617         90.o31        77        49   \n",
      "1  SAMN03753619         90.o38        84        54   \n",
      "2  SAMN03753618         90.o50        83        55   \n",
      "3  SAMN03753613         90.o53        74        46   \n",
      "4  SAMN03753601           A056        70        44   \n",
      "\n",
      "                        Organism_s       Run_s SRA_Sample_s  \\\n",
      "0                Vidua hypocherina  SRR2053224    SRS954896   \n",
      "1                      Vidua regia  SRR2053225    SRS954895   \n",
      "2                   Vidua fischeri  SRR2053226    SRS954894   \n",
      "3                 Vidua orientalis  SRR2053227    SRS954893   \n",
      "4  Lagonosticta senegala rhodopsis  SRR2053228    SRS954892   \n",
      "\n",
      "                                  Sample_Name_s isolate_s Assay_Type_s  \\\n",
      "0              Vidua hypocherina, sample 90.o31    90.o31        OTHER   \n",
      "1                    Vidua regia, sample 90.o38    90.o38        OTHER   \n",
      "2                 Vidua fischeri, sample 90.o50    90.o50        OTHER   \n",
      "3               Vidua orientalis, sample 90.o53    90.o53        OTHER   \n",
      "4  Lagonosticta senegala rhodopsis, sample A056      A056        OTHER   \n",
      "\n",
      "    ...       LoadDate_s Platform_s ReleaseDate_s SRA_Study_s  \\\n",
      "0   ...     Jun 05, 2015   ILLUMINA  Jun 24, 2015   SRP059199   \n",
      "1   ...     Jun 05, 2015   ILLUMINA  Jun 24, 2015   SRP059199   \n",
      "2   ...     Jun 05, 2015   ILLUMINA  Jun 24, 2015   SRP059199   \n",
      "3   ...     Jun 05, 2015   ILLUMINA  Jun 24, 2015   SRP059199   \n",
      "4   ...     Jun 05, 2015   ILLUMINA  Jun 24, 2015   SRP059199   \n",
      "\n",
      "              age_s  g1k_analysis_group_s  g1k_pop_code_s sex_s  \\\n",
      "0  after hatch year        <not provided>  <not provided>  male   \n",
      "1  after hatch year        <not provided>  <not provided>  male   \n",
      "2  after hatch year        <not provided>  <not provided>  male   \n",
      "3  after hatch year        <not provided>  <not provided>  male   \n",
      "4  after hatch year        <not provided>  <not provided>  male   \n",
      "\n",
      "         source_s tissue_s  \n",
      "0  <not provided>   muscle  \n",
      "1  <not provided>   muscle  \n",
      "2  <not provided>   muscle  \n",
      "3  <not provided>   muscle  \n",
      "4  <not provided>   muscle  \n",
      "\n",
      "[5 rows x 27 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_6_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": 4,
   "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": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for ID, SRR in zip(indata.Library_Name_s, indata.Run_s):\n",
    "    wget_download(SRR, \"empirical_6/fastq/\", ID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read 889818 spots for empirical_6/fastq/90.o31.sra\n",
      "Written 889818 spots for empirical_6/fastq/90.o31.sra\n",
      "Read 1012887 spots for empirical_6/fastq/90.o38.sra\n",
      "Written 1012887 spots for empirical_6/fastq/90.o38.sra\n",
      "Read 998161 spots for empirical_6/fastq/90.o50.sra\n",
      "Written 998161 spots for empirical_6/fastq/90.o50.sra\n",
      "Read 862619 spots for empirical_6/fastq/90.o53.sra\n",
      "Written 862619 spots for empirical_6/fastq/90.o53.sra\n",
      "Read 792799 spots for empirical_6/fastq/A056.sra\n",
      "Written 792799 spots for empirical_6/fastq/A056.sra\n",
      "Read 833981 spots for empirical_6/fastq/A081.sra\n",
      "Written 833981 spots for empirical_6/fastq/A081.sra\n",
      "Read 809186 spots for empirical_6/fastq/A082.sra\n",
      "Written 809186 spots for empirical_6/fastq/A082.sra\n",
      "Read 939649 spots for empirical_6/fastq/A089.sra\n",
      "Written 939649 spots for empirical_6/fastq/A089.sra\n",
      "Read 1020850 spots for empirical_6/fastq/A107.sra\n",
      "Written 1020850 spots for empirical_6/fastq/A107.sra\n",
      "Read 1001743 spots for empirical_6/fastq/A145.sra\n",
      "Written 1001743 spots for empirical_6/fastq/A145.sra\n",
      "Read 729322 spots for empirical_6/fastq/A147.sra\n",
      "Written 729322 spots for empirical_6/fastq/A147.sra\n",
      "Read 834688 spots for empirical_6/fastq/A167.sra\n",
      "Written 834688 spots for empirical_6/fastq/A167.sra\n",
      "Read 709824 spots for empirical_6/fastq/A177.sra\n",
      "Written 709824 spots for empirical_6/fastq/A177.sra\n",
      "Read 935241 spots for empirical_6/fastq/A1794.sra\n",
      "Written 935241 spots for empirical_6/fastq/A1794.sra\n",
      "Read 898117 spots for empirical_6/fastq/A204.sra\n",
      "Written 898117 spots for empirical_6/fastq/A204.sra\n",
      "Read 992534 spots for empirical_6/fastq/A248.sra\n",
      "Written 992534 spots for empirical_6/fastq/A248.sra\n",
      "Read 1034739 spots for empirical_6/fastq/A328.sra\n",
      "Written 1034739 spots for empirical_6/fastq/A328.sra\n",
      "Read 927116 spots for empirical_6/fastq/JMD1249.sra\n",
      "Written 927116 spots for empirical_6/fastq/JMD1249.sra\n",
      "Read 1031674 spots for empirical_6/fastq/JMD1285.sra\n",
      "Written 1031674 spots for empirical_6/fastq/JMD1285.sra\n",
      "Read 956686 spots for empirical_6/fastq/MDS004.sra\n",
      "Written 956686 spots for empirical_6/fastq/MDS004.sra\n",
      "Read 1220879 spots for empirical_6/fastq/MDS054.sra\n",
      "Written 1220879 spots for empirical_6/fastq/MDS054.sra\n",
      "Read 1064587 spots for empirical_6/fastq/MDS065.sra\n",
      "Written 1064587 spots for empirical_6/fastq/MDS065.sra\n",
      "Read 741675 spots for empirical_6/fastq/NKK511.sra\n",
      "Written 741675 spots for empirical_6/fastq/NKK511.sra\n",
      "Read 1079701 spots for empirical_6/fastq/SAR6894.sra\n",
      "Written 1079701 spots for empirical_6/fastq/SAR6894.sra\n",
      "Read 22318476 spots total\n",
      "Written 22318476 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_6/fastq/ empirical_6/fastq/*.sra\n",
    "\n",
    "## remove .sra files\n",
    "rm empirical_6/fastq/*.sra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 1668628\n",
      "-rw-rw-r-- 1 deren deren 69481962 Nov 23 13:13 90.o31.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 75367262 Nov 23 13:13 90.o38.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 76680615 Nov 23 13:13 90.o50.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 66086600 Nov 23 13:13 90.o53.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 62055035 Nov 23 13:13 A056.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 64828556 Nov 23 13:13 A081.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 63527984 Nov 23 13:13 A082.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 73811356 Nov 23 13:13 A089.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 79426948 Nov 23 13:13 A107.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 78226977 Nov 23 13:13 A145.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 56867082 Nov 23 13:13 A147.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 67136950 Nov 23 13:13 A167.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 54625508 Nov 23 13:13 A177.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 68861314 Nov 23 13:13 A1794.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 70348288 Nov 23 13:13 A204.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 72142052 Nov 23 13:14 A248.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 73305922 Nov 23 13:14 A328.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 70371041 Nov 23 13:14 JMD1249.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 77220081 Nov 23 13:14 JMD1285.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 76259313 Nov 23 13:14 MDS004.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 87175055 Nov 23 13:14 MDS054.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 83991953 Nov 23 13:14 MDS065.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 57656728 Nov 23 13:14 NKK511.fastq.gz\n",
      "-rw-rw-r-- 1 deren deren 83176262 Nov 23 13:14 SAR6894.fastq.gz\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "ls -l empirical_6/fastq/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Merge technical replicates\n",
    "This study includes several technical replicates per sequenced individuals, which we combine into a single file for each individual here. "
   ]
  },
  {
   "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": 11,
   "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_6/           ## 1. working directory ' params.txt\n",
    "sed -i '/## 6. /c\\CCTGCAGG,AATTC         ## 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_6_m4          ## 14. output name      ' params.txt\n",
    "sed -i '/## 18./c\\empirical_6/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": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==** parameter inputs for pyRAD version 3.0.63  **======================== affected step ==\r\n",
      "empirical_6/           ## 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",
      "CCTGCAGG,AATTC         ## 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_6_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_6/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": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "sed -i '/## 12./c\\2                    ## 12. MinCov           ' params.txt\n",
    "sed -i '/## 14./c\\empirical_6_m2       ## 14. output name      ' params.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count        24.000000\n",
      "mean     609459.416667\n",
      "std      114564.902701\n",
      "min      472250.000000\n",
      "25%      531561.250000\n",
      "50%      593659.500000\n",
      "75%      648491.250000\n",
      "max      961219.000000\n",
      "Name: passed.total, dtype: float64\n",
      "\n",
      "most raw data in sample:\n",
      "18    JMD1285\n",
      "Name: sample , dtype: object\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "## read in the data\n",
    "s2dat = pd.read_table(\"empirical_6/stats/s2.rawedit.txt\", header=0, nrows=25)\n",
    "\n",
    "## print summary stats\n",
    "print s2dat[\"passed.total\"].describe()\n",
    "\n",
    "## find which sample has the most raw data\n",
    "maxraw = s2dat[\"passed.total\"].max()\n",
    "print \"\\nmost raw data in sample:\"\n",
    "print s2dat['sample '][s2dat['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    24.000000\n",
      "mean     38.139250\n",
      "std       5.881965\n",
      "min      28.454000\n",
      "25%      33.512000\n",
      "50%      37.240000\n",
      "75%      41.636500\n",
      "max      53.215000\n",
      "Name: dpt.me, dtype: float64\n",
      "\n",
      "summary of std\n",
      "==================\n",
      "count     24.000000\n",
      "mean     112.613125\n",
      "std       46.878437\n",
      "min       70.168000\n",
      "25%       85.451500\n",
      "50%       96.007500\n",
      "75%      121.462750\n",
      "max      241.616000\n",
      "Name: dpt.sd, dtype: float64\n",
      "\n",
      "summary of proportion lowdepth\n",
      "==================\n",
      "count    24.000000\n",
      "mean      0.458336\n",
      "std       0.074790\n",
      "min       0.308587\n",
      "25%       0.416855\n",
      "50%       0.463149\n",
      "75%       0.526700\n",
      "max       0.572036\n",
      "dtype: float64\n",
      "\n",
      "highest coverage in sample:\n",
      "11    A167\n",
      "Name: taxa, dtype: object\n"
     ]
    }
   ],
   "source": [
    "## read in the s3 results\n",
    "s6dat = pd.read_table(\"empirical_6/stats/s3.clusters.txt\", header=0, nrows=25)\n",
    "\n",
    "## print summary stats\n",
    "print \"summary of means\\n==================\"\n",
    "print s6dat['dpt.me'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of std\\n==================\"\n",
    "print s6dat['dpt.sd'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of proportion lowdepth\\n==================\"\n",
    "print pd.Series(1-s6dat['d>5.tot']/s6dat[\"total\"]).describe()\n",
    "\n",
    "## find which sample has the greatest depth of retained loci\n",
    "max_hiprop = (s6dat[\"d>5.tot\"]/s6dat[\"total\"]).max()\n",
    "print \"\\nhighest coverage in sample:\"\n",
    "print s6dat['taxa'][s6dat['d>5.tot']/s6dat[\"total\"]==max_hiprop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "41.6207345663 129.154402626\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "## print mean and std of coverage for the highest coverage sample\n",
    "with open(\"empirical_6/clust.85/A167.depths\", 'rb') as indat:\n",
    "    depths = np.array(indat.read().strip().split(\",\"), dtype=int)\n",
    "    \n",
    "print depths.mean(), depths.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 0x7ffa54764090>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<div align=\"center\" class=\"toyplot\" id=\"t8f815267e1a6455b98768789bde1239b\"><svg height=\"300.0px\" id=\"t79044aedf6e04cec92ecfbf13d9f259d\" 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=\"tf9d5e7238f5149509ba5315128fd73ae\"><clipPath id=\"t970ab55cdc724f8da83e1693449fe809\"><rect height=\"200.0\" width=\"250.0\" x=\"50.0\" y=\"50.0\"></rect></clipPath><g class=\"toyplot-coordinate-events\" clip-path=\"url(#t970ab55cdc724f8da83e1693449fe809)\" 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=\"t675875093dee4a2d9d32170c1807f2a4\" 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=\"171.63\" 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=\"68.370000000000005\"></rect><rect class=\"toyplot-Datum\" height=\"69.389999999999986\" 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=\"170.61000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"39.150000000000006\" 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=\"200.84999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"23.309999999999974\" 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=\"216.69000000000003\"></rect><rect class=\"toyplot-Datum\" height=\"16.379999999999995\" 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=\"223.62\"></rect><rect class=\"toyplot-Datum\" height=\"14.489999999999981\" 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=\"225.51000000000002\"></rect><rect class=\"toyplot-Datum\" height=\"11.969999999999999\" 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=\"228.03\"></rect><rect class=\"toyplot-Datum\" height=\"11.609999999999985\" 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=\"228.39000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"11.339999999999975\" 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=\"228.66000000000003\"></rect><rect class=\"toyplot-Datum\" height=\"8.3700000000000045\" 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.63\"></rect><rect class=\"toyplot-Datum\" height=\"10.530000000000001\" 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=\"229.47\"></rect><rect class=\"toyplot-Datum\" height=\"9.6299999999999955\" 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=\"230.37\"></rect><rect class=\"toyplot-Datum\" height=\"11.069999999999993\" 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=\"228.93000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"9.3600000000000136\" 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=\"230.63999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"11.25\" 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=\"228.75\"></rect><rect class=\"toyplot-Datum\" height=\"11.699999999999989\" 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=\"228.30000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"8.0100000000000193\" 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=\"231.98999999999998\"></rect><rect class=\"toyplot-Datum\" height=\"9.8099999999999739\" 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=\"230.19000000000003\"></rect><rect class=\"toyplot-Datum\" height=\"10.97999999999999\" 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=\"229.02000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"10.97999999999999\" 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=\"229.02000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"10.890000000000015\" 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=\"229.10999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"11.160000000000025\" 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=\"228.83999999999997\"></rect><rect class=\"toyplot-Datum\" height=\"11.879999999999995\" 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=\"228.12\"></rect><rect class=\"toyplot-Datum\" height=\"10.800000000000011\" 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=\"229.19999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"11.160000000000025\" 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=\"228.83999999999997\"></rect><rect class=\"toyplot-Datum\" height=\"8.6400000000000148\" 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=\"231.35999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"10.079999999999984\" 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=\"229.92000000000002\"></rect><rect class=\"toyplot-Datum\" height=\"20.610000000000014\" 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=\"219.38999999999999\"></rect></g></g><g class=\"toyplot-mark-BarMagnitudes\" id=\"tff5b2babd5bc42c5b2e553142071ee10\" 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=\"14.489999999999981\" 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=\"225.51000000000002\"></rect><rect class=\"toyplot-Datum\" height=\"11.969999999999999\" 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=\"228.03\"></rect><rect class=\"toyplot-Datum\" height=\"11.609999999999985\" 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=\"228.39000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"11.339999999999975\" 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=\"228.66000000000003\"></rect><rect class=\"toyplot-Datum\" height=\"8.3700000000000045\" 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.63\"></rect><rect class=\"toyplot-Datum\" height=\"10.530000000000001\" 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=\"229.47\"></rect><rect class=\"toyplot-Datum\" height=\"9.6299999999999955\" 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=\"230.37\"></rect><rect class=\"toyplot-Datum\" height=\"11.069999999999993\" 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=\"228.93000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"9.3600000000000136\" 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=\"230.63999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"11.25\" 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=\"228.75\"></rect><rect class=\"toyplot-Datum\" height=\"11.699999999999989\" 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=\"228.30000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"8.0100000000000193\" 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=\"231.98999999999998\"></rect><rect class=\"toyplot-Datum\" height=\"9.8099999999999739\" 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=\"230.19000000000003\"></rect><rect class=\"toyplot-Datum\" height=\"10.97999999999999\" 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=\"229.02000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"10.97999999999999\" 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=\"229.02000000000001\"></rect><rect class=\"toyplot-Datum\" height=\"10.890000000000015\" 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=\"229.10999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"11.160000000000025\" 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=\"228.83999999999997\"></rect><rect class=\"toyplot-Datum\" height=\"11.879999999999995\" 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=\"228.12\"></rect><rect class=\"toyplot-Datum\" height=\"10.800000000000011\" 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=\"229.19999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"11.160000000000025\" 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=\"228.83999999999997\"></rect><rect class=\"toyplot-Datum\" height=\"8.6400000000000148\" 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=\"231.35999999999999\"></rect><rect class=\"toyplot-Datum\" height=\"10.079999999999984\" 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=\"229.92000000000002\"></rect><rect class=\"toyplot-Datum\" height=\"20.610000000000014\" 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=\"219.38999999999999\"></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=\"68.37\" 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, 195.0)\" x=\"50.0\" y=\"195.0\">500</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, 150.0)\" x=\"50.0\" y=\"150.0\">1000</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, 105.0)\" x=\"50.0\" y=\"105.0\">1500</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\">2000</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\">dataset6/sample=A167</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(\"#t8f815267e1a6455b98768789bde1239b 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(\"#t8f815267e1a6455b98768789bde1239b 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, 1907.0, 771.0, 435.0, 259.0, 182.0, 161.0, 133.0, 129.0, 126.0, 93.0, 117.0, 107.0, 123.0, 104.0, 125.0, 130.0, 89.0, 109.0, 122.0, 122.0, 121.0, 124.0, 132.0, 120.0, 124.0, 96.0, 112.0, 229.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"t675875093dee4a2d9d32170c1807f2a4\", \"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, 161.0, 133.0, 129.0, 126.0, 93.0, 117.0, 107.0, 123.0, 104.0, 125.0, 130.0, 89.0, 109.0, 122.0, 122.0, 121.0, 124.0, 132.0, 120.0, 124.0, 96.0, 112.0, 229.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"tff5b2babd5bc42c5b2e553142071ee10\", \"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(\"#t8f815267e1a6455b98768789bde1239b .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 = {\"tf9d5e7238f5149509ba5315128fd73ae\": {\"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\": 2000.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 0x7ffa547d04d0>"
      ]
     },
     "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_6/clust.85/A167.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=\"dataset6/sample=A167\")\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_6_depthplot.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print final stats table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "\r\n",
      "13093       ## loci with > minsp containing data\r\n",
      "13093       ## loci with > minsp containing data & paralogs removed\r\n",
      "12864       ## 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",
      "90.o31 \t6706\r\n",
      "90.o38 \t6694\r\n",
      "90.o50 \t6490\r\n",
      "90.o53 \t7132\r\n",
      "A056   \t5731\r\n",
      "A081   \t7098\r\n",
      "A082   \t7229\r\n",
      "A089   \t6918\r\n",
      "A107   \t6255\r\n",
      "A145   \t6608\r\n",
      "A147   \t6797\r\n",
      "A167   \t6015\r\n",
      "A177   \t7009\r\n",
      "A1794  \t4790\r\n",
      "A204   \t6040\r\n",
      "A248   \t6077\r\n",
      "A328   \t5912\r\n",
      "JMD1249\t7054\r\n",
      "JMD1285\t6707\r\n",
      "MDS004 \t7271\r\n",
      "MDS054 \t5489\r\n",
      "MDS065 \t6500\r\n",
      "NKK511 \t7033\r\n",
      "SAR6894\t4845\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\t1345\t*\t12864\r\n",
      "5\t844\t*\t11519\r\n",
      "6\t886\t*\t10675\r\n",
      "7\t658\t*\t9789\r\n",
      "8\t695\t*\t9131\r\n",
      "9\t878\t*\t8436\r\n",
      "10\t998\t*\t7558\r\n",
      "11\t455\t*\t6560\r\n",
      "12\t567\t*\t6105\r\n",
      "13\t1175\t*\t5538\r\n",
      "14\t847\t*\t4363\r\n",
      "15\t268\t*\t3516\r\n",
      "16\t205\t*\t3248\r\n",
      "17\t209\t*\t3043\r\n",
      "18\t223\t*\t2834\r\n",
      "19\t232\t*\t2611\r\n",
      "20\t281\t*\t2379\r\n",
      "21\t361\t*\t2098\r\n",
      "22\t498\t*\t1737\r\n",
      "23\t674\t*\t1239\r\n",
      "24\t565\t*\t565\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\t623\t0\t2713\t0\r\n",
      "1\t903\t903\t2414\t2414\r\n",
      "2\t1025\t2953\t1969\t6352\r\n",
      "3\t1095\t6238\t1452\t10708\r\n",
      "4\t1069\t10514\t1073\t15000\r\n",
      "5\t1013\t15579\t757\t18785\r\n",
      "6\t1002\t21591\t628\t22553\r\n",
      "7\t900\t27891\t499\t26046\r\n",
      "8\t745\t33851\t424\t29438\r\n",
      "9\t697\t40124\t284\t31994\r\n",
      "10\t610\t46224\t222\t34214\r\n",
      "11\t545\t52219\t160\t35974\r\n",
      "12\t510\t58339\t110\t37294\r\n",
      "13\t438\t64033\t73\t38243\r\n",
      "14\t372\t69241\t36\t38747\r\n",
      "15\t271\t73306\t22\t39077\r\n",
      "16\t240\t77146\t14\t39301\r\n",
      "17\t200\t80546\t9\t39454\r\n",
      "18\t173\t83660\t2\t39490\r\n",
      "19\t113\t85807\t0\t39490\r\n",
      "20\t103\t87867\t1\t39510\r\n",
      "21\t73\t89400\t0\t39510\r\n",
      "22\t55\t90610\t0\t39510\r\n",
      "23\t29\t91277\t1\t39533\r\n",
      "24\t22\t91805\t1\t39557\r\n",
      "25\t15\t92180\t0\t39557\r\n",
      "26\t8\t92388\t0\t39557\r\n",
      "27\t8\t92604\t0\t39557\r\n",
      "28\t4\t92716\t0\t39557\r\n",
      "29\t0\t92716\t0\t39557\r\n",
      "30\t2\t92776\t0\t39557\r\n",
      "31\t1\t92807\t0\t39557\r\n",
      "total var= 92807\r\n",
      "total pis= 39557\r\n",
      "sampled unlinked SNPs= 12241\r\n",
      "sampled unlinked bi-allelic SNPs= 10885\r\n"
     ]
    }
   ],
   "source": [
    "cat empirical_6/stats/empirical_6_m4.stats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "17797       ## loci with > minsp containing data\n",
      "17797       ## loci with > minsp containing data & paralogs removed\n",
      "17477       ## 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",
      "90.o31 \t6934\n",
      "90.o38 \t7055\n",
      "90.o50 \t6786\n",
      "90.o53 \t7440\n",
      "A056   \t6719\n",
      "A081   \t7315\n",
      "A082   \t7503\n",
      "A089   \t7430\n",
      "A107   \t6794\n",
      "A145   \t7235\n",
      "A147   \t7216\n",
      "A167   \t7098\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head -n 20 empirical_6/stats/empirical_6_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_6/ \\\n",
    "                      -n empirical_6_m4 -s empirical_6/outfiles/empirical_6_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_6/ \\\n",
    "                      -n empirical_6_m2 -s empirical_6/outfiles/empirical_6_m2.phy\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "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 83996 distinct alignment patterns\n",
      "\n",
      "Proportion of gaps and completely undetermined characters in this alignment: 50.46%\n",
      "\n",
      "RAxML rapid bootstrapping and subsequent ML search\n"
     ]
    }
   ],
   "source": [
    "%%bash \n",
    "head -n 20 empirical_6/RAxML_info.empirical_6_m4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash \n",
    "head -n 20 empirical_6/RAxML_info.empirical_6_m2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the tree in R using `ape`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAMgCAIAAABUEpE/AAAgAElEQVR4nOzdZ1wU1/4G8LOFXkQR\nkLoUEYkYFFBCSRRRREHFRJAoIhLEAkJsNxoDxkKiIShee0GjiQ0N5oIoAoqCBGmKShSVILgKiEiv\ny+7O/8Xc/15ijw4ugef7yYvlN2cPZ/gYeHbmzDksiqIIAAAAADCHLe0BAAAAAPQ0CFgAAAAADEPA\nAgAAAGAYAhYAAAAAwxCwAAAAABiGgAUAAADAMAQsAAAAAIYhYAEAAAAwDAELAAAAgGEIWAAAAAAM\nQ8ACAAAAYBgCFgAAAADDELAAAAAAGIaABQAAAMAwBCwAAAAAhiFgAQAAADAMAQsAAACAYQhYAAAA\nAAxDwAIAAABgGAIWAAAAAMMQsAAAAAAYhoAFAAAAwDAELAAAAACGIWABAAAAMAwBCwAAAIBhCFgA\nAAAADEPAAgAAAGAYAhYAAAAAwxCwAAAAABiGgAUAAADAMAQsAAAAAIYhYAEAAAAwDAELAAAAgGEI\nWAAAAAAMQ8ACAAAAYBgCFgAAAADDELAAAAAAGIaABQAAAMAwBCwAAAAAhiFgAQAAADAMAQsAAACA\nYQhYAAAAAAxDwAIAAABgGAIWAAAAAMMQsAAAAAAYhoAFAAAAwDAELAAAAACGIWABAAAAMAwBCwAA\nAIBhCFgAAAAADEPAAgAAAGAYAhYAAAAAwxCwAAAAABiGgAUAAADAMAQsAAAAAIYhYAEAAAAwDAEL\nAAAAgGEIWAAAAAAMQ8ACAAAAYBgCFgAAAADDELAAAAAAGIaABQAAAMAwBCwAAAAAhiFgAQAAADAM\nAQsAAACAYQhYAAAAAAxDwAIAAABgGAIWAAAAAMMQsAAAAAAYhoAFAAAAwDAELAAAAACGIWABAAAA\nMAwBCwAAAIBhCFgAAAAADEPAAgAAAGAYAhYAAAAAwxCwAAAAABiGgAUAAADAMAQsgO6Coih7e3t5\neXmxWCztsQAAwDtBwALoLk6dOkUIoSiqsrJSUszJyRk+fLiGhsahQ4foSn19fWhoKI/HU1NTW716\nNV28e/eum5ubpqamiYlJfHz8+x88AAB0xqIoStpjAADS0dFhaWn5888/e3t7HzhwwNHRkRDS2Ng4\nZMiQ+Pj4qqoqHx+fqqoqQoi7u7udnd28efOqq6s//PDDtra2uro6GxubHTt2jBs3LjY2NigoqLq6\nms3GxycAAKnBr2CAbiEmJsbe3t7a2prH45WUlNDF2NjYSZMmDRs2zNHRsaamhr51mJmZOWTIEA6H\nc+LECRsbGzabHRcXN2rUKFdXVw6H4+Tk1N7eLtVTAQAAwpX2AJg3bdo0fHaHf5CwsDAjI6Po6OhL\nly4RQgwNDTsHrODgYEJIZWWlgYEB/Q/b399/6tSpdIPz588TQgQCQUFBwcOHDwUCQXBwsLm5Of4X\nAACQrh4YsCorKyWzVQC6OTabraen99133925c2fAgAF00dfXl35x7969oUOHEkJu3LhhampKCNmx\nY0dBQUFRUZG6uvovv/wSFBR0+/btWbNmJSQkmJmZWVpaGhoaampqSut0AACA1gMDloyMjLGxsbRH\nAfCmqqqqYmJiGhoaVFRUCCEHDx7ct28fIUQsFvP5fC0tLUJIUlKSl5cXIWTbtm0JCQkmJiaEEA8P\nj+joaEJIfX19YmIim80WCoUWFhZ0EQAApAj3EQCkbO3atXPmzKHTFSFER0eHvkXIYrFUVVUzMzOv\nXr16+vRpOmA9ffo0IyNDLBaXlpYGBgYuX76cEGJhYZGQkNDe3r5ixQojIyMXFxcpng4AAJAeeQUL\n4B+kuLj4yJEjxcXFkoqurm55eXlra6uCgkJUVJS3t7eenl5cXBydwPbu3bt06dKlS5caGRktW7Zs\n+vTphJDo6Oi5c+fKyMhMmTLl+PHjmIAFACB1PXCZBicnp7S0NGmPAgAAAHovfNIFAAAAYBgCFgAA\nAADDELAAAAAAGIaABQAAAMAwBCwAAAAAhiFgAQAAADAMAQsAAACAYQhYAAAAAAxDwAIAAABgGAIW\nAAAAAMMQsAAAAAAYhoAFAAAAwDAELAAAAACGIWABAAAAMAwBCwAAAIBhCFgAAAAADEPAAgAAAGAY\nAhZAd0FRlL29vby8vFgslvZYAADgnSBgAXQXp06dIoRQFFVZWdm5LhAI7O3ti4qKnml/8eJFHR0d\nkUhECNm6dSvrr55vDwAA7w0CFkC30NHR8c0332zdutXAwKCkpIQuCoXClStX2traFhUVDRo0qHP7\n9PT06dOnW1pacjgcQkhgYGDj/9u4ceP8+fPNzMykcBoAAEAIIYQr7QEAACGExMTE2NvbW1tb83i8\nkpISR0dHQkhHR4eDg4NQKCwsLGSz//dxKDMzMzIy0sPDQ19fn67IycnJyckRQpKTk9PS0uLj41ks\nllROBAAASI8MWMXFxTY2NtIeBcAbYbPZu3fvNjU1jY6OvnTpEiHE0NBQcgVLQUHB3d09PT3dzs5O\n8pYrV65ERkYeO3ZszJgx3t7enXu7fft2eHj42bNnZWRk3udZAADAM3pgwBo4cGBaWpq0RwHwN6xd\nu/bOnTsDBgygv/T19e18NCsrKywsjH6dk5MTGRl55MgRFot18+bNkSNHSprV1tZ6eXn9+uuvffv2\nfW8jBwCAF8IcLAApq6qqiomJaWhooCiKoqiffvpJcgWLECIQCPLz8yVBatWqVXFxcUpKSoqKiq2t\nrcrKynV1dfShNWvWzJo165mpWgAAIBU98AoWwD/L2rVr58yZo6KiQn+po6PTOWAVFBQYGRmpqanR\nX6akpNAvtmzZcuvWrd27d9NfFhUVnTt37saNG+9x4AAA8FIIWADSVFxcfOTIkeLiYklFV1e3vLy8\ntbVVQUGBEPL77793noAlkZOTM2bMGMmXy5YtCw0NxdQrAIBugkVRlLTHwDAnJyfMwQIAAAApwhws\nAAAAAIYhYAEAAAAwDAELAAAAgGEIWAAAAAAMQ8ACAAAAYBgCFgAAAADDELAAAAAAGIaABQAAAMAw\nBCwAAAAAhiFgAQAAADAMAQsAAACAYQhYAAAAAAxDwAIAAABgGAIWAAAAAMMQsAAAAAAYhoAFAAAA\nwDAELAAAAACGcaU9ACbl5eXV1dU1NTVJeyAAAADQq/WoK1huUz/bEftbZV2DtAcC8PfcvXvXzc1N\nU1PTxMQkPj5eUs/JyRk+fLiGhsahQ4fepJ+GhoaAgAAtLS0dHZ21a9d2PiQQCOzt7YuKihgeOgAA\nvAiLoihpj4ExFrYOaw+dDPedVpidKe2xALypmpoaGxubHTt2jBs3LjY2NigoqLq6ms1mNzY2Dhky\nJD4+vqqqysfHp6qq6rVdTZ48WV1d/YcffmhsbBw8eHBbWxubzRYKhWFhYUlJSWVlZXTP7+GkAAB6\nOfyqBZCyuLi4UaNGubq6cjgcJyen9vZ2uh4bGztp0qRhw4Y5OjrW1NSIxeLO78rJyXF2dtbQ0DAz\nM0tOTiaECIVCf3//7du3i0SixMTEESNG0Fmqo6PDwcFh7Nixtra2SFcAAO9Hj5qD9e4sLCzk5eWl\nPQroRdhstrOzc0FBwcOHDwUCQXBwsLm5OR2DYmNjg4ODCSGVlZUGBgads1F+fr6np+ehQ4ccHBwO\nHDjg6+tbUVHB5XI9PDyuXr1qbW1NCElKSqIbKygouLu7p6en29nZSeMUAQB6IwSsv9DQ0EhLS5P2\nKKB3aWxsLCgoMDMzs7S0NDQ01NTUpOv37t0bOnQoIeTGjRumpqad3xIeHr5x48ZRo0YRQjw9PRcv\nXkxRFIvFIoQMHz68trb2wIEDAQEBfD5f8pasrKywsLD3d1YAAL0b7hcASFl9fX1iYmJzc3N6evrV\nq1e9vb0JIWKxmM/na2lpEUKSkpK8vLw6vyU7O3v06NH064yMDGtrazabPXPmzNbWVhaLpaamNmXK\nFBkZGUl7gUCQn58/cuTI93dWAAC9GwIWgJRZWFgkJCS0t7evWLHCyMjIxcWFEMJisVRVVTMzM69e\nvXr69Gk6YD158iQnJ4cQIi8vHxcXJxKJsrOzQ0NDIyIiCCHJyclHjx4ViUSPHj0KCgoKCQmRfIuC\nggIjIyM1NTUpnSIAQK+DgAUgZdHR0XPnzjU2Nm5paTl+/Dg914rFYkVFRXl7e/v7+8fFxamoqBBC\n9u/fv2TJEkJITEzM5s2bNTU1ly9ffvDgQUdHR0LIzp07161bp6Gh4ebm5unpGRoaKvkWv//+OyZg\nAQC8T1im4S+cnJwwBwsAAADeEa5gAQAAADAMAQsAAACAYQhYAAAAAAzrUQGrovTPqC/n1TyukPZA\nAAAAoFfrUQuN3rx2ta2tzdfXV9oDAQAAgF6tRwUsHR0dQkjn9RUBAAAA3r8edYsQAAAAoDtAwAIA\nAABgGAIWAAAAAMMQsAAAAAAYhoAFAAAAwDAELAAAAACGIWABAAAAMAwBCwAAAIBhCFgAAAAADEPA\nAgAAAGAYAhYAAAAAwxCwAKTv7t27bm5umpqaJiYm8fHxknpOTs7w4cM1NDQOHTr0Jv00NDQEBARo\naWnp6OisXbuWLtbX14eGhvJ4PDU1tdWrV3fJCQAAwF8hYAFIWU1Njaur66JFiyoqKtavX+/n5ycW\niwkhjY2N06ZNO3DgwOHDh5ctW/YmXfn4+IhEosLCwsuXL69fv57uZ+bMmZqamvn5+VeuXPn+++/p\nIgAAdCkELAApi4uLGzVqlKurK4fDcXJyam9vp+uxsbGTJk0aNmyYo6NjTU3NM8EoJyfH2dlZQ0PD\nzMwsOTmZECIUCv39/bdv3y4SiRITE0eMGMFmswkhmZmZQ4YM4XA4J06csLGxoYsAANCluNIeAEBv\nJxAICgoKHj58KBAIgoODzc3N6QwUGxsbHBxMCKmsrDQwMOgcjPLz8z09PQ8dOuTg4HDgwAFfX9+K\nigoul+vh4XH16lVra2tCSFJSEt3Y399/6tSp9Ovz58+/79MDAOiVELAAutamTZv279+vra39sgZf\nf/31gAEDzMzMLC0tDQ0NNTU16fq9e/eGDh1KCLlx44apqWnnt4SHh2/cuHHUqFGEEE9Pz8WLF1MU\nxWKxCCHDhw+vra09cOBAQEAAn8/fsWNHQUFBUVGRurr6L7/8EhQUdPv27S48WwAAIIQgYAF0tQED\nBoSEhAQGBr6swcOHDxMTE9lstlAotLCwiI6OJoSIxWI+n6+lpUUISUpK8vLy6vyW7OzsmJgY+nVG\nRoa1tTWbzZ45c+a+ffsUFBTU1NSmTJmydetWQsi2bdsSEhJMTEwIIR4eHnTnAADQ1TAbA0DKLCws\nEhIS2tvbV6xYYWRk5OLiQghhsViqqqqZmZlXr149ffo0HbCePHmSk5NDCJGXl4+LixOJRNnZ2aGh\noREREYSQ5OTko0ePikSiR48eBQUFhYSEEEKePn2akZEhFotLS0sDAwOXL18u1XMFAOgtELAApCw6\nOnru3LnGxsYtLS3Hjx+n51qxWKyoqChvb29/f/+4uDgVFRVCyP79+5csWUIIiYmJ2bx5s6am5vLl\nyw8ePOjo6EgI2blz57p16zQ0NNzc3Dw9PUNDQwkhe/fujYiI0NDQmDZtmr+//8KFC6V6rgAAvQWL\noihpj4FhTk5OaWlp7/+9AC905MiRpqamV9wiBACAngdXsAAAAAAYhknuAAyorq4uKCh44aHCwsLW\n1tbU1NRn6h999JGysnLXDw0AAKQAAQuAAVu3bs2+W6qhq/+ig2zCUTqYdLFzqbSocLZHaUBAwHsZ\nHQAAvG8IWAAMoCjqY/ep5ja2b9g+4/QpbFkDANCDYQ4WAAAAAMNwBQuAYdczL/20YU3VI76OobH/\nqnXm1iMJIV9/PvnOtTy6wSeTPrUa5SzVMQIAQNdCwAJgEiUWb1qyYN63G0aOnZB64nD0sqDdabkU\nRZXf/3PvpauKKiqEEA6HeyXljLRHCgAAXQi3CAEYxWIpKqs0Nza2tba0tbQoqaoSQuprqjsEgg0L\n/b5wHPZjaGB9TbW0RwkAAF0LV7Cgt8jJyVm3bp2Ojk5XdJ6Xl/epmRUhhMViBX23afVsz13hyzlc\nmeiEC4SQmseVvEGDvwiL0NTV/2nDmq0rvhzrOaMrhgEAAN0EAhb0FnV1dQMHDly0aFFXdL5582b6\nRVND/eYlC5ds2vXBCNvjW6P2fxf+zd7Dxh8M/e5YAt3AZ+nKuZ9YjfbwxP99AAA9GH7FQ2/B5XL7\n9u1rbGzcFZ337duXfnG3IF/LgOcwcTIhZPzns8N8phJC7t24JmhvGzLCjhAiIyvH5nDpDQcBAKCn\nwm95ACYZmJo9Kim+/nt6e2vLpfiTPDNzQkhHe/uPoYH84jsiYcfJndE2TmPZHI60RwoAAF0IV7AA\nmNRfWzf4++gD361+Uv7QyHxIUMQmQsgHIz6aOjfo29leIpFwqN3H89f8cP33dGmPFAAAuhACFgDD\nRoxxGTHG5Zni5DnzJ8+ZL5XxAADA+4eABcAANpt9bOuPqn37vWH7JxWPLBYt7NIhAQCAFCFgATDg\n66+/Dg0NfeGhkydPNjc3z549+5m6iopK148LAACkAwELgAGysrKysrIvPKSkpERRlOQxQwAA6A3w\nFCEAAAAAwxCwAKTv7t27bm5umpqaJiYm8fHxknpOTs7w4cM1NDQOHTr0Jv1QFGVhYcH6f3PmzHm7\nfgAA4B3hFiGAlNXU1Li6uu7YsWPcuHGxsbF+fn7V1dVsNruxsXHatGnx8fFVVVU+Pj6+vr6v7Soz\nM1NdXb2xsZH+UkZGhhDyFv0AAMA7whUsACmLi4sbNWqUq6srh8NxcnJqb2+n67GxsZMmTRo2bJij\no2NNTY1YLO78rpycHGdnZw0NDTMzs+TkZLq4e/fu5uZmXV3doUOHZmRkyMnJvbYfAADoCriCBdC1\nRCJRdXV1SUnJyxq0tLQUFBQ8fPhQIBAEBwebm5vTG+nExsYGBwcTQiorKw0MDDrvrpOfn+/p6Xno\n0CEHB4cDBw74+vpWVFS0trYWFxevWrVq1KhRW7ZsCQgIePTo0av7AQCALoKABdC15OXlc3Nzy8rK\nXtZg1qxZZ8+eNTMzs7S0NDQ01NTUpOv37t0bOnQoIeTGjRumpqad3xIeHr5x48ZRo0YRQjw9PRcv\nXkxRlKKiYlZWFt3gyy+/jIyMFIvFbDb7Ff0AAEAXQcAC6Fqenp6enp6vaPDw4cPExEQ2my0UCi0s\nLKKjowkhYrGYz+draWkRQpKSkry8vDq/JTs7OyYmhn6dkZFhbW3NZrNPnTplaWlJ72ZdXFxsZWXF\nZrNf3Q8AAHQR3CwAkDILC4uEhIT29vYVK1YYGRm5uLgQQlgslqqqamZm5tWrV0+fPk0HoydPnuTk\n5BBC5OXl4+LiRCJRdnZ2aGhoREQEIeTYsWPh4eHNzc0PHjwICQlZt27dy/oBAICuhitYAFIWHR09\nd+5cGRmZKVOmHD9+nJ4jxWKxoqKivL299fT04uLi6GXf9+/fn5CQcPny5ZiYmODg4LCwsCFDhhw8\neNDR0ZEQsm7dOl9fX11d3YEDB0ZERDg5Ob2sHwAA6GosiqKkPQaGOTk5paWlvf/3Qjd34cKFy5cv\nh4eHS3sgAADQ8+EWIQAAAADDELAAAAAAGIaABQAAAMAwBCwAAAAAhiFgAQAAADAMAQsAAACAYQhY\nAAAAAAxDwAIAAABgGAIWAAAAAMMQsAAAAAAYhoAFAAAAwDAELAAAAACGIWABAAAAMAwBCwAAAIBh\nCFgAAAAADONKewAADBg0zFpBSelNWv56LvVlh2oeV/KL7zI3KAAA6L0QsKAnkJWTX73/+Dt2Eu47\njZHBAAAA4BYhAAAAAMMQsAC6kerqag6H09jYKBAI5OTk3N3dOx/18/NjsViPHz9uaWnhcDhcLpfD\n4SgqKrq5uZWXl0trzAAA8DwELOiBOgSCBc62ddVVkkrd0yffzvGaPdJ8w0K/lsaGVxSlKy8vb/Dg\nwSoqKoWFhWw2u6ysTHLo1q1bJ0+e1NfX19LSKigo0NPTEwqFIpGorKxMRUUlKCiocz8CgcDe3r6o\nqEhSycrKsre379evn7GxcXx8fOfGFy9e1NHREYlEhJC7d++6ublpamqamJg80wwAAN4cAhb0NImH\n9oX5TK16xO9cPLhhDW+Q+b7L1+UUFE/u3PKKonTl5uaOGDGCEJKXl+fi4vLgwQOKouhDYWFhY8eO\ntbGxoY/SLwghGhoamzZtSk5OplsKhcKVK1fa2toWFRUNGjSIblNRUeHu7h4eHl5dXR0WFjZv3jzJ\nd0xPT58+fbqlpSWHw6mpqXF1dV20aFFFRcX69ev9/PzEYvH7PH0AgB4Dk9z/hvj4eEVFRWmPAl6D\nZ2auZcD7fv5sSYWiqJzz577Ze1hGVtZ23ISjW37w/VfY88W+GlpdMR4+n5+Xl6eiovKyBjY2Nmpq\navTr3Nzc8ePHE0Ly8/MdHR2zs7Pr6ur69u2bl5fH5/NHjhypo6ND/hqwCCFKSkosFot+3dHR4eDg\nIBQK6WtgdDEjI8PPz8/V1ZUQMmzYMA6HQ9czMzMjIyM9PDz09fUJIXFxcaNGjaKbOTk5tbe3M/3D\nAADoLRCw/oZ//etfc+bMkfYo4DUsbB2eqbS1NLe1NOsaDySE6JmY1lY9fmGxiwLWr7/+mpGRMXLk\nyJc1MDAwoAMWRVG5ubmrVq0ihOTl5Xl7e/N4vLKysr59+3799dcRERHh4eGTJ0+mj/r4+Eh6KC0t\nNTU1pTOWgoKCu7t7enq6nZ2dpIGXl5eXlxchpLi4OCgoKDIykhBy5cqVyMjIY8eOjRkzxtvbmxAi\nEAgKCgoePnwoEAiCg4PNzc0lEQ0AAP4WBKy/QVtb+6uvvpL2KOAFfo57/WwhOn9QFCUWiV5dZJam\npub48eMDAwNf1qC9vV0sFrPZ7EePHlVXV1taWra3t//xxx9WVlZ0wKqrqxMKhaNHjy4oKLC2tm5s\nbLxz5461tbWkh+Tk5I8++qhzn1lZWWFhYZ0rNTU14eHhV69e3bVr14cffpiTkxMZGXnkyBEWi3Xz\n5k06/82aNSshIcHMzMzS0tLQ0FBTU5PpHwYAQG+Bj6fQ88krKskrKpXf/5MQUlFaoj5A+2XF94+i\nKG1t7cePHxNCEhMTrays5OXlb968yePx+vTpQwcs+vLVrVu3tLW11dXVr127xuPx1NXV6R7q6up2\n7ty5fPlySZ8CgSA/P7/zNbPExEQPDw9nZ+fMzMwPP/yQELJq1aq4uDglJSVFRcXW1lZlZeW6urr6\n+vrExMTm5ub09PSrV6/Sl7UAAOAt9JYrWN9///32ffvlFV4/g2qgheXLDjXW1f72228eHh6MDg26\nHIvFGjHGJf/S+UHDrK9lpI1wdn1hsfhmwfsfG0VRzc3NxcXFDQ0NUVFRy5YtI4Tk5+fTU6x4PN6u\nXbtMTEzs7OxiYmKemeHe0dGRk5MTEhKycOFCY2NjSZ8FBQVGRkaSeV1paWnffvvt2bNn5eXlm5ub\nWSyWkpJSSkoKfXTLli23bt3avXs3IcTQ0PDgwYOurq6rVq0yMjJycXF5vz8MAICeo7cErPr6+kXf\nR5taWr1LJ2mnjtfV1TE1JHif/Fas3rw0KHC09UALS7+vVr+w+P1CKUywY7PZERER7u7uSkpKCxYs\nmDt3LiEkLy+PvgNoaGh469atI0eOkE65Ki8v78SJE1wul8vljhw5cvny5c9cavr99987T8D64Ycf\n8vLyNDQ06C8dHR0zMjIkR3NycsaMGUO/jo6Onjt3royMzJQpU44fP44JWNAz3L17d/Hixbm5uSoq\nKps3b6YnMhJCcnJyvvjii8LCwt27d/v5+amoqIwbN+706dP00erqak1NTYqi+Hy+iYlJR0cHm80W\n/XUuQWVlpYqKioqKCkVRFEXJysqOHTt279699MMoFEU5ODhcvXq1paUF/zf1QizJQ+A9hpOTU1pa\n2jPFFStW9LX6+N0DlvUANT8/v3fpBLqC8ZChnwYuesdOYrdtenCv6PXt/qYjR440NTW9Yg4WAHSd\nmpoaGxubHTt2jBs3LjY2NigoqLq6ms1mNzY2DhkyZMWKFRs3bmxtbU1KSnJwcBg4cODNmzfpN+7Z\ns2fBggW6urq//fabnZ0dRVEXLlxYtWrVwYMH586dm5qaqqCg0NLSkpqa6urqmpWVNW7cuPr6+ilT\npojF4sTEREKIZJUTQ0PDs2fPDh48mBDS2to6ceLEixcvslgsOTm5L7/8Mj4+/tatW53HfPv27ZSU\nlJCQkM7FqVOnxsXFvY8fGTCkt1zBgp7tu7BVDQ2vWSm0qKiouLj4mbXRO1v3zcq/+32///77+vr6\nV7e5deuWQCAoKSl5WQMDA4OFCxf+3W8NAJ297DJVVFRUfX39rFmzoqKiXFxcJIuPxMbGTpo06enT\npx9//PGxY8dyc3NdXFwuXrxYX1+/dOnSuLg4+n5FbW3toUOHTExM7ty54+rq2tbWRuckQoiMjAwh\nZN++fSoqKr6+vvSvoG3btpmammpra7e2tiopKdnb28fGxpaWltKL0gmFQgsLi5KSEg6HU1paGhIS\nsmHDhvj4+G3btgkEgoyMDFlZ2ba2NlNTUyMjo8GDBy9fvlwkEpWUlIwYMYK+ki1x8eLFGTNm8Pl8\nyaor0N300oDVIRCETPj4++MJav3/+5xU3dMn0cuC7v9x09zGNmTjvxVVVJ8vSnXI8CpvMh37woUL\nly9fZvBKkkgk2vvLEd9lYa9uZmX2muum0dEbELAA3gW9Rq7kMpWfn5/kMtXOnTs1NTXDw8MXLVpk\na2srWXwkNjY2ODh479691tbWBgYGV69epVee8/b2lpGRkZGRUVBQaG1tdXZ23rNnj6ysLJfLNTU1\n/eKLLxYvXqykpNTe3t7e3k5RVHp6upWVFZ/P79evX2trq4qKilgs1tXVVVVV1dPT++yzz9LS0qqq\nqugHli9dulRaWqqurt6vXz89Pb1BgwZxuahPE/kAACAASURBVNxJkybNmjWLw+Gw2ex+/fpVVVVx\nOBwOhxMTE3P79m02m01RVEFBgby8vOR86fWBrayskK66s94YsBIP7cs4feqFK32v2v3LtpVfnty5\nxfdfYc8X9U0HSWnI0E0p91H70P7jd+zk5K5usYg8wD/Xy9bIjY2NnTZtGp/Pnzt3bmtra58+fSSL\nj9y7d08kEp07d+7ChQv07T9vb28DAwMXF5dz586pq6s3NTUpKCj4+vomJCS0t7dra2tv2bLl0qVL\nSkpK5ubmV65cUVRUrKurq66uHjVq1M6dO/v160cIKS0tVVZWnjhx4tGjR3/55ZfGxsampiY60l25\ncmXr1q1GRkaPHj2ys7O7d+9ebGysrq7uw4cPP/zwQ3pm5KNHj9hstlgsbmtry83NtbW1bWlp+eST\nTzZv3kwv5kKeWx8Yuq0eGLAaGhpOnDjxTLGoqMjO6r9/CN96pW8ELACAbuhla+TGxsZ6e3vv2rWr\ntLTU2dn52rVr0dHR7e3tHA7nwYMHQUFBIpHIx8dHIBAcPnx4+PDhBgYGRkZGly9fbm5uJoTIyso2\nNTURQiiKevTo0ahRowghw4YN+/PPPwkhRkZGhYWFHR0d58+f79+//6NHj/z8/OhNq0pKSu7cuTNg\nwAB6eCoqKrm5ufTKc/fv3x86dGhSUlJKSsoPP/xw9uxZExMTsVisoKAwZsyY4cOHb9y4kc1mKyoq\n8ni8kpKS5OTk5ubm6Ohouqvn1weGbqsHPtcwYcKEkud0fvrPwtbBZvS4zm95w5W+3+95wD/PP3eT\naYB/tFmzZg0YMMDMzMzHx0dNTc3CwoKu37t3LzQ0NCEhIT8/nxBiZGQ0btw4bW3tqqoqNpttaWk5\nePDgc+fOff755ywWq62tLSEhgc5MTk5OAwYMEIlEK1euVFJSkpeXNzIyGjZsGI/He/z48ZMnT/r3\n729iYpKens5ms+vq6u7du1dXVzdnzpzt27fX19efOXPG3t4+Li6ura2Nw+EoKChIVp6zsLCgKEoo\nFC5atGjDhg2nT582MTGRl5dXUFAIDg7W19fncrmEkFOnTuXm5n7++eeDBg06f/58nz592Gy2ZH1g\nNpstWR8Yuq0eeAVr/fr1zxdra2tf+0ZprfQNPcNb33qW0ngBeg56jVw2m03PIqev94jFYj6fv337\n9rlz57a0tIwYMeL48eMsFqu5ufnPP/+UkZFJTk7W1taOi4s7ePBg//79f/rpJ4FAEBMT88EHH1y+\nfJn+5T9jxoxt27Y5OTnZ2Nj8+OOPAoGAvtrk5OTE4/Hy8vIUFRVLSkooilJXV/fz85s4ceLBgwdl\nZWX/85//yMvLp6SksFgsFouVkpIiFotlZGQ0NDTa29v5fH51dXVUVFRCQsLjx4+bmpqcnZ0HDRo0\nd+7cTz/9lBCyZ8+epqam0NDQBw8e/PjjjwEBAYSQVatWpaamSp4lVFZWrq2tlax4B91ND7yC9Ra6\n80rf8E/BMzOftvDLzhX6LvNHLm70XebctOSXFQHgXVhYWNAzpVasWCFZI5fFYqmqqhobGyclJamp\nqcXHx6uqqj59+nT+/Pnu7u7t7e0TJky4f/8+RVE///yzl5fXihUr5OXlS0tLHz58GBoaqqWlNXHi\nRHqjqrNnz0ZERNAPBtrZ2enp6dnY2PB4vLS0tKamppaWFqFQSAhZvny5mppaW1tbeXm5hoaGiorK\npEmT5OTkHj9+3Nraymaz6dd1dXUqKipGRkbq6uqBgYGNjY0sFis7O9vd3b2lpWXHjh2EkLa2Njab\nPXDgQFdXV6FQGB4eTghJSUmhF9yKjo4ODAykKArpqjtDwCKk06LeFEU9v9J35yL0Tr/++uvAgQNt\n/uqZ6/O49QwgLfQaucbGxi0tLZI1clksVlRUlLe3t7+/f1xcnIqKCiFk//79+fn59BWvoqIiLS2t\n5cuXJyYmbt26lRDy888/GxoaikSilJSUb7/99j//+c/evXvp//3V1NSsra137979wQcfzJ8//1//\n+pehoWFDQ8P169dFIlFgYCC95G9paalYLOZwOHJych9//PHRo0ebmpooilJQUCCEpKSkDBs2TFFR\n0czMbNOmTZWVlVVVVQKBYP/+/XJycrKysmfPnlVVVSWEpKWlxcTEyMvLy8rKnjlzhh68RE5ODu4P\ndn89cKHRF3p+odHPBuvEXC743zIN1VWblwaVl/450MIyZOO/FZRVni9eSTmDhUb/uehlGugPgn9X\nfHx8cXHxkiVLOhdFItHwj0d/eyC2c7Hzv6vW5iYf60E/XflDRa1v2d3bK73cjxT8+Xxx4IfDC7Mz\n3+XUAACgu+mBc7BeJvt8Uumd/62WO3/tD7kX/nJ3xtFtCv3i8pn/vLB451q+9WRcx4I3JbnLbDbc\n5vlbz52LAADQw/SWgDVz5sysrKzXNouKilq6dOnLjqbdKhi3NYrRcUFP1p03mQaAt0NR1KRJkyRb\nFnZWVFRkamqKxT+B1lsC1tChQ4cOHfraZkePHn3FSt9Hjx7V1dVldFzQw3XbTaYB4O2wWCxDQ8Pf\nf//d3t7+mUMTJkz46aefRo0aVVFRsWHDhi1b8Ixwr9ZbAhbA+/FrUXnnL9X6a645+Oyyty8sAkA3\nt3nz5gkTJgwePDgoKOi77757PmC1tLQoKSnNmzfvjz/+2Llzp1QGCd0HAhbAW2KxWI21Nf/67F2n\n5fVR78/IeACgS02bNm369Olnz541Nzd//Pjx48ePtbS0OjcYO3Ysh8NJSkpavXr1m9wzgZ4NAQvg\nLbHZ7Pu3/3htsyNHjjQ1NTG4yTQASIW+vn5YWJifn19cXJy/v//+/ftXrlz5TJtvvvmGx+P5+/tL\nZYTQrWAdLAAAgDcyYcIECwuLLVu2fPrpp6dOnRL9dYePjo6O69evnzp1SlrDg24FAQsAAKTs7t27\nbm5umpqaJiYm8fHxknpOTs7w4cM1NDQOHTr0Jv00NDQEBARoaWnp6OisXbu28yGBQGBvb19UVPSO\nQ12zZk1qamp2draLi8uZM2c6H8rOzg4LC1NXV3/HbwE9AwIWQDdSXV3N4XAaGxsFAoGcnJy7u3vn\no35+fiwW6/Hjxy0tLRwOh8vlcjgcRUVFNze38vLyl/UJ0M3V1NS4urouWrSooqJi/fr1fn5+YrGY\nENLY2Dht2rQDBw4cPnx42bJlb9KVj4+PSCQqLCy8fPny+vXr6X6EQuHKlSttbW2LiooGDRr0jqNl\ns9kHDhxYsmTJ1KlT9+zZ0/mQjo4OJgOABAIWQDeSl5c3ePBgFRWVwsJCNptdVlYmOXTr1q2TJ0/q\n6+traWkVFBTo6ekJhUKRSFRWVqaiohIUFNS5n+c/rGdlZdnb2/fr18/Y2Ji+QlBfXx8aGsrj8dTU\n1Fav/u8SEhRFWVhYsP7fnDlYQgK6XFxc3KhRo1xdXTkcjpOTU3t7O12PjY2dNGnSsGHDHB0da2pq\n6LQkkZOT4+zsrKGhYWZmlpycTAgRCoX+/v7bt28XiUSJiYkjRoyg98zp6OhwcHAYO3asra0tXXlH\nGhoa27ZtW7VqFZvNvnfvHl0sKSkxNjamt4gGIAhYAN1Kbm7uiBEjCCF5eXkuLi4PHjyQbGYVFhY2\nduxYGxsb+ij9ghCioaGxadOm5ORkuuULP6xXVFS4u7uHh4dXV1eHhYXNmzePEDJz5kxNTc38/Pwr\nV658//339F+vzMxMdXX1xv+3a9eu9/4zgF5HIBAUFBQ8fPiwpKTE39/f3NycjkGxsbGurq6EkMrK\nSgMDg87ZKD8/39PTMzw8vKKiYtmyZb6+vhRFcblcDw+PoqIibW3tkJAQyb5YCgoK7u7uLBbLzs6O\nqTHb2tq6uroqKSnt3buXrqSmpo4dO5ap/qEHwFOEAF2rtrb2ypUrr2gwceJEPT09+nVubu748eMJ\nIfn5+Y6OjtnZ2XV1dX379s3Ly+Pz+SNHjtTR0SF/DViEECUlJcnnZvrDulAopK+B0cWMjAw/Pz/6\nb9WwYcPolaYzMzMDAgI4HM6JEydsbGzoxrt3725ubtbV1e3Xr9+OHTsmTJjA+A8E4BmzZs1KSEgw\nMzOztLQ0NDTU1PzvFrH37t2jFzu4ceOGqalp57eEh4dv3Lhx1KhRhBBPT8/FixdTFEX/XzB8+PDa\n2toDBw4EBATw+XzJW7KyssLCwhgc9pdffjl9+vT4+Pi1a9fKy8unpqZu3ryZwf7hnw5XsAC6loOD\ng4ODQ9+Xk2QjiqJyc3M7X6Pi8Xj0XcKvv/46IiIiPz9fctTa2lryLUpLS01NTel+Xvhh3cvLKyoq\nihBSXFwcFBQUGRlJCPH39586dWq/fv3Cw8PXr19PCGlpaSkuLl61atX9+/d9fX0DAgLe208JerP6\n+vrExMTm5ub09PSrV696e3sTQsRiMZ/PpxeaSkpK8vLy6vyW7Ozs0aNH068zMjKsra3ZbPbMmTNb\nW1tZLJaamtqUKVNkZGQk7QUCQX5+/siRI5kdeUxMTFtb27Zt20Qi0cOHD7HVB/wFBZ2MHj36rY9C\nN3f+/Pk1a9a83Xv/85//REVFMTseiba2NpFIRFEUn8/ncrmtra1tbW1ycnJ1dXWenp6//fZbWlqa\nk5OTQCCQl5evrq5uaGhgs9nV1dWSHn744Yf58+d37tPR0fHcuXOdK0+fPg0KCrKzs7t+/TpFUdu3\nbx8zZkxRUdGTJ082b948ePDgZ0ZVU1OjoKBADwygS/Xp0+e3335ra2tbunSpq6sr/a9OLBb369cv\nJSUlPz9fV1e3oaGBoqiqqqrs7GyKonR1dbdv3y4UCq9cuWJkZJSRkUFRVP/+/WNiYoRC4cOHD11d\nXTdv3iz5FtnZ2R988EFXDP78+fP9+vXLzs5etGhRV/QP/1y4ggUgTRRFaWtrP378mBCSmJhoZWUl\nLy9/8+ZNHo/Xp08f+goWffnq1q1b2tra6urq165d4/F4kkfB6+rqdu7cuXz5ckmfz39YT0xM9PDw\ncHZ2zszM/PDDDwkh27Zt27Nnj5mZWf/+/T08PFpbWwkhp06dKikpod9SXFxsZWXFyIxggFeLjo6e\nO3eusbFxS0vL8ePH6X91LBYrKirK29vb398/Li5ORUWFELJ///4lS5YQQmJiYjZv3qypqbl8+fKD\nBw86OjoSQnbu3Llu3ToNDQ03NzdPT8/Q0FDJt/j999+fmYBVU1NTV1f37oMfM2ZMaGhoZGTkuHHj\n3r036FGknfC6F1zB6sG65xUskUgkKyubnp5eVFRkamq6e/duiqJ27do1Y8YMiqK2bt06ZMiQyZMn\nUxS1b98+T09PiqKioqLoFwKB4PLly1ZWVpGRkZ37fObD+oULF2xsbJ48eULPW29qaqIoSlNT88CB\nAyKR6P79++PGjdu2bRtFUV5eXjNnzmxqaiorK/voo48uXLjQFacM0B0UFBR89NFH+fn5795VY2Oj\nsrJyXV3du3cFPQkmuQNIE5vNjoiIcHd3V1JSWrBgwdy5c0mnKVaGhoa3bt06cuQI6TSxPS8v78SJ\nE1wul8vljhw5cvny5fScFYlnPqz/8MMPeXl5Ghoa9JeOjo4ZGRl79+5dunTp0qVLjYyMli1bNn36\ndELIunXrfH19dXV1Bw4cGBER4eTk9L5+DADvm6WlZUJCwuzZsz/55JOvvvrqXbpSVla2sbHp06cP\nU2ODnoFF/f9D4EAIcXJySktLe7uj0M1duHDh8uXLkie3/5b4+Pji4mL63gQASAtFUQ4ODlevXm1p\naWHk/rVIJFq3bt3Nmzf379//LgkJfx3gebiCBb1dcHBwVVXVq9uUl5c3NTW9YrUFbW3tLVu2MD00\nAPgfiqLMzc3v3bvH5XIrKyvpJUtoAoFg9OjR+/fvHzx4cOe3XLx4ccaMGXw+n8PhbN26NSQkpPPR\n27dvDx48+Ntvvz1z5oyrq+uePXvoVSG6MzpiZmdny8jIMJUyoYsgYEFvl3L5ypc/bnvHTrYsD2Zk\nMADwMidPnnzw4AGHw9HR0SkpKaEDVlZW1tSpU58+fSonJ0evrFtfXx8eHv7bb789ffqUEPLxxx9z\nOJy7d++eOXOmf//+ysrKGzduLC0tvX//vpmZGd3zxIkThw4dOnv27JkzZ37xxRdSPMfXOnXqVE1N\nDUVRYrFYkjJzcnJmzZpVUlKioKCwbds2X1/f+vr6FStWHD58uLm5WUZG5quvvlqzZg0hJCsra+nS\npUVFRS0tLZs2bVq4cKG0T6gnQ/iF3k5GTk5Ln/eO/3Fl5aR9HgA9WUdHR0hIiIuLC/0ILf24a2Nj\no5eX19dff+3h4dHR0UFfzqG3KNi5c+eIESNaW1vt7e3pvQ5DQ0MrKyu/++67L7744sKFC//+9787\nb2ujr69/9uzZGzduzJ49m36othvq6OhYtWqVQCDQ09PT0tKS/BA+++yzhoaGf//732w2m960cebM\nmWlpaYsWLbp8+bJIJKK3aqB3dNDX1+fxeFwud926ddI+oR4OAQsAALq7HTt2tLa27t69m8fjKSsr\n09kiNjZ28uTJISEh9Nacku2eZGRktm/f7ujoqKio6Ojo2HmvQz09vZaWlsOHD3dehpQmJye3ZcuW\n8ePHjx8/XrJeSbcSExOjrq4+ZsyYQYMGqaqqSn4Ipqamn3766ezZs5uamuhNGzMzMx8+fGhubp6Q\nkKCtrU1v1ZCRkeHj4zNr1qyxY8daWlrSOzpA18EtQoC/6BAIQiZ8/P3xBLX+/92vo+7pk+hlQff/\nuGluYxuy8d+KKqrPF6U6ZIAerqmpaf369fX19QMGDKArPB6PEBIbGxscHEwIuXz5sqamJn0Fa+LE\nifRTgWfOnJGTkxs5cuTt27fpvQ6fPHni5uZmbm4uWUbueTNmzLC2tp41a9bSpUs//fTT93F6b6au\nrm7hwoWqqqq//vrrqlWrampqEhMTf/nlFxaLVVtb+9NPP02ZMkVHR4fL5f75558sFqu5uXnWrFn0\ne3/66SdCyPDhw8PCwgYPHlxXV8fn8+kdHaDr4AoWwP8kHtoX5jO16hG/c/HghjW8Qeb7Ll+XU1A8\nuXPLK4oA0BXWrVtXXV0t+dLBwYG+eENvVigQCG7cuPHBBx8QQnbs2JGUlCRp2d7erqysPHny5AED\nBpiZmbm5uZmYmHTex/OFzMzMzp07d/LkydDQ0I6Ojq45p78tKCiIoig6ZcbExFy7du327ds2NjaF\nhYWFhYWtra2nT59WV1fX1NS0s7NTVFT86KOPYmNjZWRkNDU1g4KCrly54uLisnbt2j/++GP//v18\nPn/GjBksFovFYs2ZM6ehoSEgIEBLS0tbW9vAwKCoqKjztxYIBPb29s8U4bUQsAD+h2dmPm3hl50r\nFEXlnD/3kYubjKys7bgJuWnJLysCQFeoqqrauXPnypUr6cUbf/rpp/r6+pKSEslmhQUFBUpKSjNm\nzCCEbNu2LScnh24ZHh6urKxML0WUmJiYn5+vqqra0tLyzLpxL6SsrHzkyJEhQ4a4ublVVFR0+Um+\nzt27d48dOyYnJ2diYhIUFKSkpGRqasrn883NzcvLy8Vi8eDBg/fu3dvR0VFSUqKqqiovL//LL790\ndHRYWlrKyMg8ffrUx8dn8eLF27dvp7d3tLW1raysrKysrK6u3rVrl4+PT0dHh7e3d79+/fh8/sCB\nA+nvKxQKV65caWtrW1RURD9DAG8Otwiht6Ao6tChQ/Hx8c8e4PxvfrqFrcMzB9tamttamnWNBxJC\n9ExMa6sev7DYV0OrSwcP0GstXbq0o6ODnrtNCNHR0amqqqqqqmpra1NVVc3MzExKSmpvb6d3g376\n9GlGRoaRkdGDBw/27NkzZcoUQoiFhcXBgwd3796tr6/P5XJdXFwknTc2NtKb8LxQYGCgtbX11KlT\nIyIinJ2du/hEX2X+/PmysrKOjo7jx4+Pi4sTCoU1NTUNDQ3l5eVycnIURWlqat6+fbuoqOjmzZuj\nR49ua2tLT0//7bffbty4IRQKCSFz5879+eefz549O3PmTBaLlZeXp6Ojo6GhceDAgXHjxvn7+3/8\n8cdnzpwpLCwsLy/ncv+bDTo6OhwcHIRCYWFhIZaE+Lvw84LewtnZubi4OO85b/Je+mkjiqLEItGr\niwDAoOLi4sTExEePHvXr14+u6OrqVlVVtbS0KCoq0psVpqampqWl0Tlp7969ERERGhoa06ZN27x5\n888//0z+f6/D69evm5qaSvY6JIRcunTJ3d193LhxiYmJ9AT551lbWycmJkZHR2/cuFFa63Jfu3Yt\nMzOTDpFOTk7Z2dn6+vpPnz4Vi8Wpqanz588XiUTp6elsNpvL5SooKOzdu1dJSSkgICA+Pl4sFtPz\n2X/88Ud6R4fk5GSRSGRiYrJ06dKWlpaAgAAul+vh4VFWVubr63vhwoVJkyZJvrWCgoK7uzuLxXpm\nJ0d4EwhYAK8ir6gkr6hUfv9PQkhFaYn6AO2XFQGAcQMHDqypqZGkK0IIvc+mgoICIcTPz6+6urqg\noECytfnkyZPv3bv39OnTvLw8b29v+lOQn59fVVXVo0ePduzYoaqqKumqvr5+9OjRCgoKCxYs6N+/\nf3Bw8G+//VZfX//MGNTV1ekr359++mltbW1Xn/IzCgoK/Pz8xowZQwgJDQ1dtGiRvr6+nJzciRMn\nxo8f/+eff4aEhCgoKEyaNMnZ2XnatGnR0dGTJ08uLy8fP368vLx8SEhIXl6eh4cHRVEikUiyTd6d\nO3dWrlzZ0dFRW1tLh8vhw4fX1tYaGRmdPXv2mTFkZWV99NFH7/nEewAELIBXYbFYI8a45F86T1HU\ntYy0Ec6uLysCwD/L5MmT16xZEx8f/+DBg5qamiVLllRVVc2dO9fKymrevHknTpxobGykW7JYrK++\n+iokJGTy5MnXr19/P8OjKMrT09PJyWnQoEETJkwwMjJqaWlJT0+nJ59lZmZaWVnx+fyGhgYWizVg\nwIDp06cvWrRo3759dArMysricrmRkZFsNtvW1lYsFtfX1586dWrXrl30IwLFxcVmZmZWVlazZs1q\nbW1lsViKiooVFRVKSkqdhyEQCPLz8yURFt4cAhbAa/itWH3nWl7gaOuGmqeeC0JfUQSAfy5jY+PA\nwMDY2Njc3NzAwMCSkhIvLy97e/sVK1akpqYKhUInJ6ejR48uWbJkz549XT0Y+qE/sVj8xx9/nDhx\nYsWKFdra2u3t7StWrBgyZIi2tvbu3bttbGxUVVWPHz8+aNCgxMRELy8vW1tbMzOzb775hhDCYrG0\ntbUpisrOzt6+ffugQYN27dp17NixDRs2fP3110VFRfPmzevo6Fi3bl1ycvLRo0dFIlFKSgqXy/3y\ny7886FNQUGBkZKSmptbVp9zzYLPnv8Bmz72Qha3D2kMn37GTcN9phdmZjIwHALqP5ubmrKys1NTU\n7OxsFRUVd3d3JyenmJiYioqKnTt3Kioq0s2Y+usgEolSU1OXLVumrKxsYGCgoKDQ0tJSW1tbVlZW\nXl7ep0+fKVOmbNiwYf/+/YsXLy4rK7tw4cKCBQtUVFROnz5NX2Ty9vY+depUQ0PD559/npaWxmaz\nhwwZEhER8eDBg2XLlqWkpPj4+BQWFhJCjI2Nt27dOn78+JMnTy5fvry+vl5BQcHExOTSpUud17iP\njo4uLCzct2/fu59db4OA9RcIWL2Qp6dnXV3dq9s8efKktbXVwMDgZQ0GDBhAT6cFgJ7qzp07qamp\nqampVVVV6urq9+/fP3jwoJWVFXnnvw737t2je66srBwyZIiysrKdnd3ly5ebmpra2tpqa2u/+eYb\ne3t75k4F3gcs0wC93YkTJ17bJj4+vri4eMmSJe9hPADQPZmZmZmZmQUFBQmFwtzc3F9++WX8+PFu\nbm4LFix4i95KS0uTk5NTU1PLysqGDRs2duzY3bt39+/fX9LA09OTubGDFCBgAQAA/A1cLtfOzs7O\nzm779u1ZWVnbt2/Pzc0VCASysrKvfmNtbS19peratWuGhoZjx47dsGGDsbHx+xk2vGeY5A4gfQKB\nQE5Ozt3dvXPRz8+PxWI9fvy4paWFw+FwuVwOh6OoqOjm5lZeXi5p1traumTJEnoH3E8++eQNV/YC\nAEbY2dkdOnRozZo1L7sW3t7enpqaumLFCkdHx88//7ykpCQwMPDKlSuxsbGBgYFIVz0YAhaA9NGr\nJJeVlUkqt27dOnnypL6+Pr0TiJ6enlAoFIlEZWVlKioqQUFBkpbLli3r6Oi4efNmcXHxlClTpk+f\nTtdzcnKGDx+uoaFx6NAhSeMXFmkXL17U0dERYdFUgL/viy++OHjwoORLsVicn5+/cePG8ePHjxs3\nLjU1dezYsSkpKUlJSV999ZW1tTVWRe8NcIsQQPry8vJcXFwuXrxIURT9/E5YWNjYsWPp38J5eXmS\n7Wk1NDQ2bdpkampKt6Qo6vDhw/TuY6qqqvPmzfvjjz8IIY2NjdOmTYuPj6+qqvLx8fH19X1ZkZae\nnj59+nQrKysOhyOF8wf4h1NTU9PT0zt9+nR5eXnnaVUBAQHq6urSHh1IBwIWgPTl5+c7OjpmZ2fX\n1dX17ds3Ly+Pz+ePHDlSR0eH/DVgEUKUlJQ6P0Stqam5YMGCpUuXjhgxQllZef/+/YSQ2NjYSZMm\nDRs2rKWlpaamRiwWs9nsFxYJIZmZmZGRkR4eHvr6+u/91AF6iGHDhi1YsGDZsmWrV68eMmSItIcD\n0oeABdC1UlJSjh07pqGh8bIGfn5+9LYePB6vrKysb9++X3/9dURERHh4+OTJkwkheXl5Pj4+kval\npaWmpqZ0xmKxWGfPng0PD3d2dlZXV9+3b9/YsWMJIbGxscHBwYSQyspKAwMDOki9sHjlypXIyMhj\nx46NGTPG29u7a38WAD1Xnz59Vq9eHRAQIO2BQHeBgAXQtfh8/oABAz799NOXNVBXV//jjz+srKzo\ngFVXVycUCkePHl1QUGBtbd3Y2HjnUMVWywAAIABJREFUzh1ra2tJ++Tk5M77gpmYmBw+fLi1tfW7\n776bM2fOgwcPWCzWvXv3hg4dSgi5ceOGqakp3fL5Yk5OTmRk5JEjR1gs1s2bN7EbBsBbu3///ief\nfCLtUUA3goAF0LXk5eV5PF7nhPSMvLw8Ho/Xp08fOmAdO3YsKirq1q1b2tra6urq6enpPB5PMo2j\nrq5u586dqamphBCxWCwrK9vU1CQvL6+goBASEhITE0PX+Xy+lpYWISQpKcnLy+tlxVWrVqWmpsbF\nxdGdKysr19bWYk8MgLdw//59Pz8/aY8CuhE8yAAgZfn5+fQUKx6Pt2fPHg0NDTs7O8m8K8mLjo6O\nzMxMZ2fnhQsX0o92s9lsAwODHTt21NfXP3z4cOHChSEhISwWi8ViqaqqZmZmXr169fTp03SWemEx\nJSWFoiiKoqKjowMDAymKQroCeDt8Ph+zGKEzXMECkLK8vDz6+pahoeGtW7eOHDlCOuWqvLy8EydO\ncLlcLpc7cuTI5cuXd54p9fPPPwcHB4eFhenr68+bNy8kJIQQwmKxoqKivL299fT04uLiVFRUXlaU\nyMnJGTNmzPs8a4AeRigU4iFc6Ax7Ef4F9iKEF3qXrXKOHDnS1NQUGBjI+KgAoJtobW318PA4d+6c\ntAcC3QhuEQIAALyT0tJSIyMjaY8CuhfcIgR4SxRFTZkypbW19dXNKisrRSLRK7aU1tfXpxevAoB/\nqJKSEgQseAYCFsBbEovFpdW13x6Ifcd+wn2nMTIeAJCW+/fvI2DBM3CLEAAA4J3cv38f2zbDMxCw\nAAAA3klJSQkCFjwDAQuAMR0CwQJn27rqKkml7umTb+d4zR5pvmGhX0tjwyuKAPDPVVNT069fP2mP\nAroXBCwAZiQe2hfmM7XqEb9z8eCGNbxB5vsuX5dTUDy5c8srigAA0JMgYAG8HkVRra2ttc/p3IZn\nZj5t4ZfPvCvn/LmPXNxkZGVtx03ITUt+WREA/rmePHnyit3codfCU4QArycnJ5eamnr9+vXOxWcW\n6bWwdXjmXW0tzW0tzbrGAwkheiamtVWPX1jsq6HVtaMHgK6ENRrghRCwAF7P1dXV1dX1maJIJBr+\n8ejXvpfFYhFCKIoSi0SvLgLAPxHWaIAXwi1CgK4ir6gkr6hUfv9PQkhFaYn6AO2XFQHgnwuPEMIL\nIWABdBUWizVijEv+pfMURV3LSBvh7PqyIgD8c2ERLHghBCyALuS3YvWda3mBo60bap56Lgh9RREA\n/qHKysoMDQ2lPQrodjAHC+DtCdraHvPLOld2pGS1t7Z2Li5c/yP9oqG2pqG25vmiSNjxXgYLAF2i\nvb1dVlZW2qOAbgcBC+Atsdnssfa2BScOvLpZcXFxR0eHubn5yxp8NsGF6aEBwHvS0dHB5eIvKbwA\n/lkAvCUWi7Vt27bXNjty5EhTU1NgYOB7GBIAvGd8Pp/H40l7FNAdYQ4WAADAW8IiWPAyCFgA3QVF\nUfb29vLy8mKxWNpjAYA3goAFL4OABdBdnDp1ihBCUVRlZWXnukAgsLe3Lyoqor/cunUr66+Kiope\nWJTCOQD0MlijAV4Gc7AAuoWOjo5vvvnm559/9vb2Likp0dHRIYQIhcKwsLCkpKSysrJBgwbRLQMD\nA+fMmUO/3rFjx/37983MzIyMjJ4vSuVEAHoVrDIKL4OABdAtxMTE2NvbW1tb83i8kpISR0dHQkhH\nR4eDg4NQKCwsLGSz/3u9WU5OTk7u/9i797ia0v0P4M/epftVEkU7kjARQlSuXWSIGBVnksLkXgwz\nzKSQW6ZhajIhNZQRU9RMckQRJdLsiIxqbIWURBfSvfb6/bHO2b99UrntWrvd5/06f+z9Xc961nd1\nBt/Wei6yhJCLFy8mJyfHxcWxWKxWgwzeDkA38fz5cy0tbCcKrUCBBdCxnj59Gh0dnZSU1FaDTZs2\nGRoaBgQEXL16lRCip6eXn59PH5KXl581a1ZKSsqECRNanJWTk+Pj43P+/PkePXq0HwSAjkNRFH6Z\ngVahwALoWOvWrXNzc2tnpRxlZeXdu3fn5eX16dOHjri4uAg3uHHjhre3t3CkoqLC0dHxzJkz6urq\n7QcBoOO8evVKVVWV6SxATKHAAuhYMjIympqa7TQoLS0NCwt7/fq1srIyISQ8PDw0NFRwtKGhITMz\nc9y4ccKnbN++fdGiRYJRWe0EAaDjFBQUYAohtAUFFgDDfH193dzc6OqKEKKtrS14RUgIycrKGjBg\ngJqamiCSm5t74cKFu3fvCnfSahAAOhTWaIB2oMACYBKPx4uMjOTxeIKIjo5OcXFxbW2tvLw8IeT6\n9estBmBt3LjR09OzxSirVoMA0KHy8/PxzBjaggILgEmDBg0qLy8XjgwbNoyiKMHXdevWtTglPj7+\n7X5aDQJAhyooKLC1tWU6CxBTWGgUAADgY2AMFrQDBRYAAMDHqK6uVlRUZDoLEFMosAAAAD4Yn88X\nLP8L8Db8xwEAAPDBiouLdXR0mM4CxBcKLAAAgA+GNRqgfSiwAAAAPhhGuEP7UGABAAB8sIKCgoED\nBzKdBYgvrIPV0okTJ6qrq1s9VFxcHBIS0iLI4XCmT5/e8XkBAIAYyc/PX7JkCdNZgPhiCS9pCBMn\nTnz6stz2X67vf0rm2eiMjIwOywgAAMTRtGnTEhMTpaSkmE4ExBSeYLXUW6e/taPz+7fPvfzvjksG\nAADEU3NzM6oraAfGYAEAAHyYmpoaBQUFprMAsYYCq3WVZS8CNq52Mxu+1GJk8JYNtdVv6HhjfX3Q\n5nWu4z/znDk599ZfzCYJAACMePTokZ6eHtNZgFhDgdW6X77/WklV7aezl/1jLrDZ7KgD++l49MEA\nfnPTkZTbCzy+uXEB2+sCAHRHWKMB3gkFVuvu/5U+02WZmoZmz95adq7L0xP/M9AqJe7MnKWresjI\nTJg+y+277cwmCQAAjMjPz8caDdA+FFj/48mTJ/QH/c9GxP166Hnh46cPH8SEHKgoLSGENDc1vih+\neu1c7OJxQzfYW/1zJ5PRZAEAgBlYBAveCQXW/7h79y79Ye3ewJclxV/PsdqzcrHOQH0lVXVCSM2b\nN4QQwmIFnLsydtr0n75ehUUuAAC6ITzBgnfCMg2tk5VX+Pbn0B6ysoSQGxfi++kbEEKUVFTZUlIz\nFy1V19Sy/dfi6OCf6mtrmM4UAAA6W2VlpZqaGtNZgFjDE6zWRf60J+CbNdWvX798Vhz9y0+2/1pM\nCGGx2SPNJydEhtdWv0mMOtFP30BOQZHpTEFyUBRlZmYmJyfH5/OZzgUA2oPXF/BOKLBa57zBi89v\nXj5t7HY3R5sFi8bbzKTjy7x33buZ9tWk0beuXv76p0PMJgkSJjY2lhBCUVRJSYkgmJGRMWrUKE1N\nzYiICDry6tUrT09PDoejpqa2detWQcsbN26YmZn17Nlz4MCBcXFxnZw8QPdRWlqqpaXFdBYg7vCK\nsHVKqmqbDvz6dlyrP2dX5J+dnw9IvMbGxi1bthw/fnzBggX5+fna2tqEkKqqqvnz58fFxZWWljo7\nO7u4uBBCvvzyywkTJmRmZr58+XLEiBFbt25ls9nPnj2bNWvWiRMnbGxswsPDly9fPnv2bKbvCUAy\nYY0GeB8osP6HtLR07m3ut1/Yvv8pvVXwlhBEICwszMzMzMTEhMPh5OfnW1hYEEKioqLs7OxGjhxZ\nU1NTXl7O5/PZbHZaWtqyZcukpKSio6PHjBnDZrMJIampqa6urra2toSQkSNHYgcPgI6Tn5+PAgve\nCQXW/1BUVHxV9rKto1OnTk1OTu7MfEAChIaGBgcHa2hotNVg586dn332WUBAwNWrVwkhenp6+fn5\n9KGoqKg1a9YQQkpKSnR1delaasmSJXPnzqUbXLp0if7g6Ojo6OhICOHxeKtXr/b39+/IewLo1vLz\n88eOHct0FiDuUGABdCwFBYUVK1a4u7u308bX1zcvL69Pnz70V/pVICHkwYMHw4cPJ4TcvXvXwMCA\nEBIcHJyVlZWbm6uhofHbb7+tXr06JyeHblxeXu7j43Pr1q1Dhw6NGDGiA28JoHsrKChwcnJiOgsQ\ndxjkDsCw0tLSsLCw169fUxRFUdSxY8foJ1h8Pr+wsJAeS5uQkEA/oDpw4EBISIihoWGvXr3s7e1r\na2vpTs6dO2dvb29paZmWlobqCqBDPXnyhMPhMJ0FiDs8wQJgmK+vr5ubm7KyMv1VW1ubLrBYLJaK\nikpaWlrPnj3j4+Ppt35lZWWpqakDBgx48uSJu7v7N998QwhJTk7etm3b+fPn5eTkqqurWSyWoiKG\nBgJ0lIaGhh49ejCdBYg7FFgATOLxeJGRkTweTxDR0dEpLi6ura2Vl5fft2/fggUL+vXrFxMTQ1dg\nR44c2bBhw4YNGwYMGLBx40b6PcUPP/zA5XI1NTXpHiwsLFJTUxm5HQCJh+oK3hMKLAAmDRo0qLy8\nXDgybNgwwRqGrq6urq6uwkdnz5799voL58+f78gcAeD/FRYW6urqMp0FdAEYgwUAAPC+sEYDvCc8\nwQL4SHw+X7Ovdr9Bg9+n8S/hv7V16FXZy0e590WXFwB0IGzzDO8JBRbAR6IoSkffYOuvv39iPz4u\n80WSDwB0goKCAnt7e6azgC4ArwgBAADeV0FBAZ5gwftAgQUAAPC+SktLe/fuzXQW0AWgwAIQmcaG\nhpWWppUvSwWRyrIX29wcF48b6rfKtabqdTtBAOgSBJN8AdqHAgtANM5FhHo7zy0tKhQOhvtt5wwe\nGnrtjqy8wumDge0EAUD8VVZWqqmpMZ0FdA0Y5A7wbo8fPz5//nyLDZv5fL7wV47hUC1dzp4ViwUR\niqIyLl3YcuREDxkZU+sZJwN/cPnW++2guqZWJ90GAHyagoICrNEA7wkFFsC7Xbly5fLly1ZWVsLB\nFgWWkal5i7Pqaqrraqp1Bg4ihPTTN6gofd5qEAUWQFeBRbDg/aHAAng3dXX18ePHu7u7Cwebm5uD\nI06881wWi0UIoSiK39zcfhAAxFxBQcGQIUOYzgK6BozBAugocgqKcgqKxQUPCSHPHuVr9OnbVhAA\nugSs0QDvDwUWQEdhsVhjp9lkXr1EUdTt1OSxlrZtBQGgS8AYLHh/KLAAOpDr5q15t7nuU0xel5c5\nrPRsJwgA4q+mpkZeXp7pLKBrwBgsgI9X/rwkMep/Nhlc4fvDX5cvCkcsZs6hP1z795+tBmveYCks\ngC6Az+ez2XgqAe8LBRbAR5KSkjqwz7+xsbH9ZteuXaurq2sxA1GY/f4fRZ0aAIheUVGRjo4O01lA\nl4ECC+Djvc+er42NjW/evHFwcOiEfACg42CEO3wQPO0EAAB4NyyCBR8EBRaAuKAoyszMTE5OrsUS\npgAgDvAECz4ICiwAcREbG0sIoSiqpKREEMzIyBg1apSmpmZERAQdefXqlaenJ4fDUVNT27p1q3AP\nDQ0NZmZmubm5nZk2QDeRn5+PAgveHwosALHQ2Ni4ZcuWoKAgXV3d/Px8OlhVVTV//vyjR4+eOHFi\n48aNdPDLL7/s3bt3ZmZmenr6nj176MddTU1N3333nampaW5u7uDBgxm7DQDJVVxcrK2tzXQW0GVg\nkDuAWAgLCzMzMzMxMeFwOPn5+RYWFoSQqKgoOzu7kSNH1tTUlJeX07PE09LSli1bJiUlFR0dPWbM\nGHreeGNjo7m5eVNT07179zCTHKAjYJkG+CAosAA61ps3b/z9/UNCQlo9ymazDx8+bGBgEBAQcPXq\nVUKInp6e4AlWVFTUmjVrCCElJSW6urr0X+5LliyZO3cu3eDSpUv0B3l5+VmzZqWkpEyYMKGj7wig\nGyovL+dyubt27Vq2bJmWFjZoh3dDgQXQsdzd3VvsEv02X1/fvLy8Pn360F9dXFzoDw8ePBg+fDgh\n5O7duwYGBoSQ4ODgrKys3NxcDQ2N3377bfXq1Tk5OYJ+bty44e3t3SG3AdC9JScne3t7Dx061MXF\nhaKotWvXzpo1i961HaBVeNoJwLDS0tKwsLDXr19TFEVR1LFjx+gnWHw+v7CwkP5dOSEhwdHRkRBy\n4MCBkJAQQ0PDXr162dvb19bWCvppaGjIzMwcN24cUzcCIMGSkpKmT59uZ2d34cKFAwcOpKWljR8/\nfu/evWVlZUynBmIKBRYAw3x9fd3c3JSVlemv2tradIHFYrFUVFTS0tJu3boVHx9PF1hlZWWpqal8\nPv/Ro0fu7u7ffPONoJ+srKwBAwaoqakxchcAki07O3vEiBH058GDB/v5+V29erVv374zZ850dHS8\nfv06s+mBGEKBBcAkHo8XGRnp4eEhiOjo6BQXF9fW1rJYrH379i1YsGDJkiUxMTF0BXbkyJFdu3Zp\namrOnz9/yZIlq1atEpx4/fp1DMAC6AgFBQW6urotXgjKycm5uLikp6dv2rQpPDx84sSJISEhNTU1\nTCUJ4oZFURTTOXQZU6dOTU5OZjoLYEBcXByPx/v666+ZTgQAGHDkyBFpaWk3N7d22pSWlh49evTU\nqVOTJk3y9PTEilmAQe7Q3RUXF9fV1bXfpqSkpKysTDC5723y8vJ9+/YVdWoAIBYSExN/+umn9tv0\n7t1706ZN33zzzeXLl7///vvKysqlS5fOnTtXWhr/znZTeIL1AfAESyJp9ecYm036xE6y0689e1wg\nknwAQKzw+fyJEyempaV90FkPHz48cuRIQkKCvb396tWrNTU1Oyg9EFuorKG709Tut8L3h0/sxMdl\nvkiSAQBxc/v27dGjR3/oWfr6+n5+ftu3b4+Li3N1dVVUVHR3d7eysuqIDEE8YZA7AABAmxITEz+6\nMJKVlXVwcDh37tymTZuio6PHjh0bGBhYXV0t2gxBPKHAAvgfjQ0NKy1NK1+WCiKVZS+2uTkuHjfU\nb5VrTdXrdoIAIHmuXr06bdq0T+zExMTk8OHDSUlJ8vLy1tbWy5cvv3fvnkjSA7GFAgvg/52LCPV2\nnltaVCgcDPfbzhk8NPTaHVl5hdMHA9sJAoCEqa2tra+vFyxT94lUVVXd3d2vX7++aNEiX1/fqVOn\nRkdHNzY2iqRzEDcosAD+H8dw6PxV64QjFEVlXLow3mZmDxkZU+sZfyVfbCsIAJLn2rVrEydOFHm3\nFhYWUVFRkZGR+fn55ubmmzdvLiwsfPdp0KVgkLska2pqWrt2rZ6eHtOJdBlGpuYtInU11XU11ToD\nBxFC+ukbVJQ+bzWoronNXwEkUGJiop2dXQd13rdv302bNq1fv/7PP//86quvZGRkPDw8LC0tscWh\nZECBJcnq6+uzsrK++OILphMRa8dj4t7Zhv77jqIofnNz+0EAkCQZGRm7d+/u0EvIyMg4ODg4ODjk\n5eUdPXrU29vb3t7e3d1dXV29Q68LHQ0FloTT0NDAxOB38NrazkE5BUU5BcXigoeGo8Y8e5Sv0adv\nW0EAkDClpaWqqqqdtlKooaGhn59fVVXVyZMn7ezsPvvss1WrVhkbG3fO1UHkMAYLoD0sFmvsNJvM\nq5coirqdmjzW0ratIABImMuXL3/6/MEPpays7O7ufu3aNXd39+DgYAsLi5CQkHfuNgFiCAUWwDu4\nbt6ad5vrPsXkdXmZw0rPdoIAIEmSkpKsra2Zujq9ssOZM2cqKiqmTJmyefPmR48eMZUMfAS8IgRo\n6UxusfBXtV69t4dHt2jTahAAJElOTs6wYcOYzUFLS0uwxeG6detqamrWrl07a9YsDIQXfyiwoLsr\nKyn+bd+uT+ykorREJMkAgJh48ODBoEGDmM7iP9hstpWVlZWVFY/HCw0N3b17t729/bJlyzQ0NJhO\nDdqEAgu6u8RzZ+vr69tvk5KSUlhY+OWXX7bVQM5zpajzAgAmJSYmMvh+sC2DBg3y8/Orr6+Pi4tz\ndHTU0NBYt26dmZkZ03lBK1BgQXdnZGT0zjZFRUUURZmYmHRCPgAgDi5duhQcHMx0Fq2jtzh0cHDI\nzMwMCQnZvHmzs7Ozs7OzgoIC06nB/8MgdwAAgP/R1NRUWlqqpSXuCwjTA+HPnj1LCLG2tvb09MzP\nz2c6KfgPFFgAzGtqapKXl1+/fr0gEhMTY2NjQwjh8/mqqqr379+n4w8ePDA0NDx48GBb8c5PHkDy\ncLncsWPHMp3F+6K3OLx27Zqdnd33338/Y8aM6OjopqYmpvPq7lBgATDv/v37zc3NERERgtVuuFzu\nmDFjCCE8Hq+pqcnQ0JAQkp6ebmNj89NPP61cubKtOEVRRkZGrP9yc3MTXKWhocHMzCw3N1f40q0G\nAbq5pKSkLrdEM4vFsrKyOnXqVFhYWH5+voWFxbZt2168eMF0Xt0XCiwA5nG53NmzZ2toaMTExAgi\n9JAvLpc7atQoKSmpP//809nZOSYm5vPPP28nnpaWpqGhUfVfhw4dIoQ0NTV99913pqamubm5gwcP\npi/RahAACCEpKSmTJ09mOouPpK2tvWnTppSUlM8++8zV1dXJySkpKYnppLojDHIH6HAVFRXtHFVX\nV8/MzBw7duz48eNDQkL+9a9/URSVmZlJP8HKzMw0MTEJDg7etWtXenp6//796bPaih8+fLi6ulpH\nR6dnz57BwcEzZswghDQ2Npqbmzc1Nd27d4/N/s+vVa0GAeDNmzeEEEVFRaYT+SSCLQ5zc3OPHTu2\nZ88eJyenL7/8sqvfVxeCAgugYwUEBISGhvbt2/p+hWw228/Pj8vlzp07d+TIkVu2bMnLy5ORkWGx\nWLq6uoQQLpdbVlaWn59fUVHRLLSrdKvxmpoaHo/n5eU1efLkwMDAZcuWFRUVEULk5eVnzZqVkpIy\nYcIEQQ+tBgEgJSVl0qRJTGchMkOGDBFscWhrazts2LC1a9e+z+xp+FQUvLcpU6YwncKHefPmzcyZ\nM5nOQhL8+eef+/bt+7hzT5w4cfjw4XYaNDQ0yMnJlZeXUxS1cOHCDRs2REdHT58+naKo5uZmJSWl\nr7/+urm5+csvv/T09KRPaSsurLy8XF5evrm5WRCxsLC4cOFCi2atBgG6s/Xr19+4cYPpLDoKl8td\ntGjR1KlTo6KiGhoamE5HkuG9AADD/v77b21tbXV1dULI8uXLjx07lpaWRg/A+ueff1gslr+/P5vN\nXrt2bWhoKP22sa14bGysYJI2j8cbPXq04N1fQ0NDZmbmuHHjhC/dahCgm+taUwg/lImJSURExMmT\nJ/Pz8ydNmrR58+bCwkKmk5JMKLAAGCYYbkUImTRpUq9evQ4fPiwYgCUokkxNTT/77DN60Hpb8VOn\nTvn4+FRXVz958sTDw2PHjh2Cq2RlZQ0YMEBNTU340q0GAbqzZ8+eaWpqSklJMZ1Ix6K3OExLS7Oy\nslq7du2cOXOSkpIoimI6L4mCAguAYYIJg4QQFou1YsWK2tpawRRC4eXjPTw8fv755/r6+rbiO3bs\n4PF4Ojo68+bN27Zt29SpUwVtrl+//vZYq1aDAN3ZpUuXLC0tmc6ik9BbHP7xxx/79u1LSkqysLDY\nu3dv+5Ny4P2xULG+v6lTpyYnJzOdxQeorq52cnKKj49nOpEuLy4ujsfjff311x9xbmRk5Js3b9zd\n3UWeFQCInKurq5eXl4GBAdOJMKCuru7s2bNhYWEcDmfVqlXGxsZMZ9S1YRYhwMe7du2aYGnQtty7\nd6+2traddWjU1NQErwgBgEEURfF4vO5ZXRFC5OTkBFscBgcH5+TkODs7u7i4yMnJMZ1al4QCC+Aj\nNTc3O3y5yMrhy3c1ZBMpxfCEK20dvnYutiDnb5GmBgAfIycnZ8iQIUxnwTx6i8PKysrw8PBp06ZN\nmjRp5cqVHA6H6by6GBRYAB9PU6ff3K9Wf2Int1O70ntnAAmWlJRkbW3NdBbiQk1NzdPTc+3atZcv\nX960adOrV69WrVo1a9YsFov10X1SFGVubn7r1q2amhqJX99Ywm8PAADgPV26dKnLbUHY0eiB8KdO\nnQoNDc3MzKQHwpeVlX1cb7GxsYQQiqJKSkqE4y02RQ0KCmL9r9zc3H/++WfmzJm9e/fW19ePi4tr\n0fOVK1e0tbWFV2NmHAosAJFpbGhYaWla+bJUEKkse7HNzXHxuKF+q1xrql63EwQAZjU2NlZWVmpo\naDCdiJjS0dHZtm1bcnLywIEDXV1dHR0dr1+//kE9NDY2btmyJSgoSFdXV7BiX6uborq7uws2VN27\nd++KFSs0NTVtbW3Xrl377NmznTt3urq68vl8Qc8pKSlOTk7GxsZitb4GCiwA0TgXEertPLe06H+W\n7Av3284ZPDT02h1ZeYXTBwPbCQIAs27evDl+/HimsxB39BaHZ8+e9fX1/f3336dMmRISElJTU/M+\n54aFhZmZmZmYmHA4HEGBRW+KamVlZWpqKnhpKCsrq6SkpKSkdP369eTk5J9//jk2Nnby5Mm2trZS\nUlJTp06tr68XdJuWlubv729vb29ubi7y+/0UKLAARINjOHT+qnXCEYqiMi5dGG8zs4eMjKn1jL+S\nL7YVBADGJSUl4f3g+xsyZEhgYCD9qs7W1tbT01NQM7XqzZs3AQEBu3btIoTo6ekJGtOborJYrLfX\n5MvJyfHx8YmMjOzRo0dDQ0NWVtbTp0/z8/OXLFkydOhQuhpLT0/39/c/depUVlaWuBVYGOQO8G7P\nnz8/ffp0Xl6ecLDFGnJGpi3/bNfVVNfVVOsMHEQI6advUFH6vNWguqZWx2YPAO/h2rVrmzdvZjqL\nLkZFRcXd3d3d3T0zM9PX17e0tNTNzW3evHlvv6rbv39/Xl5enz596K8uLi7CR2/cuOHt7S0cqaio\ncHR0PHPmDL2N2KJFi86ePWtoaGhsbKynp9e7d29CSEZGhr+/f2RkJIvFys7OFrddv1BgAbzbwoUL\nJ0yYoKCgIBxsbm7+YvGSd55Lz7ihKIovNPqy1SAAMOXVq1c9evTAgk8fzcTE5NixYyUlJeHh4RMn\nTrSxsVm9erWmpiZ9tLS0NCymVl4IAAAgAElEQVQs7PXr18rKyoSQ8PDw0NBQwbmtboq6ffv2RYsW\nCUZlvXr16ty5c2w2u6mpycjIKCAggBDi5eWVlJQUExNDt1FSUqqoqBCfvb/wihDg3ZSUlIyMjAa+\npf2z5BQU5RQUiwseEkKePcrX6NO3rSAAMOvq1atTpkxhOosur0+fPps2bbp27ZqJiYmbm5uTkxO9\nxrKvr6+bmxtdXRFCtLW1hd8nvr0pam5u7oULF9avXy+IGBkZnT17tr6+fvPmzQMGDLCxsSGEJCYm\nUhRFUVRAQIC7uztFUeJTXREUWAAdh8VijZ1mk3n1EkVRt1OTx1rathUEAGYlJiZiBSxRYbPZdnZ2\n8fHx3333XXR0tIWFRXh4uIeHh6CBjo5OcXFxbW0t/fXtTVE3btzo6enZo0cPQSQgIOCrr74aOHBg\nTU3N77//3mINrYyMDHF7P0iwF+EHwV6EIKy5uXnUxCnbjkYJB78Yoh12LUutV2/6a+XL0p82rC5+\n9HCQkbHH3p/llZTfDu5Z5XbvZhoDNwAA/zVp0qQrV65I/NKXjHj16tXTp08/++wzphPpbBiDBSBK\nZ3KLhb+q9eq9PTy6RZtWgwDAlMLCwr59+6K66iCqqqqqqqpMZ8EAFFgAH4nFYhU++GffuuWf2M+r\nshciyQcAPg4WaICOgAIL4COx2ezHDx+8c2eG06dPV1dXL168uK0GwuMMAKDzJSUl7d69m+ksQNKg\nwAL4eCoqKu9so6ioSFEUvZQLAIgbiqIeP37M4XCYTgQkDV45AwBA95WdnT1ixAimswAJhAILgHlt\n7RKfkZExatQoTU3NiIiI9+mHoigjIyPB/vNubm6EkNevXy9btkxLS0tbW9vX17ej7gGga0pMTMQA\nLOgIeEUIwLDy8nJbW9vg4GBra+uoqChXV9eXL1+y2eyqqqr58+fHxcWVlpY6Ozu32FmiVWlpaRoa\nGlVVVfRXenSXs7OzhobGvXv3qqqqhgwZsmXLFsyWAhC4cuXK0qVLmc4CJBD+ngVgWExMTKu7xEdF\nRdnZ2Y0cOdLCwqK8vJzP5wuflZGRYWlpqampaWhoePHif3aMPnz4cHV1tY6OzvDhw1NTU2VlZZua\nmpYsWfLLL780NzefO3du7NixqK4ABOrr69+8eSNWy3+DxMATLICO1dzc/Pjx48zMzLYa1NTU0LvE\nNzQ0rFmzRrBLfFRU1Jo1awghJSUlurq6woVRZmamg4NDRESEubn50aNHXVxcnj17Vltby+PxvLy8\nJk+eHBgYuGzZsqKiImlpaXt7+1u3bpmYmBBCEhISOv6OAbqM69evm5u33KYdQCRQYAF0LFVV1dzc\n3JCQkFaPstlsNze38+fPt9glnhDy4MGD4cOHE0Lu3r1rYGAgfJaPj8/evXsnT55MCHFwcFi/fj1F\nUQoKCjdu3KAbrFu3zt/fn8/n02XZqFGjKioqjh49umzZssLCwo67WYCuJSkpCTvkQAdBgQXQsWbP\nnj179ux2Gjx9+vTtXeL5fH5hYaGWlhYhJCEhwdHRUfiUmzdvhoWF0Z9TU1NNTEzYbHZsbKyxsTG9\nBTWPxxs9ejSbzf7yyy9DQ0Pl5eXV1NTmzJkTFBTUUfcJ0AXduHFj69atTGcBkgmjMQAY1uou8SwW\nS0VFJS0t7datW/Hx8XSB9eLFi4yMDEKInJxcTExMc3PzzZs3PT09d+3aRQg5deqUj49PdXX1kydP\nPDw8duzYQQi5ePHiyZMnm5ubi4qKVq9eLbzfKkA3V15erqCgICMjw3QiIKEoeG9TpkxhOoUP8+bN\nm5kzZzKdBbzD0aNHNTU1tbW1V65c+erVK+G4hoaGsbHxzZs36Yifn5+5uTlFUQkJCYMGDerZs+fE\niRNTUlLoo3l5eaampqqqqiYmJgkJCXQwOjpaT09PXV3d2Ng4LCyMz+d37s0BiK/Tp0/v27eP6SxA\nYrEoimK6xusypk6dmpyczHQWH6C6utrJySk+Pp7pRAAAxM7KlStXrVpFj3QEEDm8IgQAgO7o3r17\nRkZGTGcBEguD3Luq77//PiMjg8VitdOmubk5Ozu7nTkydXV14eHh9LBoAIDuo6CggMPhtP9XKMCn\nQIHVVWVnZy/xPySnoPgpnUT47ywrK0OBBQDdTVJSEnbIgQ6FV4QAANDtYAUs6GgosAAAoHvh8/nF\nxcU6OjpMJwKSDAWWhGhsaFhpaVr5slQQqSx7sc3NcfG4oX6rXGuqXrcTBADoVm7fvj169GimswAJ\nhwJLEpyLCPV2nlta9D9boIT7becMHhp67Y6svMLpg4HtBAEAuhUMwIJOgEHuH6ChoSE/P5/pLP6j\npqZG8JljOFRLl7NnxWJBhKKojEsXthw50UNGxtR6xsnAH1y+9X47OHbadCZyBwBg0pUrV1atWsV0\nFiDhUGB9AHV19b179zKdxX/weDzBZyPTlrvB19VU19VU6wwcRAjpp29QUfq8rSAAQLdSW1vb2Nio\nrKzMdCIg4VBgfQCxWhLdzs7unW3oJV4oiuI3N7cfBADoJq5duzZx4kSmswDJhzFYkklOQVFOQbG4\n4CEh5NmjfI0+fdsKAgB0KxiABZ0DBZZkYrFYY6fZZF69RFHU7dTksZa2bQUBALqVjIwMU1NTprMA\nyYcCS2K5bt6ad5vrPsXkdXmZw0rPdoIAAN3EixcvVFVVpaUxPAY6HP4j68IK7t+TkZcXfP3hTELZ\n85Ky5yWCiMu33vSH4scFrQaF180CAJB4ly9ftrS0ZDoL6BZQYHVVCxcuvHs3vf02jY2Nf/755/z5\n89tqMH6IPofDEXVqAABiKikpacOGDUxnAd0Ci6IopnOAjlJdXe3k5CRWkx8BABhkYWFx7do1prOA\nbgFjsAAAoFt48OCBgYEB01lAd4ECCwAAuoXExEQs0ACdBgUWgLigKMrMzExOTo7P5zOdC4AEunTp\nEgos6DQosADERWxsLCGEoqiSkhLheENDg5mZWW5uLv01KCiI9b9yc3P/+eefmTNn9u7dW19fPy4u\nrkXPV65c0dbWbsba/dCNNTc3v3z5UktLi+lEoLtAgQUgFhobG7ds2RIUFKSrqyvYU7ypqem7774z\nNTXNzc0dPHgwHXR3d6/6r717965YsUJTU9PW1nbt2rXPnj3buXOnq6ur8DOwlJQUJycnY2NjKSkp\nBm4MQDxwudyxY8cynQV0I1imAUSvoaEhNjZWQ0OD6US6hjFjxqipqYWFhZmZmZmYmHA4nPz8fAsL\nC0JIY2Ojubl5U1PTvXv32Oz//DokKysrKytLCLl48WJycnJcXFx4ePjkyZNtbW0JIVOnTq2vrxd0\nnpaW5u/vb29v379/fyZuDkBcYAAWdDIUWCB6PB5v//798+bNYzqRLoDNZuvr60tLSwcEBFy9epUQ\noqenJ3iCJS8vP2vWrJSUlAkTJrQ4MScnx8fH5/z58z169GhoaMjKynr69GlDQ8OaNWuGDh1KV2Pp\n6en+/v6nTp2aNm3aggULOvnWAMRKamrq119/zXQW0I2gwALRk5WVHT58+KZNm5hOpMvw9fXNy8vr\n06cP/dXFxUX46I0bN7y9vYUjFRUVjo6OZ86cUVdXJ4QsWrTo7NmzhoaGxsbGenp6vXv3JoRkZGT4\n+/tHRkayWKzs7Oxx48Z11t0AiJ03b96wWCwFBQWmE4FuBGOwABhWWloaFhb2+vVriqIoijp27Jjg\nCRYhpKGhITMzs0V5tH379kWLFglGZb169ercuXPV1dUpKSm3bt2iH1Z5eXnFxMQoKioqKCjU1tYq\nKSlVVlZ25n0BiI+UlJRJkyYxnQV0LyiwABjm6+vr5uamrKxMf9XW1hYusLKysgYMGKCmpiaI5Obm\nXrhwYf369YKIkZHR2bNn6+vrN2/ePGDAABsbG0JIYmIiXbEFBAS4u7tTFCXcCUC3kpSUhAFY0MlQ\nYAEwicfjRUZGenh4CCI6OjrFxcW1tbX01+vXr7cYgLVx40ZPT88ePXoIIgEBAV999dXAgQNramp+\n//13wXB4WkZGBt4PQjd369YtExMTprOA7gVjsACYNGjQoPLycuHIsGHDhHcIXbduXYtT3t5c0tXV\n1dXVta1LnDhx4lOzBOjKnj17pqmpiWVKoJPhCRYAAEiyS5cuWVpaMp0FdDsosAAAQJIlJSVZW1sz\nnQV0OyiwAABAYlEU9fDhQ319faYTgW4HBRYAAEisnJycoUOHMp0FdEcosAAAQGLh/SAwBQUWAABI\nrMuXL0+bNo3pLKA7QoEFAACSqbGx8dWrV9h4HhiBAgsAACTTzZs3x48fz3QW0E2hwAIAAMmEHXKA\nQSiwAABAMqWlpZmbmzOdBXRTKLAAAEACvXr1SkZGRk5OjulEoJtCgQUAABLo6tWrU6ZMYToL6L5Q\nYAEAgATCACxgFgosAACQQHfu3Bk5ciTTWUD3hQILAAAkTWFhoba2NovFYjoR6L5QYAEAgKS5dOkS\n3g8Cs1BgATCvoaFBVlZ21qxZwkFXV1cWi/X8+fOamhopKSlpaWkpKSkFBYWZM2cWFxcLmtXW1n79\n9dccDkdJSWnSpElcLrfT0wcQO9iCEBiHAguAeffu3WOz2Y8fPxZE7t+/f/r06f79+2tpaWVlZfXr\n16+pqam5ufnx48fKysqrV68WtNy4cWNjY2N2djaPx5szZ46TkxMdz8jIGDVqlKamZkREBB2hKMrI\nyIj1X25ubsI5XLlyRVtbu7m5ueNvF6BjURT15MkTXV1dphOBbk2a6QQAgHC5XBsbmytXrlAURY8a\n8fb2trKyYrPZ9NExY8bQLTU1Nffv329gYEC3pCjqxIkT+fn5KioqKioqy5cv//vvvwkhVVVV8+fP\nj4uLKy0tdXZ2dnFxIYSkpaVpaGhUVVXRXfXo0UOQQEpKipOT0+jRo6WkpDr53gFELjs7e/jw4Uxn\nAd0dnmABMC8zM9PCwkJeXr6yspIQwuVyCwsL+/XrR9dVwgUWIURRUVF46G7v3r1XrlyZkZFBUZSS\nktKvv/5KCImKirKzsxs5cqSFhUV5eTmfzyeEHD58uLq6WkdHZ/jw4ampqbKysnQPaWlp/v7+9vb2\nWPMaJENiYiLeDwLj8AQLOlt5eXlAQEC/fv2YTkRczJgxg8vlLliwgMPhPH78WF1d/fvvv9+1a5eP\nj8/s2bMJIVwu19nZWdD+0aNHBgYGdI3FYrHOnz/v4+NjaWmpoaERGhpKD+yNiopas2YNIaSkpERX\nV5fNZtfU1PB4PC8vr8mTJwcGBi5btqyoqIgQkp6e7u/vf+rUqWnTpi1YsICZHwGASF25cmXZsmVM\nZwHdHQos6GxPnz69efMmHuALUBT1999/jx49mi6wKisrm5qapkyZkpWVZWJiUlVVlZeXZ2JiImh/\n8eLF8ePHC77q6+ufOHGitrZ29+7dbm5uT548YbFYDx48oH/Cd+/eNTAwIIQoKCjcuHGDPmXdunX+\n/v58Pp/L5fr7+0dGRrJYrOzs7HHjxnXurQOIXn19fU1NjaqqKtOJQHeHAgs6W48ePfT19R0cHJhO\nRFxwuVwOh6OqqkoXWKdOndq3b9/9+/f79u2roaGRkpLC4XA0NDToxpWVlQcPHkxKSiKE8Pl8GRmZ\nN2/eyMnJycvLe3h4hIWF0fHCwkItLS1CSEJCgqOjIyEkNjbW2Nh44MCBhBAejzd69Gg2m+3l5ZWU\nlBQTE0N3rqSkVFFRoaamxsjPAUAkbty4gZfdIA4wBguAYZmZmfQQKw6HExISoqmpOWHCBMG4K8GH\nxsbGtLQ0S0vLVatW0XUSm83W1dUNDg5+9erV06dPV61a5eHhQc8QVFFRSUtLu3XrVnx8PF1gnTp1\nysfHp7q6+smTJx4eHjt27CCEJCYmUhRFUVRAQIC7uztFUaiuoKvDDjkgJvAEC4BhXC6XfgOop6d3\n//79yMhIIlRXcbnc6OhoaWlpaWnpcePGffPNN8IjpY4fP75mzRpvb+/+/fsvX77cw8ODEMJisfbt\n27dgwYJ+/frFxMQoKysTQnbs2OHi4qKjozNo0KBdu3ZNnTpVOIeMjIxp06Z15l0DdJD09HQfHx+m\nswAgLIqimM4BOkp1dbWTk1N8fHwnX/fhw4d79uwJDQ1t9WhOTk5QUFBwcHAnZwUAEq+8vNzV1TUu\nLo7pRADwBAs+UHR0dH5+fvttysrK7t69u3fv3laPlpaWZmVl1dfXC5YJAAAQieTkZDyLBTGBAgs+\nTFBQ0DS3Ne23Ueg94POhY2rbOKrce0B1+l+lpaX9+/cXeXoA0J0lJSUJ73MAwCAUWPBhpKSkRphN\n/MROrvwZLZJkAACEnT9//pdffmE6CwBCMIsQAAAkw/Xr13V0dBwcHMrLy5nOBQAFFnyaxoaGlZam\nlS9LBZHKshfb3BwXjxvqt8q1pup1O0EAABE6cOBAeHj4rl27FixYkJaWxnQ60N2hwIKPdy4i1Nt5\nbmlRoXAw3G87Z/DQ0Gt3ZOUVTh8MbCcIACAqxcXFFRUVgwYNGjJkyJkzZw4ePLh3715MkwcGocCC\nj8cxHDp/1TrhCEVRGZcujLeZ2UNGxtR6xl/JF9sKAgCI0K+//vrVV1/Rn5WVlX/77Td1dXVHR0d6\nA3WAzodB7pJMSkrq3r17ot1V/s6dO4LPRqYt96Ooq6muq6nWGTiIENJP36Ci9HlbQQAAUWlqavr3\nv/+9efNm4aC7u/vo0aMdHBx+/PFHY2NjpnKDbgsFliSTk5N79OiRaPtssQJ4q1gsFiGEoih+c3P7\nQQCAT3f27FlbW1tp6Zb/oo0ZM+bkyZPLli2bOXOm4PkWQOfAK0IQJTkFRTkFxeKCh4SQZ4/yNfr0\nbSsIACAqoaGh7u7urR7q1atXbGxseXm5m5tbTU1NJycG3RkKLBAlFos1dppN5tVLFEXdTk0ea2nb\nVhAAQCTy8vJUVFT69OnTVgMWi7Vp06ZFixbNnTuXx+N1Zm7QnaHAAhFz3bw17zbXfYrJ6/Iyh5We\n7QQBAD7d4cOHV6xY8c5m06ZNCw0NXb169enTpzshKwBs9gwfZurUqWsPRX5iJz9v8jgetB9b5QDA\nJ6qurp4xY0ZKSsp7tq+vr9+0aRMhxN/fv0ePHh2ZGnR3KLDgwzg5OeU+ffaJnbwuL/sr9WqvXr1E\nkhIAdFthYWF1dXUfuv/gyZMnIyIiwsLCtLW1OygxABRYIHoPHz7cs2dPaGhoq0dzcnKCgoKCg4M7\nOSsAkDyTJ0+Oj49XVlb+0BPz8vJWrlz53XffiXYhGwABjMECAIAuKT09/bPPPvuI6ooQYmhoGBcX\nFxERsW3bNj6fL/LcAFBgAQBAl3Tw4ME1a9Z89OlKSkrHjx/X1taeP38+9ocGkUOBBSAuKIoyMzOT\nk5PD79MA71RaWvrs2bNhw4Z9Yj/u7u5eXl5z585NT08XSWIANBRYAOIiNjaWEEJRVElJiSB448YN\nMzOznj17Dhw4MC4ujg5mZGSMGjVKU1MzIiKiRSdXrlzR1tZuxnL5IOmOHTu2ZMkSkXRlYmISExPj\n5+cXGIit6EFkUGABiIXGxsYtW7YEBQXp6urm5+fTwWfPns2aNcvHx+fly5fe3t7Lly8nhFRVVc2f\nP//o0aMnTpzYuHGjcCcpKSlOTk7GxsZSUlIM3ANAZ+Hz+WfPnv3iiy9E1aGGhkZsbGxdXd3ixYur\nq6tF1S10ZyiwAMRCWFiYmZmZiYkJh8MRFFipqamurq62trZsNnvkyJF02RQVFWVnZzdy5EgLC4vy\n8nLB+8S0tDR/f397e3tz85abcANImISEhClTpoh2ISt6wXdXV9eZM2f+/fffIuwZuids9gzwqa5c\nuRIREdG7d++PO33x4sX9+/cPCAi4evUqIURPT09QYDk6Ojo6OhJCeDze6tWr/f39CSFRUVH0wN6S\nkhJdXV02m00ISU9P9/f3P3Xq1LRp0xYsWCCS+wIQW3v37g0PD++InqdOnWpgYLB06VI3Nzf8UYJP\ngQIL4FM9ffpUU1PTwcHh407X0dHZv39/Xl6eYDM1FxcXwdHy8nIfH59bt24dOnRoxIgRhJAHDx4M\nHz6cEHL37l0DAwNCSEZGhr+/f2RkJIvFys7OHjdu3KfeEoAY8/Pzy8nJUVVV7aD++/Xrd+7cuS1b\ntiQnJwcFBcnIyHTQhUCyocAC+FRycnL9+/c3MTH5uNNLS0vDwsJev35NL+cTHh4uWKP13Llze/fu\nXb9+fVBQEIvFIoTw+fzCwkItLS1CSEJCAv18y8vLKykpKSYmhj5LSUmpoqJCTU3t028NQNxs3769\ntrY2MDDwyJEj3377bQddRVpa2s/P748//vj888+PHDkyYMCADroQSDCMwQJgmK+vr5ubm2CxRG1t\nbfoVYXJy8rZt22JiYqytraurq+mBtywWS0VFJS0t7datW/Hx8XSBlZiYSFEURVEBAQHu7u4URaG6\nAolEV1d+fn5ffPHFH3/80dELmtjb2wcHBy9ZsiQhIaFDLwQSCQUWAJN4PF5kZKSHh4cgoqOjU1xc\nXFtb+8MPP3C5XE1NTWVlZWVlZVtbW0IIi8Xat2/fggULlixZEhMT02IN64yMDLwfBEm1ffv2ioqK\nPXv2EEJkZGSmTZvWCXXP4MGD//3vf585c2bz5s1Y/QQ+CAosACYNGjSovLy8Z8+egsiwYcMoipKX\nlz9//jwlJDU1lW7g6ur68uXLrKyst2upEydOLF26tPOyB+gsdHX1008/0e/KCSHLly8PCQnphEvL\ny8sfOXJk2LBhs2fPfv78eSdcESQDCiwAABBrvr6+5eXlwtUVIaR///5SUlIFBQWdk4OLi8uOHTsc\nHByuX7/eOVeErg4FFgAAiC9fX9+ysrKAgADh6orWaQ+xaKNHjz579mxAQMDevXspiuq060IXhQIL\nAADEVDvVFSHE2to6JSWlrq6u0/JRVVX9/fffCSHz58+vrKzstOtCV4QCCwAAxJGvr+/Lly/bqq4I\nISwW64svvjh9+nRnZkUv+O7h4TFnzpy7d+925qWha0GBBQAAYmfHjh0vX74MDAxsq7qiLVmypIOW\ndG/f5MmTIyMjN2zYIFi1DqAFFFgAACBeduzY8eLFi3dWV4QQNTU1XV3dW7dudU5iwnR0dM6fP8/j\n8RYvXlxTU9P5CYCYQ4EFAABi5P2rK9qaNWsOHz7c0Vm1il7w/Ysvvpg1a9bDhw8ZyQHEFrbKAWjP\nmzdvGhsb229TXV1dU1NTUVHRVgMZGRlFRUVRpwYggXbu3FlaWvrzzz+/Z3VFCBk1atQ///xTUVGh\nrq7eobm1Zfbs2UOHDnV1dV23bt0XX3zBSA4ghlBgAbRnwGDDwcbvtcngn8nX2jr0+J+cpw8fiC4p\nAMm0c+fO58+ff1B1RVu8ePHx48eFd0ToZAYGBomJiR4eHikpKT/++GOPHj2YygTEBwosgPZo9dfb\nEPCpbx98XOaLJBkACfbR1RUhxMnJydraeu3atR9xrqjIycmFhIRERETMnj37119/7du3L1OZgJjA\nGCwAAGDYzp07S0pKPq66IoTIy8tPmDAhOTlZ5Il9KBcXl927d8+bNy8pKYnpXIBhKLAAAIBJu3bt\nKikpCQoK+pTnTytWrDh06JAIs/poo0aNunDhwpEjR7Zt28bn85lOBxiDAgvgAzQ2NKy0NK18WSqI\nVJa92ObmuHjcUL9VrjVVr9sJAsDbdu3a9ezZs0+srggh+vr6tbW1T58+FVVin0JFReXUqVPq6urz\n589vZ/oLSDYWNlQCEZo4ceJr/qeOgXj68AEv5z5TE4JaMDI19434zzrR5yJCU+NjH9y9HXYtS61X\nbzoY+M0alZ4azhu8Dny3TkOrr8u33m8Heffu3LuZxtg9AIgrUVVXtLi4uFu3bm3btu3TuxKVlJSU\nb7/99ueffx43bhzTuUBnwyB3ECVpaemthyI/sZP9X68ICwtTUVERSUoixDEcqqXL2bNisSBCUVTG\npQtbjpzoISNjaj3jZOAPLt96vx1U19RiMG0A8bR79+7i4uIDBw6IamT6rFmz/Pz8vLy8xGcS36RJ\nk+Li4pYsWWJtbe3p6cl0OtCpUGCBOFJTU1NVVWU6i5aMTM1bROpqqutqqnUGDiKE9NM3qCh93moQ\nBRZAC7t37y4qKhJhdUUIYbPZdnZ2f/zxh4ODg6j6/HS9e/f+448/du7cuXjx4uDg4HaWxGtoaFBW\nVra2to6PjxcEXV1dw8PDS0pKlJWVlZWVWSwWRVGysrJTp049cuSItrY23ay2ttbLy+vMmTNlZWWj\nR4/ev3//mDFjOvzeoF0osEAc2dnZaWmJRVGy/ceAd7ah/4WgKIrf3Nx+EABIx1RXtK+++mrRokVi\nVWARQqSlpbdt23b27Flra+vQ0NBhw4a12uzevXtsNvvx48eCyP3790+fPt2/f38tLa3r16/369eP\nPvrixYu1a9euXr06NjaWbrlx40Y2m52dnV1TU3PixAknJyd6ZfmMjIzly5c/ffp03759Li4uhBCK\nooYPH/7333/TJ7q6uh49elRwxStXrvzrX/8qLCyUkpLqmB9GN4JB7tCBJH5IuJyCopyCYnHBQ0LI\ns0f5Gn36thV8HxRFmZmZycnJYeYRSKqamppJkyZdvnz5o1dkaF+vXr00NDTu3bsn8p4/nZ2d3cmT\nJz08PKKiolptwOVybWxsnjx5Ihgb7e3tbWVlRT+L4nK5godSmpqa+/fvv3jxIt2SoqgTJ05s375d\nRUWlT58+y5cvnzx5MiGkqqpq/vz5R48ePXHixMaNG+lz09LSNDQ0qv5LeOplSkqKk5OTsbExqiuR\nQIEFHeVcRKi389zSokLhYLjfds7goaHX7sjKK5w+GNhOsEtgsVhjp9lkXr1EUdTt1OSxlrZtBd8H\n/csoRVElJSXC8YaGBjMzs9zcXEGksrLSycmpb9++48eP//XXXwkh9DBhYcLtARjX3NwcGBg4c+bM\nxsbG33//veP+CV+5cg+iXnkAACAASURBVGVISEgHdf6JOBxOfHz8pUuXli9f3tDQ0OJoZmamhYWF\nvLx8ZWUlIYTL5RYWFvbr1+/tAosQoqioKFyh9u7de+XKlRkZGRRFKSkp0X8tREVF2dnZjRw50sLC\nory8nP7N7fDhw9XV1To6OsOHD09NTZWVlaV7SEtL8/f3t7e3NzdvORYCPg4KLOgoHMOh81etE47Q\no7/H28ykR3//lXyxrWAX4rp5a95trvsUk9flZQ4rPdsJtq+xsXHLli1BQUG6urr5+fl0sKmp6bvv\nvjM1Nc3NzR08eDAd5PP5s2fPnj59elFR0Y8//kgP13B3dxf8Srp3794VK1YYGhp2wO0CfIzU1NTP\nP/9cXl4+Pj6+R48eGhoaHXctc3Pzu3fvVlVVddwlPoWcnNzhw4fNzc1tbGwePXokfIguoTgcDv0e\n8Pvvv9+1a1dmZqagwDIx+f9tux49emRgYEDXWCwW6/z589LS0paWlgMGDBCscRoVFWVra0sIKSkp\n0dXVZbPZNTU1PB7Py8uroKDAxcVl2bJldMv09HR/f/9Tp05lZWWhwBIVjMGCjvLRQ8IHDP2s87Nt\nS1lJceyRX4Qjzhu8kmOjhSMjLaaMtJhCCEk4GdFq8HV52TsvFBYWZmZmZmJiwuFw8vPzLSwsCCGN\njY3m5uZNTU304Ay65eXLl2VlZZcsWUIIsbCwoFvKysrSv4levHgxOTk5Li6OwT1DAAQePXq0fv16\nXV3dkydP9uzZMyUlxczMrKMvunDhwsjIyOXLl3f0hT6ai4vLiBEjFi5cuG3btunTpxNC6uvr//77\n79GjR9MFVmVlZVNT05QpU7KyskxMTKqqqvLy8oQLrIsXL44fP17wVV9f/8SJE7W1tbt373Zzc3vy\n5AmLxXrw4MHw4cMJIXfv3jUwMCCEKCgo3Lhxgz5l3bp1/v7+fD6fy+X6+/tHRkayWKzs7GysKCEq\nKLCgs3WtIeG/Hw9/+0l+CykpKRUVFXPmzGn1aHJycm9VJUdHx7ZO37Rpk6GhYUBAwNWrVwkhenp6\ngidY8vLys2bNSklJmTBhgqB9eHi4np6ekZFRSUmJv7+/m5ub4FBOTo6Pj8/58+fFZ5o6dFtv3rzZ\nuXPngwcP9uzZM2TIEDqYlJRkaWnZ0ZdeuHDhgAEDXF1dBe+/xNDIkSMvXLjw1VdfJScn7969Ozs7\nm8PhqKqq0gXWqVOn9u3bd//+/b59+2poaKSkpHA4HMGTv8rKyoMHD9JPqvh8voyMzJs3b+Tk5OTl\n5T08PMLCwuh4YWEhPVsoISGB/isoNjbW2Nh44MCBhBAejzd69Gg2m+3l5ZWUlBQTE0N3rqSkVFFR\noaamxsiPRZKgwILOIxj9bThqzNtDwoWD4mPSpEnvbFNZWVlSUmJlZdXq0ZKSEk1NzcWLF7d6lBCi\nrKy8e/fuvLy8Pn360BF6so/AjRs3vL29BV9v3bplbGycnJycmJi4cuXKxYsX0w+3KioqHB0dz5w5\nIyZrtEI7JHtCPkVRx48fDwsL27hxo5+fn/ChtLS077//vqMTUFFR0dPTc3Z2jo6Ofndr5tALvv/8\n88+ff/65paUl/f8jh8M5dOiQvr7+hAkTwsLCWgzAamxszMjI8PDwWLVqFV0nsdlsXV3d4ODgpUuX\nVlVVrV+/3sPDg/6VVUVFJS0trWfPnvHx8f7+/oSQU6dOnTlz5vDhw2VlZR4eHrt37yaEJCYm0vkE\nBgbev3//8OFP3d4eaBiDBZ1HtEPCuwo2m62goKDetvLy8rCwsNevX1MURVHUsWPHBE+wCCENDQ2Z\nmZmCh/Z8Pj8/Pz8gIEBTU9PGxkZWVlbwKnD79u2LFi0SDNUCcdb+hPysrKx+/fo1NTU1Nzc/fvxY\nWVl59erVgpYbN25sbGzMzs7m8Xhz5sxxcnKi4xkZGaNGjdLU1IyI+P9X1W/Ph6AoysjISDAZQvgJ\nqEjcuHHDxsamoqIiKSnJzs5O+NCrV69kZGTk5OREe8VWOTs719TU0M9yxBmLxfL09Ny4ceOxY8fo\nN4B6enr379/fsWMHEaqruFxudHS0tLS0srLyd99998033whmBRJCjh8/fvz4cW1tbSsrKzMzs2++\n+Ybued++fQsWLFiyZElMTIyysjIhZMeOHTweT0dHZ968edu2bZs6dapwMhkZGXg/KEoUgOhMmTLl\nTG6x8P8IIWHXsgRfw65lGZma99TqM85y+m/cvFaD5p/PLikpYfpWPkB0dHRQUFBbR0+cOHH48OF2\nTl+9evXWrVsFXy9evKitrS34evPmzWHDhgm+8vl8dXX1uLi4+vr6lStXfvvtt3Q8JydnyJAhDQ0N\nH30X0JkOHz48e/ZsFRUVPp9PR+bNmzdnzpy5c+dSFBUYGDhv3jxB46KiIgUFBboln89XVVUtKyuj\nD1VVVbm5uVEU9fr16/79+9++ffvChQuampr00ebm5okTJ4aFhTU3N6emptKdp6amTpo0STAloq6u\nTlQ39ezZM1dX1+XLlxcXF7faIDY21s/PT1SXa9/du3eXLl1qZWWVmZnZOVcEaAGvCEGU6uvrE6N+\nE46s8P3hr8v/MzHQYuZ/xipd+/efrQZLnjwm3QaPx4uMjOTxeIKIjo5OcXFxbW2tvLw8IeT69evC\nA7BYLNYvv/zi7u4uLS29cOFCX19fOr5x40ZPT08Mveoq6An5N2/erKysVFdXpyfkjxs3jn4P+D4T\n8jds2DB27Ni3J+TX1NTQE/LZbHar8yEEs/R79uwZHBw8Y8aMT7+duro6Pz+/9PT0Xbt2CQ/EbiEp\nKUnkD8zaYmRk9Pfff0dHRzs4OJw/fx4jiqDzocACUdq9e/eLFy/ab1NUVHT+/HnB9OAWHj9+/ERX\np1evXh2QnTgaNGhQeXm5cIR+XiX4um7duhanLFy4cOHChS2CwkN5QPxxudwFCxbQw5nV1dXpCfk+\nPj6zZ8+mjzo7Owsavz0h38fHx9LSUkNDIzQ0lB78FxUVtWbNGiI0IZ+0Nh9CMEt/8uTJgYGBy5Yt\nKyoq+pQbKS4u9vHxuX///tKlS318fARzXVt16dKlwMBOWuiOxWINGjSorq5u165dbm5uMTExmFcL\nnQwFFojSlClT3tkmJyfnn3/+aWsvizt37jx9+hTrCINk4HA4bw9Fr6ury8rKmjFjhp2dnWBC/vff\nf8/lcj///HOKovh8fnp6uo2NDd3+4ybkk9bmQ7Q6S7/9quhtVVVVV65cSUpK4nK5bDa7qKgoKSmJ\nHnDdjrS0tPLy8s78o21lZUXfeGpq6r59+4QHLQF0AhRYAB+Jz+cbGRkpKChUVFQ0NjYeO3as1Wb0\n+5p2lpZWVFSkF2gAyfP23nCEkKCgIBkZGT6f36tXL3pCvp+fn5WVVf/+/Z88eXL16tUFCxb8+OOP\nPj4+hJDKyspffvnFyMiob9++HA7nr7/+qq6urq+v9/HxiYmJKSoq2rZt29atWwsLC1+/fj1z5szL\nly8rKyvHxcXNmjUrPz8/OTmZng/BZrP79esXFBQ0atSoFrP03+dGqqqqEhISkpKSMjIyDAwMrKys\nPD0979+//8svv9y7d09BQeGdPfj7+0+bNu0TfpYfzNra2sPDY+XKld7e3nPmzBk3btz7TAoGEBUU\nWAAfiaIoaTWN74+2vq3Y+/NxmS+SfEAMqaioqKioLF++XLC3bmNj4/79+6dMmcLj8eTk5EJCQvT1\n9ZWUlBQVFenHVJmZmSNGjHj48KFgQj6bzZ43b15cXNz169enT58eHBx88eJFY2PjUaNGOTo67tmz\nx8fHR1lZefr06V5eXnfu3NmxY4erq+uLFy/k5eVv3rw5ffr0ZcuWNTQ0jB8/PioqKjY2tsUs/bfV\n1dUVFBTk5+fn5+enpKSkpaXp6+ubmprOnTv3p59+osup+Pj44ODgM2fOvE91VVZW9ujRo0WLFonu\nR/tu2traz549a25ulpKSCg8PnzlzZkxMTN++4rUQDEgwFFgAAB0lIyNDeCg6ISQsLExVVdXa2rqx\nsZHFYt2/fz8yMjI9PV1dXZ1+h3jjxo2LFy+yWCxlZeVx48bZ2tpmZGQIBqpfvHhxzZo1d+7cuXfv\n3qpVqxobG8eMGSMlJWVvbx8ZGRkaGhoTE6Orq7tmzRrBfAg+n6+srLxw4UI9PT1HR0cXFxcdHZ1B\ngwbt2rVLMEu/vLz8wIEDBQUFT58+bW5ulpWV1dHRKS8v/+eff8aMGXPhwgX6/aMAXV2dPn36faor\nQsixY8eGDx/+zteIIjd27Fgul2tqatqzZ8/AwEA3N7f4+HhpafzDB50B62ABAHSUFnvDvXnzhl6y\n/+uvv9bT06OXaRgxYgSXyy0oKNi0aROLxTp9+vTw4cObmprq6upSUlKePHlCD1Tv1avX0aNHzc3N\nb9++vX79+qKiIi8vr23btu3cuZMQMmbMmKFDh8bHx/fq1WvJkiVDhw5ls9kLFy6MjY2dMGHCnTt3\nsrOzzc3NBw8enJ6eXllZyeVy6R1aaLKyspMnT966dWtCQsKOHTtUVFRKSkrc3Nzu3LlD10bCNxUf\nH3/w4MH3r674fH5MTIyqquqAAQNE96N9L1ZWVoIf/rhx4z7//POtW7d2cg7QbaHAAhCZxoaGlZam\nlS9LBZHKshfb3BwXjxvqt8q1pup1O0GQSKWlpYsWLaKXqiKE7N+/n16yn8VihYWFFRQU0M24XG5C\nQgK9ds5ff/3F4/H4fD596NatW9XV1cnJyT///PO6dev4fH5wcHBWVlZubu6LFy9++ukneg3SRYsW\n9enTx9DQ0NnZWU1NzcjIiBCSkZFB7zHHZrPb32NOUVFx7NixFy9etLW1PXv27I4dO+Lj4+3s7N4e\nkx4fH3/o0KH3r64IIYmJiZMnT378+LG+vv4H/vw+1eTJk4UHOHp4eBQWFv7xxx+dnAZ0TyiwAETj\nXESot/Pc0qJC4WC433bO4KGh1+7IyiucPhjYThAkEr03XHNzMyGktLS01SX76+rqsrOzR40aRZ8y\nZsyY2trapqYm0sbC/QcOHAgJCTE0NOzVq5e9vX1tbS0h5NWrV+fOnauurk5JSbl169aCBQsIIf/X\n3p3HRVXvfxz/DouAIiIIIiggWuaWC67ANXft55L2MNSrEnANA0zK8FZXwcRskTL8Ua4gaoqmhWX6\nc0tJERcCwzTE6wSigoiyiCzKMuf3x7l37lwUND0whK/nXzOf+czpe+Zh+e6c7/l+Fy5cGBcX16JF\ni+bNm5eXl5ubmxcVFd0/yLS0tDlz5sgLYslrgda2H0BsbOzatWt37twpL9L2iNatW/f6668XFxfL\ni4k3pJYtW2o0mpKSEm1l1apV4eHhupslAPWEgAUow6lL1ykB/7VmlSRJSYcPDBo9zrhZs4GjXvw5\n/mBtRTRV165dCwgIkPeGCwsL8/Hx0YYMe3t7+a/5X3/91c7OztbW9v6vq1QqeaJ6RUVFaGioj4+P\nSqXKz89PSEjQaDSZmZm9e/fOzs6WH2j94Ycf7t279+6773bs2FFe4uHQoUNymIuIiPDz85MkSXe9\nzaqqqp07d06cODEyMtLPz+/o0aN+fn7m5ua1ncvChQv9/f23b9/+h9LV5cuXq6urO3TooK9lqIYM\nGZKQkKB9a25uvm7dulmzZsnBFKg/BCxAGT0GuvcbOkq3cres9G5ZqYNLZyFE+07PFObdqK2Ipkq7\nN5y8ZP+8efO0H2mX7E9OTu7bt2+NL2qXFZUnqstPGsoL969fv37ZsmU2NjajRo2ytbU1MDDIzc2N\niIh47bXXXFxcysrKvvrqKw8Pj/T0dO3R9u3bt2/fPt1tCqdNm2ZsbOzp6fnDDz+sWbOmX79+uv01\nFBYWjh8/fsuWLefPn2/RosUf+gWioqL8/Pxyc3O1e5k3sFGjRmk3M5Z179597ty5b7/9tl7Gg6cH\nAQuoX/LflJIkaaqr6y6i6UlPT3/rrbcMDQ3lJfutrKy0H8lL9puZmQUEBHz//fe635IkycTERH49\nffr069evX716dfny5XJx4sSJly5dys3Nbdas2bZt2xwdHTMyMry9vfPy8rKysuRHFNPT07W3+e7c\nuZOWlrZ79+6tW7dqF9sMCgrKz8+XtyP85JNPXn/99S5dujzwFM6ePTtx4sTbt2/HxMR06NDhD51+\nRUXF0aNHX3zxxczMzIZ/hFDWv3//n3/+uUZx+vTpkiRt2rRJL0PCU4KnVYH6Ytq8hWnzFjmZv3fp\n0+/65Qxru3a1FYE/Kjo62s3NzdXV1cnJKSMjQ95nsLKy0t3dvaqq6vz589oVRB+4TaF2g8uDBw/G\nx8fv3r37gbfw1q1b9/XXXw8dOtTIyOgxlgmNi4ubMGGCSqXKyMho+EcIZUZGRm3atLl+/XqNFbAi\nIiLGjBnTr1+/7t2762VgaPK4ggXUF5VK1X/46JSjhyVJ+iUhvv+IsbUVgT9EXu5h2bJlQghnZ2ft\nlG0zM7Px48erVCrdDcJ37NgxduxY8d/bFMouXLgQGhoaGxt7/zbhZWVlr776akZGRkhISGpqakhI\nyGOMMyYm5m9/+5sQIiMjQ19XsIQQI0eOPHz4cI2iiYnJ5s2bX3vtteJinuRFvSBgAfXI+93FF39J\n9hvqWlyQ/4p/UB1F4NHVttyD7OTJk7p7Fz5wm0IhRGFhoaen5+bNm1u3bl3j+Gq1+n/+538mT54c\nFBT03nvvbdq06Y/uVyiEOHfunJ2dnbW1tRAiMzNTX1ewhBCjRo3Sroaly9HRceHChXPmzGn4IeFp\nwC1C4PHdvnXrxP4fdCtvR6xNSz6tWxkz3Ut+8cvxnx5YvFtaIoBHpl3uQX4gcdOmTVFRUdpPKyoq\nUlJStEteaTSaq1evtm3bVgixf/9+T09PbeeSJUtmzZp1/4oM33zzTURExIYNGzp16jRu3LiVK1fq\nTh17dGvWrPH395dfX7161dHR8TEOoohnn3320qVLkiTdfxt03LhxiYmJkZGRb7zxhl7GhiaMgAU8\nJkNDw7fn+peXF9bdlpqaWlFRUccaj4vfCVZ6aGjKalvuQZaamtqxY0ftcgwqlcrCwiIxMdHKymrP\nnj3h4eFyPT09/cCBA7/++qvukSsrK4ODg8vKyn788UdTU9NFixa9+OKLdfzRrUNxcXFaWpr2QlpV\nVZV+N6h57rnnLly40K1bt/s/+uCDDyZMmNC3b193d/eGHxiaMAIW8Ph0n7qvTWxsbElJiZ+fXwOM\nB02evNyDWq3WVrTLPcjLU504cUJ3ApZKpfrss8+mTZvWvn37uLg4bSwLDg4OCgrSnXpVXFw8ffr0\nUaNGBQUFqVSq77//PisrS96H5zFs2bJlxowZ8ut79+41a9bs8Y6jFHmxhgcGLAMDg5iYmIkTJ8ob\nDTX82NBUMQcLAP406ljuQQghSdKOHTu2bNmi3WlHCOHt7X3r1q3U1FTda1F79ux5/fXXdY+sUqk+\n/PDDN998U6VSZWZmLl++fPXq1Y83SEmSYmNjtQErKyvL2dn58Q6llJEjRx45cqS2T21tbT///PNZ\ns2ZVs2wKlEPAAhoLSZLc3NxMTU11/3YEHt2uXbuEEJIk5ebmypXi4uLZs2e3bdvW3t5eXqdUdvLk\nSTc3NysrKxcXl927dwshNBrNhg0bnJycLC0thw0btnbt2jpWda9bQkJCv379tAu+63GNBq02bdoU\nFhZWVlbW1jB48ODhw4fLD2YCiiBgAY3F/X87CiGSkpL69Omjuwa3JEk9evRQ/ZuPj49cLyoqmjp1\nart27QYNGrRhw4aGHz/0q7KyctGiRZGRkfLSo3Jx5syZ1dXV58+fP378+AcffCBn9+vXr48fPz40\nNPTWrVshISHyY3QzZsywtbVNSUnx9va+du3aA++mPaI1a9boXh7T4yqjugYNGnTq1Kk6GoKDg8+d\nO3fgwIEGGxKaNgIW0Cg88G/HO3fuTJkyJSYmRncN7sTERGtr6zv/tmbNGiGERqOZOHHimDFjsrOz\nP/300z179ujtTKAnNZYeFUJUVVX5+vp++eWX1dXVe/fu7d+/v7zaQkJCgre399ixYw0MDHr37m1o\naCiESExM7N69u6GhYatWrQYMGPAY6zLIbty4UVRU9Nxzz2krjSRgjRw58oGLNWipVKr169eHhIRk\nZWU12KjQhDHJHWgUHrgw9wPX4F67dm1paamDg4OVldWqVatefPFFIcSRI0dMTEx8fX2FEB4eHvLX\n8fSQlx49evSo0Fl61MjIaNKkSWfOnHF1dRVC7N+/X2729PSU12tQq9WBgYHyo4W+vr6TJ0+WGw4f\nPixJkru7+5kzZ8rKymoLWw/siYqKkhcX1WoMtwiFEEOGDPnkk0/q7rG0tIyKivLx8dm/f7+xsfFD\nfwGgDvyhAerXjRs3IiMjR9Xu9OnTtS3Mff8a3GVlZWq1euHChZmZmV5eXrNnz5Y7N23a5Ozs3KNH\njzZt2sTExOjrZKEvdSw92qdPn8LCwhUrVmj/tAghCgoK5s6d6+XltWbNmunTp69atSo1NTU9Pf3m\nzZuff/55YGDgA29Y15i5pdvz008/2dvbV1RU7N27183Nzd7e3t7e3tLScvHixfn5+fLTeXKPPJG8\ntjvd9cfU1NTY2Pj27dt1tz3//PNTpkwJDw9/lF9ALt5/H19L95Tx1JGAhpWWlubv71/bp6mpqUFB\nQQ05nsZgyZIluv9Wenl5yfWOHTtmZmZKkrRr167Ro0fX+FZBQYGZmVl1dbUkSd26dZs+fXpeXt7W\nrVstLCzkIp4SN27ccHR0LC4ult9u3LjRw8NDkqS//vWvZWVlcvH333/v2LGj/HrPnj1/+ctf4uLi\nNBqNXOnatatarZZfZ2ZmOjo6du3aNTk5uXPnzgkJCXI9JyfHyspq37591dXVGzZssLOz0/Z88cUX\ntra2Y8eOPXPmzLBhwywsLMzMzG7cuHHhwgVjY+MXX3xRkqSjR4/KPfLREhIShgwZor3Tfffu3Qb4\noT766KPvvvvuUTrv3LnzKL+AJEnFxcUdOnT45ZdfDhw4YGNjo3uQGqeMpw1XsAA90y7MLf87uXHj\nRvkK1gPX4N61a5f2+pZare7bt6+BgYFGo8nIyIiIiLCxsRk9erSJickDN+5FU1Xb0qMHDx7ctm1b\ndXV1dnZ2YGCgvGxbfHz8+++/HxcXN2rUqNLS0tLSUiFEfn5+QkKCRqO5fPmyn5+fm5tbjelc4r6Z\nW3fv3pV7WrVqFRUVNWnSJHd39+eee27+/Pn37t177rnnjI2Nd+7c2a9fv//7v/9LTEwMDw+Xe+Sj\nae909+zZMyEhwcTEpAF+KHk1rEfp3LJly0N/AXnumvY+voeHh3wfX26+/5Tx1NFvvsNTiCtYNQQG\nBi5evFj79uDBg/b29pIkaTQaKyurQ4cOpaSkODg4yAnM09NzxowZJSUlWVlZgwYNOnLkiNzZunXr\n3bt337t3z9/f/+9//7ueTgV6cOnSpdatW+fn52srv/32mxCirKxs586dzs7OrVu37tWrV3R0tHy9\nSr7prCVf6/r+++87d+5sZWXl6uoaExPTpUuX3NxcSZL+9re/6f7h1P4TBwwY0K5du9zc3JMnTzo5\nOf3jH/8YMGCA/KdRkqR+/fppj3/48OGTJ0++9NJLJSUl2p7S0tJBgwbFxcXl5+eHhobKf+AbQHV1\ndZ8+fR7adufOnYf+AoMHD46NjZUkafTo0bt375b++xrh/aeMpxABCw2NgKWrjr8dJUmKiYmxtrbu\n1avX6dOn5U8vXrw4cODAVq1aubq67t+/X/ut2NhYOzu79u3bL1iwoGHutqCpqu2GtSRJ+fn5gYGB\ngwcPDggI0O2ZMWOGmZlZSUmJJElffvmlpaVldHS0PJ3L2dn55ZdfLi0tLSsr0/bo0r3T3QDatGmT\nnJxcd8+j/AJnz56VK/ffxz99+nTdp4ynBAELDY2ABTRatU3nknRmbsmPXBQXF48cObLGLZHCwsLn\nnnvO1NS0sLBQkqTMzExTU9P7e+Li4n7//Xf5sElJSe7u7g12gq+88oqnp2cdDY/yC2jnrlVXVxsZ\nGcn/OzRnzpyoqChJkh74s9TjKaGxYg4WAOBfapvOpTtzKzQ0dObMmS1btjx06JAkSQcPHmzVqpWf\nn58kSZaWljdu3LCysrKwsJCnc3366afyXzYRERHanu3bt4eGhpaWll65cmXevHlLly5tsBMcNGhQ\nWlpaHc8SPsovoJ27pt1L+8yZM3v27JEnSso/S41TbpCTQ+PCOlgAACHq3El6+fLlycnJNjY2cn3Q\noEG6Pbdv3+7du7f89uWXX965c6eNjU3Hjh2Dg4OnTp0q15OSkoYPHy6/Xrp0qZeXl4ODQ+fOnZct\nWzZs2LCGOD0hJElat27dpUuXNm/e/MYbb9zf8Oi/gIeHR0JCQm17act0TxlPIZUkSfoeA54uFy5c\niIyMXLVq1QM/PXv2bExMTERERAOPCkCTFxcXFxYWdu7cOVdX19OnT2sftq2oqBg6dOiGDRu0C9AX\nFRXNmTPn2LFjTk5Ofn5+vr6+kZGR8mOYWhcuXNBdsB6ogVuEAICmRrpv63TtblRmZmYeHh7yXcKq\nqqr33ntv4MCB6enpzz77rNz5wI2n/Pz8tKt2ffLJJ6+//nqXLl30dXb4UyBgAQCamvvXYY+Ojn7m\nmWfmzZtXVlZWWloqz4uqqKiIjY1NTU0tLCw0NDSUF5Q/cuSIoaHhgQMHHBwcgoODx48fL4QwMTEx\nNzc3Nzc/ceJEfHz8//7v/7LaHOpGwAIANCn3b51eUlKyYsWK5OTkmJiYnj17btmyRe48c+aMs7Nz\nUFDQwoULtVunb9y4MT09/fjx45WVlcOHD9fdOv3ChQuhoaGxsbHGxsZ6OTX8iTDJHQDQpNy/dfqK\nFSsuXbokhOjTp48QQqVS6W6d/sUXX7Rp08bd3V3eOv348eNCiNTU1EOHDvn7+xcWFsqHLSws9PT0\n/Pbbb1u3bq2/j9ntCAAAIABJREFUk8OfBlewAABNx/1bp8u7UQ0fPlxecj08PNzExES7dfo777xj\nbGz86quvypthazSaa9eujRgxYtiwYW+88YYkSdpbgUuWLJk1a5Z2qhZQNwIWAKDpWLFixcWLF+3s\n7FQqVXR0dGZmpry0VWZmZs+ePYUQpaWl8uPzzZs3l7f6cXFxeffddwsLCzUajRynrl69euDAgb59\n+967d09uTk9PP3DgwFtvvaXfs8OfCAELANBE3L91elpaWmxs7Ny5c7Vbp6elpd27d6+8vFzeOv3E\niRODBw/Wbp0uSZKBgUF6evqgQYOee+45CwsLOXIFBwcHBQUx9QqPjjlYAIAm4v512HNycgoKCiRJ\nkpdct7KySkxMLC4uNjMz2759+7fffrt27dr8/PypU6d++OGHQgiVSmVubh4VFTVmzJg333zT19dX\nDli6U92BR8FCo2g4o0ePvlF69wkPcisn+5/nf23RooUiQwLQZKjV6gEDBqjVaisrK7mSlpbWvXt3\nedPljRs3BgcHt2/fft26dQMGDBBC/POf//Ty8kpPT5cXlB8zZoz8rW3bts2fP9/IyGj69OlLly41\nMTHR2ynhz4yAhYYzbNiwN9bEPuFBlr/xt31fb+UpHgBAY8YcLAAAAIURsAAAABRGwILeVFZU+I8Y\nWHQrT1spyr/5vo/nqwO6fhzgXXanuI4iAACNGQEL+rF3c1TIzMl52Vd1i5s+XuL0bNeo42dNzJp/\ns3plHUUAABozAhb0w6lL1ykBb+pWJElKOnxg0Ohxxs2aDRz14s/xB2srNlWSJLm5uZmammo0Gn2P\nBQDwRAhY0I8eA937DR2lW7lbVnq3rNTBpbMQon2nZwrzbtRWbKp27dolhJAkKTc3V7deUVHh5uaW\nnp5eo/+nn36yt7evrq6Wv9WjRw/Vv/n4+DTYsAEA92OhUTQu8pp+kiRpqqvrLjYxlZWVixYt+uqr\nr6ZNm5aRkWFvby+EqKqqCgkJ2b9/f1ZWVo0d0I4dOzZ16tS+ffsaGhoKIRITE62tre/cuSN/ynrT\nAKBfXMFCY2HavIVp8xY5mb8LIa5fzrC2a1dbsUmKjo52c3NzdXV1cnLKyMiQi5WVle7u7iNHjhw4\ncKCBwX/+bU1MTAwPD580aZK7u7tcWbt2bWlpqYODQ8+ePRMSElgaEQD0iytYaCxUKlX/4aNTjh5+\ntrfrLwnx/UeMfWAx90qmvkf6x0RGRq5fv17eBO2BPvjgg+7du0dERBw9elQI4ezsrA1YZmZm48eP\nP3bs2ODBg7X9p06dCg8P3759+/Dhw6dNmyaEKCsrU6vVCxcufOGFF1auXDl79uzs7Ox6Pi0AQF0I\nWGhEvN9d/PnbgX5DXTv36OX9zuIHFiPfe7PugzQ21tbWc+fO9fPzq6MnLCzs4sWLdnZ28lsvLy/d\nT0+ePBkSEiK/TkpKCg8Pj42NValU586dk3f8aN68+cmTJ+WGN998Mzw8XKPR6F7xAgA0MLbKQcMZ\nPHiw+7QnnXy9e8Oa0z8d/hNtlRMbG1tSUlJHwMrLy+vfv//58+flHWo3bdoUFRWVkJAgf1pRUWFp\naZmTk2NpaSmEGDVq1I8//qj79cLCwvj4+F69erm4uAghfv7557feeuv48eP1eEoAgIfhChYazjvv\nvJOXl1d3z7Vr106cOOHp6Vlbw9/fCLCwsFB6aPoUFhbm4+MjpyshhL29vfYWoRAiNTW1Y8eOcroS\nQhw6dEh+sXLlyrS0tLVr1wohtm/f/u23365duzY/P3/evHkffvhhw54BAKAmAhYazqRJkx7ac/bs\n2eLi4rpvqDUlarU6NjZWrVZrKw4ODjk5OeXl5WZmZkKIEydO6E7A0kpKSho+fLj8eunSpV5eXg4O\nDp07d162bNmwYcMaZvAAgNoQsAB96ty5c0FBgW6lW7duujfu33zzwXPOtm7dqn397LPPnjp1qp5G\nCAB4DAQs4PFlZ2ffu3ev7p68vLyysjLdu341tGjRoo5nDAEAf0YELOAxaTSaXv369x0y4lGafzq7\nuLaP0s/8fOVSzVXaAQB/agQs4DFJkmTfsdPrYcuf8DihXlMUGQ8AoPFgpRwAAACFEbAAxVRWVPiP\nGFh06z9LURTl33zfx/PVAV0/DvAuu1NcRxEA0JQQsABl7N0cFTJzcl72Vd3ipo+XOD3bNer4WROz\n5t+sXllHEQDQlBCwAGU4dek6JeC/llSQJCnp8IFBo8cZN2s2cNSLP8cfrK0IAGhimOSOp0V6evrq\n1avt7e2VOqBGo9F922Oge42Gu2Wld8tKHVw6CyHad3qmMO/GA4utbVijAQCaGgIWnhaXL1++ffv2\nhAkTlDqgRqPZ9sP/PbRNpVIJISRJ0lRX110EADQZBCw8LZo1a+bi4jJy5EilDlhdXS3eX1pHg2nz\nFqbNW+Rk/t6lT7/rlzOs7drVVgQANDHMwQLqi0ql6j98dMrRw5Ik/ZIQ33/E2NqKAIAmhoAF1CPv\ndxdf/CXZb6hrcUH+K/5BdRQBAE0JtwgBJX2bnqP71rKN7ZJNO2v0PLAIAGhKCFjA47txJWvLZ8ue\n8CAlRYWKDAYA0HgQsIDHZGho+OO+vRUVFXW37d+/v6ys7OWXX66twfztN5QeGgBAzwhYwOPr2bPn\nQ3suXrxYUlLi6uraAOMBADQSTHIHAABQGAELaCwkSXJzczM1Na2xRjwA4E+HgAU0Frt27RJCSJKU\nm5urLZ48edLNzc3KysrFxWX37t26/T/99JO9vX11dbUQ4vbt20FBQU5OTpaWlosXL27gkQMAaiBg\nAY1CZWXlokWLIiMjHR0dMzIy5OL169fHjx8fGhp669atkJCQOXPmaPuPHTs2derUXr16GRoaCiFm\nzJhha2ubkpJy6tSpjz76iGtgAKBfTHIH6pckSWVlZYWFta7F0KpVKwMDg+joaDc3N1dXVycnp4yM\nDA8PDyFEQkKCt7f32LFjhRC9e/eWs5QQIjExMTw8fNKkSR06dNBWZs+ebWhouHPnzn79+hkY8P9O\nAKBPBCygfmk0mu+///7EiRO1NSxYsKBr164RERFHjx4VQjg7O2uvYHl6enp6egoh1Gp1YGBgeHi4\nEOLUqVPh4eHbt28fPnz4tGnT5E5fX9/JkyfLrw8fPlyvZwQAeCgCFlC/Zs2aNWvWrLp7wsLCLl68\naGdnJ7/18vLSflRQUBAaGnrmzJk1a9Y8//zzSUlJ4eHhsbGxKpXq3LlzAwYMEEKsWrUqNTU1PT3d\n2tp6y5YtgYGBFy5cqL8zAgA8FPcRAD3Ly8uLjo4uLi6WJEmSpI0bN2qvYO3du3fSpEkjRoxITEx8\n/vnnhRALFy6Mi4tr0aJF8+bNy8vLzc3Ni4qKvvjii3Xr1nXp0qVNmzaTJk0qLy/X6wkBALiCBehb\nWFiYj49Py5Yt5bf29vZywIqPj3///ff37dtnampaWlqqUqlatGhx6NAhuW3lypVpaWlr164VQuTn\n5yckJHTs2PHKlSt+fn4LFizQ17kAAGQELECf1Gp1bGysWq3WVhwcHHJycsrLy5cvX56cnGxjYyPX\nPTw8EhIStG1JSUnDhw+XX69fv/7tt99+++23O3bsGBwcPHXq1IY8BQDA/QhYgD517ty5oKBAt9Kt\nWzdJkoQQ+/btq+OLW7du1b6eOHHixIkT62mEAIDHwBwsAAAAhRGwAAAAFEbAAgAAUBgBCwAAQGEE\nLAAAAIURsAAAABRGwAIAAFAYAQsAAEBhBCwAAACFEbAAAAAURsACAABQGAELAABAYQQsAAAAhRGw\nAAAAFGak7wEAChg7duytW7fq7rlz505JScnu3btrazA2Nj558qTSQwMAPI0IWGgKrhXeCdv6wxMe\nJNRriiKDAQCAW4QAAAAKI2ABjYUkSW5ubqamphqNRt9jAQA8EQIWmqDKigr/EQOLbuVpK0X5N9/3\n8Xx1QNePA7zL7hTXUdSjXbt2CSEkScrNzZUrxcXFs2fPbtu2rb29fVhYmLbz5MmTbm5uVlZWLi4u\n8qyy27dvBwUFOTk5WVpaLl68WC/jBwBoEbDQ1OzdHBUyc3Je9lXd4qaPlzg92zXq+FkTs+bfrF5Z\nR1FfKisrFy1aFBkZ6ejomJGRIRdnzpxZXV19/vz548ePf/DBB/KVrevXr48fPz40NPTWrVshISFz\n5swRQsyYMcPW1jYlJeXUqVMfffQR18AAQL+Y5I7GRZKkvLy8lJSUxz6CU5eubR2dPnr9Vd1jJh0+\nsGj9VuNmzQaOenHbyuVefw+5v9japq0SZ1BTUVHRb7/9ZmpqWltDt27dzMzMoqOj3dzcXF1dnZyc\nMjIyPDw8qqqqfH19R48eXVxcvHfv3v79+xsYGAghEhISvL29x44dK4To3bu3oaGhECIxMXH27NmG\nhoY7d+7s16+f3AkA0BcCFhoXc3PzqqqqnTt3PvYRegx0r1G5W1Z6t6zUwaWzEKJ9p2cK8248sFhP\nASsmJubQoUPPP/98bQ2zZs1ycnKKiIg4evSoEMLZ2Vm+gmVkZDRp0qQzZ864uroKIfbv3y/3e3p6\nenp6CiHUanVgYGB4eLgQwtfXd/LkyXLD4cOH6+NEAACPjoCFxqVz5847duz4o9/ac1+oup9KpRJC\nSJKkqa6uu6istm3bTpo0yc/Pr46esLCwixcv2tnZyW+9vLy0H/Xp06ewsDAmJmb27NlXr/7rvmdB\nQUFoaOiZM2fWrFnz/PPPr1q1KjU1NT093draesuWLYGBgRcuXKin0wEAPAruI6DpM23ewrR5i5zM\n34UQ1y9nWNu1q62oF3l5edHR0cXFxZIkSZK0ceNG+QrWjBkzysvLVSqVpaXlSy+9ZGxsLPfv3bt3\n0qRJI0aMSExMlC+MffHFF+vWrevSpUubNm0mTZpUXl6ur3MBAMgIWGj6VCpV/+GjU44eliTpl4T4\n/iPG1lbUi7CwMB8fn5YtW8pv7e3t5YB18ODBbdu2VVdXZ2dnBwYGzps3TwgRHx///vvvx8XFjRo1\nqrS0tLS0VAiRn5+fkJCg0WguX77s5+e3YMECfZ0LAEBGwMJTwfvdxRd/SfYb6lpckP+Kf1AdxQam\nVqtjY2Pl8CRzcHDIyckpLy9fvXr10qVLbWxsxo0b98orrwQFBQkhli9fnpycbGNj07Jly5YtW8pT\n3devX79s2TIbG5spU6b4+voGBATo5VwAAFoqSZL0PQbgSXXu2Xvm/H884UE2LV+SeeE3RcajKzY2\ntqSkpO45WACAJoZJ7mgK/L1n3bx4pu6erKysa9euubvXOh3ez2vmH/3nfvXVVw+d8JSUlHTv3r06\nGtq0afPyyy//0X80AKAx4woWnhZHjhw5fvx4aGioUgesrq526drjlcC3nvA4363/Un3+rCJDAgA0\nElzBAh5fqzZt3MZOeMKD7I/dpMhgAACNB5PcAQAAFEbAAhTzJ91kGgCgOAIWoIw/6SbTAID6QMAC\nlOHUpeuUgDd1K/J+0oNGj5P3k/45/mBtRQBAE8Mkd+Dhdu3aFRYW1qZNG91ijSdwG9sm0wAAPSJg\nAQ9naGg4a9as+fPn6xarq6v7/GXoQ7+rr02mAQB6xC1CoL408k2mAQD1h4AF1JdGvsk0AKD+ELCA\netRoN5kGANQr5mABSvo2PUf3rWUb2yWbdtboeWARANCUELCAx6RSqW5mX1viO/UJj1Nx764i4wEA\nNB4ELOAxGRgYXM/KfGhbbGxsSUmJn59fAwwJANBIMAcLAABAYQQsAAAAhRGwgMZCkiQ3NzdTU1ON\nRqPvsQAAnggBC2gsdu3aJYSQJCk3N1e3XlFR4ebmlp6erq0kJSX16dPHxsZm8+bNciUyMlL133T7\nAQANjIAFNAqVlZWLFi2KjIx0dHTMyMiQi1VVVe+9997AgQPT09OfffZZuXjnzp0pU6bExMRs3bo1\nODhYLvr5+d35t08++eT111/v0qWLfs4EAMBThEB9KysrS0lJ2bmz1oWvRo0aZWlpGR0d7ebm5urq\n6uTklJGR4eHhIYSorKx0d3evqqo6f/68gcG//ndox44dEyZM6N27d1lZWUFBgUajMTAwMDExMTEx\nEUIcPHgwPj5+9+7d8naHAAC9IGAB9at79+43btwoLCysraGwsNDIyCgiIuLo0aNCCGdnZ+0VLDMz\ns/Hjxx87dmzw4MHa/h07dsydO1cIkZub6+joqA1eQogLFy6Ehobu27fP2Ni4vs4HAPAICFhA/Ro8\neLBuPHqgsLCwixcv2tnZyW+9vLx0Pz158mRISIj27aVLl3r27CmE+PXXX5955hltvbCw0NPT89tv\nv23durViowcAPBbmYAF6lpeXFx0dXVxcLEmSJEkbN27UXsESQlRUVKSkpAwYMEB+q9Forl692rZt\nWyHE/v37PT09tZ1LliyZNWuWdqoWAECPuIIF6FlYWJiPj0/Lli3lt/b29roBKzU1tWPHjpaWlvJb\nlUplYWGRmJhoZWW1Z8+e8PBwuZ6enn7gwIFff/21gQcPAHggAhagT2q1OjY2Vq1WaysODg45OTnl\n5eVmZmZCiBMnTujeYVSpVJ999tm0adPat28fFxenjWXBwcFBQUFMvQKARoKABehT586dCwoKdCvd\nunWTJEn79s0336zxFW9vb29v7xrFPXv21M8AAQCPg4CFp91HH310+/btunvUanVhYWFeXl5tDVZW\nVn//+9+VHhoA4M+KgIWn3YZtO/ze/7juHusBw+puWL/kPQIWAECLgIWnnYlZ807dn3/CgzQzNVNk\nMACApoFlGgAAABRGwAL+S2VFhf+IgUW3/jPdqij/5vs+nq8O6PpxgHfZneI6igAAyAhYwH/s3RwV\nMnNyXvZV3eKmj5c4Pds16vhZE7Pm36xeWUcRAAAZAQv4D6cuXacE/NeyCJIkJR0+MGj0OONmzQaO\nevHn+IO1FQEA0GKSO54WpaWl27dvT0hIqKOnx0D3GpW7ZaV3y0odXDoLIdp3eqYw78YDi61t2tbb\nwAEAfz4ELDwtJkyYMGHChPvr94eq+6lUKiGEJEma6uq6iwAACG4RAnUzbd7CtHmLnMzfhRDXL2dY\n27WrrQgAgBYBC6iLSqXqP3x0ytHDkiT9khDff8TY2ooAAGgRsICH8H538cVfkv2GuhYX5L/iH1RH\nEQAAGXOwgJq+Tc/RfWvZxnbJpp01eh5YBABARsDC086updnquV5199y8ebO8vNzR0bG2hvatWyo9\nLgDAn5hKkiR9jwFo7Hbv3q1Wq+fPn6/vgQAA/hyYgwUAAKAwAhbQKJSXl8+fP9/Jycnc3HzIkCHJ\nycnajyRJcnNzMzU11Wg02uKAAQMMDQ2NjIwMDQ0tLS0//fRTfYwaAPBgBCygUQgODq6srDx37pxa\nrX7ppZemTp2q/WjXrl1CCEmScnNz5UpFRcXZs2ezsrKqqqqqqqq2bdsWFham7S8qKpo6dWq7du0G\nDRq0YcMGIcTt27eDgoKcnJwsLS0XL14st/3zn/8cN26cra1tp06ddu/eXWM8P/30k729fTVrqALA\nYyFgAfonSdLWrVuXLFliYWFhZ2c3Z86cF154Qf6osrJy0aJFkZGRjo6OGRkZcvH8+fNt2rRp3769\nEEKlUnXo0MHW1lb+SKPRTJw4ccyYMdnZ2Z9++umePXuEEDNmzLC1tU1JSTl16tRHH32k0WgKCgrG\njh37xhtvXL9+/YMPPvD29ta9PHbs2LGpU6f26tXL0NCwQX8IAGgqCFhAo2Bra+vv75+UlCRJkrm5\nuXzlSQgRHR3t5ubm6urq5OSkDVjJyckDBgyQXxcXF3/++eczZsyQ3x45csTExMTX19fAwMDDwyMu\nLk4IkZiY2L17d0NDw507d/br18/AwCAuLu6FF14YO3asoaHhsGHD7t27px1JYmJieHj4pEmT3N0f\nvokQAOCBWKYBqF+//PLL119/3bp169oa/vrXv3bo0GHfvn2hoaEjRoywtraOiooaOXKkEKKkpCQi\nIuLo0aNCCGdnZ92A9d1338mbIQohevbsuX79evn1pk2bnJ2de/TokZubGx4e7uPjI4Tw9fWdPHmy\n3HD48GEhREVFRWpq6rVr1yoqKubOndu1a1cDAwMhxKlTp8LDw7dv3z58+PBp06bV028CAE0eAQuo\nX0lJSYWFha6urg/81MDAwMzMTAjRqVOnrVu3lpeXf/jhhz4+PleuXFGpVCtWrLh48aKdnZ3c7OX1\nr/W6kpOT9+/fP2bMGPn1kCFDNBqNnJDOnDnTq1ev+Pj4Q4cO+fv7v/rqq2vWrElNTU1PT7e2tt6y\nZUtgYOCFCxdmzZr1ww8/dOnSpVevXs7OzvIdxqSkpPDw8NjYWJVKde7cOe1FMgDAH8U6WMDDPck6\nWLGxsSUlJX5+frU1aDSaZs2alZSUmJqaCiFu3rzZq1ev7Ozsmzdv9u/f//z58y1bthRCbNq0KSoq\nKiEh4e7duy1btszOztbOu1KpVOXl5fJjhi1atMjKyrK1tb1161a3bt1u3LjRvXv3H374oVOnTkKI\ny5cvDx069PLly9euXbO3tzcwMKiqqurRo0dERMTYsWNHjRr1448/6o6tsLDQ0tLyMc4aAJ5yzMEC\n9MzAwMDR0XHVqlW3b9++du1aQEDAvHnzVCpVWFiYj4+PnK6EEPb29vItwl9//dXOzk6brnSpVCoz\nM7PTp09XVFSEhob6+PioVKr8/PyEhASNRnP58mU/P78FCxYIIXr06PHDDz/cu3fv3Xff7dix4+jR\no4UQhw4dkiRJkqSIiAg/Pz9JkkhXAPB4CFiA/n311VdfffWVvb39yJEj3dzcFixYoFarY2Nj582b\np+1xcHDIyckpLy9PTk7u27dvjSPI87FUKtWXX37p5+fXqVMnc3Nzee2G9evXL1u2zMbGZsqUKb6+\nvgEBAUKIiIiI1157zcXFpays7Ouvv5ZvL2olJSVxfxAAngS3CIGHq9dbhACApocrWAAAAArjKULg\nMUmS5O3tXV5eXnfblStXqqqqakwe1+Xg4PD5558rPToAgD4RsIDHpNFoki/8850vY57wOB8HeCsx\nHABAI0LAAh6foZGRuUWrJzyIgQHb0QBAU8McLAAAAIURsADFVFZU+I8YWHQrT1spyr/5vo/nqwO6\nfhzgXXanuI4iAKApIWAByti7OSpk5uS87Ku6xU0fL3F6tmvU8bMmZs2/Wb2yjiIAoCkhYAHKcOrS\ndUrAm7oVSZKSDh8YNHqccbNmA0e9+HP8wdqKAIAmhknuwMPdu3dv8+bN+/bt0y3WWKS3x0D3Gt+6\nW1Z6t6zUwaWzEKJ9p2cK8248sNjapm39jh4A0OAIWMDDvfLKK6+88kqNYnV1dZ+/DH3od+VNbCRJ\n0lRX110EADQZ3CIE6otp8xamzVvkZP4uhLh+OcParl1tRQBAE0PAAuqLSqXqP3x0ytHDkiT9khDf\nf8TY2ooAgCaGgAXUI+93F1/8JdlvqGtxQf4r/kF1FAEATYmqxkRdAI+ourr6+cEeimyVk/bzKUWG\nBABoJJjkDjwmAwODXp077gkPrbtN3uzZxcWltobhA/spPTQAgJ5xBQuoX7GxsSUlJX5+fvoeCACg\n4TAHCwAAQGEELED/KioqTExMxo8fr1v09vZWqVQ3btwoKyszNDQ0MjIyNDRs3rz5uHHjcnJytG3l\n5eXz5893cnIyNzcfMmRIcnJygw8fAFATAQvQv/PnzxsYGGRlZWkraWlp33zzTYcOHdq2bZuamtq+\nffuqqqrq6uqsrKyWLVsGBgZqO4ODgysrK8+dO6dWq1966aWpU6fK9aSkpD59+tjY2GzevFnbXFRU\nNHXq1Hbt2g0aNGjDhg1CCEmSevToofo3Hx+fhjppAGjKmOQO6F9ycvLo0aN/+uknSZLkRd5DQkJG\njhxpYGAgf9qv378mwtvY2KxYseKZZ56ROyVJ2rp1a0ZGhoWFhYWFxZw5c3777TchxJ07d6ZMmbJ7\n9+68vLyZM2d6eXkJITQazcSJE729vbdt23bixIkVK1b4+vomJiZaW1vfuXNHPr6xsbF+fgIAaFq4\nggXoX0pKioeHh5mZWVFRkRAiOTn56tWr7du3l3OVbsASQrRo0UIOYTJbW1t/f/+kpCRJkszNzeXr\nUjt27JgwYULv3r09PDwKCgo0Go0Q4siRIyYmJr6+vgYGBh4eHnFxcUKItWvXlpaWOjg49OzZMyEh\nwcTEpIHPHQCaJK5gAfUrLy9v165dKSkptTW89dZbycnJ06ZNc3JyysrKat269T/+8Y9ly5aFhoZO\nnDhRCJGcnDxz5kxt/+XLl5955hk5Y6lUqn379oWGho4YMcLa2joqKmrkyJFCiB07dsydO1cIkZub\n6+joKF8J27Rpk7Ozc48ePXJzc8PDw318fMrKytRq9cKFC1944YWVK1fOnj07Ozu7vn8QAHgasEwD\nUL8KCgouXLhgamr6wE8NDAw6duxoZ2d348aN1157bcaMGa1atQoLCztw4ICFhcW1a9eaNWtmaWmZ\nl5dnbW0tfyU8PDwjI2P16tW6xykvL//www83btx45coVlUrl4uJy5MgRZ2fn7777bvXq1QcOHBBC\ndO/evVevXitXrjx06JC/v39hYaEcvGSFhYUODg4lJSW6RQDA4+EKFlC/rKys3N3d62hITk52cnJq\n1aqVfAVr+/btn332WVpaWrt27aytrY8dO+bk5KRNV0VFRatXr/7xxx+FEBqNplmzZiUlJaampmZm\nZvPmzYuOjpbrV69ebdu2rRBi//79np6ecjEjIyM+Pt7Gxmb06NEmJiYqlWrXrl29evWSF0FVq9V9\n+/YlXQGAIviPKaBnKSkp8hQrJyendevW2djYDB48WDvvSvuisrIyMTFxxIgRAQEBciQyMDBwdHRc\ntWrV7du3r127FhAQMG/ePPlhQAsLi8TExDNnzuzZs0cOWCqVyszM7PTp0xUVFaGhoT4+PiqVavv2\n7aGhoaWqgNHYAAAC7klEQVSlpVeuXJk3b97SpUv1+ksAQNPBFSxAz5KTk11dXYUQzs7OaWlpsbGx\nQidXJScn79y508jIyMjIaMCAAQsWLJg2bZr2u1999dXcuXNDQkI6dOgwZ86cefPmCSFUKtVnn302\nbdq09u3bx8XFtWzZUi5++eWXfn5+RkZG06dPDwsLE0IsXbrUy8vLwcGhc+fOy5YtGzZsmF5+AQBo\nepiDBQAAoDBuEQIAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAA\nCiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAo\njIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAw\nAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMII\nWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAACiNg\nAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAF\nAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYA\nAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAACiNgAQAAKIyABQAAoDACFgAAgMIIWAAA\nAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACgMAIWAACAwghYAAAACiNgAQAA\nKIyABQAAoDACFgAAgMIIWAAAAAojYAEAACiMgAUAAKAwAhYAAIDCCFgAAAAKI2ABAAAojIAFAACg\nMAIWAACAwghYAAAACiNgAQAAKOz/AelUAMU9rUh/AAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -h 800 -w 800\n",
    "library(ape)\n",
    "tre <- read.tree(\"empirical_6/RAxML_bipartitions.empirical_6\")\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": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 0.04078121\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "mean(cophenetic.phylo(ltre))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Translation to taxon names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Library_Name_s                        Organism_s\n",
      "0          90.o31                 Vidua hypocherina\n",
      "1          90.o38                       Vidua regia\n",
      "2          90.o50                    Vidua fischeri\n",
      "3          90.o53                  Vidua orientalis\n",
      "4            A056   Lagonosticta senegala rhodopsis\n",
      "5            A081                  Vidua paradisaea\n",
      "6            A082                      Vidua obtusa\n",
      "7            A089                    Vidua macroura\n",
      "8            A107          Lagonosticta rhodopareia\n",
      "9            A145              Lagonosticta larvata\n",
      "10           A147                    Vidua macroura\n",
      "11           A167    Lagonosticta senegala rendalli\n",
      "12           A177                  Vidua chalybeata\n",
      "13          A1794             Anomalospiza imberbis\n",
      "14           A204            Lagonosticta rufopicta\n",
      "15           A248                 Lagonosticta rara\n",
      "16           A328     Lagonosticta sanguinodorsalis\n",
      "17        JMD1249                Vidua purpurascens\n",
      "18        JMD1285                  Vidua chalybeata\n",
      "19         MDS004                    Vidua raricola\n",
      "20         MDS054               Clytospiza monteiri\n",
      "21         MDS065    Lagonosticta rubricata congica\n",
      "22         NKK511                  Vidua interjecta\n",
      "23        SAR6894  Lagonosticta rubricata rubricata\n"
     ]
    }
   ],
   "source": [
    "print pd.DataFrame([indata.Library_Name_s, indata.Organism_s]).T"
   ]
  }
 ],
 "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
}