20170803_RBP2b_diversity_for_Wai-Hong.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "This notebook is used to create a figure on RBP2b nucleotide diversity and Tajima's D for Wai-Hong Tham at WEHI, for a manuscript under preparation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import allel\n",
    "import vcf\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_context('paper')\n",
    "sns.set_style('white')\n",
    "sns.set_style('ticks')\n",
    "rcParams = plt.rcParams\n",
    "rcParams['font.size'] = 9\n",
    "rcParams['axes.labelsize'] = 9\n",
    "rcParams['xtick.labelsize'] = 9\n",
    "rcParams['ytick.labelsize'] = 9\n",
    "rcParams['legend.fontsize'] = 9\n",
    "rcParams['axes.linewidth'] = .5\n",
    "rcParams['lines.linewidth'] = .5\n",
    "rcParams['patch.linewidth'] = .5\n",
    "rcParams['font.family'] = 'arial'\n",
    "rcParams['ytick.direction'] = 'out'\n",
    "rcParams['xtick.direction'] = 'out'\n",
    "rcParams['savefig.jpeg_quality'] = 100\n",
    "rcParams['savefig.dpi'] = 80\n",
    "rcParams['legend.frameon'] = True\n",
    "from matplotlib.ticker import MultipleLocator, FormatStrFormatter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "output_dir = '20170803_RBP2b_diversity_for_Wai-Hong'\n",
    "!mkdir -p {output_dir}\n",
    "\n",
    "manuscript_vcf_ftp = 'ftp://ngs.sanger.ac.uk/production/malaria/pvgv/May2016_release/pv_1_0_NG_manuscript.vcf.gz'\n",
    "manuscript_vcf = \"%s/%s\" % (output_dir, os.path.basename(manuscript_vcf_ftp))\n",
    "rbp2b_vcf = \"%s/PvRBP2b.vcf.gz\" % output_dir\n",
    "diploid_vcf = \"%s/PvRBP2b.diploid.vcf.gz\" % output_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(148,)\n"
     ]
    }
   ],
   "source": [
    "# These are the positions (in bp) of the two exons of RBP2b, using Sal1 coordinates\n",
    "RBP2b_exons = [[54438, 54492], [54781, 63146]]\n",
    "    \n",
    "# The following are the 148 samples used for most analyses in the PvGv manuscript (http://www.nature.com/ng/journal/v48/n8/full/ng.3599.html)\n",
    "analysis_samples = np.array(\n",
    "    [\n",
    "        'PD0165-C', 'PD0166-C', 'PD0167-C', 'PD0168-C', 'PD0169-C', 'PD0170-C', 'PD0171-C', 'PD0172-C',\n",
    "        'PD0173-C', 'PD0174-C', 'PD0175-C', 'PD0176-C', 'PD0177-C', 'PD0178-C', 'PD0179-C', 'PD0180-C',\n",
    "        'PD0181-C', 'PD0182-C', 'PD0183-C', 'PD0184-C', 'PD0593-C', 'PD0594-C', 'PD0595-C', 'PD0596-C',\n",
    "        'PD0597-C', 'PD0598-C', 'PD0599-C', 'PD0600-C', 'PD0601-C', 'PD0602-C', 'PD0603-C', 'PD0604-C',\n",
    "        'PD0605-C', 'PD0606-C', 'PD0607-C', 'PD0608-C', 'PD0609-C', 'PD0610-C', 'PD0611-C', 'PD0612-C',\n",
    "        'PD0613-C', 'PD0614-C', 'PD0615-C', 'PD0617-C', 'PD0618-C', 'PD0619-C', 'PD0620-C', 'PD0624-C',\n",
    "        'PD0625-C', 'PD0674-C', 'PD0675-C', 'PD0676-C', 'PD0678-C', 'PD0679-C', 'PD0680-C', 'PD0681-C',\n",
    "        'PD0684-C', 'PD0685-C', 'PD0686-C', 'PD0687-C', 'PD0688-C', 'PD0689-C', 'PD0690-C', 'PD0691-C',\n",
    "        'PD0692-C', 'PD0693-C', 'PD0694-C', 'PD0695-C', 'PD0696-C', 'PD0697-C', 'PD0698-C', 'PD0699-C',\n",
    "        'PD0700-C', 'PD0701-C', 'PD0702-C', 'PD0703-C', 'PD0704-C', 'PD0705-C', 'PD0706-C', 'PD0707-C',\n",
    "        'PD0708-C', 'PD0709-C', 'PD0710-C', 'PD0711-C', 'PD0712-C', 'PD0713-C', 'PD0714-C', 'PD0715-C',\n",
    "        'PH0177-C', 'PH0178-C', 'PH0180-C', 'PH0182-C', 'PH0187-C', 'PH0188-C', 'PH0190-C', 'PH0310-C',\n",
    "        'PH0312-C', 'PH0320-C', 'PH0893-C', 'PH0914-Cx', 'PH1078-C', 'PH1109-C', 'PH1112-C', 'PH1113-C',\n",
    "        'PH1124-C', 'PH1125-C', 'PH1133-C', 'PJ0006-CW', 'PJ0007-C', 'PJ0008-C', 'PJ0009-C', 'PJ0011-C',\n",
    "        'PJ0012-Cx', 'PJ0016-C', 'PJ0018-C', 'PJ0030-C', 'PJ0032-C', 'PJ0033-C', 'PJ0035-C', 'PJ0037-C',\n",
    "        'PJ0041-C', 'PJ0042-C', 'PJ0043-C', 'PJ0047-C', 'PJ0049-C', 'PJ0055-C', 'PJ0057-C', 'PJ0059-C',\n",
    "        'PJ0063-C', 'PJ0065-C', 'PJ0068-C', 'PJ0069-Cx', 'PJ0071-C', 'PJ0072-C', 'PJ0074-C', 'PJ0079-C',\n",
    "        'PJ0083-C', 'PJ0085-C', 'PJ0086-C', 'PJ0088-Cx', 'PJ0090-Cx', 'PJ0097-C', 'PJ0099-C', 'PJ0106-Cx',\n",
    "        'PJ0117-Cx', 'PJ0123-C', 'PJ0127-Cx', 'PJ0167-C'\n",
    "    ]\n",
    ")\n",
    "print(analysis_samples.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Download files from MalariaGEN website"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2017-08-03 09:45:12--  ftp://ngs.sanger.ac.uk/production/malaria/pvgv/May2016_release/pv_1_0_NG_manuscript.vcf.gz\n",
      "           => `/nfs/team112_internal/rp7/data/pv/analysis/20170803_RBP2b_diversity_for_Wai-Hong/pv_1_0_NG_manuscript.vcf.gz'\n",
      "Resolving ngs.sanger.ac.uk (ngs.sanger.ac.uk)... 193.62.203.121\n",
      "Connecting to ngs.sanger.ac.uk (ngs.sanger.ac.uk)|193.62.203.121|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /production/malaria/pvgv/May2016_release ... done.\n",
      "==> SIZE pv_1_0_NG_manuscript.vcf.gz ... 276257852\n",
      "==> PASV ... done.    ==> RETR pv_1_0_NG_manuscript.vcf.gz ... done.\n",
      "Length: 276257852 (263M) (unauthoritative)\n",
      "\n",
      "100%[======================================>] 276,257,852  321M/s   in 0.8s    \n",
      "\n",
      "2017-08-03 09:45:13 (321 MB/s) - `/nfs/team112_internal/rp7/data/pv/analysis/20170803_RBP2b_diversity_for_Wai-Hong/pv_1_0_NG_manuscript.vcf.gz' saved [276257852]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O {manuscript_vcf} {manuscript_vcf_ftp}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2017-08-03 09:45:14--  ftp://ngs.sanger.ac.uk/production/malaria/pvgv/May2016_release/pv_1_0_NG_manuscript.vcf.gz.tbi\r\n",
      "           => `/nfs/team112_internal/rp7/data/pv/analysis/20170803_RBP2b_diversity_for_Wai-Hong/pv_1_0_NG_manuscript.vcf.gz.tbi'\r\n",
      "Resolving ngs.sanger.ac.uk (ngs.sanger.ac.uk)... 193.62.203.121\r\n",
      "Connecting to ngs.sanger.ac.uk (ngs.sanger.ac.uk)|193.62.203.121|:21... connected.\r\n",
      "Logging in as anonymous ... Logged in!\r\n",
      "==> SYST ... done.    ==> PWD ... done.\r\n",
      "==> TYPE I ... done.  ==> CWD (1) /production/malaria/pvgv/May2016_release ... done.\r\n",
      "==> SIZE pv_1_0_NG_manuscript.vcf.gz.tbi ... 15050\r\n",
      "==> PASV ... done.    ==> RETR pv_1_0_NG_manuscript.vcf.gz.tbi ... done.\r\n",
      "Length: 15050 (15K) (unauthoritative)\r\n",
      "\r\n",
      "\r",
      " 0% [                                       ] 0           --.-K/s              \r",
      "100%[======================================>] 15,050      --.-K/s   in 0.001s  \r\n",
      "\r\n",
      "2017-08-03 09:45:14 (24.9 MB/s) - `/nfs/team112_internal/rp7/data/pv/analysis/20170803_RBP2b_diversity_for_Wai-Hong/pv_1_0_NG_manuscript.vcf.gz.tbi' saved [15050]\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "!wget -O {manuscript_vcf + '.tbi'} {manuscript_vcf_ftp + '.tbi'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create VCF for RBP2b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!bcftools view --regions Pv_Sal1_chr08:54438-54492,Pv_Sal1_chr08:54781-63146 --output-file {rbp2b_vcf} --output-type z {manuscript_vcf}\n",
    "!bcftools index -f --tbi {rbp2b_vcf}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create truly diploid version of vcf\n",
    "!zcat {rbp2b_vcf} | \\\n",
    "sed 's/\\t0:/\\t0\\/0:/g' | \\\n",
    "sed 's/\\t1:/\\t1\\/1:/g' | \\\n",
    "bgzip -c > \\\n",
    "{diploid_vcf}\n",
    "!bcftools index -f --tbi {diploid_vcf}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read VCF data into numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vcf_data = allel.read_vcf(diploid_vcf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([False,  True], dtype=bool), array([ 80, 148]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vcf_samples_in_148 = np.in1d(vcf_data['samples'], analysis_samples)\n",
    "np.unique(vcf_samples_in_148, return_counts=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(159, 228, 2)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# How many variants and samples, and what is ploidy?\n",
    "vcf_data['calldata/GT'].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert genomic coordinate to position within gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4275\n",
      "1\n",
      "55\n",
      "nan\n",
      "56\n",
      "8421\n"
     ]
    }
   ],
   "source": [
    "def position_in_gene(pos=59000, exons=RBP2b_exons):\n",
    "    offset = 0\n",
    "    for exon in exons:\n",
    "        if pos >= exon[0] and pos <= exon[1]:\n",
    "            return(pos-exon[0]+offset+1)\n",
    "        else:\n",
    "            offset += (exon[1] - exon[0] + 1)\n",
    "    return(np.nan)\n",
    "print(position_in_gene())\n",
    "print(position_in_gene(54438))\n",
    "print(position_in_gene(54492))\n",
    "print(position_in_gene(54493))\n",
    "print(position_in_gene(54781))\n",
    "print(position_in_gene(63146))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_tajd(vcf_data=vcf_data, samples_to_use=vcf_samples_in_148,\n",
    "              size=300, start=54438, stop=63146 , step=3, despine_bottom=False,\n",
    "              highlight_regions=[[169, 470]], ax=None):\n",
    "    pos_in_gene = [position_in_gene(x) for x in vcf_data['variants/POS']]\n",
    "    ac = allel.GenotypeArray(vcf_data['calldata/GT'][: ,samples_to_use, :]).count_alleles()\n",
    "    D, windows, counts = allel.stats.diversity.windowed_tajima_d(\n",
    "        pos=pos_in_gene, ac=ac, size=size, start=position_in_gene(start), stop=position_in_gene(stop), step=step,\n",
    "        fill=0.0\n",
    "    )\n",
    "    D[np.isnan(D)] = 0.0\n",
    "    \n",
    "    # set up axes\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(figsize=(8, 2))\n",
    "#         fig.tight_layout()\n",
    "    \n",
    "    x = (windows[:, 0] + windows[:, 1]) / 6\n",
    "    \n",
    "    # plot\n",
    "    sns.despine(ax=ax, offset=5, bottom=despine_bottom)\n",
    "    ax.plot(x, D)\n",
    "    if despine_bottom:\n",
    "        ax.set_xlabel('')        \n",
    "        ax.set_xticks([])        \n",
    "    else:\n",
    "        ax.set_xlabel('PvRBP2b amino acid position')\n",
    "    ax.set_xlim((position_in_gene(start)/3, position_in_gene(stop)/3))\n",
    "    for region in highlight_regions:\n",
    "        ax.axvspan(region[1], region[0], facecolor='k', alpha=0.1)\n",
    "    ax.axhline(0)\n",
    "    ax.set_ylabel('Tajima\\'s $D$')\n",
    "    \n",
    "    minorLocator = MultipleLocator(100)\n",
    "    ax.xaxis.set_minor_locator(minorLocator)\n",
    "\n",
    "    return(ax, x, D)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nfs/users/nfs_r/rp7/miniconda3/envs/biipy/lib/python3.5/site-packages/allel/stats/diversity.py:1040: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  wD = d / d_stdev\n",
      "/nfs/users/nfs_r/rp7/miniconda3/envs/biipy/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['arial'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAC5CAYAAADKxC6cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAIABJREFUeJzt3Xl4U2XaP/DvyUnSpvu+70Apa0Eoi2yyC6VUUUZxpYPL\nOAI6jj8FnUGH4R1GR718B5fBrQI6KqBQC6KvoAioiCyyUwq0tKV039KmSc72+yM0NE3aJjRplt6f\n6/K6zDkn59w5pMmd57mf52EkSZJACCGEEOJhZM4OgBBCCCHEESjJIYQQQohHoiSHEEIIIR6JkhxC\nCCGEeCRKcgghhBDikSjJIYQQQohHoiSHEEIIIR6JkhxCCCGEeCRKcgghhBDikSjJIYQQQohHoiSH\nEEIIIR6JkhxCCCGEeCS5swNwlA0bNuDBBx90dhgAgPLycmeH4BAxMTHODoEQQgjplMe25HhqYkEI\nIYQQ63hskkMIIYSQvo2SHEIIIYR4JEpyCCGEEOKRKMnpQ4quNuNoYZ2zwyCEEEJ6BSU5fcjhgjoU\nlKidHQYhhBDSKyjJ6UMkCWBZxtlhEEIIIb2CkhxCCCGEeCRKcvoSasQhhBDSh1CS05dIzg6AEEII\n6T2U5BBCCCHEI1GS05dQdxUhhJA+hJKcPqKxWY8Gtd7ZYRBCCCG9hpKcPqKkSoNJw8OdHQYhhBDS\na1wuydHpdLjzzjsxf/58ZGZm4t///rezQ3J71Q1anLjYAD+VAhFB3tj8fYmzQyKEEEIcTu7sADpS\nKpXYsGEDfH19wXEc7rnnHkyePBkjRoxwdmhu69eCOmSkhSAuXIX4CB9KcgghhPQJLteSwzAMfH19\nAQA8z4PneTAMVcz2hChKiA33oftICCGkT3G5lhwAEAQBCxYsQElJCe655x6kp6dbPC43Nxe5ubkW\n9yUkJDgyRLciiBJYGSU4hBBC+haXTHJYlkVeXh6amprw+OOP4/z580hNTTU7LicnBzk5ORbPsXbt\nWkeH6TZESnIIIYT0QS7XXdVeQEAAxo4di/379zs7FLcmCJTkEEII6XtcLsmpq6tDU1MTAECr1eKn\nn35CSkqKk6NybyIAGSU5hBBC+hiX666qqqrCihUrIAgCJEnCrbfeiqlTpzo7LEIIIYS4GZdLctLS\n0rB9+3Znh+FZOizMybIMtuwtgYxhkHVzLJQKl2vQI4QQQnrM5ZIc4ngLJsUBAPYcrUSThkNYoJeT\nIyKEEELsj5IcD/fbhXpU1mtNtrXNl6OUy8ALojPCIoQQQhyO+ik8XEFJEx6e18/iPjkrA8dLFvcR\nQggh7o6SHA/HyBh4K1mL+xRyhlpyCCGEeCxKcvowBSsDJ1BLDiGEEM9ESU4fFuinwP4TVfj422LU\nNumcHQ4hhBBiV1R43IclR/thSbQfjl+oR02DDqEBNMqKEEKI56CWHA+m1nAQxe67o5QKGfQ81eYQ\nQgjxLJTkeLAPvrqE4SlB3R6nlMug5yjJIYQQ4lkoyfFQHC8iPsIHg5MCuz1WIaeWHEIIIZ6HkhwP\nJEkSvj1SAX8fhVXHhwZ64cSlBhSUNjk4MkIIIaT3uFySc/XqVdx///2YM2cOMjMzsWHDBmeH5HZE\nUUJjM4dpIyOsOt7XW45Zo6PQ2Mw5ODJCCCGk99g0uqqiogKNjY2Ij4+Hj4+PQwJiWRYrVqzAkCFD\n0NzcjDvuuAMTJkxA//79HXI9T8QJEny8WbCs9TmsyotFq15wYFSEEEJI77IqySkrK8OyZctQU1MD\nLy8v1NTUYNy4cXj22WeRnJxs14AiIiIQEWFogfDz80NKSgoqKyspybEBz4tQ2JDgAECgrwJFV5vR\n1MIh6+ZYB0VGCCGE9B6rvglfeeUV3HXXXdi/fz92796Nw4cPY+rUqXj44YdRXFzssODKyspw9uxZ\npKenO+wanogTJMhZxqbnqLzkWHxrCrXmEEII8RhWteQUFxfj9ddfv/4kuRx33XUXYmNj8eabb+Jf\n//qX3QNraWnB8uXL8dxzz8HPz8/iMbm5ucjNzbW4LyEhwe4xuQteECG3sSXHiFZ5IIQQ4iGsSnIY\nxnKrwMSJE/Haa6/ZNSAA4DgOy5cvR1ZWFmbNmtXpcTk5OcjJybG4b+3atXaPy9UVXW1Gk4ZHg1oP\nudy2lhxCCCHE01j1c7+6uhpbtmzB8ePH0dLSYrKvswToRkmShOeffx4pKSmdJjDEsr2/VcFLIUNk\niDdGp4Y4OxxCCCHEqaxqyVm2bBnOnj2LvLw8FBYWwtfXFwMGDMCAAQNQU1Nj14COHDmCvLw8pKam\nIjs7GwDw1FNPYcqUKXa9jifyVcmRlhDQo3O06gRU1WsREextp6gIIYQQ57AqybnrrrtMHldUVKCg\noAAFBQXIyMiwa0CjR49GQUGBXc/ZZ9ihnia9fxDqm/WU5BBCCHF7N7QKeVRUFKKioqh1xQPJWQaC\nQNXHhBBC3J/LzXhMnIuVycBTkkMIIcQDUJJDTMhZBrxIi3USQghxfzYlObt27UJzczMA4K233sLS\npUtx+vRphwRGboAdBrpRdxUhhBBPYVOS89Zbb8HPzw+HDx/GgQMHcNttt+HFF190UGjEGVgZQ91V\nhBBCPIJNSQ7LsgCAH374AYsWLcKMGTPAcbRytSs4fqEegb6KHp9HzspQWq3B2cuNOF+mhiRRwkMI\nIcQ92ZTkREZGYtWqVdi1axemTJkCvV4Pkeo3XEKzlsfgxMAenycyxBuDEwPACxL2n6iCRkdrWRFC\nCHFPNg0hf/3117F//378/ve/R0BAAKqqqvDMM884KjZiA0GUwNq4KKclrIzB0OQgAEBptQaCSC05\nhBBC3JNNSY5KpTJZSyoiIgIRERF2D4rYThQNCYo9sTIqQiaEEOK+rEpypk2bBoZhEBISgi1btjg6\nJnIDBFG0e5IjZxlqySGEEOK2rEpyvvvuO0fHQXrIMS05MkpyCCGEuC2bl3VobGzE5cuXodPpjNvs\nvX4VsZ0gSJA5oruKkhxCCCFuyqYkZ8uWLdi4cSMqKiqQlpaG48ePY8SIEdi4caPdAlq5ciX27t2L\n0NBQ7Nixw27n9XSCJFFNDiGEENKOTUPIN27ciK1btyImJgabNm3Ctm3bEBISYteAFixYgPfee8+u\n5+wLBMEBSQ7V5BBCCHFjNrXkKJVKeHl5AQD0ej369euHoqIiuwaUkZGBsrIyu56zLxAlCaydVyJT\nyBkcPFODwjI1bkoNRmiAl30vQAghpM/asPMMOF6EQi7DfbemgbX3lxhsTHKioqLQ1NSEGTNmICcn\nBwEBAU4dQp6bm4vc3FyL+xISEno5Guf48OtL8PGSQ5QkyO38BhnRLxjJUX64UqNBYZkaoYMpySGE\nEGIfel7Aw9nDHHoNm5KcN998EwCwbNkyjB07Fmq1GpMmTXJIYNbIyclBTk6OxX1r167t5Wicw8db\njt/d4piETiZjEOyvRIuWR1ML75BrEEIIIY5i8+iqNmPGjLFnHMSFyWQMRKrNIYQQYoMDx6/gbHEd\nGBjqRRkGCPS73iPA845fFsqmJOfkyZP4z3/+g/LycvD89V/2+fn5dg+MuA6WYSDQQp2EEEJscLao\nDg9lDwXDGJKcR/6xGw/OG4zByYYBSyqvG25nsZpNV3j66afxzDPPIDU1FTKZ/QuEAOCpp57CoUOH\nUF9fj8mTJ2PZsmVYuHChQ65FrCOTgVpyCCGEWI3jRVyuaDImOIZtAkL8vRHs791rcdiU5ISEhGD6\n9OmOigUA8Nprrzn0/MR2NzIpYFW9Bhotj4hgFXy8FWb7q+tbUVLZBABIig5AaKDKLrESQghxPp2e\nx9B+YSbbOEGEQuGYBpLO2JTkLF++HM8//zzGjx8PpVJp3N5+0U7ieW6kJufDHWcwrF8ofjyuw723\nppntzz9wCekDwsDzIr76qRj3zxlkr3AJIYQ4GSeYr6eo5wzDxXuTTUnO559/jkuXLoHneZPuKkpy\nPJshybHtOSEB3pg1Lgkf7jhtcb8oShiVFgleEHGmqM4OURJCCHEVgmA+rQnHC1DK2V6Nw6Ykp6Cg\ngIqM+yCWYVBS1YLfCusBBhiSFGhVNs7KGIjdFCyzMga8rRkUIYQQl8YLolmSkxgd0CvFxu3Z1G6U\nnp6OCxcuOCoWj/bujov4Yl+ps8O4ITIZgyPn6+GjkuNyZQuu1rZ2+xwJhuSmoUkHqYtEh2EY4/BC\nQgghnsGQ5Jh+tr/+p1sQ5N+7k8ralFIdOXIE27ZtQ1xcnElNDrXudC/QTwHeTiOU9JwAXpDg7dU7\nzX5tyXhqnD/UGg6cDXMbyOUy6DgBDMNg/RcnMGtsItKS7LveGSGEENciCJJDlmmwlU1Jzvvvvw9J\nkkyGhJHe98GuIvh4sRjRP7hXrieXy/BoVj8AgIKVQc93n6y1tc5Eh/lCz4ngeAFeChbnLtdTkkMI\nIR7OUneVM1iV5CxatAiffPIJMjMzTRKctoTn6NGjDguQmAsJUGJUagga1PpeuZ6MYXBTqiExUcgZ\n8EL3LTlt3VVKOQs9J0DPCwgLUqGovBGnL9VCrbkeuyCK+CD/tMkILkmScPesgfD3UZqdmxBCiGuz\n1F3lDFYlOZ988gkA4NixYw4Npq/QcyKUNswV8J8vL4AB8Oj8/sZtrAzYe7wKyVG+Doiwcwq5DK06\nzurjvRQy6DkBOk5AoJ8Sk6Jj0disw6yxicZjgvy9kN4/3KSF58fj5SitVGNwcqhd4yeEEOJ4vDt2\nV5EbI4oSZADa2j+W/fsI1v85w+rnh/gr0Vab23YuGcOgpLIFT9yRau9wuyRnZeAE62uLlAoWxwur\nwfEiggO8MXpQpNkxd80YaLYtMtQHn/5fASJDfGiiQEIIcTOCKELhLklOW3fVyJEjwTCMyWgZhmEQ\nFBSEJUuW4N5773VYoO5MECWwLAPRyuTg60NX0aThoNZwmH5TJOQsA/7ac7V6AV5KFizLgGGYXu/z\nZG2cGHD8sGicL2kAAKQlWV9D1D8uCFNHx6OyTkNJDiGEuBmeN3zvOZtduqvq6+tx9913U5LTCUGU\nIJMxgJVJTlMLh99NTcCZ4kYcPV+PsEAvVDXoIAgitHoB3krWsGimDS0q9iKTwaYlHvx8lLgpLeKG\nrsXKnPMaSd+26+diXKlqNj4elBSCCekxzguIEDfEiyK8enniP0ts7q5qbGzE5cuXodPpjNsyMjKw\nadMmuwW1b98+/M///A9EUcTChQvxyCOP2O3cziCIktn01tZQymVobOEwNDkQDc0cfi2og0ZrSHJk\nVky05wgs031LTlfz4thCzsog0ESBpJddqWrGQ9lDjY/fzTtpkuQIooRXPz4ChVwGPx8FhqaEYvww\nSoIIaU8QRLC9PPGfJTZFsGXLFmzcuBEVFRVIS0vD8ePHMWLECGzcuBERETf2a70jQRCwevVq5Obm\nIjIyEnfeeSemTZuG/v37d/9kFyW2teRY+91/LR9SKmRQa3goFSymjYxARb0WCATiwn3Q3MrbvGim\nPRiSq66P+Wz3efipej4qSnYDC4MSYm8Xyxrx+qdH4adSoqWVw5L5QxAV6oMH5g6GKErI3XGakhxC\nOuB5CfJeXqfKEpsi2LhxI7Zu3YqYmBhs2rQJ27ZtQ0iIfec8OXHiBBITExEfHw+lUonMzEzs2bPH\nrtfoTVV1Gmz4pgj9YvwAxrZWDoVchnq1Ht4Kwy/G/rH+6B/rD2+loQnQGb2d1izW2azhsGiWeTGx\nrdp3V9m6QGhn6tVaXK1p6fa/xmZd9ycjfUJ1vQYzMhLwUPZQ+KoU4HgRPt4KAIa/Byc0qBLi8nhR\nhFzm/CTHppYcpVIJLy/DlMx6vR79+vVDUVGRXQOqrKxEVFSU8XFkZCROnDhh8djc3Fzk5uZa3JeQ\nkGBx+94jpbhQ1tjzQK3U2KIDL0g4frEBtY06vPTfs1B5sVizyfLClQDg48Vi894SiKKEID8F9h6v\nMrQEtdOq5SFKwOa9JY5+CSZ0egGVdVrUNOng52t5Yc2TF2vwXt6pHl+rtrEVek7AiQs12HHgEuZN\nTOnxOX89U4GwQO/ur92kszgSjHi+Ux3ev43Neuz6uRgHT1Xg1MUaaLQcrta2oEGts3g8IQQoqWhC\noL8XAn3tt4xD+25ka9mU5ERFRaGpqQkzZsxATk4OAgIC7NZN1cZSS0dnMyzn5OQgJyfH4r61a9da\n3H7LqHjcMir+xgO00YWyBpw4V4Jxg8MAAJu/L0FKrB9+d4vlJMxaDc16VNbrenweWzU263HwTC1m\nj4lGTIx5E70kSXjvy1M39Gbs6FxxHWoaWzExPRZ5+y7a5ZwSJDycPazb497Ls89rIO6n47/97l9L\ncOe0AUiOCURu/mlMz4jHqUu1mHtzssXjCSHAR7vOYnpGAqLDencut46sSnJ4nodcLsebb74JAFi2\nbBnGjh0LtVqNSZMm2TWgqKgoVFRUGB9XVlbaPZHqqW9/uYwLZQ1QdFE5rucFpMYHIyHK/4aKjrsT\n5KfEM4sG2f283ZHJGAhdtM/zgv3mRpDJrg+dJ8QRTl6swa6fihEScL11LyU2wOQYjhOgVBj+1hUK\nGTRaHsp2tQa8IOKTb85h0ey03gmaEDeg0fFOT3AAK5OchQsXYtu2bSbbxowZ45CAhg0bhuLiYpSW\nliIyMhI7d+7Eq6++6pBr3aiahlY8MHcwfFWKTo9pbNbh64PFiA33M+1qcv60AT1iqSbn8NlKnLpY\nAwDQ8yJUdqqol7MyiDS6ijhQs4bDnPFJGNY/rNNj9LwIxbWkRilncaGsAX7tlhv5w4Lh1F1FSAeu\nssSlVd9G9hoSbA25XI5Vq1bhoYcegiAIuOOOOzBgwIBeu741OEHstmpcIZeB5yXwogi2/b+2BLdO\ndAyTAZpuO3a+CjnzhgAAWlo5u01QyMoYqDUcSivVdjkfIR1JkmRW79bRbVP6IeBaUjN5ZCzOFtdh\nxIDw3giPENJDViU5dXV1nRb4Aui0LuZGTZkyBVOmTLHrOe2J47vvkpGzMvCCCFHo/kPUncgYWJyf\npy2xCfSzX5FZbIQfcnecxtmiOqQPCLs+FJ+4LEmSUFKhRoCfEsH+3Rd4O5skdf+Lc8n86/U2UaG+\niAp1fhM8IcQ6ViU5oiiipaXF0bG4DdGKX3/stSSnbUmHNlq9YGz6dke9OXeNnJXh+ZwxEAQJm/ec\nx96jZfD1lmPs0OgbPifjzs1oLuzy1SaovOUQRQlbvy+El4LF0oUjnB1WtxyROFfUtkDPCcbHgX5e\ndk3+CXEHOr3Q/UG9wKokJzw8HEuXLnV0LB6FvTYjsSCKkLX7qXjz0DDjPDfuqH1NTm1jK7R6ARzn\nuLoZhZyFQg7cNqU/rtY0I39/EcYMiep0xB1xjk27zmJ4/zCMHBiBIcmhKGu3LIIrEyXJ5O/THt7+\n4gRubjc54MWyBvzxznS7XoMQV8cLrlFP6XI1Oa5Kq+dR32SYF6NVy1v1HAbMtSUdrm/rH+vviPB6\njYy5PvnZW1tPYPSgCIwdGtX1k+wgwFeJAN8Q/HzyKgRRgtzCwm8cL3Q54s0WktXTUxPAsNSBXG5o\nvZSzjNvcP3snOZIkIT7CH7PHJRq3vZvXZLfzE+IOBFGCv0/PZ723B6v6TT788EMHh+H68vdfwo8n\nynH0XCXSrSw6rGlsRWFpg+fVkVx7OVFhPphzczJGpfXepHlyVgaeN/+FcOlKIx76n2/R0sp1+Xx3\n+fJ1N4IggpW1JTnu0x0rSZJdR4GIEjzv750QG+n0PLxcpMfCqpacoKAgR8fh8nhBQtbEJJOho925\nfUo/tGh5hKr0DozMCZyYJ7CsDLyFmiCNlkNEsA90nNDl0H5iX7+dr8KRc1W4XNGE8cNjIAiGFp34\nCH+TYdVV9Ro8t9gx0070hL1rckRRNJsXy5o6MI2Ww+5D12cvZ1kZbh2XCNaNEsa+prJOg6YWHQbE\nBzs7FJfTquOhUjp/cU7gBlYh76tu5MNwYKJhXa/y8nJHhORUomj/WgZryFnGYkuOnhfhq1KYFHxa\nQoXH9vXr2UosyRqKzAnJOHSmAty1Fp1bxyeZHOeq88iIEuzyPm5rIRQEySzJ0Wg5HC+sxvD+YcZa\nsuZWDlodjyB/L8hZGUoq1GBlDDIGG7p+d/5YBLWGQ5A/FSy7qq9/LkZdkxasrBj3zRlkMqGkrTRa\nzmT+MW8vuVu1iHb02/lqRIT4ODsMAJTkWE0QRWqGbocXROMssL1JzsogWJggkLcyySH2xYCBTMYY\n54USBNGtRg9aM0+OLQw1eKbny57SDx9/fQ4XShtwxzTDnF9vbz2O8GAV4iL8MGNMIrR6HpGhvsYv\nBn9fJdVCujiOF+Hno8D3h0uRNSmlR0nO2g9/xZB+oQBg7PK9e2bPFzl2ljNFdfjjHcOdHQYAG1ch\n78tEUQLrAiuquoo9RyvtNrOxLViWAWexJcfQTaW3sI84TlsLhkLO4mJZA2obtRaLwgHXHMAgivat\nyRFECbIOrz8xKgC+3gp8uPOMcVtQgBeyJqWgqcXQld2qE0xGXcoYxuJ8VMT1yGQMhB4uP5MYHYC7\nZw7E3TMH4p5Zad3WFro6by/WZbpaqSXHSgJNRGeiXs1h8fx+vX5dBStDbaPWrLWgvkmH0EBvbNt7\nwWwSOkmScPst/REWpLL6Or3VrdWq4/HypsP46+/HuvX7K8BXieljElBaqTZ2ubSnuDbyyl6j3+zF\n3oXCgiha/DFkKWEJ8FXi1KVa1DXpUN2gwQNzBxv3yWTotfmoyI0TryU3fA+WnxE6JNoyGSW49kRJ\njpUMQ02dHYVr4AQRPl5yp2Tq6anhOHiyAgWX6022MwyQNTEFPt7mRcfHCqpQfLXJpiSnt0ZhtbRy\n+O18NThBhJfMtRKArjSodSitVKNZc/0X500DI3DTQMuL6SrkMug525OcxmadyUR61fWtaG61XMiv\nkMsQF2HbFA32qi1rS4pFC91VwPVWrC17zl/bYGj9WrVknMXztZ+qgfScIEp4+n9/wOpHb7Y4tLmq\nToNtey+AZWVgGNNZrgFg109FuFLdAo2Wwx8WDIdSwUKChKmj4+GlZC225DQ266DW6BEe7AOvLrr2\neQtdvPVNOhw5V4l+sUFUl9VDlOTYgCagM9BzIoL8nNMUGRPmhwVT+9v0HF+VAldrWuw+XNge9LwA\nH285OF7s8oPQ1ez48RISowIw5+Ykq45PjgnAR1+fNbZyiJKElJhAzBiT0Olz1Bo97nvha+S/mm3c\ntmHnGYwcaHkKhx9PlGPlg2NsqgmS7DRPTleFxwDw6O3DsWT+UGPi4uPd9UevpYVwyXV6TsCL7x4E\nK2Pw9z/c3O3xvCCioVmPo+eqkBQTgMQow0rzOk6AgpXh4OmrmDIqDmmJIdiw8wz07VaeB4Cy6mY8\nnD0Mn39XiOZWDiHX9qUmBKO+SWvx3yp3x2nEhvtBqWCRPbnzVm9LywTdNqUfrta04LvDpTZ/3hFT\nLpXk7Nq1C2+88QYuXryILVu2YNiwYc4OiVig0QrwUrhGf6s1AnyV+OT/ClBU3mh1IiFjGGxoV0Nh\nydybkxEebH3rkCUcL8LXWwGOEwA3Gvqu0wu4eXiMxS90S8YPi8H4drMA6zkBn35b0OnxgiCivLrZ\nbHbwIH8vTM+wnBhVN7RCp+ehkFs/zYMoSXb98WKp8BiAzdMaUJdF11p1PNJTw9DUor/246Xrf0OO\nFzFuaBQEUcIn3xRg6cJ0+Pko8df//ISZYxJQcLkec29OBmD4t2pp5SwOrFB5y7HvWBnUGs44KSzL\nylBe04wTF2oQGaLCjDGJxvPMuTkZX+672GVsPG++4HNqQjCiQn2x6+ciq+8JscylkpzU1FSsW7cO\nL7zwgrNDIV0YNzgUcRGuMTzQGlGhvnhucQYAWF08vnjeEIhd9LOfL2nAueI6hAfHdnme/35zDppO\nZsiWIGFSeix8VHKLxdSuTJQsf5lbi2VllgvIOQHVDa3Y9VMxfL3lZhNvdtWN6K2Uo1UnwM+Gt6Yo\nGupfeiohMgDv5Z0CL4iYMDym+yd0Q8ZYbslpbuWQm3/aWPTfquPxwNxBfW5tLK1egEopR6uWt6rW\ni+MFhAR4Y9roeOg5AfVqHfx8lAgPVqGpRW8cyg8AAxOC8dnu8yZrDkZfW5TV11uBH46W4W+PjDce\nL5MxOHWpFrPGJuKnE+WY0e66Ki85ahpaUV7d+TIn9WqdxdZHhVzm0CVzHMmVpupwqSSnX7/eL2Tt\n6P9+uYySCjU0Wg6LZqX1+Je6Jxo1MMTZIdjM1loQVsaA7aJGJiLYB4dOV2DSyK6THI2Wx0PZQy3u\n++L7QtQ2ag0tOS6yzos1OF6EVtezofqszHLNyZ7DpWhQ6zAgPghTborDu3knrT6nykuO7w6XIKhd\n4Xn6gLAuVw23V3dV+2Uc7IHpZHRVXWMrUhOCMHtcEgBgz68lqG3UulySI0kSJDsVdWv1vFnNS71a\nC28vObyUcug6qfUqr25GQ7MOg5JCDF1C147x91Eaa8l8vRUorVKbTFw3rH8YhvUPsxjLhPQYjEgN\nh3e741kZg1YtjwBfpdmPKFbGYFj/MPxWWN3laxwzxLxYXym//kOgqk6Drw8Wm+zvFxuECek9T6g9\nnUslObbKzc1Fbm6uxX0JCZ339Xdl1ljDh9X+366gpqGVkhxiUXiwqscL0Km85LhU3oiwIBU27z5v\nUhApQULmhGTEhPn1NFS7++FoGRKjHLMGm54TMD0jHhHBhuYYQwHu9e6Irn4hTh0Vh6u1LcbHTS16\n/Hi8HHdMG4CvfipCVZ0GqQnBECUJ54oNhetF5Y12T1DsgZXBYkuOVi+YfMF6K+XQ6q1bS6835e44\nA62Ot8vCpK9+fAQpMYFm26eMisPBk1dRXt2MyBAfs0Rv857zkCRDSy7Pi1Bca5nx9mLRqjPcM4VC\nhpvSIhBp5cR1clZmdh1WxqBVx8NLwZrU/LW9V6eOirf6tZqcl5UZR22dKarF6EGR6B93ffWBDTvP\nuGyS40rL5/R6krN48WLU1NSYbX/yyScxY8YMC8/oXE5ODnJycizuW7t27Q3F10Ypl4ETDL9WW3U8\nKmo0PTrA86nUAAAaaUlEQVQf8Tzd1QF0N+2At5cchSX1eC5njMkXFwCcL6lHYUmDSyY5Wj1v8Zen\nrSx9EOo5Acp2v8rlrAxHzlUZ597pqk5FqWCNBaWAoaj0WEEVAKCsqhk584bgw52nIWOYTlvXXEVn\nhcdaPW8yP5WXkoVO73oTYIqiZLfJQiNDfLFodprFfWOGROHouSps3n0ef/n9WOP2F9/9GRPTY6DV\nC9DqecP8TdfeV95KORqaDYstM2Awbmh0j+JjZTJodLzZ67XHF31lrQY/HC3DsfPVuGd2msk1WnU8\nfjhaBm8li7E9fA2erNeTHHdZ7FMhZ6G/1h9aWqnulZW2iWfhOAHKLkb69I8LQkmF2izBAQBvJeuS\nv9ABQ9Gxo0aCcbwIZbui9rk3J6P4aiM43vCFkTUxxepzKeUy6K7NgM0whhoHBu4xNLvjEPJWHY81\nH/yCkko1Xnp8onG7l5LFL6cqMCgpBN5OmJyzMxIkq+syquo1xpY7W8VF+CMuwh9V9a3GbRwvIik6\nADPGJGL3ocvQ6QXk7buIR243DGQJCfBG3r6LOHWhxqzg90b4quTQtHLG5PPrn4sNcdihnuaBuYOg\n1QvIikhBRIdehUWz0lCv1iJ//yWMHhzVoxq5ntr6XSEa1Drj49ZO6hCdwXX+KlyMQiEz9ttq9TyC\nXKzPmzifRsvh5AXzVsn2+7v6NRsf6Y8HMwdb3OftZSiidUU6TnDYCsN6TjCprwgPVt1wl7GlljZb\nvnydqWNNjkbLITrMF2eK6hATfr11b3ByKC6WNeJCWQOG9rNcR+LKBFHCkjXfmkwT0B7HC53OoN0Z\njZYzzpflpZSjur4VMeF+xkQqOswXzy0egw/yTxtGNfZQQlQA1q809ELcPXMg6pq0AICb7VCAntCu\nZbKjtr+NqFBf6PS8xTnCekuDWueyraMuleR8++23+Pvf/466ujo8+uijGDRoEN5//32nxKKUy4wL\nQWp1Ary93GcOE9I75k5IRm1Da5fHTM+4sf74QD8vVNS24L28UxBEEbHhfphnQyuGI+n0gl26Iipr\nNfji+0KTbYWlDTZ/qXmijt1Vek5EgK/SrA6MlTFIjgmARmf9L2etjsfFK43GxwMTgyFnZTh+vhqH\nz1Xintlp3S7ZUlalRvHVJkxMNy+8/+lEORqadAi2Yi0n3bXWyrYC86gQX2RNuv4+12h5qLqZUwgA\n6pq0WLvhEMKDfCAIorE7dVBSCL4/UorJI8zjTIoOQG1j13+/tvJVKWyeLqCnVF4stHrBaUmOKy7X\n0p5LJTkzZ87EzJkznR0GAEMF/pY9hThw4goEQcIDcwc5OyTiYvrHBZkUAtqTl4LFHxYYFrgTRAkf\n7jjd7XNqG1tR33S9yTg8WOWQUTeCKNllheQn7x5pttbYjDGJdp23RqcX8O72k8ZZY0VBcqmiyM50\nnCdHzwnwU1me/8fHW46Cy/XGBSJZGYOk6IBO7+PBU1dR16RDVKgPLl1pBMcLOHWpFjKGudb1ozGp\nbbLkaEEVTl6oMSY5umtdswzD4GhBFR65fRi27Cns8hyA4d/ngbmDsHB6KgBg/RcnTFaslyAhKqTz\n0XFtahpaMXdCMm65Kc5ke1iQynjujqaNvrEfIK7G30eJ/35zDn4qBSJDfTFnfFKvXt+RLbv24FJJ\njiuJCffDX34/Fll/zsOHq2YhNJBGWRHnYK2cGO7jr89heLuhr98fKcXDt7nuhJp+FqbXt7elC0eY\nPH50gWusjNydjvPk6DgB/j6Wf6knRgXgSlUzSirUAAxF60H+XmZruLWprtfgllFxiIvwR0iAN6ob\nWtGq4/Fw9jD8eLzcuDjkzyfLcaaozuI5LpY1IjbierfZs2/sx8CEYIQH+0DlJbc6udZxpqPFOv77\nPPvGfiRHd51wAcCy341AcB9d/mDGmARMvZawfZDf/Y8he9Noefi4UD1YR64bmQuxVBhKSG9iwJj8\nwm3ToNZh0eyBiA33g8pbjlvaDVctLGvozRCJHcWG+2Lnj0U4eKoCgKE25dbxSbjDwhT/SgVr8u8+\neWRsp5PP1TRo8dnuAsydYJjdNzxYhfwDl4xJSVSoD/L2XcSB4+Xw9VYgZ96QTmN8P//6+3FISiju\nmDoADGOYe8ZahiHxnbcCPHH3yE5bsNqLj3TMlAbugGEYYxcvz4tWzQBtrXq1Frt+KgYAjB8WjWQL\nQ/kNNVCu+x3pupG5kK7+CAnpDZ0V9f16pgL7jl1BkL+X2VDiukat2QKXxD0kRAXgsTvM55ix9CXT\nkZyVdVqwGh3mixEDwo2tQqGBKvy/+0Yb9/eLC8JT94yyKkY9J2LngUsAgIoajbG7rL3uvnC1er7L\nrg5XnELBlcVF+uHSlUb0s1M3+oXSBvSLDURidAD2/3bF+P5rbNahqNxQ13WlugV+LrwkDSU53Ugf\nEOaU1bYJsYaPtwK7firCP/44wWw5gekZCTh1qda4vbmVg1bH27Qae5u1Gw4ZnxfoS0mTu1LIWSy/\na6RdzvVg5mCoWwwrwo+3MJIoY1Ak3vvylNloNh0n4PFrkwTqOkxuSHomOtQXao3eZJsoSqjppMDa\nT6XosmC5RcsjJswXoYHexm5MwLAyQGSID3y8FYgM8cGgJNedBZ/eXd1Y84cJzg6BkE5FhvjgprQI\nxEWYN9dHhvjg4Kmrxsf/d/AyCkvr8ewDGTZdQxQlhAf5uOwQUeIcfipFl7/g01PDkZ5qvmJ8+25X\nnd61i1bdTUy4HzbvPo9fTldgaL8wTBgeg6MFVTh8thLRYeYF3GVVzcaEs70TF6rx2/lqFJU34dHb\nh0EhZ1FY2oBNu84CAM4V1+Hvj95sl2U7HI2SHELcWFiQCk/efZPFfdFhvjhWUG1sYr54peGGRkWd\nvlTr0n3uxL2EBXkbh4y3ann8bobl0U/EdrHhfvjTopvA8SI+/vosJgyPQVOLDjPHJFjswupsbbhf\nTlfgwbmDr02iaUhCX3x4HIRrxfAyhnGLBAegJIcQjyVnZTh5sQaZE5Lh4y3H5BGx+HL/JVTVaRDR\nYa2ei2UN0HMiBiWbNzsfLajCXfRFROzktinmxdPEvhRyGS5eacQ720+iqk6DkakRFo9r35W4bvNv\nxlXTvZSs2VxYCjkL16286RwlOYR4sPmTUjB+WLTxV5eXksXliiazJOf7I2UQRPMkp6xKjer6Vpda\nMoAQ0r2//H6sce6irmj1PCTJ8NnwUPZQbNlzHr+eqeylKB2PPrkI8WAd58kJ9vfGJ98WoPhq07UW\nHgWq6jWoqtcg0M/LbJh6Y7MO2VNcY6ZlQoj1rFlfLj7SD5u+OguVtxwjBhjqp2aMScDIgZZbftwR\nJTmE9CGJ0QH44x3pOHS6AkXlTRiSEorickPCkz7AvEiUEOK5Zo9LMtsW7O/d6USS7oiSHEL6mABf\nJfrFBRqWLTl+Bc2tXJeTvhFCiLtyqSTnpZdewvfffw+FQoGEhASsXbsWAQHdT+lNCLFNckwgnrl/\ndPcHEkKIG3OpWe4mTJiAHTt2ID8/H0lJSVi/fr2zQyKEEEKIm3KpJGfixImQyw2NSyNGjEBFRYWT\nIyKEEEKIu3KpJKe9zz//HJMnT3Z2GIQQQghxU71ek7N48WLU1NSYbX/yyScxY8YMAMDbb78NlmUx\nf/78Ls+Vm5uL3Nxci/sSEhJ6HiwhhBBC3BYjSZLk7CDa27ZtGz799FN8+OGHUKlsX0iwzYYNG/Dg\ngw/aMTJCCCGEuBOXSnL27duHf/7zn/joo48QEuK6q5oSQgghxPW5VJIzc+ZM6PV6BAUZFhJLT0/H\n6tWrnRwVIYQQQtyRSyU5hBBCCCH24rKjqwghhBBCeoKSHEIIIYR4JJda1sHRjh07hldeeQVDhw69\n4XNUVVUBACIierZKqyudx56xFBUVYezYsT0+T0/jcaX7a6/zuNL9tdd5XCmWtvP09B674mvq6Xlc\n6f7aKx5Xur/2Oo8r3V97naf9OWJiYmweNd2nanKOHTuGxx57DAcPHuzROQBg5MiRPY7FVc5jz1h6\nen/tFY8r3V97nceV7q+9zuNKsbSdhz4jHBuLq7yHXen+2us8rnR/7XWenp6DuqsIIYQQ4pH6VEsO\nAEyePBn79u1zdhgei+6vY9H9dTy6x45F99ex6P6aopYcQgghhHgkSnIIIYQQ4pEoySGEEEKIR2Jf\nfPHFF50dRG/racU46RrdX8ei++t4dI8di+6vY9H9va7PFR4TQgghpG+g7ipCCCGEeCRKcgghhBDi\nkSjJIYQQQohHoiSHEEIIIR6pzyQ5+/btw+zZszFz5ky88847zg7HbU2bNg1ZWVnIzs7GggULAAAN\nDQ3IycnBrFmzkJOTg8bGRgCAJElYs2YNZs6ciaysLJw+fdqZobuslStXYvz48Zg3b55x243c023b\ntmHWrFmYNWsWtm3b1uuvw1VZur/r1q3DpEmTkJ2djezsbPzwww/GfevXr8fMmTMxe/Zs7N+/37id\nPkMsu3r1Ku6//37MmTMHmZmZ2LBhAwB6D9tTZ/eY3sdWkPoAnuel6dOnSyUlJZJOp5OysrKkwsJC\nZ4fllqZOnSrV1taabHvppZek9evXS5IkSevXr5defvllSZIkae/evdKSJUskURSlY8eOSXfeeWev\nx+sODh06JJ06dUrKzMw0brP1ntbX10vTpk2T6uvrpYaGBmnatGlSQ0ND778YF2Tp/v773/+W3nvv\nPbNjCwsLpaysLEmn00klJSXS9OnTJZ7n6TOkC5WVldKpU6ckSZIktVotzZo1SyosLKT3sB11do/p\nfdy9PtGSc+LECSQmJiI+Ph5KpRKZmZnYs2ePs8PyGHv27MFtt90GALjtttuwe/duk+0Mw2DEiBFo\nampCVVWVM0N1SRkZGQgMDDTZZus9PXDgACZMmICgoCAEBgZiwoQJJr/e+jJL97cze/bsQWZmJpRK\nJeLj45GYmIgTJ07QZ0gXIiIiMGTIEACAn58fUlJSUFlZSe9hO+rsHneG3sfX9Ykkp7KyElFRUcbH\nkZGRXb5BSNeWLFmCBQsW4LPPPgMA1NbWIiIiAoDhj7Gurg6A+X2Pioqi+24lW+8pvcdt9/HHHyMr\nKwsrV640dqV0dh/p/lqnrKwMZ8+eRXp6Or2HHaT9PQbofdydPpHkSBbmO2QYxgmRuL9PPvkE27Zt\nw7vvvouPP/4Yv/76a6fH0n23v87uKd1r2yxatAjffvst8vLyEBERgX/+858A6P72REtLC5YvX47n\nnnsOfn5+nR5H9/jGdbzH9D7uXp9IcqKiolBRUWF8XFlZafyFQWwTGRkJAAgNDcXMmTNx4sQJhIaG\nGruhqqqqEBISAsD8vldUVNB9t5Kt95Te47YJCwsDy7KQyWRYuHAhTp48CaDzzwq6v13jOA7Lly9H\nVlYWZs2aBYDew/Zm6R7T+7h7fSLJGTZsGIqLi1FaWgq9Xo+dO3di2rRpzg7L7Wg0GjQ3Nxv//8cf\nf8SAAQMwbdo0bN++HQCwfft2TJ8+HQCM2yVJwm+//QZ/f3+P/4OyF1vv6cSJE3HgwAE0NjaisbER\nBw4cwMSJE535Elxa+9qw3bt3Y8CAAQAM93fnzp3Q6/UoLS1FcXExhg8fTp8hXZAkCc8//zxSUlKQ\nk5Nj3E7vYfvp7B7T+7h7cmcH0BvkcjlWrVqFhx56CIIg4I477jC+GYj1amtr8fjjjwMABEHAvHnz\nMHnyZAwbNgxPPvkktm7diujoaPzv//4vAGDKlCn44YcfMHPmTKhUKvzjH/9wZvgu66mnnsKhQ4dQ\nX1+PyZMnY9myZXjkkUdsuqdBQUH44x//iDvvvBMA8PjjjyMoKMhpr8mVWLq/hw4dwrlz5wAAsbGx\nWL16NQBgwIABmDNnDubOnQuWZbFq1SqwLAsA9BnSiSNHjiAvLw+pqanIzs4GYLjn9B62n87u8Y4d\nO+h93A1aoJMQQgghHqlPdFcRQgghpO+hJIcQQgghHomSHEIIIYR4JEpyCCGEEOKRKMkhhBBCiEei\nJIcQQgghHomSHELcyKBBg5CdnY158+Zh+fLlaG1t7fTYsrIyDB8+HNnZ2Zg7dy6eeeYZcBwHAPjl\nl18watQoZGdnIysrC4sXL0ZtbS0A4IsvvsC4ceOMz9u8eTMA4Msvv0RWVhaysrJw9913G+fnKCsr\nw7x58xz8yoGHH34YTU1NDr+Otdddt24d3n//fYdd9+677wZguL/5+fnG7SdPnsSaNWscdl1CPAkl\nOYS4EW9vb+Tl5WHHjh1QKBT49NNPuzw+ISEBeXl5yM/PR0VFBXbt2mXcN3r0aOO+YcOG4eOPPzbu\nmzt3LvLy8rBp0ya89tprqKmpQVxcHD766CPk5+fjsccew1//+leHvU5L3n33XQQEBPTqNZ153bZ/\n2ytXrmDHjh3G7cOGDcNf/vKXXo+HEHfUJ2Y8JsQTjR49GgUFBfjXv/6FmJgY3HvvvQAMLQy+vr7G\n9W0AgGVZDB8+3OKKw5IkoaWlBQkJCWb7QkNDkZCQgPLyctx0003G7SNGjDBZA4fneTz77LM4c+YM\nkpOT8dJLL0GlUpmca/Pmzfjss8/AcRwSExPx8ssvQ6VSYcWKFfDy8sKlS5dQXl6OtWvXYtu2bfjt\nt9+Qnp5uXHRw2rRp2Lp1KzQaDR5++GGMGjUKx44dQ2RkJN566y14e3vj7NmzeOGFF9Da2oqEhAT8\n4x//QGBgoEkc3333Hd5++21wHIegoCC88sorCAsLQ0tLC9asWYNTp04BAJYuXYrZs2cbrxsSEoK3\n334b27dvR3R0NEJCQjBkyBCze7ZixQoolUpcuHABtbW1WLFiBaZOnQqdTocXX3wRp06dAsuyWLFi\nBcaNG4fCwkKsXLkSHMdBFEWsW7cOSUlJGDlyJI4dO4ZXX30VFy9eRHZ2Nm6//XYMGjQIH3zwAdav\nX4+GhgY899xzKC0thUqlwurVq5GWloZ169ahvLwcZWVlKC8vx4MPPogHHnig2/cUIZ6GWnIIcUM8\nz2Pfvn1ITU1FZmamSQvNrl27cOutt5ocr9PpcPz4cUyaNMm47fDhw8jOzsYtt9yCn376yTidfnul\npaUoLS01S4C2bt2KyZMnGx8XFRXhd7/7HfLz8+Hr64v//ve/ZueaOXMmPv/8c3z55ZdISUnB1q1b\njfuampqwceNGrFy5En/4wx+wePFi7Ny5E+fPn8fZs2fNznX58mXce++92LlzJ/z9/fHNN98AAJ55\n5hk8/fTTyM/PR2pqKt544w2z544aNQqbN2/G9u3bkZmZiffeew8A8NZbb8HPzw/5+fnIz8/HuHHj\nTJ536tQpfPXVV9i+fTveeOMN42KIlly5cgUfffQR1q9fjxdeeAE6nc7YUpafn49XX30VK1asgE6n\nw6effooHHngAeXl5+PzzzxEVFWVyrj//+c/GVrfFixeb7Fu3bh0GDx6M/Px8/OlPf8Kzzz5r8m/y\n/vvvY8uWLXjzzTeNXZWE9CXUkkOIG9Fqtca1a0aPHo0777wTSqUStbW1qKysRH19PQICAhATE4Oy\nsjKUlJQgOzsbly9fxuzZs5GWlmY81+jRo7F+/XoAwDvvvIOXX37ZuPbNV199hSNHjkCpVGL16tUm\nawgdPHgQW7duNUlkoqOjMWrUKADA/PnzsWnTJixZssQk9sLCQrz++utQq9VoaWkxWXxx6tSpYBgG\nAwcORFhYGAYOHAgA6N+/P65cuYJBgwaZnCsuLs64bciQIbhy5QrUajXUajXGjBkDALj99tvxxBNP\nmN3DiooK/OlPf0J1dTX0ej3i4uIAAD///DNee+0143EdW4AOHz6MGTNmGFuoulrYcM6cOZDJZEhK\nSkJ8fDwuXbqEI0eO4L777gMA9OvXDzExMSgqKsKIESPwn//8BxUVFZg1axaSkpI6PW9HR44cwbp1\n6wAA48ePR0NDA9RqNQDDGlFKpRIhISEICQlBbW2tWQJFiKejJIcQN9JWk9PR7Nmz8c0336CmpgaZ\nmZnG7W01OVVVVbj//vuxZ88e42rQ7U2fPh3Lli0zPp47dy5WrVpldty5c+fwl7/8Be+++y6Cg4ON\n2xmGMTmu42PA0I3z1ltvIS0tDV988QUOHTpk3KdUKo3Pa/t/AJDJZOB53uxc7Y9hWRY6nc7smM6s\nWbMGixcvxvTp0/HLL78YW3skSbIYd3evy5rjGIZBZ8sEZmVlIT09HXv37sWSJUuwZs0ajB8/3qrr\nWDpn27U73iNL95EQT0fdVYR4gMzMTHz11Vf45ptvMHv2bLP9ERERePrpp/HOO+9YfP6RI0cs1uS0\nV15ejmXLluHll19GcnKy2b5jx44BAHbu3Gls1WmvpaUF4eHh4DjOZLSQvfj7+yMgIACHDx8GAOTl\n5SEjI8PsOLVajcjISADA9u3bjdsnTJiAjz76yPi4sbHR5HkZGRn49ttvodVq0dzcjO+//77TWL7+\n+muIooiSkhKUlpYiOTkZGRkZxtddVFSEq1evIiUlBaWlpYiPj8cDDzyAadOmoaCgwORcvr6+aGlp\nsXidjIwMfPnllwAMI+aCg4Ph5+fXaVyE9DXUkkOIBxgwYABaWloQERGBiIgIi8fMmDED69atMyYB\nbTU5kiTB39+/22HJb775JhoaGvC3v/0NgKF14IsvvgBg6H7Ztm0bVq1ahaSkJCxatMjs+U888QQW\nLlyI2NhYpKamdvrF3RMvvfSSsfA4Pj4ea9euNTtm6dKleOKJJxAZGYn09HSUlZUBAB577DGsXr0a\n8+bNg0wmw9KlS02Kt4cMGYK5c+ciOzsbsbGxFhO5NsnJybjvvvtQW1uLv/3tb/Dy8sI999yDF154\nAVlZWWBZFmvXroVSqcRXX32FL7/8EnK5HGFhYXj88cdNzjVw4ECwLIv58+djwYIFJl13S5cuxcqV\nK5GVlQWVSmUs0iaEGDBSZ22ohBBCbLZixQrccsstZsXfhJDeR91VhBBCCPFI1JJDCCGEEI9ELTmE\nEEII8UiU5BBCCCHEI1GSQwghhBCPREkOIYQQQjwSJTmEEEII8UiU5BBCCCHEI/1/G02+v8sSmREA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc090dcd68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax, x, D = plot_tajd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_diversity(vcf_data=vcf_data, samples_to_use=vcf_samples_in_148,\n",
    "              size=300, start=54438, stop=63146 , step=3, despine_bottom=False,\n",
    "              highlight_regions=[[169, 470]], ax=None):\n",
    "    pos_in_gene = [position_in_gene(x) for x in vcf_data['variants/POS']]\n",
    "    ac = allel.GenotypeArray(vcf_data['calldata/GT'][: ,samples_to_use, :]).count_alleles()\n",
    "    pi, windows, number_of_bases, counts = allel.stats.diversity.windowed_diversity(\n",
    "        pos=pos_in_gene, ac=ac, size=size, start=position_in_gene(start), stop=position_in_gene(stop), step=step,\n",
    "        fill=0.0\n",
    "    )\n",
    "    pi[np.isnan(pi)] = 0.0\n",
    "    \n",
    "    # set up axes\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(figsize=(8, 2))\n",
    "#         fig.tight_layout()\n",
    "    \n",
    "    x = (windows[:, 0] + windows[:, 1]) / 6\n",
    "    \n",
    "    # plot\n",
    "    sns.despine(ax=ax, offset=5, bottom=despine_bottom)\n",
    "    ax.plot(x, pi)\n",
    "    if despine_bottom:\n",
    "        ax.set_xlabel('')        \n",
    "        ax.set_xticks([])        \n",
    "    else:\n",
    "        ax.set_xlabel('PvRBP2b amino acid position')\n",
    "    ax.set_xlim((position_in_gene(start)/3, position_in_gene(stop)/3))\n",
    "    for region in highlight_regions:\n",
    "        ax.axvspan(region[1], region[0], facecolor='k', alpha=0.1)\n",
    "    ax.set_ylabel('Nucleotide diversity (${\\pi}$)')\n",
    "    \n",
    "    return(ax, x, pi)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nfs/users/nfs_r/rp7/miniconda3/envs/biipy/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['arial'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAC5CAYAAAA1Zgt+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAIABJREFUeJzs3XtYVHX+wPH3mRkGkDsq4AVQDO+XNC39lbZiSIqEF6y1\nsiSr3dKym5a7LRVPu7bdzLW1LDe2q6UmXprKFCu6aKXpoqWmKYEigyJ3gWFmzu+PydGR2wwOcvHz\neh6fhznne875zPEw8+F7VVRVVRFCCCGEEA3StHQAQgghhBBtgSRNQgghhBBOkKRJCCGEEMIJkjQJ\nIYQQQjhBkiYhhBBCCCdI0iSEEEII4QRJmoQQQgghnCBJkxBCCCGEEyRpEkIIIYRwgiRNQgghhBBO\nkKRJCCGEEMIJkjQJIYQQQjhB19IBtBVvvvkmt99+e0uHYZeXl9fSIbhd165dWzoEIYQQol5S0+Sk\n9pikCCGEEMJ5kjQJIYQQQjhBkiYhhBBCCCc0KWk6ffo0FovF3bEIIYQQQrRaTnUEt1qtGAwGNm7c\nyJ49e9Dr9ZhMJoKDgxkzZgw33XQTPXr0aOZQRWtXZbLw+S4jKnBl3450CvBs6ZCEEEIIt3Eqabrt\nttsYNWoUDz30EL1790ajsVVQFRcX89133/HCCy9w3XXXkZiY2KzBitbt2MlK/Dp40NHfk1/zyiVp\nEkII0a44lTSlpaXh4eFBRUWFPWECCAwMJC4ujri4OGpqapy6YGZmJn//+9+xWq1Mnz6du+++22G/\nyWRiwYIF/PTTTwQGBrJ48WK6d+9OUVER999/P3v37mXKlCmkpKQAUFlZybx588jJyUGr1TJ27Fge\neeQRANauXcuzzz5LaGgoALfeeivTp093Kk7hnJJyE1mHSzheWElxuYnxI7qg1SgUl5taOjQhhBDC\nrZzq0+Th4QHAlClTOHToEGaz2b7vyJEjDmUaYrFYSE1NZcWKFRgMBj766CMOHTrkUGb16tX4+/uz\nefNmZs2axfPPPw+Ap6cn8+bNY8GCBbXOe8cdd/Dpp5+Snp7Ojz/+yJdffmnfN3HiRNavX8/69esl\nYXKz/x0qYv03x+jor+fGsRHcnXAZPcJ80GjAalVbOjwhhBDCrVya3LKkpISnnnqKI0eOEBQURFRU\nFD///DObN2926visrCwiIyMJDw8HID4+noyMDC677DJ7ma1btzJ37lwA4uLiSE1NRVVVOnTowPDh\nw8nJyXE4p7e3NyNHjgRAr9fTv39/jEajK29LNFFReQ1xI7oQGuzlsF2jKEjOJIQQor1xKWnq0qUL\nb7/9NgD5+fn8+uuveHo632/FaDQSFhZmfx0aGkpWVlatMl26dLEFp9Ph5+dHUVERwcHBjZ6/tLSU\nzz//3GHm7s8++4wffviBnj17snDhQvu5xYWzWlU0ddRVahRFapqEEEK0Oy4lTWVlZfz4449ER0cT\nFhbmkAA5Q1Vrf5EqiuJymbqYzWYeeughZs6caa/JGjt2LJMmTUKv17Ny5UoeffRR3nrrrXrPkZaW\nRlpaWp37IiIiGo3hUmNVVTSa2v83Go1tnxBCCNGeuJQ0nT59mjfeeIODBw9SU1NDVFQU0dHRPPro\no04dHxYWRn5+vv210WgkJCSkVpnjx48TFhaG2WymrKyMwMDARs/9t7/9jR49ejBr1iz7tqCgIPvP\nN954o71/VH2Sk5NJTk6uc9+iRYsajeFSY7WqaOpIaDUaBYvUNAkhhGhnnEqaVFVFURQ++OADe41L\ndXU1hw4d4uDBgw5lGjJo0CCys7PJzc0lNDQUg8HACy+84FAmJiaG9PR0hg4dyqZNmxg5cmSj5128\neDHl5eX8/e9/d9heUFBgT8q2bt1Kr169nHm7wklWlbprmhQFqWgSQgjR3jg9T9P48eMZN26cfZun\npyfR0dGUlpby6KOPctVVVzF16tSGL6bTkZKSwp133onFYmHatGlER0ezZMkSBg4cyLhx40hKSmL+\n/PnExsYSEBDA4sWL7cfHxMRQXl5OTU0NW7Zs4Y033sDX15dXX32VqKgopkyZApydWuDtt99m69at\naLVaAgICpLbIzWw1TbW3n6lp+mjbMfy8Pbj28pDahYQQQog2RlHr6kR0nurqatasWcPGjRs5evQo\n/v7+VFVVoaoqV199Nbfccgv9+vW7GPG2mEWLFrFw4cKWDsMuLy+vpUPgo23HGDs0FB8vx9y7/HQN\nX2adoKLSNjXFjWOd6w/WtWtXt8cohBBCuItTNU2enp7ccsst3HLLLdTU1FBUVISXlxf+/v7NHZ9o\nxaxW6uzTpGgUqqoteOo1VJusLRCZEEII4X4udQQH2ySW53feFpcmq6qirWfKgSqThY5enpI0CSGE\naDecmhFciLpYrHV3/tdqbOvQeerl8RJCCNF+yLeaaDLVqqKtoye43kPLreN7MGpApxaISgghhGge\nLjfPCXGGVYX6ZoPo2tH74gYjhBBCNDOXapreeecdSkpKmisW0QY5M1u7EEII0R64lDSdOHGCpKQk\n5s2bR2ZmZp1LngghhBBCtEcuJU0PPvggn332GUlJSaSnpzN+/HhefPFFcnJymis+IYQQQohWweWO\n4Iqi0LlzZzp16oRWq6WkpIT777+fZ599tjniE0IIIYRoFVzqCP7WW2+xbt06goKCSEpKYsGCBXh4\neGC1Whk/fjwLFixorjiFEEIIIVqUS0mT0Whk6dKldOvWzb7tueeeY/78+SxfvtztwYnWo6LSjIdO\nYV9OGftzSukS7EVFlbmlwxJCCCEuGpea57799luHhAkgMzMTgF69erkvKtHqvPXZETZ8c4wDuaVM\nuaYbWo3CkeMVjR6nAjnGxssJIYQQrZ1TSdN7771HQkICR44cISEhwf4vJiaGPn36uHTBzMxM4uLi\niI2N5bXXXqu132Qy8cADDxAbG8v06dM5evQoAEVFRcycOZOhQ4eSmprqcMzevXtJSEggNjaWp59+\n2j6qr7i4mOTkZMaPH09ycrJMl3ABOgZ4cmZBFL2Hls5BXpgtjS+RMrJ/R46fqmre4IQQQoiLwKmk\nKSEhgVdeeYWYmBheffVV+7+1a9fy/PPPO30xi8VCamoqK1aswGAw8NFHH3Ho0CGHMqtXr8bf35/N\nmzcza9Ys+/k9PT2ZN29enf2mnnzySVJTU/nss8/Izs6213699tprjBo1is8++4xRo0bVmaSJptFp\nFCzWxqecUECmphBCCNEuOJU0+fn50b17d1588UW6detm/xcYGOjSxbKysoiMjCQ8PBy9Xk98fDwZ\nGRkOZbZu3cqUKVMAiIuLY9u2baiqSocOHRg+fDienp4O5QsKCigvL2fo0KEoisLkyZPt58zIyGDy\n5MkATJ48mS1btrgUrzjPObmPTqtgsTiRNCkKkjMJIYRoD5zqCD5jxgxWrlxpT0zgbO2Boij8+OOP\nTl3MaDQSFhZmfx0aGkpWVlatMl26dLEFp9Ph5+dHUVERwcHBTp0zLCwMo9EIQGFhISEhIQCEhIRw\n6tQpp+IUjdNqFczO1DQpSNIkhBCiXXAqaVq5ciUAu3btuqCL1dVMc/4yHM6UuZDyDUlLSyMtLa3O\nfREREU06Z3ul02gwO1PThDTPCSGEaB9cmnLgk08+YfTo0fj6+rJs2TJ+/vln7r33Xvr37+/U8WFh\nYeTn59tfG41Ge03QuWWOHz9OWFgYZrOZsrKyBpsBzz9nfn6+/ZwdO3akoKCAkJAQCgoK6q2tOiM5\nOZnk5OQ69y1atKjR93cp0WoVpzqCK4qCpExCCCHaA5emHFi2bBm+vr7s2LGDr7/+msmTJ/PEE084\nffygQYPIzs4mNzcXk8mEwWAgJibGoUxMTAzp6ekAbNq0iZEjRzZYcxQSEoKPjw+7d+9GVVXWrVvH\nuHHj7Odat24dgMN20UTn/Dd08NQysGeAU8dIRZMQQoj2wKWkSavVAvDll18yY8YMrrvuOmpqapw+\nXqfTkZKSwp133snEiROZMGEC0dHRLFmyxN55OykpieLiYmJjY0lLS+ORRx6xHx8TE8MzzzxDeno6\nY8aMsY+8e/LJJ3n88ceJjY0lIiKCMWPGAHD33XfzzTffMH78eL755hvuvvtuV96uON85yU8HLx1z\np/Ru9BAN0jwnhBCifVBUF77R/vSnPxEaGsq3337L2rVr8fLyIikpiQ0bNjRnjK3CokWLWLhwYUuH\nYZeXl3dRr7fqixxb0qTAjX9wvn+X8VQV+3NLuXZISKNlu3btegERCiGEEM3LpZqml156iWuuuYYV\nK1bg7+9PcXGxrDd3CbCqqusrO58hzXNCCCHaCac7gquqSnFxMePHj7dvCwkJqdWRW7Q/VquKolHA\nidFy59Mo0jwnhBCifXC6AkFRFObMmdOcsYhWymoFraZp0ziA1DQJIYRoH1yacmDIkCFkZWUxePDg\n5opHtCIZP+ZTWGLCYlXp1smb/KIqXE2dbFMOSNYkhBCi7XMpafruu+94//336datG97e3vbtGzdu\ndHtgouUVlpi4ceyFTeqpkT5NQggh2gmXkqbXX3+9ueIQ7ZgkTUIIIdoDlwZFde3alR07dpCenk63\nbt1QFIXCwsLmik20tKZ3Yzp7CkWRjuBCCCHaBZeSpieffJLdu3djMBgA8PHx4amnnmqWwEQr4IZc\nR6O45TRCCCFEi3MpacrKyuKJJ57A09MTgICAAJdmBBeXIEWR5jkhhBDtgktJk06nw2Kx2NeCO3Xq\nFBpNk6c9FK2dW5rnkKomIYQQ7YJLHcFnzpzJnDlzKCwsZPHixXz66ac88MADzRWbaAcUwCpZkxBC\niHbApaTphhtuYMCAAWzfvh1VVVm2bBm9evVqrthEO6BI85wQQoh2wqW2tf/+97/4+Phwyy23cOut\ntzYpYcrMzCQuLo7Y2Fhee+21WvtNJhMPPPAAsbGxTJ8+naNHj9r3LV++nNjYWOLi4vjqq68AOHz4\nMImJifZ/w4YN47///S8AS5cuZfTo0fZ9X375pcvxXtLckOxI85wQQoj2wqWapvLycmbPnk1AQADx\n8fHExcXRqVMnp4+3WCykpqaSlpZGaGgoSUlJxMTEcNlll9nLrF69Gn9/fzZv3ozBYOD555/npZde\n4tChQxgMBgwGA0ajkeTkZDZt2kRUVBTr16+3n3/MmDHExsbazzdr1ixmz57tytsUbqRgW/BXCCGE\naOtcqmmaO3cuBoOBlJQUCgoKuPXWW5k1a5bTx2dlZREZGUl4eDh6vZ74+HgyMjIcymzdupUpU6YA\nEBcXx7Zt21BVlYyMDOLj49Hr9YSHhxMZGUlWVpbDsdu2bSM8PJxu3bq58rZEfdw2T9OFn0cIIYRo\naU0a+taxY0c6depEYGCgS5NbGo1GwsLC7K9DQ0MxGo21ynTp0gWwjdbz8/OjqKjIqWMNBgOTJk1y\n2Pbuu++SkJDAwoULKSkpcTpWgfua54QQQoh2wKXmuffee49PPvmEwsJCrr/+ep5++mmHprXG1DUz\ntHLet2p9ZRo71mQysXXrVh5++GH7thkzZnDvvfeiKApLlizhmWeeYdGiRfXGl5aWRlpaWp37IiIu\nbA22S5WiSPOcEEKI9sGlpOn48eP89a9/pW/fvk26WFhYGPn5+fbXRqORkJCQWmWOHz9OWFgYZrOZ\nsrIyAgMDGz02MzOTAQMGOPSxOvfn6dOn8+c//7nB+JKTk0lOTq5zX0PJVrvljuY5pHlOCCFE++BU\n89yMGTMAeOedd7j55psZNmwYw4YNY+jQoQwbNszpiw0aNIjs7Gxyc3MxmUwYDAZiYmIcysTExJCe\nng7Apk2bGDlyJIqiEBMTg8FgwGQykZubS3Z2NoMHD7YfZzAYiI+PdzhXQUGB/ectW7YQHR3tdKzC\nPbRahVNl1az6Ioc3Pj5MUZmppUMSQgghmsSpmqaVK1cCsGvXrgu7mE5HSkoKd955JxaLhWnTphEd\nHc2SJUsYOHAg48aNIykpifnz5xMbG0tAQACLFy8GIDo6mgkTJjBx4kS0Wi0pKSlotVoAKisr+fbb\nb0lNTXW43nPPPcf+/fsB6NatW639ovlpNQqzJ9qmptj+80nKTtcQ5Kdv4aiEEEII1ymqLEHvlEWL\nFrFw4cKWDsMuLy+v2a+x6vMcbhzrvr5c3+8rJCTIkx5hvnXu79q1q9uuJYQQQribUzVN9XWOPqO+\nfkBCnEtRwGpt6SiEEEKIpnEqaaqoqADgyJEj7Nmzx94P6fPPP2f48OHNF51oWW6eLkCrUWQknRBC\niDbLqaRp7ty5ANxxxx2sXbsWX19f+/Z58+Y1X3SiXdEoClarJE1CCCHaJpcmt8zLy0OvP9uJV6/X\nc+zYMbcHJdonRaMgOZMQQoi2yqV5mhITE0lKSiI2NhZFUdi8ebN9yRMhGqPVIDVNQggh2iyXkqZ7\n7rmHMWPGsGPHDsA2oqx///7NEphofzSK9GkSQgjRdrmUNAEMGDCAAQMGNEcs4iJ54+PDnCoz8WBS\nb7TaJi0/2CSKomC1SNIkhBDCdRaryumqGvvr4rJqOgd546V3OZVpsot3JdFq+HrrCPbXY7aq/D4/\n6EWh1UC11DQJIYRogs935HA4r5QAH1vf6hqLlbIKEx462xeZxWLlT1MHN3SKCyZJ0yVKq1GwWFTw\nuHjXlNFzQgghmqq6xsqka3rStVPdEyS/vn5Ps8fgUtuMqqqsX7+el19+GbCNpsvKymqWwETz0moV\nLBc5gZHRc0IIIZrKYrWiUdw8gaCLXEqannzySXbv3o3BYADAx8eHp556qlkCE81IAa3ScNJ0vLCS\napN7p+/WahQO5JTy7d6TFBRVufXcQggh2jerFbSa+tMWHy8PVqzfy4r1e6morKm33IVwqXkuKyuL\n9PR0Jk+eDEBAQAA1Nc0TmGhejdU0lZ2uYVDPALdeMyK0A1pNCJUmCz8cOEX8SFlrTgghhHOsVisN\n5EzcHNe32WNwqaZJp9NhsVhQfq8eO3XqFJqG3kEdMjMziYuLIzY2ltdee63WfpPJxAMPPEBsbCzT\np0/n6NGj9n3Lly8nNjaWuLg4vvrqK/v2mJgYEhISSExMZOrUqfbtxcXFJCcnM378eJKTkykpKXEp\n1vaixmylqMzEqdJq9h4podpkPdunqR4WK2g07q0G1SgK4SEd6N7JW/o2CSGEcInFqjZY03QxuHT1\nmTNnMmfOHE6ePMnixYuZMWMGf/rTn5w+3mKxkJqayooVKzAYDHz00UccOnTIoczq1avx9/dn8+bN\nzJo1i+effx6AQ4cOYTAYMBgMrFixgqeeegqLxWI/7s0332T9+vWsXbvWvu21115j1KhRfPbZZ4wa\nNarOJO1SsHlnPt/uPcnOX4rIO3maQ8fKbElTA4mLqqoNZvQXQqNRkEF0QgghXGG1qm7/Y95VLn0t\n3nDDDcyfP58///nPhISEsGzZMiZMmOD08VlZWURGRhIeHo5eryc+Pp6MjAyHMlu3brXPMh4XF8e2\nbdtQVZWMjAzi4+PR6/WEh4cTGRnZaCf0jIwMe1Pi5MmT2bJliytvt90wm1VihoUSOzyM2OFhzP9j\nv0aTJotVbbYOd5pG+lMJIYQQ52sNSZNTfZrS0tLq3J6ZmUlmZibJyclOXcxoNBIWFmZ/HRoaWivx\nMRqNdOnSxRacToefnx9FRUUYjUaGDBnicKzRaLS/nj17NoqicNNNN3HTTTcBUFhYSEhICAAhISGc\nOnXKqTjbG+s5tUaKouDrrSM0yIsvdxfgodMQFuzFmCEh5x0DHs30cGo0yMzgQgghXGJRVbRtIWmq\nqKgA4MiRI+zZs4eYmBgAPv/8c4YPH+70xdQ6viiV82oz6ivT0LErV64kNDSUwsJCkpOTiYqKYsSI\nEU7HdUZaWlq9CWJERITL52strHXUGl3W3Y/LuvsB8MHnOXUf01xJk6Kg/l7TpKoqNWYVvUfLtlML\nIYRo3dpMTdPcuXMBuOOOO1i7di2+vr727fPmzXP6YmFhYeTn59tfG41Ge03QuWWOHz9OWFgYZrOZ\nsrIyAgMDGzw2NDQUgI4dOxIbG0tWVhYjRoygY8eOFBQUEBISQkFBAcHBwQ3Gl5ycXG+t2aJFi5x+\nn62NVW24U3ddrXBWVaW5nk1FgZyC03z5vwJOFFezP6eUx2fK0jxCCCHqZ7W2fE2TS3/e5+Xlodfr\n7a/1ej3Hjh1z+vhBgwaRnZ1Nbm4uJpMJg8Fgr7U6IyYmhvT0dAA2bdrEyJEjURSFmJgYDAYDJpOJ\n3NxcsrOzGTx4MKdPn6a8vByA06dP88033xAdHW0/17p16wBYt24d48aNc+XtthtNSYCa8+FUFIUD\nuWX0DfdneJ8gBjZxaoOCotMO6xAJIYRov5qzr62zXJqnKTExkaSkJGJjY1EUhc2bN9s7bTt1MZ2O\nlJQU7rzzTiwWC9OmTSM6OpolS5YwcOBAxo0bR1JSEvPnzyc2NpaAgAAWL14MQHR0NBMmTGDixIlo\ntVpSUlLQarUUFhYyZ84cwDY6b9KkSYwZMwaAu+++mwceeIA1a9bQpUsXlixZ4srbbVfObwY9l06j\n8P7W3xjSK5B+kbYEprHaKXcIDfaiymRBUYqbdPzSVbvpExnErdf3c3NkQgghWhuzxdo2mufOuOee\nexg9ejQ7d+4EbE1W/fv3d+mC1157Lddee63DtnOb+Dw9PfnXv/5V7/Xvueceh23h4eFs2LChzvJB\nQUG8+eabLsXXLjXS53rqmHCKykzsOHDqbNJ0kTJ6rabp69F16+xLTU3ds5av/fwghaXnzDquwg1j\nehEa3KFJ1xJCCNGycvLLWjoE1xfsHThwIAMHDmyOWEQL0iiOI9psHe4u0nWbuFqLVqtgrufgU6XV\n3JU4yP4669AJsvNKJGkSQog2qmdX965S0RROJU0zZsxg5cqVDB061KGZR1VVFEXhxx9/bLYAxcWh\nKArqOfmHRb04NU2KRsHShOkHzrRtO3uoh1ZLqcXk8nWEEEKIM5xKmlauXAnArl27mjUY0XI0ClhR\n+WDrbyi/Tz55Re+GRxu657pKo82HdTFbrHjoNJjqaZ5Tzzuph05Djdm9CxALIYS4tLjUAPPcc885\ntU20Mk5UGCkaW02TolG4cWwEM8ZF4uPlcuut03y9L+zcZrMVnbbux7euOb0kaRJCiLarrs/1luBS\n0vTtt9/W2paZmem2YETLOb9PU3N7/p7LL+h4s8WKVqtQVFbFjn1Gsg6dsO8rr6zB47yESpImIYRo\nu1rDHE3gZPPce++9x8qVK8nNzSUhIcG+vaKigmHDhjVbcMI5X2ed4FhhZb0VStX1NGGdS1GUi5o0\nNTQFgjPMFiseWg1Trr2MstMmtmUd54efjXjoNOTkl/HH2D4O5XU6DcZTpzHVWNB7aC/o2kIIIS4u\ns1VFq20jSVNCQgJjxozhxRdf5OGHH7Zv9/HxITAwsNmCE87JO1XJTWMvbJkXVzpVu50LvwcHfjtF\n+pe/0tHfi/8b3JXLwm3P35DozpgttuRQq9XU+ouko78XHbx0bNtznGuHdXdb6EII0ZDyyhr2/nrS\n/vqy7oF0CvRuwYjaJovFivZiDOluhFNJk5+fH35+frz44ovs37+fHTt2ADB8+HBJmtoJRaHJ8yVd\nTDn5ZUz9w2X0jghy2K7RKOg19dcgabUaRvQL5bdWMM+HEOc7WlCGRqPQtZNvS4ci3KC4rJolH+yi\nW2dfLFYrI/qF4eGhodpkYfP3OcwY36fxkwgHlrZU03TGW2+9xapVq4iNjQVg/vz53HjjjcycObNZ\nghMXj0aj0AZyJn7YZ2T05d2adKxWq8Ha1EmhhGhGH3+bjU6r4Y4EWYOxPSivNHFl/1Am/F9Ph+0W\nq8r/Dp6o5yjREIulDfVpOmP16tWsWrWKDh1sEwTedddd3HTTTZI0tQMapQVHJ7hw2c5B3nh5Nm3k\nnVZjm0pBiNZGUcAiCX27UW2y4KmvXfOt1SicKK7kl5yiWrXlomEWaxtqnjuXVqut82fRtimK0pTp\nktoUSZqEEBdDdY0FT4+6v14nj+nFN//L48tdR6msMjNzQj+C/L0ucoRtj8WiomtrzXNTp05l+vTp\n9ua5LVu2MG3atGYJTIhzvfXxz5hqrESE+jX5HBqNgsUiSZNonZTfR0S89fHPlFfWcO+0IS0ckXDV\niaJKNn//G7/kFDHpmqg6y/TtEUzfHraJgz/fmUthSZUkTU6wtKUpB85ITk7myiuvZOfOnaiq2qQF\nezMzM/n73/+O1Wpl+vTp3H333Q77TSYTCxYs4KeffiIwMJDFixfTvbtttNPy5ctZs2YNGo2Gxx9/\nnNGjR3P8+HEWLFjAyZMn0Wg03Hjjjdx+++0ALF26lFWrVhEcbHtAH3rooVqLBYu2wVRj5c7EC1vz\nUKvRSE2TaHFbd+QwJLozHQO8+elwIZm7jvK/gycID/XjvU37KT9dI9NitFHvfLqPaWMv4/pRPfD3\n0TdaXuaPc57ZYkVTz4TGF5NLSZOqqhw6dIjS0lLmzp1LXl4eWVlZDB482KnjLRYLqamppKWlERoa\nSlJSEjExMVx22WX2MqtXr8bf35/NmzdjMBh4/vnneemllzh06BAGgwGDwYDRaCQ5OZlNmzah1Wp5\n7LHHGDBgAOXl5UybNo2rr77afs5Zs2Yxe/ZsV96maKe0WoWtO3I4UXSahNFRhHX0aemQxCXox/0n\n+PVYCQoKigKzbxhIRWUNHbw9sFiseOi0vGn4udHznFn7U7Qevh08iAjzd7q8XqfFZLY0Y0Tth9Wq\nomtrNU1PPvkkGo2G7du3M3fuXHx8fLjvvvv48MMPnTo+KyuLyMhIwsPDAYiPjycjI8Mhadq6dStz\n584FIC4ujtTUVFRVJSMjg/j4ePR6PeHh4URGRpKVlcXQoUMJCQkBwNfXl6ioKIxGo8M526s9h4v5\n+bdS5HPTOVqNwpG8Uu6ZOoRfj5ZI0iRaRKCfZ61aU3vNUh01TCs/O0B55dnFpquqLZwqrWJE/1Am\nnjc6S7QtOqlpctqZVSBamktJU1ZWFunp6UyePBmAgIAAampqnD7eaDQSFhZmfx0aGkpWVlatMl26\ndLEFp9NWcwLnAAAgAElEQVTh5+dHUVERRqORIUOGOBxrNBodjj169Cj79u1zKPfuu++ybt06Bg4c\nyGOPPUZAQIDzb7iVMxZVUVRmIti/8WpgZ6gqdHTTuVqjM+3hHbx0HC1w/rkVwlkl5dW8u2k/Hrqz\nzQiBvp5MH9fb/vr8xaQbU1FZw12Jgxy25Z0oZ9cvMnS9rdPrNJSdls8iZ9j6NLWx5jmdTofFYrFX\nCZ86dQqNC2+iriHt51cv11emsWMrKiq4//77+ctf/oKvr22CuBkzZnDvvfeiKApLlizhmWeeYdGi\nRfXGl5aWRlpaWp37IiIubMbt5lBlsuBVx7DWprrQWcVbO83vSZO3l47s/NIWjka0R4UlVfSJCGLc\niLO/S8vXZvHV7mP211XVjTfHNJZYeeq1VJvMTQ9UNAvFleUNgABfTzZ/n8OoQV2aKaL2o8bcBmua\nZs6cyZw5cygsLGTx4sV8+umnPPDAA04fHxYWRn5+vv210Wi0N62dW+b48eOEhYVhNpspKysjMDCw\nwWNramq4//77SUhIYPz48fYynTp1sv88ffp0/vznPzcYX3JyMsnJyXXuayjZailV1RY8PVo+824r\nPPU65t10OZ0DvSktNzV+gBAuslpVe3J+xo3X9cZYdNr++qbY3ucf1qC6EigvvY60j36msLQKgPLT\nNSSO6UVocAd8vD2aELloCeGhfvjK/1ejrFaVx/79NX+ZdWVLh+Ja0nTDDTcwYMAAtm/fjqqqLFu2\njF69ejl9/KBBg8jOziY3N5fQ0FAMBgMvvPCCQ5mYmBjS09MZOnQomzZtYuTIkSiKQkxMDA8//DDJ\nyckYjUays7MZPHgwqqry17/+laioqFoJT0FBgT2x2rJlC9HR0a683Vapxmzh6bTv6eBhRaMoeEjS\n5DStRuG6KyMByC+s4PV1ezhZUsm904YQ4OvZwtGJ9sCqqmjOqz0P8vdyeUh5YzUWXnot14/qYW+2\nO/DbKX7JKWL73uPcHNfXtaBFi/LUa1n7+UE6B3Vo8moH7Z2pxlY72+ZqmgB69erlUqLkcDGdjpSU\nFO68804sFgvTpk0jOjqaJUuWMHDgQMaNG0dSUhLz588nNjaWgIAAFi9eDEB0dDQTJkxg4sSJaLVa\nUlJS0Gq17Nixg/Xr19O7d28SExOBs1MLPPfcc+zfvx+Abt26kZqa2qS4WxNTjZXe4UGMHWwbobHq\ni5wWjqhtembuaLQaBcPXhzldZZakSbiF1U3rY51bu1RXAqXVapiTdLbvZp/IYKK6BfLOJ/su+Nri\n4vpjbB/KTpv4YPMvDOsTQpXJTMcAWdD3XFUmW9LUoYmrQbiTUxEMHTq0zqGtZ4a8/vjjj05f8Npr\nr601V9K8efPsP3t6evKvf/2rzmPvuece7rnnHodtw4cP58CBA3WWf+6555yOq60wW6ytYlbUtu5M\np3CtVoPZIqNXhHtYrLVrmi4WnVbBLEuxtChXO/mDbeRkxwBvvD11LF75IzqdhkdnDpfpJM5RXWPh\nrsSBDOzVqfHCzcyppGnXrl3NHYdwkm2l5/bVJGexqi027FanVSRpEm5TV5+mptAoisvnki/ZlrX3\n15N4XsCkpDdf35eKyhq+2JnLjwcKuKJvqBuja9s+35lLn1ayVl/L13UJl5y//k5YkJe96rKt6t8j\ngHc2ZxPZtYye3QLw9fZgSHRn+/4Ptx6sc/FLd9BpNbK0SguzWKzt5g8BtyZNqooGxaXai/LTNRzK\nLQYgIsxPZha/CAqKTvP+Zwcw1Vi5a3LTVy3w9NDi6aGlT2QwJeXVboyw7ausMreKWiYAlz6pHn30\nUUpLzw7VLikpYeHChW4PStTv/JWexwwJYfyItj1cdUivQJInRHHtsO4E+Xmy5XvHflqnyqqYOaFf\ns1xbmuda3ooNe/nrK9/w+c5ctu3Jq3N6kbbCorpnfSzN74tLu3ovxgztxtGCMnYfPEHmrmONHyAu\n2JFjJYy9IpxHbr3CLX0jtVpZI/N8Zmvr6ZbiUk3TgQMH8Pc/O0V8QEAA+/ZJx8OLqT33aQr/fTHe\n3b+c4IPNB9BoFNt6cc34AaKTpKnFdQrwZtI1UVRU1nD8ZAUr1u/laEE5D84YRqBf2+qgb3XToqIa\nja15zqri0oR+Z5p0qqrNpH30E0eOlwC2hWRvm9iP7iFNX/BaOCosqeS9TQeorDYza5Jra7A2RKfR\nSN+0OrSW5meXkiar1UpJSYl9Vu3i4mIslrbdNNTWWCztr0/T+aaP601FpW2W3E++PUJxM1ZVS5+m\nlme2WPH30dOtsy+9I4K4dlh3tu7IJXPXUW4Y07SRui3FalUdZgNvqr6RQbaRcAoMjHK9WcLLU8c9\n086Orvtq9zHKK2Xm6fzCCiqrzfTseuErQ+SdqODK/qFcNdC9Nf22mib5TGqtXEqa7rjjDm666Sau\nv/56AD799NNGJ4wUdbNYrJwqtSUDOq3i9Dwu7bmm6QwPncZewxDo58kvv/fRaA4Bvp6s+/JXduwr\ncNhusVoZc3l3+vUMbrZrC5sasxXdeX8IDO3dmUVv/kDC6KhW8xemMyxu6tM0on8YI/qHNV7QSVrN\npdPks2n7b+QaywDbfFa3ntO0v2n7b1RU1tSe306Frp18GDW4K8dOlKPXaegT2fDv/ra9x/nDsO5u\nj99W+31p/F+1RS4lTZMnT2bgwIFs374dgJdffvmSWBi3OXybdZxfcovoGODN3l9PMnf65Y02Rbxp\n+JmS8mrGXhEOXBozWnvqdVRWN99yEb0jglgwc3it7dU1Fpa8vwtvLx09uji/armw1bZkHTpBgK8n\nBadO19rft0ewQ9+PGrO1Vu1MkL8Xgy/rRJXJgrenjm17jrP38En7/uLSavx99Q5NVxqNwh0JA5rh\nHTnP6qY+Te6m1ShYLoEmH1VV2b73OE/cORKAbXvyeH39Hvt+jaJw7znzW51rw1e/smlbNhFd/Mnc\ndZSHb76iwY70qqrSuxlGdGl/788mznJ1eZrm5FLS9PLLLzu8/vTTTwGYO3eu+yJqZ37JKeL4yQoA\nLu/d2f5lUV1jIWZ4OD27BhAdHsh7m/Y3OtKlW2cfbo+3tZ3n5eU1b+CthJdeS2XVxV9jy9NDy/8N\n7kJVMyZs7dXJkkrWZx6mT2QQ/c+rqasyWVjywS4evW2EfXi2pZ5+QB28bAmzt6eOvYdP1lq09nwr\n1u9135toIneNnnM3bTOOEj1dVYO3pw5FUTDVWPjPhr3otBruSBhA3skKNn51mCsHhDG8X/MPoT9d\nZXZ4lkYN6sqoQV2dOvaG0WebgnOOl2IyWxv8TG6uxWNt/1ftP8Ftq1xKmjp06GD/ubq6mi+++IKo\nqCi3B9WeBPnZmt1y8svIOniS0UNt0+Tbmtlsv3QDojoyIKpji8XYmg3rE0Kv7oEtcm29Tkt1jfTZ\nc1VltZmRA7sQNzKyzv0dPHX8e/Vu+vawJVS5BWV1NsH169GR1Rm/UFphwtqK//L+fGeubZi/AsVl\n1dx6ffOM9LwQzVV7oaoqT7/xPRFhfrbpO6xW/jAsnPxTFby48kciwvyYPq43G7769aIkTSazhT6R\nF177o/fQ2pbuaGBduKZMZOkMWz/L1ve8nxnJ2RzN5bnGMj76+jAeOluS2q9nMFcPdi7Zvdhc7tN0\nrtmzZ9eaoVs46hzkTecgbzx0Gn49WmLffm7SJOrn5anDq4WmzvfQaVps0s22rLLa3OByBwN7daJb\niK+9BvHKevru9OsZTL+ewU5/2Tf0JVZjtpJfWEHnQG+3P0+/Hi3hzsSBrbrvlVarNPlZLiyp5E3D\nz/j56Gvts1pVxl8VwR+uCHfY3q9n8O/dCGyOnShn645cYoaHn38KtzLVWC9ogskzPDw09vXO6tKc\n02J4e+rY8+tJvL10dArwcmvfNleVVph4f/MBNIrCwdwipo/r3SzJb3ZeKTHDw+39yN76+Gf2HTll\n33+6qvUMYrigT4/Kykpyc3PdFUu7ZvsCPvtLaLaokjS1cnoPrUwy5yKrVeW9T/dz1+SGm9KC/LwI\ncnL0uzv6CG3bk8ehoyV467XMcHJB25Ly6lpzHa3L/JX//DW2VtnWnDCBrSnJbGnaF8+xE+VcPbjr\nBY0Sm3fTUN7+ZB9jr+jerPfKVGPBww1Jk+eZmqZzHMkrYdeBExw+VkKnQC90zdQ818HLg78mX8nh\nYyVs/PrwRU+aLFaVUyVVAKzP/JWJ/9eD7iF+HDtRzhc7j9I50JtuIb51fn/t/qUA46lKh1rm/MIK\nXl2bxZN3jar3mqera+jgdbbv6G0T3TeFg7u5lDQlJCTYf7ZarZw6dYp7773X7UG1R3oPLaZz/tIz\nW6zodK37g/ZS59fBg7WfH+KnI4UO21UVrhsRQVS3Cx+27E5V1Wa++ymfay7v1miikXeynNLy2oMJ\nzBYru3854TDi6FzpXxwi9soIfDvUrnUor6yhtLyaAb062ufcai3KK2sYf1UEy9Zkwab99u2hHTvU\nO6T/h31GOgZ4Ofw/Zx8vrbNsa9fU5rmVm/ZTWmFi6tjoC7p+gK8nXTv58vr6vXWuzWdVVWaM74Pf\n789V2WkTq7b8gqrCZeGBDY5Ss81nZXtvldVmPM8fGdcEnQM7cOC3IiLCzn6Rf/bdb4wbEUHimCie\ne2cnkWHN94zrtBp6RwQR6Ovpln5yRWVVVFTWsGOfkeIy2x+Ck6+9zGHw0Y8HCsg6eIL8wtOEh/rh\n4+1BRJiffW6v0OAOdOnUgaxDJ9nw1WG866ixtVitmC0qI8vD0Hto8fbU8dXuYxSXV7Puy1/RKHDl\ngDD8ffT88LORkgpbLHt/Lbwozbfu4FLS9Oqrr549UKejY8eO6HSuVVZlZmby97//HavVyvTp07n7\n7rsd9ptMJhYsWMBPP/1EYGAgixcvpnt32y/M8uXLWbNmDRqNhscff5zRo0c3eM7c3FweeughSkpK\n6N+/P88++yx6fe0P+4tBr7P95WIbXq1I81wbEBHmz99mX1Vr+7ET5WQdPOH2pGnrjhwOH2v4Szkk\nyLvePl4Hc4v4/icjfSKDCOvoU2v/R18fto9E/PVoCddcXnefgZMllaza8gudg7wdmli+/t8xtvyQ\nw2/5pfh666kxWwgJ6oDx1GkUxVbbEujnyahBLTNDfWWVmV0HHKeOUBTb/+Op0ip8vDyIvSqCEf1C\n7bUdPx8pZPve43WeT1EUrugX6tDcU9cXRVug1TYtaaqoMvOnqYPdEsPUsfWPtN6+9zirtvyC7+99\niPZlnyL2ykiuHtKVl1fvpmsnH0KDO5CdV/v3Y+vOXILO+fKPvaruvnSu6NsjiOfeOUzF781C/XoE\no9EoXPb77959N17ulvm4GnOmxuvcJuXTVTWNNrV28NLZ+wcBvLHhJ/pHdaRvj2B6dQsk69AJDuYW\nOdRi7dxnZNYk2+jTut6bTqshZnhEozF/ui2bTdt/45ecInp1C2D/b0U8d99oikqrqa6x8PmOXKwq\nRHXz56oBts+KUQO70jHAu9FztwZOfQKcP2ruXIqiMGfOHKcuZrFYSE1NJS0tjdDQUJKSkoiJiXGY\ntmD16tX4+/uzefNmDAYDzz//PC+99BKHDh3CYDBgMBgwGo0kJyezadMmgHrP+fzzzzNr1izi4+NJ\nSUlhzZo13HzzzU7F6m4+3jpOFlfyr1W78NBq2P9bEYltbOI+YeOl17p9GgSLVWXHvoI6pz84189H\nCu1/KZ4vJKgD46+K4Nus47WaFS1WlQ5eOib+X08Abhijq7fvx1UDunC6qoZPt2c7jEjTahT+9dAf\n7JOrlp82sfGrw9wW39/+ZdeSEkZHUfh7s8IZ1TUWtnyfg4+XB/4+eockEHC56aNHFz9WrN9LUWkV\nIcEdKCg63epq1eqi1Wg4WVzJ0YIy+7Zgfy/MFhX/8/oqVVabeX3dHjp4eVy0UVzD+oQQ6OdJt86+\ngK0P2plakEnXRPH5jlxOlVYxemi3WoPP46/u6fah/x28PLj/psvtIw7/s2EvunMSCZ+L9Lx76m2D\nUc5Nmhav/JHo8Ibf75G8Eq65vJv9dQcvHRNG9bC/7tk1gPc27ed/B89O4+Hr7eGWRPD6369TY7Zi\nqrEwRavgodMSEmwbSOZs83hr5VTSdO6ouTMqKytZs2YNxcXFTidNWVlZREZGEh5u++CKj48nIyPD\nIWnaunWrfQqDuLg4UlNTUVWVjIwM4uPj0ev1hIeHExkZSVZWFkCd5+zVqxfbt2/nhRdeAGDKlCm8\n/PLLLZY0eei09hl6LRYrX+0+5pYOi+Li8/bUUXa6hiqTGZ1W41SNocWqOvRpO19JuYmeXRufD6p/\nz4ZHWVosVha9+QN3TR5EpwDHCVM1GsWp/iQeOg0Bvp7cdF2fBsv5dtC3qg/Anl0D3DLTc0Ouu9JW\ni2GxWG0L6ipKm5ihv2tnH37J8WDPr783Nasqx05UUFJRzeBeney1M0VlVSxP38PE/+vB4Ms6N3BG\n99J7aOl7zmSS5zYP9+jiT3ILzL91ZuQzwM1xfVtkKonILv6s/OwAOq2GKpNtJvOwjj7ceF3vBo/L\nyS+l7PTZPmxJMY7lg/29mDv98maJ+QwPneai1MZdbE4lTeeOmisvL+ett97iww8/ZOLEibVG1DXE\naDQSFnb2L7vQ0FB74nNumS5dbFV2Op0OPz8/ioqKMBqNDBkyxOFYo9EIUOc5i4qK8Pf3tzcfhoWF\n2cu3NK1WU2u0iWg7vPQ6vPRaVmccJL+wwuHDtT6nSqvo0ql2k9m5Rp/zl2FTabUaHr+jdpOicC+t\nVkNb+pPHS6+z1wCcy2Kx8sraLLLzbc1ephorN13Xu9mTz7ampWoTr+wfZh9dWl5ZQ2lFNYG+jc9C\nfm5fLOFeTjfQFxcXk5aWxsaNG5kyZQrp6en2NeicVdcwzfP/8q2vTH3brXXMclvfX9ON/ZWdlpZG\nWlpanfsiIhpvyxWXBo1G4abYhmthhGgLtFpNs9c4CPfw9fZoFc3glzqnkqZ//vOfbN68mRtvvJGN\nGzfi49PwX8z1CQsLIz8/3/7aaDQSEhJSq8zx48cJCwvDbDZTVlZGYGBgg8fWtT0oKIjS0lLMZjM6\nnY78/Pxa1zpfcnIyycnJde5btGiRy+9XCCGEEO2HUw2OaWlpFBQU8MorrzB69GiGDRvGsGHDGDp0\nKMOGDXP6YoMGDSI7O5vc3FxMJhMGg4GYmBiHMjExMaSnpwOwadMmRo4ciaIoxMTEYDAYMJlM5Obm\nkp2dzeDBg+s9p6IoXHXVVfbO4unp6bWuJYQQQgjhLKdqmvbv3994IWcuptORkpLCnXfeicViYdq0\naURHR7NkyRIGDhzIuHHjSEpKYv78+cTGxhIQEMDixYsBiI6OZsKECUycOBGtVktKSgpara1XQV3n\nBJg/fz4PPvggL730Ev369WP69OlueR9CCCGEuPQoanPOB9+OLFq0iIULF7Z0GHbtccHerl1b51pD\nQgghBFzgMiqXktb2hd7a4hFCCCHaO6lpEkIIIYRwQvubeUoIIYQQohlI0iSEEEII4QRJmoQQQggh\nnCBJkxBCCCGEE2T03AW699577YsFC/f76aefGDDg4i/WeamQ+9v85B43L7m/za+93uOuXbty++23\nu3SMjJ67QGPGjCEzM7Olw2i35P42L7m/zU/ucfOS+9v85B6fJc1zQgghhBBOkKRJCCGEEMIJkjQJ\nIYQQQjhBkiYhhBBCCCdI0iSEEEII4QTtk08++WRLB9HWDR06tKVDaNfk/jYvub/NT+5x85L72/zk\nHtvIlANCCCGEEE6Q5jkhhBBCCCdI0iSEEEII4QRJmoQQQgghnCBJkxBCCCGEEyRpaqLMzEzi4uKI\njY3ltddea+lw2qyYmBgSEhJITExk6tSpABQXF5OcnMz48eNJTk6mpKQEAFVVefrpp4mNjSUhIYGf\nfvqpJUNvtRYuXMioUaOYNGmSfVtT7ml6ejrjx49n/PjxpKenX/T30VrVdX+XLl3K6NGjSUxMJDEx\nkS+//NK+b/ny5cTGxhIXF8dXX31l3y6fIXU7fvw4M2fOZMKECcTHx/Pmm28C8gy7U333WJ5jJ6jC\nZWazWR03bpyak5OjVldXqwkJCerBgwdbOqw2aezYsWphYaHDtn/+85/q8uXLVVVV1eXLl6vPPvus\nqqqq+sUXX6izZ89WrVarumvXLjUpKemix9sWfP/99+revXvV+Ph4+zZX72lRUZEaExOjFhUVqcXF\nxWpMTIxaXFx88d9MK1TX/f3Xv/6lrlixolbZgwcPqgkJCWp1dbWak5Ojjhs3TjWbzfIZ0gCj0aju\n3btXVVVVLSsrU8ePH68ePHhQnmE3qu8ey3PcOKlpaoKsrCwiIyMJDw9Hr9cTHx9PRkZGS4fVbmRk\nZDB58mQAJk+ezJYtWxy2K4rC5ZdfTmlpKQUFBS0Zaqs0YsQIAgICHLa5ek+//vprrr76agIDAwkI\nCODqq692+OvyUlbX/a1PRkYG8fHx6PV6wsPDiYyMJCsrSz5DGhASEsKAAQMA8PX1JSoqCqPRKM+w\nG9V3j+sjz/FZkjQ1gdFoJCwszP46NDS0wQdONGz27NlMnTqVDz74AIDCwkJCQkIA2y/3qVOngNr3\nPSwsTO67k1y9p/KMu+7dd98lISGBhQsX2puO6ruPcn+dc/ToUfbt28eQIUPkGW4m595jkOe4MZI0\nNYFax3ygiqK0QCRt38qVK0lPT+f111/n3Xff5Ycffqi3rNx396vvnsq9ds2MGTPYvHkz69evJyQk\nhGeeeQaQ+3shKioquP/++/nLX/6Cr69vveXkHjfd+fdYnuPGSdLUBGFhYeTn59tfG41G+19AwjWh\noaEAdOzYkdjYWLKysujYsaO92a2goIDg4GCg9n3Pz8+X++4kV++pPOOu6dSpE1qtFo1Gw/Tp09mz\nZw9Q/2eF3N+G1dTUcP/995OQkMD48eMBeYbdra57LM9x4yRpaoJBgwaRnZ1Nbm4uJpMJg8FATExM\nS4fV5pw+fZry8nL7z9988w3R0dHExMSwbt06ANatW8e4ceMA7NtVVWX37t34+fm1+19Qd3H1nl5z\nzTV8/fXXlJSUUFJSwtdff80111zTkm+hVTu3b92WLVuIjo4GbPfXYDBgMpnIzc0lOzubwYMHy2dI\nA1RV5a9//StRUVEkJyfbt8sz7D713WN5jhuna+kA2iKdTkdKSgp33nknFouFadOm2R8u4bzCwkLm\nzJkDgMViYdKkSYwZM4ZBgwbxwAMPsGbNGrp06cKSJUsAuPbaa/nyyy+JjY3F29ubf/zjHy0Zfqv1\n0EMP8f3331NUVMSYMWO47777uPvuu126p4GBgdx7770kJSUBMGfOHAIDA1vsPbUmdd3f77//nv37\n9wPQrVs3UlNTAYiOjmbChAlMnDgRrVZLSkoKWq0WQD5D6rFz507Wr19P7969SUxMBGz3XJ5h96nv\nHn/00UfyHDdCFuwVQgghhHCCNM8JIYQQQjhBkiYhhBBCCCdI0iSEEEII4QRJmoQQQgghnCBJkxBC\nCCGEEyRpEkIIIYRwgiRNQlzC+vXrR2JiIpMmTeL++++nsrKy3rJHjx5l8ODBJCYmMnHiRBYsWEBN\nTQ0A3333HVdccQWJiYkkJCQwa9YsCgsLAVi7di0jR460H7dq1SoANmzYQEJCAgkJCfzxj3+0zw9z\n9OhRJk2a1MzvHO666y5KS0ub/TrOXnfp0qX85z//abbr/vGPfwRs93fjxo327Xv27OHpp59utusK\n0Z5I0iTEJczLy4v169fz0Ucf4eHhwfvvv99g+YiICNavX8/GjRvJz8/nk08+se8bPny4fd+gQYN4\n99137fsmTpzI+vXrefvtt3nxxRc5efIk3bt355133mHjxo3cc889/O1vf2u291mX119/HX9//4t6\nzZa87pn/22PHjvHRRx/Ztw8aNIjHH3/8oscjRFskM4ILIQBb0nPgwAGee+45unbtyi233ALYakB8\nfHzs61MBaLVaBg8eXOeK5qqqUlFRQURERK19HTt2JCIigry8PIYNG2bffvnllzusYWU2m3n00Uf5\n+eef6dmzJ//85z/x9vZ2ONeqVav44IMPqKmpITIykmeffRZvb28ee+wxPD09OXz4MHl5eSxatIj0\n9HR2797NkCFD7IuQxsTEsGbNGk6fPs1dd93FFVdcwa5duwgNDWXZsmV4eXmxb98+nnjiCSorK4mI\niOAf//gHAQEBDnFs3bqVV155hZqaGgIDA3n++efp1KkTFRUVPP300+zduxeAuXPnEhcXZ79ucHAw\nr7zyCuvWraNLly4EBwczYMCAWvfsscceQ6/Xc+jQIQoLC3nssccYO3Ys1dXVPPnkk+zduxetVstj\njz3GyJEjOXjwIAsXLqSmpgar1crSpUvp0aMHQ4cOZdeuXbzwwgv8+uuvJCYmMmXKFPr168cbb7zB\n8uXLKS4u5i9/+Qu5ubl4e3uTmppK3759Wbp0KXl5eRw9epS8vDxuv/12brvttkafKSHaG6lpEkJg\nNpvJzMykd+/exMfHO9QgffLJJ1x//fUO5aurq/nf//7H6NGj7dt27NhBYmIif/jDH/j222/ty1ec\nKzc3l9zc3FoJ1Zo1axgzZoz99ZEjR7jxxhvZuHEjPj4+vPfee7XOFRsby4cffsiGDRuIiopizZo1\n9n2lpaW89dZbLFy4kD//+c/MmjULg8HAL7/8wr59+2qd67fffuOWW27BYDDg5+fHpk2bAFiwYAGP\nPPIIGzdupHfv3rz88su1jr3iiitYtWoV69atIz4+nhUrVgCwbNkyfH192bhxIxs3bmTkyJEOx+3d\nu5ePP/6YdevW8fLLL9sXR63LsWPHeOedd1i+fDlPPPEE1dXV9pq8jRs38sILL/DYY49RXV3N+++/\nz2233cb69ev58MMPCQsLczjXww8/bK8VnDVrlsO+pUuX0r9/fzZu3MiDDz7Io48+6vB/8p///IfV\nqxddWG4AAASbSURBVFfz73//2940K8SlRGqahLiEVVVV2deeGj58OElJSej1egoLCzEajRQVFeHv\n70/Xrl05evQoOTk5JCYm8ttvvxEXF0ffvn3t5xo+fDjLly8H4LXXXuPZZ5+1r1318ccfs3PnTvR6\nPampqQ5rgG3fvp01a9Y4JEZdunThiiuuAOCGG27g7bffZvbs2Q6xHzx4kJdeeomysjIqKiocFmMd\nO3YsiqLQp08fOnXqRJ8+fQC47LLLOHbsGP369XM4V/fu3e3bBgwYwLFjxygrK6OsrIwrr7wSgClT\npjBv3rxa9zA/P58HH3yQEydOYDKZ6N69OwDbtm3jxRdftJc7v4Zqx44dXHfddfYatIYWOp0wYQIa\njYYePXoQHh7O4cOH2blzJ7feeisAvXr1omvXrhw5coTLL7+cV199lfz8fMaPH0+PHj3qPe/5du7c\nydKlSwEYNWoUxcXFlJWVAbY13vR6PcHBwQQHB1NYWFgrIROivZOkSYhL2Jk+TeeLi4tj06ZNnDx5\nkvj4ePv2M32aCgoKmDlzJhkZGfbV5s81btw47rvvPvvriRMnkpKSUqvc/v37efzxx3n99dcJCgqy\nb1cUxaHc+a/B1my1bNky+vbty9q1a/n+++/t+/R6vf24Mz8DaDQazGZzrXOdW0ar1VJdXV2rTH2e\nfvppZs2axbhx4/juu+/stVGqqtYZd2Pvy5lyiqJQ37KhCQkJDBkyhC+++ILZs2fz9NNPM2rUKKeu\nU9c5z1z7/HtU130Uor2T5jkhRC3x8fF8/PHHbNq0ibi4uFr7Q0JCeOSRR3jttdfqPH7nzp119mk6\nV15eHvfddx/PPvssPXv2rLVv165dABgMBnut07kqKiro3LkzNTU1DqPB3MXPzw9/f3927NgBwPr1\n6xkxYkStcmVlZYSGhgKwbt06+/arr76ad955x/66pKTE4bgRI0awefNmqqqqKC8v5/PPP683lk8/\n/RSr1UpOTg65ubn07NmTESNG2N/3kSNHOH78OFFRUeTm5hIeHs5tt91GTEwMBw4ccDiXj48PFRUV\ndV5nxIgRbNiwAbCNiAwKCsLX17feuIS41EhNkxCilujoaCoqKggJCSEkJKTOMtdddx1Lly61JxVn\n+jSpqoqfn1+jw9j//e9/U1xczFNPPQXYai/Wrl0L2Jqb0tPTSUlJoUePHsyYMaPW8fPmzWP69Ol0\n69aN3r1715sIXIh//vOf9o7g4eHhLFq0qFaZuXPnMm/ePEJDQxkyZAhHjx4F4J577iE1NZVJkyah\n0WiYO3euQ2f6AQMGMHHiRBITE+nWrVudieEZPXv25NZbb6WwsJCnnnoKT09Pbr75Zp544gkSEhLQ\narUsWrQIvV7Pxx9/zIYNG9DpdHTq1Ik5c+Y4nKtPnz5otVpuuOEGpk6d6tBUOXfuXBYuXEhCQgLe\n3t72TvNCCBtFra+OVwghRIt77LHH+MMf/lCrM74Q4uKT5jkhhBBCCCdITZMQQgghhBOkpkkIIYQQ\nwgmSNAkhhBBCOEGSJiGEEEL8f7t1IAAAAAAgyN96hQGKIgZpAgAYpAkAYJAmAIAh6oc93qnADCIA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc6469ffd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plot_diversity()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_pi_d(vcf_data=vcf_data, samples_to_use=vcf_samples_in_148,\n",
    "              size=300, start=54438, stop=63146 , step=3, ax=None, filename=None):\n",
    "    fig, ax = plt.subplots(2, figsize=(8, 4))\n",
    "    \n",
    "    _, x, pi = plot_diversity(vcf_data, samples_to_use, size, start, stop , step, despine_bottom=True, ax=ax[0])\n",
    "    _, x, D = plot_tajd(vcf_data, samples_to_use, size, start, stop , step, despine_bottom=False, ax=ax[1])\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    \n",
    "    if filename is not None:\n",
    "        fig.savefig(filename)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create final plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nfs/users/nfs_r/rp7/miniconda3/envs/biipy/lib/python3.5/site-packages/allel/stats/diversity.py:1040: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  wD = d / d_stdev\n",
      "/nfs/users/nfs_r/rp7/miniconda3/envs/biipy/lib/python3.5/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['arial'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAAE5CAYAAAADe11uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAIABJREFUeJzs3Xl8VPW9P/7XOWeW7AlJSEJICGRhC/uiuAAaQBSKoGzl\nWpeotdeKdruivfZHLY/6oPfWit5avrXapkUrVUREGgURRLCKls2ERSSQQCDJZN8z2znn98fAkCHb\nTDJnZjJ5PR8PHw/mzFneE2d5n/dnE1RVVUFEREREQUP0dwBERERE5F1M8IiIiIiCDBM8IiIioiDD\nBI+IiIgoyDDBIyIiIgoyTPCIiIiIggwTPCIiIqIgwwSPiIiIKMgwwSMiIiIKMkzwiIiIiIIMEzwi\nIiKiIMMEj4iIiCjI6HpzUGtrK4xGIyRJ8nY8fve3v/0N999/v7/D6KCsrMzfIfhMcnKyv0MgIiLq\n19xK8BRFQX5+Pnbs2IHCwkIYDAZYrVbExsZi1qxZWLlyJYYPH65xqL4xkBIpIiIiCk5uJXj33Xcf\nbrjhBvz0pz/FyJEjIYqOlt36+np8+eWX+N3vfoe5c+di8eLFmgZLRERERD1zK8HLy8uDXq9HS0uL\nM7kDgJiYGMyfPx/z58+HzWbTLEgiIiIicp9bgyz0ej0A4K677kJRURHsdrvzueLiYpd9iIiIiMi/\nPBpk0dDQgF/96lcoLi7GoEGDkJ6ejpMnT2L37t1axUf9zN92FiM0REJjiw0PL8zwdzhEREQDkkcJ\n3pAhQ/D6668DACoqKnD27FkYjUZNAqP+w2qTUVlvAQCEGCWsuGUY3t53wc9RERERDVweJXhNTU04\ncuQIsrKykJSUhKSkJK3ion6gtLIV+7+uhKICI5LCIYjAzePi/R0WERHRgOdRgtfa2oq//OUvOHPm\nDGw2G9LT05GVlYWnnnpKq/gogBWcrcft1w9BXBSruERERIHErQRPVVUIgoC33noLw4YNAwBYLBYU\nFRXhzJkzLvvQwCErKsKMvZorm4iIiDTk1ija++67D6+//jp0uqs/5kajEVlZWUhMTMRTTz2Fbdu2\naRYkBSZZUSFxsTsiIqKA41b55bXXXsM777yDn/70p7h48SKioqJgNpuhqipuuukmPPDAAxgzZozW\nsVKAURQVksiqLRERUaBxK8EzGo245557cM8998Bms6Gurg4hISGIiorSOj4KYLKqQugmwWOzPRER\nkX943IFKr9cjISFBi1iov1EBsYsEThQElNeasf/rSoxMicSUkbE+Do6IiGjgYg8q0oQkAsVlzZgx\nNg5FZc3+DoeIiGhAYYJHmhBFAa0WGSF6yd+hEBERDThM8EgToiCgxWyH0cAEj4iIyNc8SvDeeOMN\nNDQ0aBULBZH4aCMuVbUhPEQCVH9HQ0RENLB4NMiiqqoKy5Ytw9ixY7F06VLMnDmToySpUzdkx+OG\nbC5bRkRE5A8eVfB+8pOf4KOPPsKyZcuwbds23HbbbXjhhRdw4QIXlh+Q3M3teQ9ARETkUx73wRME\nAYMHD0Z8fDwkSUJDQwOeeOIJ/O///q8W8RERERGRhzxqot20aRPee+89DBo0CMuWLcOaNWug1+uh\nKApuu+02rFmzRqs4iYiIiMhNHiV4JpMJv//97zF06FDntt/+9rd48skn8corr3g9OApwHDxBREQU\nkDxqov38889dkjsA2L9/PwAgIyPDe1FRcGEiSERE5FNuVfDefPNNbN68GaWlpVi0aJFze0tLC6ZM\nmaJZcBQYFFWFoqjQSSL2HjFBEIDBMUYoKjM3IiKiQORWgrdo0SLMmjULL7zwAn72s585t4eHhyMm\nJkaz4Cgw7PyqHIVn6/HTFaNRXtOG2ZMScKCgCvXNVreOr6q3oKnVhsgwvcaREhEREeBmE21kZCRS\nUlLwwgsvYOjQoc7/epPc7d+/H/Pnz8e8efPwpz/9qcPzVqsVP/7xjzFv3jwsX74cFy9eBADU1dXh\n3nvvxeTJk7Fu3TqXY44fP45FixZh3rx5+PWvfw2VlSWvslgVZCRHQFZU6PUiUgaHITk+FHbZvb/z\n1FGDUNfkXjJIREREfedWgrdq1SoAwOTJkzFlyhRMmTIFkydPdj52lyzLWLduHV577TXk5+fjn//8\nJ4qKilz22bJlC6KiorB792488MADeP755wEARqMRP/rRjzodqfvss89i3bp1+Oijj1BSUuLsF0je\noSgqdDoRdlmFdHlia50kuJ3gCYIAhTk3ERGRz7jVRLt582YAwNGjR/t0sYKCAqSlpSE1NRUAsHDh\nQuzZsweZmZnOffbu3YvVq1cDAObPn49169ZBVVWEhYVh2rRpHSZVrqysRHNzMyZPngwAWLJkCfbs\n2YPZs2f3KVa6SlZVhEgSbHYFonglwRNhlxW3jhcAVlWJiIh8yKNRtB9++CGam5sBABs3bsTq1atx\n8uRJt483mUxISkpyPk5MTITJZOqwz5AhQwAAOp0OkZGRqKurc/ucSUlJHc5JfaQCkiTALiuQxN5U\n8ADmd0RERL7j0Tx4GzduxB133IFDhw7hs88+w4MPPohf/vKX2LJli1vHd1bFuXYtW3f26cv+AJCX\nl4e8vLxOnxs2bFi3xw5UkijAZlchXr4l8KSCJwoCK3hEREQ+5FGCJ0kSAODTTz/FqlWrMHfuXLz8\n8stuH5+UlISKigrnY5PJhISEhA77lJeXIykpCXa7HU1NTd0O5rj2nBUVFR3Oea3c3Fzk5uZ2+tz6\n9evdeSkDjigCNlmBeDl51kuOhM8drOARERH5lkdNtImJiVi7di0+/PBDzJ49G1arFYriXhUHAMaP\nH4+SkhKUlpbCarUiPz8fOTk5Lvvk5ORg27ZtAIBdu3ZhxowZ3VbkEhISEB4ejmPHjkFVVbz33nuY\nM2eOJy+L3CAKgksfvPTkCDy4YIRbxwqCAPffJURERNRXHlXwXnzxRRw4cAAPPvggoqKiUFlZ6dH6\nszqdDmvXrsXDDz8MWZaxdOlSZGVl4aWXXsK4ceMwZ84cLFu2DE8++STmzZuH6OhobNiwwXl8Tk4O\nmpubYbPZ8PHHH+Mvf/kLMjMz8eyzz+LnP/85zGYzZs2ahVmzZnnysqgngqOJ1i6rzj54IQYJw5Mi\n3Dtc4CALIiIiXxJUN395VVVFRUWFcwBEsFq/fj1+/vOf+zuMDsrKyvx27bf3XUByXCiiw/WorDdj\nzpSkng9q5+uiOoQYJIwaFuXW/snJyb0Jk4iIiC5zu4lWEAQ89thjWsZCAUwSBZc+eJ7gPHhERES+\n5VEfvIkTJ6KgoECrWCiAiaIAu/1qE60n2ERLRETkWx71wfvyyy/xj3/8A0OHDkVoaKhz+44dO7we\nGPnX10V1OF3aBACobrBgbFoU9h2rxA1j4z0+lyAIYHpHRETkOx4leK+++qpWcVCAOV3ahBW3us4J\nOG6E52sPA4AoAAoreERERD7jURNtcnIyDh06hG3btmHo0KEQBAE1NTVaxUZBgvPgERER+ZZHCd6z\nzz6LY8eOIT8/HwAQHh6OX/3qV5oERsGDgyyIiIh8y6MEr6CgAL/85S9hNBoBANHR0bDZbJoERsGD\ngyyIiIh8y6MET6fTQZZl58oStbW1EEWPTkEDkCAIbKIlIiLyIY8GWdx777147LHHUFNTgw0bNmDn\nzp348Y9/rFVs5E+ez4bSJRGs4BEREfmSRwnenXfeiezsbBw8eBCqqmLjxo3IyMjQKjYKEhxkQURE\n5FseJXh//etfcfvtt+Oee+7RKh4KQpwHj4iIyLc8SvCam5vx0EMPITo6GgsXLsT8+fMRH+/5xLc0\nsHCQBRERkW95NEJi9erVyM/Px9q1a1FZWYnvfe97eOCBBzQKjYIFB1kQERH5Vq+GwMbFxSE+Ph4x\nMTGc6Jh6JHAlCyIiIp/yqIn2zTffxIcffoiamhrcfvvt+PWvf43MzEytYiN/8mI+JoCDLIiIiHzJ\nowSvvLwczzzzDEaPHq1VPBSE2ERLRETkW24leKtWrcLmzZvxxhtv4O9//7tzu6qqEAQBR44c0SxA\n6v/CQiQUnqvHxepW2OwKbpmUgKHxYf4Oi4iIKGi5leBt3rwZAHD06FFNg6EA4sWJjuOijPjBnY6m\n/G8vNqGqzsIEj4iISENcZ4x8ShIAme21REREmnKrgpeXl9ft87m5uV4JhoKfKApQFCZ4REREWnIr\nwWtpaQEAFBcXo7CwEDk5OQCATz75BNOmTdMuOgo6kihAZoJHRESkKbcSvNWrVwMAHnzwQbz77ruI\niIhwbv/Rj36kXXQUdBwVPH9HQUREFNw86oNXVlYGg8HgfGwwGHDp0iWvB0XBi020RERE2vNoHrzF\nixdj2bJlmDdvHgRBwO7du3HXXXdpFRsFIUkQOMiCiIhIYx4leI8++ihmzZqFQ4cOAQDWr1+PsWPH\nahIYed/RM3UorWxFzuQERITp/RKDKIIVPCIiIo15lOABQHZ2NrKzs7WIhTR2urQRI4ZEoLzWjKye\nEjyNcjAOsiAiomDV3GbDpvyT0Ouv9oCrrG1FVuog5+MhceGYOXmo5rF4nOBR/yUKAkINEuyy/xIs\n9sEjIqJg1dBswYjkKNxx4wjnNptdQavZ5ny8+aPTTPDI+3SSmxU0L65k0R5H0RIRUbCSZQWi6Dp+\nVa8TER1hdD6WJI1+YK/h0ShaVVWxfft2vPzyywAco2oLCgo0CYy0IUkCZLn7DKuhxYaaBos21xcE\n1DRaUFzejMo6sybXICIi8gdZUSGJ3SdwAgS89fFpvPXxaRSXNWgWi0cJ3rPPPotjx44hPz8fABAe\nHo5f/epXmgRG2tCJAuw9VPBqGy2YNipWk+uHhkgYMzwK1Q0WvP85p9ghIqLgIctqjxW6+xeOwe0z\nhuP2GcORkhCpWSweNdEWFBRg27ZtWLJkCQAgOjoaNputh6Nc7d+/H8899xwURcHy5cvxyCOPuDxv\ntVqxZs0anDhxAjExMdiwYQNSUlIAAK+88greeecdiKKIX/ziF5g5cyYAICcnB+Hh4RBFEZIk4d13\n3/UopmDVZrHDbFWgqioqas0wW2VIkgi5hz54iqJC7OEOpLdEQcCULEfyWFzRosk1iIiI/EFWlB4r\neHqdhOgISfNYPErwdDodZFmGIDiCr62t7dDW3B1ZlrFu3Trk5eUhMTERy5YtQ05ODjIzM537bNmy\nBVFRUdi9ezfy8/Px/PPP48UXX0RRURHy8/ORn58Pk8mE3Nxc7Nq1C5Lk+CP97W9/Q2ysNlWn/urN\nj89jWGI4AEcfgNOlTciZktjjIAtZUSEJvukjQEREFCxkRYUkedQ4qhmPorj33nvx2GOPobq6Ghs2\nbMCqVavwgx/8wO3jCwoKkJaWhtTUVBgMBixcuBB79uxx2Wfv3r3OyZPnz5+PL774AqqqYs+ePVi4\ncCEMBgNSU1ORlpbG/n89iAzXY960JMybloRbJiXg/7sv+3ITbfd98GQFmlXwiIiIgpU7ffB8xaMK\n3p133ons7GwcPHgQALBx40ZkZGS4fbzJZEJSUpLzcWJiYockzWQyYciQIY7gdDpERkairq4OJpMJ\nEydOdDnWZDI5Hz/00EMQBAErV67EypUrPXlZA0aIQUJ0hAEVNWa8ve8CTLVmPLIoA0a9a6lYUXvu\nQ0BERESuFLmfJXh5eXmdbt+/fz/279+P3Nxcty6mdrJElXBNU2BX+3R37ObNm5GYmIiamhrk5uYi\nPT0d06dP7zKOvLy8Ll/TsGHDun0N/V2oUcKKWx2v8bOCKtQ0WJEcH+qyj6Ko0Pm4xFxc3owhcaEI\nMWjfL4GIiEgLjgpeYDTRupXgtbQ4OsMXFxejsLAQOTk5AIBPPvkE06ZNc/tiSUlJqKiocD42mUxI\nSEjosE95eTmSkpJgt9vR1NSEmJiYbo9NTEwEAMTFxWHevHkoKCjoNsHLzc3tMildv369268n4PUw\n3V1oiIRWi73DdllRYfTBSmbNrXb8dWcxwkIk7DlswrMPjMOQuNCeDyQiIgpAsqLAoAuMQoVbaebq\n1auxevVq1NXV4d1338XTTz+Np59+Glu3bnVJunoyfvx4lJSUoLS0FFarFfn5+c5k8YqcnBxs27YN\nALBr1y7MmDEDgiAgJycH+fn5sFqtKC0tRUlJCSZMmIDW1lY0NzcDAFpbW/Gvf/0LWVlZbsc0kIUZ\nJVTUmlHT6DrnnZajaNtrbLEhKyUCK24ZhnvmpfVpbuV3PzmDA0c57QoREfmPrKgQA6SLk0d98MrK\nymAwGJyPDQYDLl1y/0dVp9Nh7dq1ePjhhyHLMpYuXYqsrCy89NJLGDduHObMmYNly5bhySefxLx5\n8xAdHY0NGzYAALKysnDHHXdgwYIFkCQJa9euhSRJqKmpwWOPPQbAMUr3O9/5DmbNmuXJyxqw0hLD\nUVFrxlt7L+C+24Yj4vL6tL7qJGqTFegvNwVLQt/WqK1pNKOm0dzl8i8NzZYOS6RFhOmhD5A7LSIi\n6v8Ki6pxy9QUf4cBwMMEb/HixVi2bBnmzZsHQRCwe/du54hXd82ePRuzZ8922fajH/3I+W+j0Yj/\n+7//6/TYRx99FI8++qjLttTUVLz//vsexTBg9JCjhYXocOvkRFjtikty5bMEz65Ap3NcRxT7luB1\np7bRjJe3HEP2iDjnNkVVYbHK+N4dYzS5JhERDTwqgKzUQf4OA4CHCd6jjz6KmTNn4vDhwwAc/dXG\njh2rSWDkO+I11TPFR9Ok2GT1agVPFKB0MpDGG9osdkwdnYiFN11d/FlVVfxlxwlNrkdERAOTGEBz\nyHqU4AHAuHHjMG7cOC1iIT8RRUBRgW0HLsImK6iut2DUMO2WT7nCLivO0bqSKKCHJXJ7JHRRsrTa\n5A5TwThGZvftekRERIHKrQRv1apV2Lx5MyZPnuwyrYmqqhAEAUeOHNEsQNKeeHkaGpusYMUtvpsm\nxqATcWU2FlEUOvSRc1dnU+i0Z7HKMHL6FSIiGkDcSvA2b94MADh69KimwZB/iIKAHha30MRPlo1C\nWIjjLdiXBM8uq9CJItqsdhSXNQAAhsSFI8ToOLfZameCR0REA4pHTbS//e1v8eSTT/a4jbS371gl\nKuvN3e4TFebeZHaihv3fuhPRLj5JRK8HWdhlBXqdiBnjUnH2YgPssoK/7/wGo9IGwWyVUWpqwvcX\nj+9wXJvFjsYWK6LCDZ2clYiIyDNqTxPQ+pBHCd7nn3/eYdv+/fuZ4PlBZb3Za82pooBeV8+8RRQ8\nTzL3HirFyeIaGA0Spo5OxJgRsRgzIhYAcPOkoZAvd+qLDDN0OmjkurGJ2PTBSaxePqnvL4CIKEDV\nNZlxvKjG+XhwbChGp8X6MSLyBbcSvDfffBObN29GaWkpFi1a5Nze0tKCKVOmaBYc+YajguffGCRR\ngM3qWTvxtxfq8MiS8TDoOza/RoT2XL28ftwQFJyt9uiaRP1FUWk97LKC0cP5Qz4Q7TpYgtPn6xAa\nooMoCLhh/BDn9Fe7v7rABE8jXQ328we3ErxFixZh1qxZeOGFF/Czn/3MuT08PBwxMTGaBUe+IQiO\nZk5/Du8WRQGyBxU8WVbQ1GrtNLkjImDfkYtQVZUJ3gBV32zB/QvHIjrC2OG5/ce46s9A4NZSZZGR\nkUhJScELL7yApqYmfPLJJ/jkk088WqaMApckCrDLqnNEqz84Blm4v3+L2Y5hiX2fyiWQ7raIvI5v\n7wHLbJERauy8hqOqwOsfnkJZVbOPoyJf8qgP3qZNm/D2229j3rx5AIAnn3wSK1aswL333qtJcOQb\noiDAZld8MrlxVyTRvX6AqqrifEUTKuta3WqGJRqo2nf2tthk1NS3ITE2DJI/7+TIZ6x2GXpd5/+v\nH1qUjeY2G7Z+UgRFUVFd34an75/u4whJax4leFu2bMHbb7+NsLAwAMD3v/99rFy5kglePyeIjqlG\n/NpE6+ZatJV1bdj6yRlMHjkYM8YP8UFkRP1Lm8WO2kazy/yP2/YVobbRjOljEjF9bJKfIyStXKho\nxFsff4voCCPqGs0u89a2J0kioiOMeHBRNgDg1e2FvgyTfMTjlSwkSer039R/iYIAu6z4ZP3Zrkhu\nzoNnlxWMGR6LnGneGUGsqioaW6zQ68QumzOIAlmpqQm/2fRvzJ3u+EwcP1uDqWMSYLXJSEuKwqvb\nC2HQSZg5aSiaW21+jpa0dPZSA+6anYnMVPaNJw8TvLvvvhvLly93NtF+/PHHWLp0qSaBUeeKLjUh\nIsS7iYgoCLDK/m2idXeQhbcT0UOnTDBbZbS02dhEQf1SQ7MFy3KyMCEzHgAw97phiAwzYMGNI1z2\nO1lc47yJUhQVZ0rrXGbsamqxoqbBDEkUMO/6NF+FT15ktyvQddEs2x32RfaOnlZV8jWPMoXc3Fxc\nd911OHz4MFRVxfr16zF27FitYqNOHCiowuCYjqOi+sJoEHHwZDXSkyO8el5PSG4OslAUFZLovT5E\nTa02LLhxOPYeLvXaOYm8paahDefLm1y2hYXqXKa4aLPYERNhRFx0aLfnar9azKWqZuz+6gImZg12\nPm/QichOj8O2fUVM8Popx/reTNb8xS6rAfX39yjBU1UVRUVFaGxsxOrVq1FWVoaCggJMmDBBq/jo\nGnqdCLvs3buEkSmRSBwUgohQ/zVRiqIA2Y0Mz9tfYK1mG8IvD9aQFdWvzdRE18r/VzFGDRsEve5q\nd5gd+89hyF3hzvdqfZMFqUk9jyiX2n3G7LKC0WmxmDlpaIf9wkI4eKm/ciQYvajgCUBxWQNGJEdr\nENXAIctKQA1i8ugX/dlnn4Uoijh48CBWr16N8PBwPP7449i6datW8dE1JFGA3e7dhWMFQUBMhH+X\n63K3gifLqlc/QLOnpCA6wojJIxPw4efF+M7N6V47N1FfKYqKqWMSXX60o8IN+ODzEudjQQBumJDc\n47kkUXQOZOruRulkcQ3+39avodOJsNsV3Do1lXPp9RNXlm301PI5I7F51zf4wd0s1vTFL1/9ApPa\nVcX9zaMEr6CgANu2bcOSJUsAANHR0bDZ2GlXa7sOnsfxM2UQBCAsRILFqgTd/FaOCp57ffC8WcH7\nySrHSiyTRg7G8p9/iYqaVgDAw4vHee0aRL0lKx1Ht2emxvSqE73U7jPW3Y2SQS/h/oVjERaih9Um\n4+87v8G7+4qweFYGstPjPH8R5DOO70fPE7yocANaLXZ8UViOwTGhHKTRSyeLazE+I97fYTh59E7Q\n6XSQZdk59Lq2thaiF/tDUedKTU1YNScN381Jw+KbUvwdjiYkN9fD9XYF7wqdJGLTs7fj4cXjAmqx\naBrYVBVeG/zU/iaquxul3zx2s7OZ1qCXkLsoG4tmpqOxxeqVOEg7dnvvmwiXz8lCRJgeH3xeDItN\nRn2TxcvRBTdVVREeovN4TXUtefROuPfee/HYY4+hpqYGGzZswKpVq/CDH/xAq9hoAHG7gqdoN51L\nZJh/m6mJtCRJVwdZeHqjFGKQYLHatQqNvMSu9L6Tf0pCJMZnxMNokPC7vx/GS28dhc3L3YGCmdWu\nYMFNI3DfgsAZeOpRE+2dd96J7OxsHDx4EKqqYuPGjcjIyNAqNhpARFFAfbMNVlv3XyiyrMKo8fqz\nAhw/hP6cNobI2xx98C4PslAU6DxofQkx6GC2ylqFRl6gKCqKyxpcBuT0xv0LxqLFbEN5dQt2HDiH\nu2/N9FKEwW3f4VKMGjbI32G48HjYZEZGBpM6P4sO12PwIO9OleJvBr2IrJQIvPXJBURE1CBhUCgM\negl3zkx3dgkoKq3Hx/++gPsXanuH5BiprMAgciLvQFfXaMagqBB/h9EviIJrBc9ocD/BMxoknCiu\nca6MMSEzvsdpWcg3zBY7XnzrKGKjQrDgxhF9buEIMeoQYtTBqJdw/FyNl6IMfrWNFtw0seOodH9y\nK8GbPHlyp0ueqKoKQRBw5MgRrwdGXZt/XfAt0SUKAmaMjceMsfGIHjQYdlnB9v1n0dRqQ1S4o+m0\nsq4Vi2dlYOhgbefr00kibHYFBo0rhdR331//Ma7PTsKgyBBERxgwMWswREFARkp0l8s0DVSSJEC+\nPMWSTVYQLrk/HcrgmFDcfUsmZFlFS5sN7+8/h3nXX11NJjrCyC4OfmKqa8XkkQmYP8O7cxfqJNHr\nMzYEs0Ccg9CtBO/o0aNax0HdGGid/q/MS5edHodNH5yE0SBBUVRU1bVh5byRml//SgWPAt/tM4Y7\nRzwXlzWg1NSEwrM1uNOQjtTEnueGG0gkUYCpthXfnK/FRVMT4qLdr3wKguCcI01VVdQ0tqHw7NXq\nTvGlBvxw2USvx0xdO3epAW/u+gYJsWFYeNOIng/wkI7fgx6x2xXoA2gOPKAXTbREvjJ1dCKmjk4E\n4Lg7uu/ZnT6pqukkEU2t1k6vpZPEXs0zRdobkRyNEcnRiAwz4PdvH8O6R25ACNcXdooMM2DSyMGo\nqm1DcnxErye1FQShw1rQr20/7o0Q+y1VVWFtV+3SSaJLU+nR05UoPFuNGeM6tr4IApA+NMbjptWK\nmhYsmZ2BcRpNyyGJAmxM8NwmB2C/bY++/Z566ik888wziIqKAgA0NDTgN7/5DdavX69JcMHiUlUz\nPvi8GAIEqFAx//o0DEuK6vG4hmYLTLWtMFvYuVkniWiz2GHoYwdid0weNRj7jlzs9LmSskbcMjUF\nNwdYX4uByNFFpOP2SSMH43xFE+qaLBhyOcFrbrW61MFLTU0ovtTgclxEmAGzpwTnNESAYyDTjW5M\niEye+/irCzh7qQFhl9cJN9W0YuqYROfzB4+XY/GsDJRVt3Q4tryqGR9+XoLQy8dabQoeWTKux8ES\n+f8qxpp7p3nxVbhiFwfPBdrfzKME7/Tp087kDnBMdHzq1CmvBxVszpTWY/bkFIwcNgi1jWb89Z8n\nMK3dh78rnxeU48YJQ5AzLRUA5ySyyyoMeu2rZ1cqQZ2RZQV/3nGCCZ4G3vjwFGobzRAEASFG1x83\nRVExfEg0brt+mPNL1GbvfNZ+QRAQFW5Ai9kxCfu5Sw14Z+8ZpLVbzisiVI8po10/g1v2fBvUCV6w\nOV/eCIuuUAivAAAgAElEQVRNRkaKo/pV09CGplYbhg/p+ebZ27bsPYM/PHmrMykrq25Gm/nqtDJZ\nqTHddhloNdsgigJCDDps2fMtWtrsiInsOsGzywoyU2IQHaHtYDsh2GbU11AgdqXyKMFTFAUNDQ2I\njnb8+NXX10OWWV0CHIuCv7uvCABQXt2CtQ/NcD5nt8vOH6LYqBDcdUumW/ML3btgjHNAQVlZmQZR\n9y/fuWmEs3+ev0iS2GnViPquzWrH/QvHIirc0OFOWFVVbPrgFLbvP4fbL3cmbzHbuqxyxEQYsX3/\nWUzMHAyLTcbtN6RhQmb3SwhdqaD4k9lqx8HjFc7HNQ1tfowmcFTUtODtj7/F4JhQ4PJoYFNtCxIG\nheHfJ01oaLHAYpUhyypyF431+Qjf68YmubwXk+M9GwjWfv1fo0GCxdb976rVJvuku0ogJi3+oKoq\nTpyrcf5uS5KAcenxAdckey2PvtEefPBBrFy5ErfffjsAYOfOnfjP//xPTQLrb+KiQ/H9xeMBAK9u\nL3SOMAYclYb2y8dwQefe4TqJwa+rioQgCLhvwRgcOHYJ7+w949x+XXZSp/tPHjUY2RlxOHG2BlX1\nbchK7Xl+KgGCy+fWH85dakBNfRsmZDn6Vd27YIzfYvGlljYbdh083+XzJeUN+M7N6UiOD3duC73c\n/H7mYj2GJUYiLESPvYdKUVRajxCDzu83g71l1Ot6nFTaYpN90pphqmnF2Yv1GBIfjhCDLmASmpqG\nNpwprQfg6N94p8azK7RZ7Nh18Lyzwl9R24ovT1QgrV1Xq8raVs2u31seJXhLlizBuHHjcPDgQQDA\nyy+/jMxMToJ4LaNegtWuOCfktfVyAWiigeKDz4t7nMBaEATMmuxeE6ogCDDqJUwZneB2DDpJgF1W\nodd1/iN2qrgWH3xRjKTYcNxz+2i3zwsAZVXNaLVc/dG2WGV8faYK/zHf9Txmi4y0IVFuJaTB5NsL\ndYiJNGJiVucDBm4RUhDbxXyHo9Ninf+ePHIwPi8sx6YPTuLRpf1zVG9XFbzKulYc+7YKJeWNSE+O\n1nzCdwD4z7snYP/Ri9j91QVkpcZgzvRhPR/kA39+/wTmz0iDKAiYOiYRp8/X9ZjgfXuhDqFGXYem\n8rwdJ/C9O0Z32+exudWGzNQYl65V5dUtsNqv/n+alNV9C4E/eJTgvfzyyy6Pd+7cCQBYvXq19yIK\nAga9BKtNdn4A7V30FaL+yWpTsOPAuU6fG5cRF/AVWptdwfmKRmQMdW+uOLPF3mWTUVVdG2RFwah2\nP7LttVnsKClrxJgRnT9/xaWqZmcF3F8Megm2dt0prlVR24LFMzOw68vz+NfXZS7NV+Mz4rusPsqK\nirx/nnAZ4KCqwKXK5g77mq12zftVBZq3dp/Gxcpm3HP76D43rQ6KCsHCm0bgd38/3OVnFACyhsW4\nJIb1TRZ8ddLRNB4fE4opo9y/MfC2QZFGfHuh3iXJb2i24NX3CjFrUgruXzgWa1/5HLdOTdU8lviY\nUNx9axYqa1vxxfFyr51XUVSoqgpTXSu+KakFAESFGzv0Ta+oacHJ4hqoKvDVyQoMjgkD4KjcT7yc\nULVZ7Hhj5ymcvVSP+iYLBkV2fiNwqaoZY0fEuiR4Ta1WfHr0IkRRwPjMeBj1ErLT4wA4bujMlyup\nptpWRFxTER7SrpocqDxK8MLCwpz/tlgs2LdvH9LT0z264P79+/Hcc89BURQsX74cjzzyiMvzVqsV\na9aswYkTJxATE4MNGzYgJcVx1/7KK6/gnXfegSiK+MUvfoGZM2e6dU5fM+olWKwyIi//ua5toqX+\n7YGFY9HU2nHhdVV1VKIeulObBO9SVTP2Hirtcb/YqBBMyOx66oRj31bh4PFyPLZ8Yqd9hRqaLfj3\nyav9wA4er8DILpbg0etEfHO+FtNGJ2JIfLjLlA2m2lb8+f3juFjZ7PzBVKFiUtZgHPu2CoqqornV\nhugII0b4oWP8tfQ6Ed9eqEPENRP26iURKQkRaGyxIiJMj8yUaBgNkvML3y4reGfvGXS1xrgKFYtn\ndZzO4uyl+g77mq0yEgwDa4LtxlYrfnbPVK+e85G7xqOlzdbl85s/Oo2Lpmbn+tf/PlmBOdNTkTw4\nAu/tO4tBkUakJUWh1NTU6eLxH3153qW5MmOo9z7z4zLi8c/PziH8cp/QUKMOyYMjMGV0ImZOdgzu\nyl2UjRSNJ3xvL8Sog9nSsdlYVlR0+cZvRxQFl5vJ1z88BVlRYdCJuGH8EAiCgO37z2Lq6ASX/T7+\n6gKmjU2EQSfh+nFDOiRZgOPv487NoXL5RuvV7YWOUc6jE1BW3YLVyyfBaHAUZUrKGvB5YRkECEiM\nDXMmcfExoRg9vPub1EDkcR+89h566CE8+uijbh8vyzLWrVuHvLw8JCYmYtmyZcjJyXFp5t2yZQui\noqKwe/du5Ofn4/nnn8eLL76IoqIi5OfnIz8/HyaTCbm5udi1axcA9HhOX0sfGo2te8+gxWyDJIq4\nVNWMRTM9S4QpcIWH6rvs3+PGd12vbdnzLe5fONbZ96grZy824GRx10sMGQ0SFt40Ake/qcQnzReh\nXhP0BVMT5s9Ic1ZTpo9N6raqVF3f5uhD9eV5lwEC1fVtePjOcUiIvXpjeOiUCZeqmvG9O8b0+Dp8\nbeakoThZXIvaRtcR681tVuw5VAqdJCAuOgTzZwzvcGxv5iILMejw6vZCNLfakDAoDLWNZqiq6tfq\nUW/Z7DIuVja5bEsYFIY2S+cVyYPHy1F4thoCBESEen8FjMgwQ7cra9x2fRokSUDY5ffg/qMXMT4j\nHhFhBiyamY6vTlTgH7tPI31oNOI7qSreNCFZ0/nn7l841jk91q4vzyO+otHlZmx0FxVzrYQaJbRd\nk+Cpqopf/ulzZAyN6fZYi012Dkq4orKuFU9+z3WKl2mjE/Ha9uMuCZ7VJnvttYqigIfudEyK3max\no6ahDeMz45EcHxEwfQu9rU/fsG1tbSgt7bmicEVBQQHS0tKQmuooLS9cuBB79uxxScb27t3rbPKd\nP38+1q1bB1VVsWfPHixcuBAGgwGpqalIS0tDQUEBAPR4Tl+bmDXYWT622mR89nUZQgyB9WNG2rDY\nZFTWtUISBY+am2oa2pzVhK6EGnRdNj+0l50e52xm6IqqqnjxH0dxxw3DO04rIcCj92t8TCjiY0Lx\nyJKe76KnjUl0a4ogf0iKC0dSnO+aXe69wzGAwmaXYbMrEAWh307MfPPEoS4rW6iqiktVzYAKxEWH\n4O5bs5zPvfXxadTUm/Ho0gl+G9By7efjuUdvcv47fWg00r1YkeuNlISrzYiJcWE4X96IzNTuEykt\n6SQRjS1WvLq9EICjm8qgSCMmj0zA0pysbo9VVRWnL9S5zByxYm7HFYlmTh7qrFBqLdSoc/kbByuP\nvk0WLVrk/LeiKKitrcUPf/hDt483mUxISro66i0xMdGZpLXfZ8gQx2zfOp0OkZGRqKurg8lkwsSJ\nE12ONZlMANDjOf3JoJcuz2NHA8GEzHh8UViOukaz485V7LlpXlYcA3IGddGJ/IobJnhvDWJBEPCT\nVVO8dj7qPb1O6nFS20A3PjMe47voFvD6h6dcVrpIigvDymWjfBVavxcZZtCsWuguQRDwxMrJzsfN\nbTbY7LJb6w8LguDziiM5eJTg/fGPf7x6oE6HuLg46HTun+LapiCg48zPXe3T1XZF6TifXE93hXl5\necjLy+v0uWHDAmOUEPVPMydxAmSi9q5UKil4OPrC9c9paAYSt7Kza0fPticIAh577DG3LpaUlISK\niqt9dEwmExISEjrsU15ejqSkJNjtdjQ1NSEmJqbbY3s657Vyc3ORm5vb6XNcdo2IiIj6O7eGdoaF\nhXX4TxAEbN26Fa+99prbFxs/fjxKSkpQWloKq9WK/Px85OTkuOyTk5ODbdu2AQB27dqFGTNmOBa3\nzslBfn4+rFYrSktLUVJSggkTJrh1TiIiIqKBxK0KXvvRs83Nzdi0aRO2bt2KBQsWdBhZ2+3FdDqs\nXbsWDz/8MGRZxtKlS5GVlYWXXnoJ48aNw5w5c7Bs2TI8+eSTmDdvHqKjo7FhwwYAQFZWFu644w4s\nWLAAkiRh7dq1kCRHv5XOzklEREQ0UAlqZ53bOlFfX4+8vDzs2LEDd911F+677z7nmrTBZP369fj5\nz3/u7zA6GEhr0SYnJ/e8ExEREXXJrQre//zP/2D37t1YsWIFduzYgfDwwJ/BubcCNbkI1LiIiIgo\n8LhVwRs9ejQMBgMkSXIZoXplYe4jR45oGiQRERERuc/tJloiIiIi6h+4QCoRERFRkGGCR0RERBRk\nmOARERERBRkmeERERERBhgkeERERUZBhgkdEREQUZJjgEREREQUZJnhEREREQYYJHhEREVGQYYJH\nREREFGSY4BEREREFGSZ4REREREGGCR4RERFRkGGCR0RERBRkmOARERERBRkmeERERERBhgkeERER\nUZBhgkdEREQUZJjgEREREQUZJnhEREREQYYJHhEREVGQ0fk7gEDzt7/9Dffff7+/w+igrKzM3yH4\nTHJysr9DICIi6tdYwbvGQEqkiIiIKDgxwSMiIiIKMkzwiIiIiIIMEzzqtbKaNhw5U+vvMIiIiOga\nTPCo1/79TQ2Onan3dxhERER0DSZ41GsWm4LwUA7EJiIiCjRM8KhPBH8HQERERB0wwSMiIiIKMkzw\niIiIiIIMEzwiIiKiIMMEj3pFVVWoKgDB8W8iIiIKHEzwqFe+udCEuCgDJFGArDDBIyIiCiRM8Mhj\niqqipsGCUalR0Esi9n9dBYVVPCIiooDBBI889s6+UlQ1mJEwyIjbrxsCU50Zdrvi77CIiIjoMs5S\nS71y18xU578jw3SwyyoMej8GRERERE6s4JHnrpndWCeJsMtsoiUiIgoUTPCoz3SSALvMJloiIqJA\nwQSP+owVPCIiosDCBI/6TCcKsHOqFCIiooDBBI/6TKcTOIqWiIgogHAULfVZbJQR+45VIsQoISU+\nDDeOi/d3SERERAMaEzzqs5EpkRiZEgkAeHvfBT9HQ0RERGyiJY/8aUcRIkJ5X0BERBTImOCR2xRF\nRWykAQuuT+56J461ICIi8jsmeOQ2m6xAp+NbhoiIKNDx15rcZrMpMPSQ4LWY7The3IDi8mYfRUVE\nRETXYoJHbrPJKvQ9JHhzpyZBEIA9h00+ioqIiIiuxd7y5DarXekxwUtNCAMAnChu8EVIRERE1AlW\n8Mgtx4sbkP9FmTOB65EAqCpHXBAREfkDK3jklso6MxbdmIzBMSFu7W/QibDZFRj0ksaRERER0bVY\nwaMelVa24kRJAyLD9G4fkzk0An/651mYrbKGkREREVFnfFLBq6ioQENDA1JTUxEW5mYTnwfKy8ux\nZs0aVFdXQxRFrFixAvfff7/XrzNQHTpdg/+Ym4YQg/vVuHEjYmCqM6O5ze7RcURERNR3miZ4Fy9e\nxOOPP47q6moYjUZUV1djxowZeOqppzBixAivXUeSJDz99NPIzs5Gc3Mzli5diptuugmZmZleu8ZA\nZpMdExx7KsyoQ6vZDkQbNYiKiIiIuqJpE+3zzz+PlStX4sCBA/j4449x6NAh3Hrrrfj+97+PkpIS\nr10nISEB2dnZAICIiAikp6fDZOI0Hd4kCILHx4wYEo4vTlTj3f2lGkREREREXdG0gldSUoIXX3zx\n6sV0OqxcuRJDhw7FH/7wB/z2t7/1+jUvXryIU6dOYeLEiV3uk5eXh7y8vE6fGzZsmNdjGqiSYkOx\nMicNb39ywd+hEBERDSiaJnhdVX1uvvlmvPDCC16/XktLC5544gn893//NyIiIrrcLzc3F7m5uZ0+\nt379eq/H1e9xthMiIqJ+RdMEr6qqClu2bMHIkSORmZmJ8PBw53O9afLrjs1mwxNPPIFFixbhtttu\n8+q5B6LPCqtQVt0GAH0fCevd/9VERETUA00TvMcffxynTp3C9u3bcebMGYSHhyMrKwtZWVmorq72\n2nVUVcUzzzyD9PT0Litz5Jmy6jasuNU7zdVmiwxZViBJnJWHiIjIFzRN8FauXOnyuKKiAqdPn8bp\n06cxffp0r13n8OHD2L59O0aOHInFixcDAH76059i9uzZXrsG9V7K4DCUVrVieFLXzeZERETkPT5d\nySIpKQlJSUleT7ymTZuG06dPe/Wc5D3REXrY7ezIR0RE5CtsMyPN6SQRNpkJHhERka8wwSPN6SUB\ndlnxdxhEREQDhk8SvA8//BDNzc0AgI0bN2L16tU4ceKELy5NAUAnCbCzgkdEROQzPknwNm7ciIiI\nCBw6dAifffYZlixZgmeffdYXl6YAoJNEVvCIiIh8yCcJniQ5Fpv/9NNPsWrVKsydOxc2m80Xl6Ze\nqKhtgzenKdTrRBRXtKDwXD0Kz9X3fV49IiIi6pZPErzExESsXbsWH374IWbPng2r1QpFYUUnULWY\nZYwY4r0pTZLjQzEhPQY6SYSpzoxvzjd67dxERETUkU+mSXnxxRdx4MABPPjgg4iKikJlZSXWrFnj\ni0tTL8iKCkn0XglPEgWMHhYFADDoRdQ1Wb12biIiIurIJwleaGioy/JhCQkJSEhI8MWlqRdkWYFO\n0mZ9MR1H1BIREWlO0wQvJycHgiAgNjYWW7Zs0fJS5EXeruC1pxMFyBxRS0REpClNE7y9e/dqeXrS\niCyr0Bu16Z6pk0TYFSZ4REREWvLZUmUNDQ04f/48LBaLc5s316Ml79GygidJrOARERFpzScJ3pYt\nW7Bp0yZUVFRg9OjR+PrrrzFp0iRs2rTJF5cnD9llDZto2QePiIhIcz6ZJmXTpk145513kJycjNdf\nfx3btm1DbGysLy5NvSArKnSSRk20ooAWsx0NLTa0WeyaXIOIiChQfVNSi/c+PYv3Pj2L8+XaTRvm\nkwqewWCA0WgEAFitVmRkZKC4uNgXlyYP7Pj8EtosMkx1ZmQO9d48eO2FhugQGabHwRPVKDG14AeL\nMjW5DhERUSA6cOwSlsx2/PZFRRg0u45PErykpCQ0NjZi7ty5yM3NRVRUFKdJCUBtVhkrbh2m6TUk\nUcDcqUkAgLf3XdD0WkRERAFHAAYPCtX8Mj5J8P7whz8AAB5//HFcf/31aGpqwsyZM31xaSIiIiKf\naDXb8M7eM1DbjSUckRzl0u2pqcU3k/37bBTtFdddd52vL0nu8vXgVg6mJSKiIFJqakJKQgRunjgU\nAGCzK/juLz7ALx+egVCjI+VaOW+UT2LxSYJXWFiIP/7xjygrK4PdfrVj/Y4dO3xxeSIiIiLN/fuU\nCddnJ8GglwAAep2jcjc+Mx7Gy9t8xScJ3n/9139hzZo1GDlyJETRJwN3qT/oxUwsiqLixLkaKIqK\njJRoRIR13kH1i8IymK0yAGDamEREdrEfERGRt7RZ7MhKHeR8LAiOHzqDzve5j0+uGBsbizlz5iA1\nNRVDhw51/keBQ1FVaDT1nVdV1rXi06MXUVHbggPHLnW6j6qq+LygHKkJkbDaZBw/W+PjKImIiK66\nkuj5kk8qeE888QSeeeYZ3HDDDTAYrlZSbrvtNl9cntygKCrEfpDhKYqK9KHRmD42CfsOl3a6j82u\nID4mFJmpMbArCqpq23wcJRERkX/5JMHbunUrzp07B7vd7tJEywTPfYqi4mJVK+KjjQgL8f7/Nn8k\neK1mGW/tPQ9BFGDUi1h8U0qPx1xZRi3UqEOrufOJkm12xVkON+gkWO2yV+MmIiJy1/wZaX65rk8S\nvNOnT3NARR8Vlzfji5M1iAzTuZUIeUpWoNnyZF2pqG3D3KlJmDYq1u058eyyAkkUEGKQcPZSg3N7\nm8UOSRRg0DsSOn27Dq42O5dGIyIi7QmddC5fvXySHyLxUR+8iRMnoqioyBeXClp2RcWYtCjYZe/N\nLdLUasP7/7qE/C/KYLUrPq/ghRokxITrATjeiIrS82uTFRWiKEIQBKQkRDiP+d3fD+O17ccBADab\n4hy5pNeJrOAREdGA45MK3uHDh7Ft2zakpKS49MFjVc99sqzCoBMhu5EEuau0shUpg0NRVtOGxhYb\nJB93Al1990hn1VCvF2G1KwgxdD+MXFFU6KTLo5L0Eqw2GSFGHZLiwp372OSrTbR6nQi7XUFlbSsU\n1fVvFxNhRIjR51NBEhERac4nv25//vOfoaqqX0aRBAu5XWLTW4dP1+JCZSvmTE1EVJgebRYZMRF6\n1DZZUdtkha9nsGnfJGy43JTaU4Inyyqky4GGGCSYrTKMl4+JjjDg1e2FsNoUzJrkGKVtNOjw6dFL\nqKprw9CEq+vrqipQXd+G3EXZ3n5ZREREfqdpgrdq1Sps3rwZCxcudEnuriR7R44c0fLyQUWWVWci\nAwD/b/sZPLo4y6NzFF1qQlZKJMqq2xA1TI9Wi4wh8aFoNcvYdqAUjyzK9HbYbjPoRVhsCiJ72M+u\nXG1KNuolWGwy7LKjSXb5nJEd9o8I1eOFH8+GKHQcpv6/rx/C+fJGpA2J8tbLICIiCgiaJnibN28G\nABw9elTLywwI9sujR684VlTf4zGfHquEqc4MVVVxy+RECIKAsBAdLDbHoINWix1hRgmKqiIrJRLJ\ncdovftwVnehe87MiX61kZqfH4f0DZwEAM7KHdHlMV4NHls/JQv6/ivHDZRN7ETEREdFViqIikBoq\n2QGpn+hNE62pzowVtw7DubJmnC5tdDRn6kU0ttgAAG0WGWFGyVFR7c2yEl4kioJbCZ5jmhRHE21G\nSgwyUmJ6fc0RydHQ67myCvUvxWUNOPxNpfOxAOCOG4cjLETvv6CIyDmNV6DwSRPt5MmTIQgC1Had\n3AVBQExMDB566CHcc889fb7W/v378dxzz0FRFCxfvhyPPPJIn88ZSK5MD4JejLEwGkTUNFgwLDEc\nRoOE0qpWDIo0oLreAr1OhKLC73cdkii4NYrWrigw6ry3np+/E1siT311ogI3TkhGRKgjofu8oAyX\nqppdlkcCgM27vkFlXRvCQ/WICNPjuz5a4JxooJIVBZIUOEUDvzbR1tXV4bvf/W6fEzxZlrFu3Trk\n5eUhMTERy5YtQ05ODjIz/denzNtkRYXO0zuDy7sb9RIaWmwYEheKxBgjMpIj0Nxmx03j4y8n3v5P\n8EQRPVbwTLWt+OjgeayY27GvHdFAYbHJGDwoFCEGx9d3dKQROw6cQ1S4EYCjj/P3l4xHVX0bfvTd\nyQDgnEKIiLSjDKQKXnsNDQ04f/48LBaLc9v06dPx+uuv9/ncBQUFSEtLQ2pqKgBg4cKF2LNnT9Ak\neH/f+Q1M1W2YPioWKuCc7sNslbtfP/ZyvmTUi6ioNWNy5iBIkohJmYOu2U/1ex3LnQre2Yv1uOuW\nTIwcNqjb/TyhtiuJ2uxX58/rK6tNxvmKRrf3T0uKgkHvvcokBS+rTYGhXRVbL4k4X96El342FcDV\nZK59k60gBF7zEVGwscuB9RkTVFX13sRqXdiyZQs2bdqEiooKjB49Gl9//TUmTZqETZs2eeX8O3fu\nxIEDB/Dcc88BAN577z0UFBRg7dq1ne6fl5eHvLy8Tp8bNmwY3njjDZdt+w6XouhiQ6f7+8Lxs9UY\nnhQCAKioMcNik2GqMyMuytjtcYMi9IiNNgIqUN1gRkyEETpdxzffmdImqABGpvY0hlU7lXVmGHQi\nYiINiAiP6HSfUlMTIsP1iIkI8dp1C89WY3xGPABg+/6zWDwrwyvnLatuRmOL1dmM1p3mVhuiIgxI\nju/8dRO11/49CwCVta0oPFeNOdOGAXB8X4xNj8Op4hpkpzv2O1VSg6zUGOgk3kQQacVstaPU1NSh\nu4Q3PLx4nMfH+KSCt2nTJrzzzjtYsWIFXn/9dZw9exa///3vvXb+znLU7ubcy83NRW5ubqfPrV+/\nvsO2W6am4papqb0P8Bp1TWbIbqxIYTRIiAwz4LXtx7FgeiwAR/PMzq/KkZYUjhW3DPNKPIdO1wIA\npo2K9cr5euProjqEGCWMSo1CcnJyp/vs/KIE2elxSE30XiL62vbjeHjxONjsMrbvP9urD1FnDhy9\nhMGxoRid1vPf9JuSWlQ3tOHmiUO9cm0KHq1mG6rq21y2NbdZXd6nhUXVKC5vcG57bftx/Mf8Udi2\n7yzuuX00AGDLnm8x97phGBTpvZsjInJ1srgGZ0rrvVYo6CufJHgGgwFGo6PaZLVakZGRgeLiYq+d\nPykpCRUVFc7HJpMJCQkJXju/N1lsMja8eQRTRvcc3zcldXj6/uku23SSCJvdu0VXfyZ2V4ii0GPS\na7HJLnMBelObxbvLmcmqCtHNjo2imwNMaOB595MiDIoKcek6cH12kss+ep0Iq+3qess6SUCbxQ5D\nuxHiE7MG44/vFmDBjSMwMWuw9oETDUA19WaMSvN+9a63fJLgJSUlobGxEXPnzkVubi6ioqK8moCN\nHz8eJSUlKC0tRWJiIvLz8/G73/3Oa+f3JptNxtj0OCyZ3XP/wNfqO3aMlkShw5JbwaCz13XolAlf\nFJYj9PJyYpeqmjHvOu9ULa/V2GLpeScPKIrq9tq+osAEjzpnlxXMu25Yt/0zdTrRZQqliDAD/r7z\nG2SlXp1CaOSwQVh40wi0mu2axks0kFlsdhgDqC+1pgme3W6HTqfDH/7wBwDA448/juuvvx5NTU2Y\nOXOm166j0+mwdu1aPPzww5BlGUuXLkVWlmerPPiK1a5AH0DDqAOFJAqwWRWXbRU1LViWk4Uh8eFd\nHNV3rWYbmtts+PCLkg6Vkb7wZDSVGKRJO/WdoqLHG4URQ6Kw9uEZzsdLb83s9AZDEkXeSBBpyGLV\nrpWpNzRN8JYvX45t27a5bLvuuus0udbs2bMxe/ZsTc7tTVab3PfJdYPwO1oUBcjXJDltFjtCjNp+\nWO66JRNv7voG8dHeXcVDYRNtQHj/wFlcNDXju7eNQmxU/+t/pig9v48kSURa0tXl9gRBgNTJpOii\nwHZIsPsAACAASURBVBsJIi1ZbPLAqeD5YIBuv2Ozu05x4A71mozOJitIiOl/P1bdkTrpg9dmsTub\nZ7WSmhiJR5aMB+Do7/Ta9uOobmjD0/dN7+HI7skeNdH2PAcgue/rb6vw8aEL+Nl/TIWpthUTRw5G\neXVLv0zwVLj/PuqJKML5GatvsmDTByddplIRRQEPLsr2yrWIBqLzFU0ID6AVZTT99aytre1yOhIA\nXY5kDWY2u+LS+bk37pk73DvBBBBHM6Xj3xu3fg2DTkJDs8Wnd0N33+roF/nq9sI+n8ujPnis4HlV\nfbMFpppW5+NQow5mK/ueSaLorOA1NFswengsbrs+zfm8N973RANZRJgeIRoXJTyhaSSKoqClpUXL\nS/QLx89WY8+/SxEeqofVJmOumwMFBMHRR2wgFELbV/D0OtFr05V4W/6/itFmsWNZTvd9PN1pWruC\nCZ532eyyy6jTUIMONY1t3RwxMLR/n9llxeO1rYmof9E0wRs8eDBWr16t5SX6hZoGM+64cbjHKzDM\nGDcEWz8pQmZKtEaRBQ5316L1BVEQupz1v6y62a1zKKpno2iv7X9Inispb0RxWQMKz9a4NO3HRYdg\nx2fn8E1JnXObXifi/oVj/RGm30ii4OwK4Hh/u7YkCBDcrjx/e6EOFtvVqYUSBoUhMTbMuwGTX1hs\nMppbrYjzcr9k8j32wfOB3i4RlJ0eh+z0OABAWVmZt8MKKFd+fDypfGlFrxNhlxVIYu+bhz2v4PX6\nUnTZR1+eR87UVKQlRWHPoQvO7QmxYVhz7zSXff/8fv9Ym1Xw4iKCHSt4rgmeQS/i7KV6pCZEujQz\nybLjzXllEXW7rOCdvWeQM+3q5O9bj53BD5dO9Fqs5D87vyjBF4XlWDFnJCRJ4LyJbqppaPP7b9e1\nNE3w/vrXv2p5+n5DURTnlyN1ThQdibDVrvh9FJJeEmHrYxyOSoh7+wZS9bI/EwQg8/Lcb3sPdf9F\n21/uPa8dYNUX7afjkeWOlbqcaak4dMqE/3vrGF766S0QRQGtZhuey/sK0RFGzJ+RholZg9FqtiNj\naDRmjBviPPb42RqvxUn+1WaxIzE2DJ99fQlx0aF9TvDssoK1r3yB9KFXW6JqG80dbrr6u48Onsfc\n6drM09pbmiZ4MTExPe80AHCR755JogCrXcEHB8swNCnOr7HodCLs9r6V1DxtouX0Fd5VXd+GkvJG\nKG4sCThQtO/nKisd++ClJERiUGQI/vz+CbRa7IgI1aO5zYbpY5MwPiMO5ysaATj6BYeGBE5HcvIu\ns8WOiDA9Gpqs8MbPltUmIzs9zrlsHhCcA3pUAPExgdWszU+pD9hlJng9kUQBDS02RIbqcKef1/HT\n60ScPl+L6Ehjh+csVhlR4Y71ga+lQsX4jHjMGDekXwyyeC7vS/zH/NEYkRx8fTzvWzAGhWerccvU\nlC73UVW12zWrg037Gwm7rHbaqnClaba51YqIUD1azXaEh+gweFAY3v2kCOcuNcJstePG8Z2vF039\nn01WEBVmQElZIxLj+t6v0mZXXAY9Ad7tehAoHN16Aut1McHzAVlRvDaXVbCSRAGtZjuSBoX4/UNy\n88ShOFVSi6q6jiMvb7s+rcvBMqqq4s/vn3AmeIG+kkV5dQuq69v6fYJ35JtKXKxqQkOT1bkteXAE\nkgdHdHmMThJgtrqOthUAt7pSfHuhDsOSIhFicHx9yoqKz45dgtxNR8oxw+M0XZHFHe1vJGRZgUHf\n8ev/ynv2H7tPIyLUgIYWC3KmpiIq3IAng6xJLVAoiooLpiaEGCQkxbn/Hnlg3S786pEbXCa5bu+L\nwjIUnKmGJIlo+P/bu/eAKOv80ePvmYGBkTsioCKKCmpewLyka1qKaIqIlbXd059rZ1u1y26ntFor\njyfTrU67dlnLtm21zVxTES9ZmmVuq6ZrGV5JUUAEBLkOMDPMPOePyYlhZgRkYAA/r390nst3vvN8\neWY+z/eqNzBvZrztbxasNbFvrDtCeEgnDCYz82b+0ody0qieDOoTxnfH850lzQ+ZlyjXG+nfM5Qu\nIVevtTKaHKcG06hVbNyTic7XmymjezXyE7dtZhcPTZ4kAV4rsFgUhw7Nwp5araLaYMZH6/nrFBas\nY2xC9yafp1KpbH2mzE1tovVADZ63twZjM5ui24IDxy5y+619ufXGHg0f/LMRN0TyyRen7LZlXSzn\npbmjGzz39X8eZsHdQ20DoCr0Rk6ev0zSyJ5Oj6821PLV4Rzundzf6f7WUncUba2TUbRXbF6R0uZ+\nqDqyzJwSvjyUQ3FZDc//z02NOsdsthAWrCPt6zNERwYwNqE7nYN0VFQZOZ1dwtC4cHYdzOHZ2SPR\nqFWk7T1DWaUR39BffvJLKgwkxHVh6q9iHJpMQwJ8CfLzYX/GRafv/+WhHMYN7c5X/83hrsS4q+bV\nZDbjXW9y//tv609ltYmPPz/l4qz2p9ZJtwdPkwCvFZilibZBXhoVFy9Xo9OGeTorzeKlVrM6LQN9\ntcnhS80VrbeGSyXVTpt96wrw8+bXE/u5I5soioKPtwZjnaku2itvr6bVfID9CPUrGrr+pRUG9mdc\nxGRW7K6b0WQmsrOfXSfyuswWxeUPZWtSq1VcKqnip5xS8i5V0jnI+coe7g7uDCazbSQuWLtANPbe\nuB7oa2oZOTCS/54qbPQ5xloLCbFdmDyqFzmFFRw5dYmJI6M5dKKAbfuy6N0tiJjugbbfnUA/LcfO\nFnG5zB+zxYJFUcjOryDI/5duKO9vycCiKAR00gLWvxdFgTfW/RcvjZpfDenGjf3CAQjopCU+tgtH\nM4sazKvJZEHrVX/EtoZQbw2dfLwwmzvOIMS21uVDAjw3qTVbKLxsnT1f662x62zZlGWrrldabw3/\n+9cD8Ne17z/J2dew1JO3l5rH7xna4HGN7Zj84bbjmK5SM5dTWMHv7oxH5+PVIQK8lvZTbil7j1wg\nt7CC22/py5yUgXZzwBlMZrRXGXGt8VATfH3+nbQM6hNGbmEFXYJ19OrqvGnP3V7++0HievzSraGw\npIon772xVd67PaiuqSUk0AeNWu10+hpnjD//zXUJ0aH1VpOZbZ3j0VRrISrCn7wiPYF+Wtvxw/pH\n8P3pQi6X16BWW1c1iezsxw0xoQD0igyk2ljL9LGO/Z/9dVqSx8Rw+GSBLcAD8NKoOXuhjA/Sj101\nr6WVBkYOjHS6T+utwWAy06mDBHhtTfv+NW1DKqtMfJ95CYADx/LtmnrMlrbXNt8WBfq1nTX82qLK\nKhM1htoGl8KpNVuuuhLI7u+yKSypws/Xm3K9kYoqo91+lUqFv659lMWug+fd+vDkbODF0cwipt0c\nQ+cgHRq1ilPnL1Nw+Zel0IwmMz6NWH6wbg2hvtpE4ogeDOoTRq3ZQq3ZglqlQuutobLaZDuu/vrM\nzaFRq66p60FzKIpCVLi/3QjKhmpK24JyvZGd+88xfVwft03b9N3xfI5nXXbYfj6/nNnTBuKr1VBj\nNOOvc/xbqqox8bf0YwT6aXlo6g12/dr8dd7oa6xL8ZlqLfTqGkj6vrMk1pmnMNBPy7ihrgccJd3k\nvHuB2WLB20uDr1aDwej4MPjiI6MbNeeQq3tU6622Bni+3mzck0lRWY3d4LScwopGdZsQzkmA5ybB\nAT5M/VUMAHmX9HY/FG1xdI1of6IjAqhqRIDXkE6+Xhw8lk9CXBglFQa2/zvLbv/J8yUsuDuB0EDn\nTXhtybmLFW5b1u5KH7X6/WiMtWZ0Pl62e1hbr2nbaLLg42TAQl1zUwfbvc7KK+NMbimD+oTx9oYf\nCAn05XhWMf+TMpCNe34i+ueO831cNPu2F82dT9JTzl0s44fMS4wZ0u2qg3Wa4vvTl5g1zfnqKd5e\nGsJDdPxj+3FUwG9SB9masc9fLGfHf84R37cLJ7OtAaJ1OT7rfo1GbevDa6o1M2ZId2bc0tcteTYY\nzQR00uLzc/AJ1v5/V+4F67/X/ttmvZesrQ3F5TXMTR1k94C1ZscJqg21divTiMaTq9YCvL3Utgl7\ni8uqOZdXftUmHCEaw9fHixpDbbPTCQn05fDJQmbc0sfpckQ795+jssrYLgI8d04EfGX+w/pNZAaj\n2S5ICQ30Zeu5Es5drACsNXhTx8Q06b10Pl5U/VyWnXy9eXDKANbsOEFJuYHxw3sw8gbnTVrthdls\nYeNXP+GlUdv6dF2hoLT5flc1BjOdg3RUu+F+s1Fx1b6HE0f2ZOLInvxr92kqq02EBFiPPXg8n5Sx\nvena2Y+T2ZdRFIWPPjtJYp01zUvKa/jueD6FJdUOU5I0R3FZDVHh/vh4a6isNlJVY7KOPm9EjXVj\n9OsZQtreM6hVKvx1Wofa86FxXVi74wQms4WYbkFtZsStxaLwU26pXdeLCr3xKmd4hgR4LUDrrcFk\nsv4oVFSZGDUoUmrwRLPpfDSN+sFpqMWkf89Q3n56gsv9nXys859db7w0akxmC/XD2vrzeAX5+7Dg\n7oRmvZfOx8uhLH21GkorDUS6Ye4xT6usNlGuN7L56zOsfznZbt+QPmH8n78d4MU23PRWbaglJMDH\nVmvVkPl/+pLXn7jFLQ/yPvWaQyurTESEdkKtVqFWWaeT6trFz64/3B3j+1JcVkP6N2d54Db3jda+\ne2Ic3cL80GjUhAXp2PTVGQC3PYD07xlK/56hLvcP6hPGoD5hGExmh1HvnpRbWMEXB7MZ0veXQYF1\nl+5rKyTAawE+P/cr8Mfad6KTb/vozyTatrjoEHZ8e44vD+dc9bjoiIBmvU+XUB3b9mWx7wfr+sd9\nooIYP6ztfXmBeydM9fZSk3GmCP96NU4lFTVuHx3XydcbfbV9gOej1VBWaWi1wQ8tqcZoto3S9dXa\nBz03DerK0TNFTZ5o2mxRMNVpGtd6a2xz++3PuEiXEB2xPZzPUVnfT7mlRIR2cqhdBOsyWufzywkJ\n9G10DV5ekZ4t35y1pde/V4jL+eka4qv1YnVahm1keF5Rpa1WOS46hDU7TjgMWugTFUzv7tbuBe6s\nwas70vzOCbFuS7eptF5qW1NuW1BRZSIhtgtj4tv2hN8S4LUAnY8XaXvP4ufrRV6Rvk1G9qL9iQoP\nYO6MwQ0f2Ez1n6pXp2U0KsAzmMx2TcheGjV+LTxYw51NtOOHRZGZU4q+ziAHgMmjnHdAbw5vLzV5\nlypZ98UpjLXWoMXP15tvj15k9OCuDZzd9hmMv/SbchbE+Wq9uFRabZt4V6XCabBV17Z9ZykoqcJX\na+2q0DXMj+EDIliz/QRDYsM4cvpSowO8D9KPMW5olK1s6wab/9p1mkF9wgj011Je2bhmt4kjoxlx\nQ4T1hQJrtp/gxv6/1LA1pR/ilVo6Z31LxyZ0dzlQRqVSse7/JnfIKWhUKhX5xXp27j8PQGyPYJfT\nErUGfbWpXQxEkwCvBUy6qSeV1SYeeOEzVj+XRHgDM30L0ZY1Noh6b/OPhIf80rx49KdLLHnkV+1m\niqDOQTqnfRJbylP3D6PWbLE16916YxTDB0TYzU3WXtXU67dY368Gd2VPnZroguIqunXxp6ur+QxV\nUFRWw4NTBuCrtU7v8/Hnp8gtrGTiyGiG9gsnu840Qvt+uMDFIr3L9/fTedvVzj298htGDowk0E9L\neZWRMfHdyMoro6C4ymUadWm9NHY1do/eOcSueXfHt+eobmS3h4G9OxMVfm0DO9rjgJbGeuT2wbYH\nyM/2n+eR7i3/sOtKZbWR0KC2X9MuAV4L0GjUti/pLsG6Njf5oRBN0cnH2+XUFiXlNTz1wDBUKhW+\nWi/unvjLrPYWRUFfY2qwZuZamZuw3m9bVH80dN3vjfamS4jObp5Gk8nCbaN7kexi8EmfqGD6RAXb\nXteaLWRml7pM/8czRRw/W2wLYK5MDr7ncA4PThkAgL/vL3+nIQE+LlcWAUgcEW2rDQLo3T2IUYOs\nNacjfu5f5qv1osZ4bX1R6z8opI7rQ5XB5OJoe13D/Dy+rF1bVPfhsaXm7zx7ocw2KfnIgZH0rfM3\nWldlldTgCVzP/yNEe3H/VTptr91xgs/2n6dcb6Cqxv4HLCTAh30/5Ll15JvFothGrlXV1NqtrSk8\nJ3Wc4wS5AL+9Y0ijzvfSqBkQ47qzfedgX/TVJruH5aceGGZ3TFOWgjObLWScKeLjz62vi8tq6FGv\n76qvj4bL5TVcLq9xON9Xq2lS32rreq3SkuMuZrOCpQUWEDhwLJ/EEdGoVSo+P3DeFuApisK3P160\nTSB/POtyu+h6Jd+OLei2NjKkW4iW4t9Jy3fH83n810PxqdeZfuLInqyut/rG6ewSAv20TV5a7IoP\nth6zG5E+qgP0VxMNCw/pdE2rxLii0ah5bvZIriwB7WwN0cBOWoL8fdhzyHFQU3ZBhazG4UFdw/yo\nqDI61HgfO1vssta1exf/Br939NUmwoJ1qFXY1biWVhg4draY8cOsk0XfM6mfw2CstkgCvBY0b2a8\np7MgRIuKjw0jKtzfadOiRq0CBYpKq21L963fdZohsWFOl0RqDIuiMGeaeyY2Fte3hmrgNBq1y5rJ\nVZuOUlKnZk9pA0vRXU+GxIaxftdpFCA8RMeMW/piNlvYuu+sy1U50vae4X/d7lijrCgKb234AV+t\nF+V6g+0BskJvtHU7qDGYGRPfrdGDeNoKCfCEENcsplsQMd1cj2abODKa5//6LcMHWEcYHjtbzLA6\nowubIr9Yj6kNTZUgrl83DYxkd52avSv990TruDLSX1EUVm+x9rss1xvpExVsNz9gXYdPFtj+X1xW\nTVZeOWDtfxka6Mt99Zr4f3+ffReA9kgCPCFEiwkL1nHhUiV/XZgIwH2T+/H/Pv4vcdEhdp3swfok\nveT9A8ydMYhuYY6jCHd9l93qa6kK4UxCXDgJcdf2oCLcR6VS4aVW8+7mH1GpYPyNjesXt+3fWdwQ\n05nL5TWs+/wUyTc3bSWa9kICPCFEiwn007L6uSTb606+3kwf24fCkiqHAK/gchWmWjNHTl2yra15\nhanWQrneaDfxqhBCNLZvpsWs8MmuU1jMCmculPHQVOu6wB9/foqbBrbvpQFdkQBPCNFiVCoVEaH2\nS2/17BrIJ1+cYn9GPqMGRTJ6sHU2+I8/P8VvUgdz4txlfvypyCGtGbf0kVHpQohrMjMxFosFOvl6\ncXudQTUf/HGSB3PVsiTAE0K0qkA/LXNnDKakooY120/go/XiUkk1Pt4aenUN7BBLdQkh2pbWnMS8\nrWj3Ad7y5cvZs2cP3t7eREdHs2zZMgID5QdCiLYu2N+HW4dFYaq1ENMtkFuGSv86IYRwF/etSuwh\nY8aMYevWraSnp9OrVy9WrVrl6SwJIRpBpVIxpG8XhvWPIC46xGFlByGEENeu3Qd4N998M15e1h+G\nhIQE8vPzPZwjIYQQQgjPavcBXl2ffvop48aN83Q2hBBCCCE8ql20icyaNYuiIsdRdU888QQTJ04E\n4J133kGj0TB9+vQG0/vggw/44IMPnO6Ljo5uXmaFEEIIITxMpXSANVY2bdrEunXr+Pvf/45O17yR\nMh9++CEPP/ywm3ImhBBCCNH62n2At3fvXl555RXWrl1LaGiop7MjhBBCCOFx7T7AS0pKwmg0Ehxs\nnRU/Pj6eJUuWeDhXQgghhBCe0+4DPCGEEEIIYa9DjaIVQgghhBAS4AkhhBBCdDjtYpoUTzpy5Aiv\nvvoqgwYNckt6hYWFAISHh7eptNydnrvTysrK4qabbmp2WlfSg7b3Od2dXkvkzV3l0NY/p7vSu17K\nwN3ptfW8SRl4Nr3r8TchISGhyTN8SB+8Bhw5coRHH32U/fv3uy09gKFDh7aptNydnrvTkjLwbFpX\n0nNXObT1z+mu9K6XMnB3em09b1IGnk1PfhMaR5pohRBCCCE6GKnBa4Rx48axd+9eT2fjuiZl0DZI\nOXielIHnSRl4npRBw6QGTwghhBCig5EATwghhBCig5EATwghhBCig9G8+OKLL3o6E+2Bu0bEiGsn\nZdA2SDl4npSB50kZeJ6UwdXJIAshhBBCiA5GmmiFEEIIIToYCfCEEEIIIToYCfCEEEIIIToYCfCE\nEEIIIToYCfAasHfvXiZPnkxSUhLvvvuup7PToU2YMIGUlBRSU1O54447ACgtLWX27NlMmjSJ2bNn\nU1ZWBoCiKCxdupSkpCRSUlI4duyYJ7Pebi1atIjRo0czbdo027ZrueabNm1i0qRJTJo0iU2bNrX6\n52jPnJXBypUrGTt2LKmpqaSmpvL111/b9q1atYqkpCQmT57MN998Y9su31XX7uLFizz44INMmTKF\n5ORkPvzwQ0DuhdbkqgzkXmgGRbhUW1urJCYmKtnZ2YrBYFBSUlKUzMxMT2erwxo/frxSXFxst235\n8uXKqlWrFEVRlFWrVikrVqxQFEVRvvrqK2XOnDmKxWJRjhw5osycObPV89sRHDx4UMnIyFCSk5Nt\n25p6zUtKSpQJEyYoJSUlSmlpqTJhwgSltLS09T9MO+WsDP7yl78oq1evdjg2MzNTSUlJUQwGg5Kd\nna0kJiYqtbW18l3VTAUFBUpGRoaiKIpSUVGhTJo0ScnMzJR7oRW5KgO5F66d1OBdxdGjR+nZsyc9\nevRAq9WSnJzM7t27PZ2t68ru3buZMWMGADNmzGDXrl1221UqFQkJCZSXl1NYWOjJrLZLI0aMICgo\nyG5bU6/5vn37GDNmDMHBwQQFBTFmzBi7p2lxdc7KwJXdu3eTnJyMVqulR48e9OzZk6NHj8p3VTOF\nh4czcOBAAPz9/enduzcFBQVyL7QiV2XgitwLDZMA7yoKCgqIjIy0vY6IiLjqH5xovjlz5nDHHXfw\nySefAFBcXEx4eDhg/QK4fPky4Fg2kZGRUjZu0tRrLvdJy/joo49ISUlh0aJFtqZBV9daysB9cnNz\nOXHiBPHx8XIveEjdMgC5F66VBHhXoTiZA1qlUnkgJ9eHjz/+mE2bNvHee+/x0Ucf8d1337k8Vsqm\n9bm65lIW7nfvvffyxRdfkJaWRnh4OK+88gogZdDS9Ho9jz32GM8++yz+/v4uj5NyaDn1y0DuhWsn\nAd5VREZGkp+fb3tdUFBge5oT7hcREQFA586dSUpK4ujRo3Tu3NnW9FpYWEhoaCjgWDb5+flSNm7S\n1Gsu94n7hYWFodFoUKvV3HXXXfz444+A6+8kKYPmM5lMPPbYY6SkpDBp0iRA7oXW5qwM5F64dhLg\nXcXgwYM5d+4cOTk5GI1Gtm3bxoQJEzydrQ6pqqqKyspK2////e9/Exsby4QJE9i8eTMAmzdvJjEx\nEcC2XVEUvv/+ewICAq7bm9jdmnrNb775Zvbt20dZWRllZWXs27ePm2++2ZMfod2r2590165dxMbG\nAtYy2LZtG0ajkZycHM6dO8eQIUPku6qZFEXhueeeo3fv3syePdu2Xe6F1uOqDOReuHZens5AW+bl\n5cXixYv5zW9+g9ls5s4777T9cQn3Ki4uZt68eQCYzWamTZvGuHHjGDx4ME888QQbNmyga9eu/PnP\nfwbglltu4euvvyYpKQmdTsfLL7/syey3W7///e85ePAgJSUljBs3jgULFvDII4806ZoHBwfzu9/9\njpkzZwIwb948goODPfaZ2htnZXDw4EFOnjwJQPfu3VmyZAkAsbGxTJkyhalTp6LRaFi8eDEajQZA\nvqua4fDhw6SlpREXF0dqaipgLRe5F1qPqzLYunWr3AvXSKU4a7AWQgghhBDtljTRCiGEEEJ0MBLg\nCSGEEEJ0MBLgCSGEEEJ0MBLgCSGEEEJ0MBLgCSGEEEJ0MBLgCSGEEEJ0MBLgCSHcYsCAAaSmpjJt\n2jQee+wxqqurXR6bm5vLkCFDSE1NZerUqTz99NOYTCYADhw4wLBhw0hNTSUlJYVZs2ZRXFwMwMaN\nGxk1apTtvPXr1wOwZcsWUlJSSElJ4Z577rHNm5Wbm8u0adNa+JPD3LlzKS8vb/H3aez7rly5kvff\nf7/F3veee+4BrNc3PT3dtv3HH39k6dKlLfa+QojGkwBPCOEWvr6+pKWlsXXrVry9vVm3bt1Vj4+O\njiYtLY309HTy8/PZsWOHbd/w4cNt+wYPHsxHH31k2zd16lTS0tJYs2YNr7/+OkVFRURFRbF27VrS\n09N59NFH+eMf/9hin9OZ9957j8DAwFZ9T0++75WyvXDhAlu3brVtHzx4MM8//3yr50cI4UhWshBC\nuN3w4cM5deoUf/rTn+jWrRv3338/YK1Z8vPzs60zCaDRaBgyZAgFBQUO6SiKgl6vJzo62mFf586d\niY6OJi8vjxtvvNG2PSEhwW4tytraWp555hmOHz9OTEwMy5cvR6fT2aW1fv16PvnkE0wmEz179mTF\nihXodDoWLlyIj48PZ8+eJS8vj2XLlrFp0ya+//574uPjbQufT5gwgQ0bNlBVVcXcuXMZNmwYR44c\nISIigrfffhtfX19OnDjBCy+8QHV1NdHR0bz88ssEBQXZ5ePLL7/knXfewWQyERwczKuvvkpYWBh6\nvZ6lS5eSkZEBwPz585k8ebLtfUNDQ3nnnXfYvHkzXbt2JTQ0lIEDBzpcs4ULF6LVavnpp58oLi5m\n4cKFjB8/HoPBwIsvvkhGRgYajYaFCxcyatQoMjMzWbRoESaTCYvFwsqVK+nVqxdDhw7lyJEjvPba\na5w5c4bU1FRuv/12BgwYwN/+9jdWrVpFaWkpzz77LDk5Oeh0OpYsWUL//v1ZuXIleXl55ObmkpeX\nx8MPP8xDDz3U4N+UEKJppAZPCOFWtbW17N27l7i4OJKTk+1q5nbs2MFtt91md7zBYOCHH35g7Nix\ntm2HDh0iNTWVW2+9lW+//da29FNdOTk55OTkOAR/GzZsYNy4cbbXWVlZ3H333aSnp+Pn58c///lP\nh7SSkpL49NNP2bJlC71792bDhg22feXl5fzjH/9g0aJF/Pa3v2XWrFls27aN06dPc+LECYe0AxxZ\nnQAABVZJREFUzp8/z/3338+2bdsICAhg586dADz99NM89dRTpKenExcXx5tvvulw7rBhw1i/fj2b\nN28mOTmZ1atXA/D222/j7+9Peno66enpjBo1yu68jIwMtm/fzubNm3nzzTdtC7I7c+HCBdauXcuq\nVat44YUXMBgMthrS9PR0XnvtNRYuXIjBYGDdunU89NBDpKWl8emnnxIZGWmX1h/+8AdbbeusWbPs\n9q1cuZIbbriB9PR0nnzySZ555hm7Mnn//ff517/+xVtvvWVrnhdCuI/U4Akh3KKmpsa2huTw4cOZ\nOXMmWq2W4uJiCgoKKCkpITAwkG7dupGbm0t2djapqamcP3+eyZMn079/f1taw4cPZ9WqVQC8++67\nrFixwrYG5fbt2zl8+DBarZYlS5bYrfW5f/9+NmzYYBfEde3alWHDhgEwffp01qxZw5w5c+zynpmZ\nyRtvvEFFRQV6vd5ugfjx48ejUqno168fYWFh9OvXD4C+ffty4cIFBgwYYJdWVFSUbdvAgQO5cOEC\nFRUVVFRUMHLkSABuv/12Hn/8cYdrmJ+fz5NPPsmlS5cwGo1ERUUB8J///IfXX3/ddlz9mr9Dhw4x\nceJEW83k1RZXnzJlCmq1ml69etGjRw/Onj3L4cOHeeCBBwDo06cP3bp1Iysri4SEBP7617+Sn5/P\npEmT6NWrl8t06zt8+DArV64EYPTo0ZSWllJRUQFY13LVarWEhoYSGhpKcXGxQ/AohGgeCfCEEG5x\npQ9efZMnT2bnzp0UFRWRnJxs236lD15hYSEPPvggu3fvJjEx0eH8xMREFixYYHs9depUFi9e7HDc\nyZMnef7553nvvfcICQmxbVepVHbH1X8N1qbLt99+m/79+7Nx40YOHjxo26fVam3nXfk/gFqtpra2\n1iGtusdoNBoMBoPDMa4sXbqUWbNmkZiYyIEDB2y1fIqiOM13Q5+rMcepVCpcLUmekpJCfHw8X331\nFXPmzGHp0qWMHj26Ue/jLM0r713/Gjm7jkKI5pEmWiFEi0pOTmb79u3s3LmTyZMnO+wPDw/nqaee\n4t1333V6/uHDh532wasrLy+PBQsWsGLFCmJiYhz2HTlyBIBt27bZavPq0uv1dOnSBZPJZDcq1F0C\nAgIIDAzk0KFDAKSlpTFixAiH4yoqKoiIiABg8+bNtu1jxoxh7dq1ttdlZWV2540YMYIvvviCmpoa\nKisr2bNnj8u8fPbZZ1gsFrKzs8nJySEmJoYRI0bYPndWVhYXL16kd+/e5OTk0KNHDx566CEmTJjA\nqVOn7NLy8/NDr9c7fZ8RI0awZcsWwDoyOiQkBH9/f5f5EkK4l9TgCSFaVGxsLHq9nvDwcMLDw50e\nM3HiRFauXGkLgK70wVMUhYCAgAan3njrrbcoLS3lpZdeAqy1Qhs3bgSsTY6bNm1i8eLF9OrVi3vv\nvdfh/Mcff5y77rqL7t27ExcX5zJoaY7ly5fbBln06NGDZcuWORwzf/58Hn/8cSIiIoiPjyc3NxeA\nRx99lCVLljBt2jTUajXz58+3G6gycOBApk6dSmpqKt27d3caxF4RExPDAw88QHFxMS+99BI+Pj7c\nd999vPDCC6SkpKDRaFi2bBlarZbt27ezZcsWvLy8CAsLY968eXZp9evXD41Gw/Tp07njjjvsmqvn\nz5/PokWLSElJQafT2QakCCFah0pxVTcvhBCiQ1m4cCG33nqrw0AXIUTHI020QgghhBAdjNTgCSGE\nEEJ0MFKDJ4QQQgjRwUiAJ4QQQgjRwUiAJ4QQQgjRwUiAJ4QQQgjRwUiAJ4QQQgjRwUiAJ4QQQgjR\nwfx/aTcqDe8P6C0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efbfefe4b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_pi_d(filename=\"%s/PvRBP2b_Pi_D_all_SNPs_20170803.pdf\" % output_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}