stories/gwss-pbs/gwss-pbs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genome-wide selection scane (GWS) using the population branching statistic (PBS)\n",
    "\n",
    "This notebook illustrates an analysis that uses genetic variation data and scans the genome for evidence of positive selection in a specific population. The analysis is illustrated using scikit-allel and the Ag1000G phase 2 dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import dask.array as da\n",
    "import allel\n",
    "import intake\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace with different setup if using a different dask scheduler\n",
    "from dask_kubernetes import KubeCluster\n",
    "from dask.distributed import Client\n",
    "cluster = KubeCluster(n_workers=10)\n",
    "client = Client(cluster)\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# open malariagen data catalog for accessing data on google cloud\n",
    "cat = intake.open_catalog('https://malariagen.github.io/intake/gcs.yml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['samples', 'snps', 'snps_pass', 'snps_pass_biallelic', 'haps']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Anopheles gambiae 1000 genomes phase 2\n",
    "list(cat.ag2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<zarr.hierarchy.Group '/' read-only>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# here are the SNP genotype data in scikit-allel-style zarr format\n",
    "callset = cat.ag2.snps_pass.to_zarr()\n",
    "callset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/\n",
      " ├── 2L\n",
      " ├── 2R\n",
      " ├── 3L\n",
      " ├── 3R\n",
      " ├── X\n",
      " └── samples (1142,) object\n"
     ]
    }
   ],
   "source": [
    "# data are grouped first by chromosome - 2L, 2R, 3L, 3R are the autosomes of Anopheles gambiae\n",
    "print(callset.tree(level=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3R\n",
      " ├── calldata\n",
      " │   └── GT (14481509, 1142, 2) int8\n",
      " ├── samples (1142,) object\n",
      " └── variants\n",
      "     ├── ABHet (14481509,) float32\n",
      "     ├── ABHom (14481509,) float32\n",
      "     ├── AC (14481509, 3) int32\n",
      "     ├── AF (14481509, 3) float32\n",
      "     ├── ALT (14481509, 3) |S1\n",
      "     ├── AN (14481509,) int32\n",
      "     ├── Accessible (14481509,) bool\n",
      "     ├── BaseCounts (14481509, 4) int32\n",
      "     ├── BaseQRankSum (14481509,) float32\n",
      "     ├── Coverage (14481509,) int32\n",
      "     ├── CoverageMQ0 (14481509,) int32\n",
      "     ├── DP (14481509,) int32\n",
      "     ├── DS (14481509,) bool\n",
      "     ├── Dels (14481509,) float32\n",
      "     ├── FILTER_BaseQRankSum (14481509,) bool\n",
      "     ├── FILTER_FS (14481509,) bool\n",
      "     ├── FILTER_HRun (14481509,) bool\n",
      "     ├── FILTER_HighCoverage (14481509,) bool\n",
      "     ├── FILTER_HighMQ0 (14481509,) bool\n",
      "     ├── FILTER_LowCoverage (14481509,) bool\n",
      "     ├── FILTER_LowMQ (14481509,) bool\n",
      "     ├── FILTER_LowQual (14481509,) bool\n",
      "     ├── FILTER_NoCoverage (14481509,) bool\n",
      "     ├── FILTER_PASS (14481509,) bool\n",
      "     ├── FILTER_QD (14481509,) bool\n",
      "     ├── FILTER_ReadPosRankSum (14481509,) bool\n",
      "     ├── FILTER_RefN (14481509,) bool\n",
      "     ├── FILTER_RepeatDUST (14481509,) bool\n",
      "     ├── FS (14481509,) float32\n",
      "     ├── HRun (14481509,) int32\n",
      "     ├── HW (14481509,) float32\n",
      "     ├── HaplotypeScore (14481509,) float32\n",
      "     ├── HighCoverage (14481509,) int32\n",
      "     ├── HighMQ0 (14481509,) int32\n",
      "     ├── InbreedingCoeff (14481509,) float32\n",
      "     ├── LowCoverage (14481509,) int32\n",
      "     ├── LowMQ (14481509,) int32\n",
      "     ├── LowPairing (14481509,) int32\n",
      "     ├── MLEAC (14481509, 3) int32\n",
      "     ├── MLEAF (14481509, 3) float32\n",
      "     ├── MQ (14481509,) float32\n",
      "     ├── MQ0 (14481509,) int32\n",
      "     ├── MQRankSum (14481509,) float32\n",
      "     ├── NDA (14481509,) int32\n",
      "     ├── NoCoverage (14481509,) int32\n",
      "     ├── OND (14481509,) float32\n",
      "     ├── POS (14481509,) int32\n",
      "     ├── QD (14481509,) float32\n",
      "     ├── QUAL (14481509,) float32\n",
      "     ├── REF (14481509,) |S1\n",
      "     ├── RPA (14481509,) int32\n",
      "     ├── RU (14481509,) object\n",
      "     ├── ReadPosRankSum (14481509,) float32\n",
      "     ├── RefMasked (14481509,) bool\n",
      "     ├── RefN (14481509,) bool\n",
      "     ├── RepeatDUST (14481509,) bool\n",
      "     ├── RepeatMasker (14481509,) bool\n",
      "     ├── RepeatTRF (14481509,) bool\n",
      "     ├── STR (14481509,) bool\n",
      "     ├── VariantType (14481509,) object\n",
      "     ├── altlen (14481509, 3) int32\n",
      "     ├── is_snp (14481509,) bool\n",
      "     └── numalt (14481509,) int32\n"
     ]
    }
   ],
   "source": [
    "# within each chromosome group (e.g., 3R), data look like this\n",
    "print(callset['3R'].tree())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"allel allel-DisplayAs2D\"><span>&lt;GenotypeDaskArray shape=(14481509, 1142, 2) dtype=int8&gt;</span><table><thead><tr><th></th><th style=\"text-align: center\">0</th><th style=\"text-align: center\">1</th><th style=\"text-align: center\">2</th><th style=\"text-align: center\">3</th><th style=\"text-align: center\">4</th><th style=\"text-align: center\">...</th><th style=\"text-align: center\">1137</th><th style=\"text-align: center\">1138</th><th style=\"text-align: center\">1139</th><th style=\"text-align: center\">1140</th><th style=\"text-align: center\">1141</th></tr></thead><tbody><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">0</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">1</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">2</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">...</th><td style=\"text-align: center\" colspan=\"12\">...</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481506</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481507</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481508</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr></tbody></table></div>"
      ],
      "text/plain": [
       "<GenotypeDaskArray shape=(14481509, 1142, 2) dtype=int8>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# inspect genotype calls from a single chromosome arm (e.g., 3R) - the genotype calls are the main data we operate on\n",
    "gt = allel.GenotypeDaskArray(callset['3R/calldata/GT'])\n",
    "gt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.7/site-packages/dask/dataframe/utils.py:14: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  import pandas.util.testing as tm\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ox_code</th>\n",
       "      <th>src_code</th>\n",
       "      <th>population</th>\n",
       "      <th>country</th>\n",
       "      <th>location</th>\n",
       "      <th>site</th>\n",
       "      <th>contributor</th>\n",
       "      <th>contact</th>\n",
       "      <th>year</th>\n",
       "      <th>m_s</th>\n",
       "      <th>sex</th>\n",
       "      <th>n_sequences</th>\n",
       "      <th>mean_coverage</th>\n",
       "      <th>ebi_sample_acc</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AA0040-C</td>\n",
       "      <td>Twifo_Praso__E2</td>\n",
       "      <td>GHcol</td>\n",
       "      <td>Ghana</td>\n",
       "      <td>Twifo Praso</td>\n",
       "      <td>Twifo Praso</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>2012</td>\n",
       "      <td>M</td>\n",
       "      <td>F</td>\n",
       "      <td>95033368</td>\n",
       "      <td>30.99</td>\n",
       "      <td>ERS311878</td>\n",
       "      <td>5.60858</td>\n",
       "      <td>-1.54926</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AA0041-C</td>\n",
       "      <td>Twifo_Praso__H3</td>\n",
       "      <td>GHcol</td>\n",
       "      <td>Ghana</td>\n",
       "      <td>Twifo Praso</td>\n",
       "      <td>Twifo Praso</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>2012</td>\n",
       "      <td>M</td>\n",
       "      <td>F</td>\n",
       "      <td>95843804</td>\n",
       "      <td>31.70</td>\n",
       "      <td>ERS311886</td>\n",
       "      <td>5.60858</td>\n",
       "      <td>-1.54926</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AA0042-C</td>\n",
       "      <td>Takoradi_C7</td>\n",
       "      <td>GHcol</td>\n",
       "      <td>Ghana</td>\n",
       "      <td>Takoradi</td>\n",
       "      <td>Takoradi</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>2012</td>\n",
       "      <td>M</td>\n",
       "      <td>F</td>\n",
       "      <td>107420666</td>\n",
       "      <td>35.65</td>\n",
       "      <td>ERS311894</td>\n",
       "      <td>4.91217</td>\n",
       "      <td>-1.77397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AA0043-C</td>\n",
       "      <td>Takoradi_H8</td>\n",
       "      <td>GHcol</td>\n",
       "      <td>Ghana</td>\n",
       "      <td>Takoradi</td>\n",
       "      <td>Takoradi</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>2012</td>\n",
       "      <td>M</td>\n",
       "      <td>F</td>\n",
       "      <td>95993752</td>\n",
       "      <td>29.46</td>\n",
       "      <td>ERS311902</td>\n",
       "      <td>4.91217</td>\n",
       "      <td>-1.77397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AA0044-C</td>\n",
       "      <td>Takoradi_D10</td>\n",
       "      <td>GHcol</td>\n",
       "      <td>Ghana</td>\n",
       "      <td>Takoradi</td>\n",
       "      <td>Takoradi</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>David Weetman</td>\n",
       "      <td>2012</td>\n",
       "      <td>M</td>\n",
       "      <td>F</td>\n",
       "      <td>103044262</td>\n",
       "      <td>33.67</td>\n",
       "      <td>ERS311910</td>\n",
       "      <td>4.91217</td>\n",
       "      <td>-1.77397</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    ox_code         src_code population country     location         site  \\\n",
       "0  AA0040-C  Twifo_Praso__E2      GHcol   Ghana  Twifo Praso  Twifo Praso   \n",
       "1  AA0041-C  Twifo_Praso__H3      GHcol   Ghana  Twifo Praso  Twifo Praso   \n",
       "2  AA0042-C      Takoradi_C7      GHcol   Ghana     Takoradi     Takoradi   \n",
       "3  AA0043-C      Takoradi_H8      GHcol   Ghana     Takoradi     Takoradi   \n",
       "4  AA0044-C     Takoradi_D10      GHcol   Ghana     Takoradi     Takoradi   \n",
       "\n",
       "     contributor        contact  year m_s sex  n_sequences  mean_coverage  \\\n",
       "0  David Weetman  David Weetman  2012   M   F     95033368          30.99   \n",
       "1  David Weetman  David Weetman  2012   M   F     95843804          31.70   \n",
       "2  David Weetman  David Weetman  2012   M   F    107420666          35.65   \n",
       "3  David Weetman  David Weetman  2012   M   F     95993752          29.46   \n",
       "4  David Weetman  David Weetman  2012   M   F    103044262          33.67   \n",
       "\n",
       "  ebi_sample_acc  latitude  longitude  \n",
       "0      ERS311878   5.60858   -1.54926  \n",
       "1      ERS311886   5.60858   -1.54926  \n",
       "2      ERS311894   4.91217   -1.77397  \n",
       "3      ERS311902   4.91217   -1.77397  \n",
       "4      ERS311910   4.91217   -1.77397  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load sample metadata - each row has data on an individual mosquito\n",
    "df_samples = cat.ag2.samples.read()\n",
    "df_samples.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Illustration of some key steps\n",
    "\n",
    "Before showing the full PBS analysis, first illustrate some key steps that prepare and compute the input data.\n",
    "\n",
    "To prepare the data, we need to be able to subset to variants that pass pre-defined filters, and to samples in a given population."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Locate pass variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True, ...,  True,  True,  True])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# e.g., for chromosome 3R\n",
    "loc_pass_variants = callset['3R/variants/FILTER_PASS'][:]\n",
    "loc_pass_variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14481509"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.count_nonzero(loc_pass_variants)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Locate samples in population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# e.g., the Anopheles gambiae mosquitoes from Burkina Faso\n",
    "pop = 'BFgam'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 67,  81,  82,  83,  91,  92,  93,  97,  98,  99, 100, 101, 102,\n",
       "       103, 104, 105, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120,\n",
       "       121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133,\n",
       "       134, 135, 136, 137, 138, 139, 140, 155, 156, 157, 158, 159, 160,\n",
       "       161, 163, 164, 165, 166, 169, 173, 174, 182, 184, 185, 186, 188,\n",
       "       189, 191, 192, 194, 197, 203, 204, 205, 206, 207, 210, 211, 214,\n",
       "       215, 218, 219, 220, 221, 222, 223, 224, 226, 227, 229, 230, 232,\n",
       "       233])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loc_pop_samples = df_samples[df_samples.population == pop].index.values\n",
    "loc_pop_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "92"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(loc_pop_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"allel allel-DisplayAs2D\"><span>&lt;GenotypeDaskArray shape=(14481509, 92, 2) dtype=int8&gt;</span><table><thead><tr><th></th><th style=\"text-align: center\">0</th><th style=\"text-align: center\">1</th><th style=\"text-align: center\">2</th><th style=\"text-align: center\">3</th><th style=\"text-align: center\">4</th><th style=\"text-align: center\">...</th><th style=\"text-align: center\">87</th><th style=\"text-align: center\">88</th><th style=\"text-align: center\">89</th><th style=\"text-align: center\">90</th><th style=\"text-align: center\">91</th></tr></thead><tbody><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">0</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">1</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">2</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">...</th><td style=\"text-align: center\" colspan=\"12\">...</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481506</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481507</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481508</th><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">...</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td><td style=\"text-align: center\">0/0</td></tr></tbody></table></div>"
      ],
      "text/plain": [
       "<GenotypeDaskArray shape=(14481509, 92, 2) dtype=int8>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# select rows and columns from the genotype calls array\n",
    "gt_pass_pop = gt.subset(loc_pass_variants, loc_pop_samples)\n",
    "gt_pass_pop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Allele count computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"allel allel-DisplayAs2D\"><span>&lt;AlleleCountsArray shape=(14481509, 4) dtype=int32&gt;</span><table><thead><tr><th></th><th style=\"text-align: center\">0</th><th style=\"text-align: center\">1</th><th style=\"text-align: center\">2</th><th style=\"text-align: center\">3</th></tr></thead><tbody><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">0</th><td style=\"text-align: center\">184</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">1</th><td style=\"text-align: center\">184</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">2</th><td style=\"text-align: center\">184</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">...</th><td style=\"text-align: center\" colspan=\"5\">...</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481506</th><td style=\"text-align: center\">184</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481507</th><td style=\"text-align: center\">180</td><td style=\"text-align: center\">  4</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td></tr><tr><th style=\"text-align: center; background-color: white; border-right: 1px solid black; \">14481508</th><td style=\"text-align: center\">184</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td><td style=\"text-align: center\">  0</td></tr></tbody></table></div>"
      ],
      "text/plain": [
       "<AlleleCountsArray shape=(14481509, 4) dtype=int32>\n",
       "184   0   0   0\n",
       "184   0   0   0\n",
       "184   0   0   0\n",
       "...\n",
       "184   0   0   0\n",
       "180   4   0   0\n",
       "184   0   0   0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The PBS test is based on comparisons of different populations. \n",
    "# To compare populations, we perform an allele count within each population. \n",
    "# This step is data-intensive because it makes a pass over the genotype calls\n",
    "# array.\n",
    "# This step is not compute-intensive.\n",
    "# This is a good step to watch the dask dashboard.\n",
    "# E.g.:\n",
    "ac_pass_pop = gt_pass_pop.count_alleles(max_allele=3).compute()\n",
    "ac_pass_pop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6022557"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# how many segregating SNPs?\n",
    "ac_pass_pop.count_segregating()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PBS multi-population genome-wide scan for signals of selection\n",
    "\n",
    "Now we illustrate the full PBS analysis, which we want to be able to run for any chromosome and any set of three populations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def population_allele_counts(chrom, pop):\n",
    "    \"\"\"Convenience function to obtain allele counts for a given chromosome and population.\"\"\"\n",
    "    gtz = callset[chrom]['calldata/GT']\n",
    "    gt = allel.GenotypeDaskArray(gtz)\n",
    "    loc_pass_variants = callset[chrom]['variants/FILTER_PASS'][:]\n",
    "    loc_pop_samples = df_samples[df_samples.population == pop].index.values\n",
    "    print(pop, 'no. samples:', len(loc_pop_samples))\n",
    "    gt_pass_pop = gt.subset(loc_pass_variants, loc_pop_samples)\n",
    "    ac_pass_pop = gt_pass_pop.count_alleles(max_allele=3)\n",
    "    return ac_pass_pop\n",
    "\n",
    "\n",
    "def pbs(chrom, pop1, pop2, pop3, window_size=100, min_maf=0.02, normed=True):\n",
    "    \"\"\"Compute the PBS statistic in windows across the genome.\"\"\"\n",
    "    \n",
    "    # load variant positions\n",
    "    loc_pass_variants = callset[chrom]['variants/FILTER_PASS'][:]\n",
    "    pos = callset[chrom]['variants/POS'][:][loc_pass_variants]\n",
    "    \n",
    "    # compute allele counts in each population\n",
    "    ac1 = population_allele_counts(chrom, pop1)\n",
    "    ac2 = population_allele_counts(chrom, pop2)\n",
    "    ac3 = population_allele_counts(chrom, pop3)\n",
    "    ac1, ac2, ac3 = da.compute(ac1, ac2, ac3)\n",
    "    ac1 = allel.AlleleCountsArray(ac1)\n",
    "    ac2 = allel.AlleleCountsArray(ac2)\n",
    "    ac3 = allel.AlleleCountsArray(ac3)\n",
    "    print('finished allele counting')\n",
    "\n",
    "    # locate and subset to segregating variants at sufficient frequency\n",
    "    ac = ac1 + ac2 + ac3\n",
    "    loc_seg = ac.is_biallelic_01() & (ac.to_frequencies()[:, :2].min(axis=1) > min_maf)\n",
    "    print('found', np.count_nonzero(loc_seg), 'segregating variants')\n",
    "    pos = pos[loc_seg]\n",
    "    ac1 = ac1[loc_seg]\n",
    "    ac2 = ac2[loc_seg]\n",
    "    ac3 = ac3[loc_seg]\n",
    "    print('finished locating segregating variants')\n",
    "    \n",
    "    # setup windows over the genome\n",
    "    starts = allel.moving_statistic(pos, statistic=lambda v: v[0], size=window_size)\n",
    "    starts[0] = 1  # fix to start of sequence\n",
    "    ends = np.append(starts[1:] - 1, [np.max(pos)])\n",
    "\n",
    "    # compute pbs in windows\n",
    "    res = allel.pbs(ac1, ac2, ac3, window_size=window_size, normed=normed)\n",
    "    print('finished computing PBS')\n",
    "\n",
    "    # Visualise the results - a peak indicates a genome location under selection\n",
    "    # in the first population.\n",
    "    fig, ax = plt.subplots(figsize=(14, 3), facecolor='w')\n",
    "    x = (starts + ends) / 2\n",
    "    ax.plot(x, res, marker='o', linestyle=' ', mfc='none', mec='k', markersize=2)\n",
    "    ax.set_title(f'Signals of selection in {pop1} (reference populations: {pop2}, {pop3})')\n",
    "    ax.set_xlabel(f'Chromosome {chrom} position (Mbp)')\n",
    "    ax.set_xticklabels(np.array(ax.get_xticks()) / 1e6)\n",
    "    ax.set_ylabel('PBS')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BFgam no. samples: 92\n",
      "UGgam no. samples: 112\n",
      "GW no. samples: 91\n",
      "finished allele counting\n",
      "found 1431465 segregating variants\n",
      "finished locating segregating variants\n",
      "finished computing PBS\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0YAAADgCAYAAAAniXZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydeXRVRbb/vzc3JAJhCFMEEsAQGRPATkBxRGRwRHFAHFqcWl1ofCBxbPuJ3SxwSARFH4jaOIv4XjeoNOAEjigQUREHZoSEeQY1kKR+f/CrY91KnXPPne/N/X7WOgtuzlSnzq46e9fetcsjhBAghBBCCCGEkCQmJdYFIIQQQgghhJBYQ8OIEEIIIYQQkvTQMCKEEEIIIYQkPTSMCCGEEEIIIUkPDSNCCCGEEEJI0kPDiBBCCCGEEJL00DAipJ7y2muvYciQIRG/z+LFi5GdnR326z744INo1aoVjj/++LBdc/z48bj22mvDdj3JxIkTcfPNN4f9utF6h/HGs88+izFjxgR9/ueff44TTzwRGRkZmDNnThhLVv8ZMGAAnn/++aDO/eWXX5CRkYGampowl4qEk+3bt6N79+6oqqqKdVEIiTtoGBGSwHz22Wc49dRT0axZM7Ro0QKnnXYali1bBgC45ppr8N5778W4hMGxefNmlJWV4YcffsC2bdtiXRwfTIbgAw88ELQy6UQo73D8+PFo0KABMjIykJGRge7du+P//u//rP2LFy9GSkqKtT8jIwMXXXRRuIoeNEeOHMGECRNw9913B32N//7v/8Ydd9yBQ4cO4ZJLLglj6YhKp06d8MEHH1i/O3TogEOHDsHr9casTCbDTm+zQgg8/fTT6NWrFxo1aoTjjz8eAwYMwKxZs6JdXFccPHgQd911Fzp16oTGjRujQ4cOuPzyy7F06VIAwK233orRo0dbxx89ehSNGzc2/u3LL79EVlYWzj77bMyYMSPqz0JIvEPDiJAE5cCBA7jwwgtRXFyMPXv2oKKiAg899BDS09NjXbSQ2bRpE1q2bIk2bdrEuigJzZVXXolDhw7h0KFDmDJlCq699lps377d2t+uXTtr/6FDh/DOO+/EsLTHmDt3Lrp164b27dsb91dXV/u9xqZNm9CzZ8+g7u/m+iSxufPOOzFlyhSUlZVh9+7dqKiowIQJE7BgwYJYF60OVVVVGDhwIFauXIl3330XBw4cwI8//oiRI0fiP//5DwDgzDPPxMcff2yds3z5cnTo0AGffPKJz98AoLCwEMCxQZdnn302ik9CSGJAw4iQBGX16tUAgKuuugperxcNGzbEkCFD0KtXLwDAiy++iNNPP906/r333kPXrl3RrFkzjB49GmeddZY1siqPLSkpQWZmJk444QTMnz/fOnfmzJno3r07mjRpgtzcXMcP6qOPPor27dujSZMm6Nq1Kz788EPjcfv378d1112H1q1bo2PHjpgwYQJqa2vxwQcfYPDgwaisrERGRgauv/76Oufu2rULF154IZo3b44WLVrgjDPOQG1tLQCgsrISl112GVq3bo0TTjgBTz31lG1Zv/zyS5x66qlo3rw5evfujcWLF1v79uzZgxtuuAHt2rVDZmYmLrnkEhw+fBjnnXeeVbaMjAxUVlbWCdF7++230bNnTzRv3hwDBgzAjz/+aO3r1KkTSktL0atXLzRr1gxXXnklfv/9d2P59Hfo8Xgwffp0nHjiicjMzMTtt98OIYTt86kMHToUTZo0wbp16/we+9tvv2HUqFHIzMxE9+7d8dhjj/mMuD/yyCPo3LkzmjRpgh49euDf//63T5lPO+00jB07Fs2bN0dubi6++OILvPjii8jJyUGbNm3w0ksv2d57/vz5OOuss6zfGzduhMfjwQsvvIAOHTpg4MCBAIB//vOf6N69OzIzMzF06FBs2rQJANC5c2esX78eF110ETIyMlBVVYX9+/fjpptuQtu2bdG+fXs8+OCDVriXWt4WLVpg/PjxqKqqQklJCTp06ICsrCzcdttt+O233wD84X0oKytDmzZt0LZtW8ycOdOn7saNG4eOHTuiWbNmOP30061zneRNp1OnTpg0aRJ69OiBzMxM3HDDDT5y8txzzyEvLw8tWrTAsGHDUFlZae3zeDx46qmnkJubi1atWuHuu++22ocuq7J+TQbhunXrMHDgQLRs2RKtWrXCNddcg3379gEA/vznP+OXX36x6vmxxx6rc63KykoMGzYMLVq0QF5eHp577jnr2uPHj8eIESNw3XXXoUmTJujZs6elvAPu+5FAWb16Nf7nf/4Hs2bNwuDBg9GwYUN4vV6cfvrpePHFF63jNmzYgDPPPBNNmjTBoEGDcPvtt/vU2xVXXIHjjz8ezZo1w5lnnolVq1ZZ+66//nqMHj0a5513HjIyMnDaaadh27ZtGDNmDDIzM9GtWzesWLHCVXlfeeUVbNmyBXPmzEF+fj68Xi8aN26Myy+/HOPHjwcAnHXWWfjxxx+xa9cuAMCnn36KkSNH4vDhwz5/69+/Pxo0aAAAOPnkk7F+/Xqr3RBCjkHDiJAEpUuXLvB6vRg1ahTmz5+PvXv32h67a9cuXH755Zg0aRJ2796Nrl274osvvvA55quvvkLXrl2xa9cu3HPPPbjpppsspbtNmzbWaOXMmTMxduxYfP3113Xu8/PPP+Ppp5/GsmXLcPDgQSxcuBCdOnUylqm4uBj79+/H+vXr8fHHH+Pll1/GzJkzMWjQIMyfP9/yZqjKiqSsrAzZ2dnYuXMntm/fjokTJ8Lj8aC2thYXXXQRevfujYqKCnz44YeYMmUKFi5cWOcaFRUVuOCCC/Dggw9iz549KC0txWWXXYadO3cCOKb4/frrr1i1ahV27NiBsWPHonHjxj5lO3ToENq1a+dz3dWrV+Oqq67ClClTsHPnTpx//vm46KKLcOTIEeuY2bNnY8GCBdiwYQO+++474zPa8e6772LZsmX49ttvMXv2bOOz6QghMG/ePBw5cgQ9evTwe/zDDz+MjRs3Yv369Xj//ffx6quv+uzv3LkzPv30U+zfvx8PPfQQrr32WmzdutXa/9VXX6FXr17YvXs3rr76aowcORLLli3D2rVr8eqrr1phbiZWrlyJrl271vn7xx9/jB9//BELFy7EnDlzMHHiRPzrX//Czp07ccYZZ+Cqq64CcEyZ79ChA9555x0cOnQI6enpGDVqFFJTU7F27VqsWLEC7733nk+41VdffYXc3Fzs2LEDf/3rX3Hvvfdi9erV+Oabb7B27VpUVFTg73//u3X8tm3bsH//flRUVOCFF17A7bffbrW/kpISlJeX44svvsCePXvw2GOPISUlxa+8mXjttdewcOFCrFu3DqtXr8aECRMAAB999BHuv/9+zJ49G1u3bkXHjh0xcuRIn3P//e9/Y/ny5fj6668xd+5c/POf/7S9jx1CCNx///2orKzEjz/+iM2bN1vK+CuvvOJTz/fcc0+d86+66ipkZ2ejsrIS//u//4sHHnjAx8B5++23MXLkSOzbtw/Dhg3DHXfcAcC5H/nss8/QvHnzgJ9F8tFHHyEnJwdFRUWOx1199dXo168fdu/ejfHjx+OVV17x2X/eeedhzZo12LFjB/70pz/hmmuu8dk/e/ZsTJgwAbt27UJ6ejr69++PP/3pT1ZffNddd7kq7wcffIChQ4eicePGtsdkZ2ejY8eO+PTTTwEAn3zyCc444wyceuqpPn8788wzrXNSU1ORl5eHb7/91lU5CEkaBCEkYfnhhx/EqFGjRPv27YXX6xUXXXSR2LZtmxBCiJkzZ4rTTjtNCCHESy+9JE455RTrvNraWpGdnS2ee+4569jOnTtb+w8fPiwAiK1btxrve/HFF4spU6YIIYRYtGiRaN++vRBCiDVr1ojWrVuL999/Xxw5csS23NXV1SItLU2sWrXK+tv06dPFWWedVeeaJv72t7+JYcOGiTVr1vj8/csvvxQ5OTk+f5s4caK4/vrrhRBCPPTQQ+Kaa64RQgjxyCOPiGuvvdbn2CFDhogXX3xRVFZWCo/HI/bs2VPn3qayqdf9+9//Lq644gprX01NjWjXrp1YtGiREEKIjh07ildeecXaf/fdd4tbb73V+JzqOxRCCADi008/tX5fccUVYtKkScZzH3roIdGgQQPRrFkz0bBhQ5GSkiIeffRRn+fweDyiWbNm1vbmm28KIYQ44YQTxIIFC6xjn3vuOcf30bt3bzFnzhyrzHl5eda+7777TgCw5FIIIVq0aCFWrFhhvFZeXp6YP3++9XvDhg0CgFi3bp31t3PPPVc8//zz1u+amhrRsGFDsXHjRiHEsTp+//33hRBCbNu2TaSlpYlff/3VOv71118XAwYMsMqrykxtba1o1KiRWLt2rfW3L774QnTq1Mmqt+OOO04cPXrU2t+6dWuxZMkSUVNTI4477jjxzTff1HkuJ3kz0bFjRzFt2jTr97x580Rubq4QQogbb7xR3H333da+gwcPitTUVLFhwwYhxDE5UevwmWeeEQMHDhRC+MqqEH/Ur3yes846y+oXdP7973+LPn36+JRR1rN+rV9++UWkpKSIAwcOWPvvu+8+MWrUKKsc55xzjrVv1apV4rjjjhNCuO9HTJjKr7bZf/zjH+Lkk0/22d++fXvRrFkzkZ6eLjZu3Cg2bdokvF6vOHz4sHXMNddc41NvKnv37hUAxL59+4QQQowaNUrcfPPN1v6nnnpKdOvWzfr93XffiWbNmrl6nnPOOUfce++91u8VK1aIZs2aiSZNmoguXbpYfx81apQYM2aMqKmpEa1btxaHDx8W06ZNs/7WvHlzsXjxYp9rn3rqqeKll15yVQ5CkgV6jAhJYLp3744XX3wRW7Zswffff4/KykpjNq/Kykrk5ORYvz0eT50EAmr2t0aNGgGANao/f/58nHLKKWjRogWaN2+O//znP1aIhkpeXh6mTJmC8ePHo02bNhg5cqRPiI9k165dOHLkCDp27Gj9rWPHjqioqHD13HfffTfy8vIwZMgQ5Obm4pFHHgFwbG5JZWUlmjdvbm0TJ070mVcj2bRpE9566y2fYz/77DNs3boVmzdvRosWLZCZmemqPCqVlZU+z5WSkoKcnByfZ9Pr2s57YiKQc0eMGIF9+/bh119/xbp16/Dyyy/7hEG2a9cO+/bts7YRI0ZYz6DKi/p/AHj55ZfRp08fq96+//57H3nIysqy/t+wYUPj3+zKnZmZiYMHD9b5u1qGTZs24b/+67+s+7do0QJCCKP8bNq0CUePHkXbtm2t42+99Vbs2LHDeO2dO3fi119/RWFhoXX8ueee6+PZadmyJVJTU63f8j3s2rULv//+Ozp37mwsh5282aGWq2PHjlZb0mUsIyMDLVu29Hl+u3MDYceOHRg5ciTat2+Ppk2b4tprrzW2exOVlZVo0aIFmjRp4lMOp3bw+++/o7q62nU/YiI1NRVHjx71+dvRo0etELKWLVvWqfMtW7Zg165dqKqqghDCKrvsBwHf+qypqcF9992Hzp07o2nTppY3y6kNuJV/Hb28ffr0wb59+/Cvf/3LJ6vcmWeeiU8++QQrV65Ebm4uGjVqhNNPP93622+//YaTTz7Z59oHDx4MyftGSH2EhhEh9YRu3brh+uuvx/fff19nX9u2bbFlyxbrtxDC57cTVVVVuOyyy1BSUoLt27dj3759OP/8823ntlx99dX47LPPsGnTJng8Htx77711jmnVqhUaNGjgE9/+yy+/2E6412nSpAnKysqwfv16vPPOO3jiiSfw4YcfIicnByeccIKPsn/w4EFrkrJKTk4O/vznP/sce/jwYdx3333IycnBnj17rPkUKh6Px7Fs7dq183kuIQQ2b97s+tkiRadOnXDeeee5SrCgy8vmzZut/2/atAl/+ctf8PTTT2P37t3Yt28f8vPzXc918kevXr2s+XMqar3n5OTg2Wef9Xl3v/32G0499dQ65+Xk5CA9PR27du2yjj1w4IDPnBD12q1atULDhg2xatUq6/j9+/e7UmRbtWqF4447zjiPy0ne7FDr/ZdffrHCNnUZO3z4MHbv3u0jY3bnNm7cGL/++qu1zynr4/333w+Px4PvvvsOBw4cwKuvvurznp3aQrt27bBnzx4fIzeQNu6mHzHRoUMHbNy40edvGzZssAzJgQMHYsuWLT7zmXTatm2LPXv2+NSTWp+vv/465s6diw8++AD79++37heuNqByzjnn4L333sPhw4cdjzvzzDPx7bffYt68eTjjjDMAAD179sTmzZsxb9489O3bF8cdd5x1fHV1NdauXYvevXuHvcyEJDI0jAhJUH766SeUlZVZCuzmzZvxxhtv4JRTTqlz7AUXXICVK1dizpw5qK6uxjPPPOM6DfaRI0dQVVWF1q1bIzU1FfPnz7dNIf3zzz/jo48+QlVVFY477jhrYrOO1+vFiBEj8Ne//hUHDx7Epk2b8MQTT7heY+jdd9/F2rVrIYRA06ZN4fV64fV60a9fPzRt2hSPPvoofvvtN9TU1OD777+3UpirXHvttXjnnXewcOFC1NTU4Pfff8fixYuxZcsWtG3bFueddx5Gjx6NvXv34ujRo1aGp6ysLOzevRv79+83lm3EiBGYN28ePvzwQxw9ehRlZWVIT083Ku3RZMuWLViwYIGrbG0jRozApEmTsHfvXlRUVODpp5+29h0+fBgejwetW7cGcCwxh8kYD5bzzz/fJ8OWidtuuw2TJk2yjJv9+/fjrbfeMh7btm1bDBkyBOPGjcOBAwdQW1uLdevW2d4jJSUFf/nLXzB27FjLq1RRUeFqLldKSgpuvPFG3HXXXaisrERNTQ2WLFmCqqoqR3mz45lnnsGWLVuwZ88eTJw4EVdeeSWAY0bDzJkz8c0336CqqgoPPPAATj75ZJ/5fI8//jj27t2LzZs348knn7TO7dOnDz755BP88ssv2L9/PyZNmmR7/4MHDyIjIwPNmzdHRUUFHn/8cZ/9WVlZWL9+vfHcnJwcnHrqqbj//vvx+++/47vvvsMLL7xQZy6OCbf9iIkrr7wSM2fOxNKlSyGEwOrVqzF58mRrDlbXrl1x6623YuTIkXj//fetfkKdc9mxY0cUFRVh/PjxOHLkCJYsWeIzoHDw4EGkp6ejZcuW+PXXX/HAAw+4Kpsd119/vTHJDABcd911aNu2LYYPH47vv//ekh3dsMvLy0NWVhaefPJJyzDyeDw4+eST8eSTT/rMLwKApUuXolOnTj6eR0IIDSNCEpYmTZrgq6++wsknn4zGjRvjlFNOQX5+PsrKyuoc26pVK7z11lu455570LJlS/zwww8oKipyldq7SZMmeOqppzBixAhkZmbi9ddfx7Bhw4zHVlVV4b777rMWZt2xYwcmTpxoPHbq1Klo3LgxcnNzcfrpp+Pqq6/GjTfe6OrZ16xZg0GDBiEjIwP9+/fH6NGjMWDAAHi9Xrzzzjv45ptvcMIJJ6BVq1a4+eabjUZMTk4O5s6di4kTJ6J169bIycnB448/bmXveuWVV9CgQQN069YNbdq0wZQpUwAc88xdddVVyM3NRfPmzeuE+HTt2hWvvvoqiouL0apVK7zzzjt45513kJaW5urZwsmbb75pZc/r27cvTjvtNDz00EN+z/vv//5vZGdn44QTTsCgQYNw+eWXW7LSo0cPjBs3Dv3790dWVhZWrlyJ0047LWxlvuiii/DTTz85hk4NHz4c9957L0aOHImmTZsiPz/fJ4uizssvv2wlnsjMzMTll1/uGML26KOPIi8vD6eccgqaNm2KQYMG4eeff3ZV/tLSUhQUFKBv375o0aIF7r33XtTW1vqVNxNXX321FS6am5uLBx98EMAxL8I//vEPXHbZZWjbti3WrVtXZw2eiy++GIWFhejTpw8uuOAC3HTTTQCAwYMH48orr0SvXr1QWFiICy+80Pb+Dz30EL7++ms0a9YMF1xwAS699FKf/ffffz8mTJiA5s2bo7S0tM75b7zxBjZu3Ih27dph+PDhePjhhzF48GC/dejUj3z66afIyMiwPXfo0KF45JFHcMMNN6BZs2Y4//zzMWrUKNxyyy3WMc888wzuvPNO3HXXXWjRogWys7Pxt7/9DW+++SY6dOgA4FjiiyVLlqBly5Z48MEHceWVV1pt4LrrrkPHjh3Rvn179OjRwzgYFQibN2+2bUPHHXccFi1ahB49euCCCy5A06ZN0bVrVyxbtgyzZ8/2OfbMM8/Ezp07fa51xhlnYMeOHXUMo9deew233XZbSOUmpD7iEZHw/RJC4pra2lpkZ2fjtddew9lnnx3r4pA4Z9q0aZg1a5ZfT064mDFjBn744QfLGE1GOnXqhOeffx6DBg0K+FyPx4M1a9YgLy8vAiVLTq688kp069YNDz/8cFive+TIEfTu3RvfffedNQ8q0uzYsQNnnXUWVqxY4RNeRwihx4iQpGHhwoXYt28fqqqqMHHiRAghQh7pJPWTrVu34vPPP0dtbS1+/vlnlJWVYfjw4VG7/y233JLURhGJPcuWLcO6detQW1uLBQsWYO7cubjkkkvCfp+0tDT8+OOPUTOKgGPLL/z44480iggxkOr/EEJIfWDJkiW4+uqrrZCiOXPmWBnDCFE5cuQIbr31VmzYsAHNmzfHyJEjMXr06FgXi5CosW3bNlx66aXYvXs3srOzMW3aNJx00kmxLhYhJMIwlI4QQgghhBCS9DCUjhBCCCGEEJL00DAihBBCCCGEJD31ao5Rq1atfNZxIIQQQgghhBCVjRs3YteuXXX+Xq8Mo06dOjmuZk0IIYQQQghJboqKiox/ZygdIYQQQgghJOmhYUQIIYQQQghJemgYEUIIIYQQQpIeGkaEEEIIIYSQpIeGESGEEBICJSUlSE9PR0lJSayLQgghJARoGBFCCCEhMHXqVCxcuBBTp06NdVEIIYSEAA0jQgghJASKi4sxdOhQFBcXx7oohBBCQsAjhBCxLkS4KCoq4jpGhBBCCCGEEFvsbAZ6jAghhJAQUecZcc4RIYQkJvQYEUIIISGSnp6OhQsXYujQoQBg/b+qqirGJSOEEKITE4/RggUL0LVrV+Tl5eGRRx6ps/+nn35C//79kZ6ejtLSUp99nTp1QkFBAfr06YOioqJIFpMQQggJCXWeEeccEUJIYhIxj1FNTQ26dOmC999/H9nZ2ejbty/eeOMN9OjRwzpmx44d2LRpE+bMmYPMzEyfsINOnTph+fLlaNWqlet70mNECCGEEEIIcSLqHqOlS5ciLy8Pubm5SEtLw8iRIzF37lyfY9q0aYO+ffuiQYMGkSoGIYQQQgghhPglYoZRRUUFcnJyrN/Z2dmoqKhwfb7H48GQIUNQWFiIGTNm2B43Y8YMFBUVoaioCDt37gypzIQQQkgwMOECIYQkPhEzjEwReh6Px/X5n3/+Ob7++mvMnz8fzzzzDD755BPjcbfccguWL1+O5cuXo3Xr1kGXlxBCCAkWLvJKCCGJT8QMo+zsbGzevNn6vWXLFrRr1871+fLYNm3aYPjw4Vi6dGnYy0gIIYSEgvQUFRQUMOECIYQkOBEzjPr27Ys1a9Zgw4YNOHLkCGbNmoVhw4a5Ovfw4cM4ePCg9f/33nsP+fn5kSoqIYQQEhTSU7Ry5UpUVVXVybBKCCEkcUiN2IVTU/H0009j6NChqKmpwY033oiePXti+vTpAIDbbrsN27ZtQ1FREQ4cOICUlBRMmTIFP/zwA3bt2oXhw4cDAKqrq3H11Vfj3HPPjVRRCSGEkKBgam5CCKk/cIFXQgghJIwUFRWhvLwchYWF/CYRQkgcEpMFXgkhhJBko7y8HJMnT0Z5eXmsi0IIISQAaBgRQgghYaSwsBBjx45FYWFhrItCCCEkAGgYEUIIISGirmO0fPlyCCEYRkcIIQkGDSNCCCEkRLiOESGEJD40jAghhJAQYXY6QghJfGgYEUIIISFSWlrKdYyIK9SwS0JIfEHDiBBCCCEkSjDskpD4hYYRIYQQQkiUYNglIfELF3glhBBCCCGEJA1c4JUQQgiJApxDQgghiQkNI0IIISSMcA4JIYQkJjSMCCGEkDDCOSSEEJKY0DAihBBCQkQNn2PqbkIISUxoGBFCCCEhwvA5Qkik4fzFyEPDiCQE7AwIIfFMQUEBzj77bBQUFMS6KISQegoHYCIPDSOSELAzIITEMytXrsSiRYuwcuXKWBeFEFJP4fzFyJMa6wIQ4gZ2BoSQeIZ9FCEk0pSWlnLuYoThAq+EEEIIIYTEISUlJZg6dSqKi4tpFIURLvBKCCGEEEJIAsGpBNGFhhEhhBBCCCFxCMN0owsNI0IIISQEmDWTEBIpuC5adKFhRAghhIQAQ10IIaR+QMOIEEIICYHi4mKcc845qKmpQVFREb1HhBCSoNAwIoQQQkKgtLQUHo8HHo8H5eXl9B4RQsIGQ3WjCw0jQgghJARKSkpQU1OD2tpaAOBEaUJI2GCobnShYUQIIYSEwNSpUzFixAjU1tZi3LhxnChNCAkbzEoXXWgYEUIIISGQmZmJ2bNnAwAWL17MsBdCiJFgwuKYlS66eIQQItaFCBd2q9gSQgghkcLj8cDr9aKmpgYAsGjRIgwdOhRVVVUxLhkhJJ5IT0/HwoUL2T/EAXY2Q0Q9RgsWLEDXrl2Rl5eHRx55pM7+n376Cf3790d6enodS9jfuYQQQkg8UFhYaBlFWVlZDHshhBhhWFz8kxqpC9fU1OD222/H+++/j+zsbPTt2xfDhg1Djx49rGNatGiBp556CnPmzAn4XEIIISQeGDBgAMrLy5GSkoKdO3daRhIhhKiUlpYyJC7OiZjHaOnSpcjLy0Nubi7S0tIwcuRIzJ071+eYNm3aoG/fvmjQoEHA5xJCCCHxwNSpU5GSkoLa2lp4PJ5YF4fEOUy/TEj8EjHDqKKiAjk5Odbv7OxsVFRURPxcQgghJJrIsBiv14sxY8bEuDQk3mH6ZRIINKSjS8QMI1NOB7cjaYGcO2PGDBQVFaGoqAg7d+4MrJCEEEIIIVGE80xIIEyZMgW1tbWYMmVKrIuSFETMMMrOzsbmzZut31u2bEG7du3Cfu4tt9yC5cuXY1G6ZRMAACAASURBVPny5WjdunVohSaEEEICZOrUqfB6vfB4PJg8eTJHd4kjTL9MAkEIAY/HY3QakPATMcOob9++WLNmDTZs2IAjR45g1qxZGDZsWMTPJYQQQqJJcXExampqLAWGYVKEkHAxduxYeDwejB07NtZFSQoiZhilpqbi6aefxtChQ9G9e3eMGDECPXv2xPTp0zF9+nQAwLZt25CdnY0nnngCEyZMQHZ2Ng4cOGB7LkkuGFdLCEk0+vTpwzCpJITfKxIp6GGMLlzglcQtXAiNEBLvlJSUoKysDCkpKdbG/ir54PeKkMTCzmaI2DpGhIQKJ6gSQuKZoqIilJeXA4CVqpv9VXLC7xUh9QN6jAghhJAg8Hg8PpOiZbpuhrwQQkh8Y2czRGyOESGEEFKfKSws9MkUVVNTg8mTJ3O+CSGEJCg0jAghhJAgWL58OdLS0pCS8sen1OPxcAFPQghJUGgYEUIIIUFSUFCA2tpa63erVq1QXV2NQYMGcb4JISRk6IGOLjSMCCGEkCBZsWKFz+/t27fjww8/hNfr5VwjknRQiQ8/9EBHFxpGhBBCSJDIBAwq9BYRJ4qKiuDxeFBUVBTrooQdf0o8DafAYcbD6ELDiBBCCAmSMWPG+BhHXbp0QU1NTYxLReKZ8vJyTJ482Ur1Xp/wp8TT+xE4XOA1utAwIoQQQkJArmEEAKtXr0ZKSgqmTJkS41KReKWwsBBjx45FYWFhrIsSVkpKSjB16lQUFxfbKvEFBQU4++yzUVBQ4Pda9CyRWMB1jAghhJAgSU9PR3V1NWpra+H1en28RfXo80qIX9LT07Fw4UIMHToUVVVVQR8TyHHxhBvDkMQPXMeIEEIICSMlJSWWUQTAxyjS5x2RxINei8BwMxfG7XyZcM2rieY7ZJhg/YAeI0IIISQI0tPTcckll2D27NlISUnxSdudkpLCuUYJTiJ6LYgv0XyH9BglFvQYEUIIIWEkMzMTs2fPBgAfo8jj8WDs2LGxKhYJE5HKBpbMnqhoP3s0M7oxSUL9gB4jQgghJAjUcDmPx4MGDRpYo9PFxcUcPSZGktkTlczPTuILeowIIYSQMJKVlWX9XwiBI0eOWBm3ON+A2FGf16Xx5xHKzMzE2WefjczMzCiXjBB30DAihBBCgmDv3r3o0qWLz988Hg/Ky8tRUFBQb5VfEhrBhlwlQgievwGB7du3Y/Lkydi+fXvY750I9UPiH0fD6Ndff8XRo0et3z///DMmT56Mf/3rXxEvGCGEEBLPFBcXY/Xq1T5/E0LA4/FgxYoVnG9AwkoieCH9ecMiuYZTqPVDw4oAfgyjc889Fxs3bgQArF27Fv3798f69evxzDPP4P77749G+RISNi5CCKn/LF682Oe3nHMkhPBJxkBIOEiEEDx/3rDly5dDCBGR+eCh1k8iGJ4mqHOGF0fDaO/evTjxxBMBAC+99BKuuuoqTJ06FfPnz8e7774blQImIonauAghhLinvLwc2dnZABCREXCSGERLMWXWM2dCrZ94NTz9yRd1zvDiaBipGXc++ugjDB48GACQlpaGlBROT7IjXhtXIsMREUJIvFFYWIgtW7YAALZs2QI1ySsNpeRh8uTJEEJg8uTJsS5KvSEW3/x4NTz9GT7UOcOLo3XTq1cvlJSUYPLkyVi7di2GDBkCANi3b19UCpeoxGvjSmQ4IkIIiTdkOJCcTM4Bw+TE4/FYc8skHMwLjVh88+P1nfkzfKhzhhfHXvy5555Dq1atsHHjRrz33nto1KgRAOCHH36IO8Eh9RuOiBBC4g35HRw7dqxP6u4BAwagvLw8LpUsEn769OmD6upq9OnTx/obB/NCIxbf/KlTp+KSSy5BWVlZXLVbGj7RJeAFXnft2oWWLVv6jIzEC1zglRBCSLTwer0AgNraWni9XowZMwZPPvkkqqurAQCpqakQQli/Sf3EtGhpSUkJF/jViPc6KSkpQVlZGUaMGIE5c+ZwAdp6TlALvH755ZcYMGAALr30UqxYsQL5+fnIz89HVlYWFixYELHCEkIIIfGOmn2upqYGZWVlPkaQDLEi9RuTd4Oj/HWZMmUKamtrMWXKlFgXxUhpaSnGjRuHOXPmRNVTFa8hfMmKo2F0xx134IEHHsBVV12FgQMH4vnnn8e2bdvwySefMF03IYQQ4kB1dTXGjh0b62KQCBOoEVRfFGF/z6HuLykpQU1NDWpra+N2sKCoqAhlZWUoKCiIqkHLsMv4wtEwqq6uxpAhQ3DFFVfg+OOPxymnnAIA6NatW1QKRwghhMQr/hQ8IQQ9BklAUVERPB4PUlJSXBk78TqXJVD8KfTq/qlTp2LEiBGora11NVgQC+OxvLwckydPRnl5edTuCbibT1VfjOmEQDhw0kknGf9v+h0PFBYWxroIhBBCkgQAjltWVlasi0iiAADh9XoFAJGWlub3+HHjxgkAYsSIEa6Oj1fGjRsn0tLSxLhx4/z+1vf5Iy0tTSxatCiq9VNYWCgAxKUuGYv6qO/YvWfH5AterxeNGzeGEAK//fablZVOCIHff/8dR48ejajRFihMvkAIISRamJIQdenSBatXr8aiRYt8JuOT+ktRURHKy8vh8Xhw1113ufISynMKCwvrjd5iSkIRLPWxfkIh3hNXJCJBJV+oqanBgQMHcPDgQVRXV+PAgQPW73gziuIJujwJISS5kGsYrV69GgBw9tlno6CgIJZFIlFi+fLlViIOt0rrypUrsWjRIqxcuTLCpYse4UyxHc76cauTxavuFmujKF7rJWI4uZl+++03MXnyZHH77beLZ599Vhw9ejQgN9X8+fNFly5dROfOncWkSZPq7K+trRXFxcWic+fOoqCgQJSXl1v7OnbsKPLz80Xv3r1duzXjxf1JlychhNR/ZOhNVlaWSElJ8Qmj4zeA2DFu3DiRkpIivF6v69CyQEPRwn1+OK4ZyPHBlNeuXt3qZPGqu8W6XLG+f6SwsxkcDaMRI0aIa665RkyfPl1cfPHF4s4773R9w+rqapGbmyvWrVsnqqqqRK9evcSqVat8jpk3b54499xzRW1trViyZIno16+fta9jx45i586dru8nRPwYRpHogAghhMQXcq4I/v+8BNUw4jeA6EjdwOPxiNTUVOH1el2fG6pyGgnl1uv1Ws/hRu+JtIKdlpYmGjRoIFJTU33u4VQ20zyowsLCuGq/sdYpY33/SBGUYZSfn2/9/+jRowElXPjiiy/EkCFDrN8TJ04UEydO9DnmlltuEa+//rr1u0uXLqKyslIIkdiGESGEkPpPWlqajzEkFUWPx1MvFQkSGtIwAGB5GBPZY5SSkiIaNGggUlJSXBk9wXjKAiGY6+vlVgc7PB5P2MsYKvXVSIkFdjaD4xyjBg0aWP9PTU0NKESvoqICOTk51u/s7GxUVFS4Psbj8WDIkCEoLCzEjBkzAro3IYQQEmn0OURjxoyx1mlp3LgxysrKUFRUFKPSkXijuLgY55xzDgCgtrYWXbp0cb12TagLxkZiwdmxY8fC4/Fg7NixxvlF+tyU0tJSpKam4oMPPojImj2lpaWoqalBdXW16+dUy11SUoKysjJrnwjjekvhmqfDNY8ij6Nh9O2336Jp06Zo2rQpmjRpgu+++876f9OmTR0vbBIoPYOP0zGff/45vv76a8yfPx/PPPMMPvnkE+N9ZsyYgaKiIhQVFWHnzp2OZSKEkGQj6SbORgmZNSszMxPAsSyuwLFBxBEjRmDv3r0YMGBA1NdEIfGLNAzURB3hSFQQK3Rjq7q6Gk888QRSU1OthAFSiZdrPaWmprpKTBKtfkt9BrnWkqSwsDBs9/Fn0Lh93nAmuCA2RMpFFWooncpDDz0kHn/8cb/3ZCgdIYT4Ul8nzsYaACIjI8Mn2YKagKFRo0ZxuyYKiQ0y1Av/P5QukDlG8Y4q+3KOjxr2BUBMnjzZdWKSYPutUELNwhmm5m+NJ51o9tMMxztGUHOMQuHo0aPihBNOEOvXr7eSL3z//fc+x7z77rs+yRf69u0rhBDi0KFD4sCBA9b/+/fvL+bPn+/3nvwAEUKIL/wIRoasrKw6yRY8Ho/P3+JtEjeJLeoco1gbzeE2IFTZ93g8dRZ1VduLvzlAocxFMhkYbpIqyHt6PJ6wzIEK1NCJZj/NwbJjRN0wEuJY1rkTTzxR5ObmigkTJgghhJg2bZqYNm2aEOJYuu7Ro0eL3NxckZ+fL5YtWyaEEGLdunWiV69eolevXqJHjx7Wuf6gYUQIISQaeL1eHyNICOHzW46gUwFJTkyKruoxGjFiRMjekFCUaVU5DvQ6JsVaPpdMPKIfZ5cxzu31/aEaP7qBoxqkdteV5VM9XqEQyVTmoRLPZYsmMTGMog0NI0IIIdFADR2S2au8Xq81cu71ei0lLVJZuEj8Yqfcp6Sk2Gak8+cp0a8Zysi/quxKmXWbhc3J6JNl11PZB+IFCkYRT0tLEyNGjPAZnJADFE4eIyeDKlrIupKyEc7rhsOgqa8eJhpGhBBCSJhQ1y3SQ4dk2JBUxOqjUkGcsVNK1bV/9GO9Xq+jVyWcHiN5vurpUaedu722ndGRlpbmM58q0l4H1biQz6IaYnb39tc+o+EtUevKNO8s2PWVwtX30GOUwNRXw6i+CiUhhCQqTqE36oh1VlYW++84JlLfVzdeFYlUYKXHIJLr/OjGi5Rhuan73KxLJMMC9TA1feBAfcZItQf1nmq7c3oWtU70Mnu9XuH1el0bF8HKkr/5TW5CAcNZnmSBhlECoy6ilqywgRNC4glViVKVGTWESFc2SfwRKY+e6bp294rE9810TV2XUD1G0sMiy+amTGr4mj+Pkdpe3ITtBRtOp3vcAgnhk21ZrY9ADDn9/Qb6DHa6XrAeo3BSH3UwGkYJjMn1nmwwHIUQEk/Y9Ul6Uob09HRmqItjIukx0hVy9V6RVjRN8mmnS6jzntQ5Qv7K53SMOkAgQ0ulESVD3ZzaRLAJGPQ6DyRUTs80GeighjT+pC7q9hlkGUz3DUROIilTpmdJdGOJhlECk+jCFw5YB4SQeEGfWK7/TTeKunTpwgx1SYiTMhlIiFYw2IXymb6jepmkkSTl1o2Sbpr/pBpD8l9pfPhrE/K4rKysoJ7ZrZdF9dKEGu4XrMdIDZUzpft3KyeRHEA2PUuiD1jTMIoDqNwTQkjiYwrZ0Ud87ZQc9v/1g2A9KmlpaZZRIJN2xAp9bouUUXVNoi5duliGgpMxpyvJ6nPaZadTk5ToAH8sCOsWNbuek5GjGk3+kjSYDD679656jALR9+R5WVlZllE6YMAAyygN1MCNlkwluk5LwygOSHTrur6S6I2bEBJdTMqdKQxHJmZI5jDo+oqbub92XhQAIjs7OyClPxAFPRBPhSmBiDrvSBoaJiPDKTRQn1OkztlRw/bC5TGSdSuNI1kej8dT5z2pHhrdSNWfw2Tw2ZVZDVUMxtMj69dkMKpGl3rtUEMekxkaRlEkXqx54g4arKQ+w34nMqjzAuTItimcLt7rXp8XQdxhN19H90aoC7mqCQkCrXOTgi7D03T5cjuvxm7tHmn4y/LbhZY53cfOc6p6o5zqwV8iAlNZTJ44dW0x3aBTPUbSMFSNFLu6sWvTapmDmRvkZCjq/YwcmPF6vbZyoNYNdZy60DCKIhTCxIKKI6nPsD+KDFIx1pU7XemL98VdgwlZigXx1k87KejSG6ErrNJwDmYhT6c5PLrnwO1cIn99gt01dcNK9ZxKo8RkFEkjS9aBPsdIvYed4WnXn9l5iPR6Mj2jKeTV35pSbt5RoDjNi7JL7KAasP4M4WDK5XauViJCwyiKxFsHTghJXtgfhQ/TSLOq4JrCh+Rx8UqieIwiZeCHu32o3haTcazP6wnlvmrZTfUj3600VnRviDqvyJ9HRL2+6f/yOaWXVF9TSJ9nZFK2/YWIqZ4Sk9Gn17mpnkyoZdfnQunzr0IJWXO7384LJOu0UaNGPu81mEQVbpFeMJOhn+jQMCIJBxU6Qkg8YVISVCXNLlTIbmI3cU846s90jXAZXLqRYmcc+0sJHUwIlp2XSFf4PR5PnYVLpdfGLsGCbujpiQXk/6Wyrnpd5H61DtQ6cmuUyX3qnCD9XN0YDURO1Hcikx5I1JBFk1dGLbOa/EHdL405Wf+yfPrz6l4gPSug6t1Vvbz+vL7Bth3VwynLLZ+lPkDDiCQcDAEihMQaVakweYO6dOli9VF2HiO7Ef1A709Cw/QOwjXarns7THNrTH+XSreq7LtdBNafTNnJo2rUq6FjurGvKsYAHMPK1AEB/VrqvRs1auTjZXJ6VtX40Q0U/Vzd4xMI6oCHer5qFOpGiukdqM+r7pcheaos6J43/R2bQgClrMpryXLaeX0DDZ00nWf3HusDNIxIwkGFgBASC0whSvqItN3otMlbIBXvYPqzcA8Q+QtHqs99runZTcqsenww3ht5XX1LS0sTmZmZdWRCHfF3CmUzzVnyN69IPU836PVrqfvl/CFd7m2VSc0ItGsr6jPapes2GU66saQatIGGvKnvTH0+1dOjl0M1bvWQWlkOk4HiFFqpP5M/D1ogRrMQvgZfIO1a7fN0A0m+2/oADaMEpb5/qAghJBoE0pea5jp4vV5rwVan0VM7z5KuMLr9XoV7DpA6gq0rWHYT3uOBcH0LdeXSFP4kcco8ZzIuTWFV0uMij5VKvDrHSJ8LZPf8qvdG9XKp+9xkybObsK/LtZv0007ybrfpCQR0j4kalmc3D0s3WtTjA0mTrXqKZJ3KMsry6IaeGoaovhNdBvQ60sP0JOr1TWVXDTDVqFIz06nlkfXrJNtO6H2OLL806qXHKhJzmqIJDaMYEuyiX0IwnIwQQsJBIH2p2k+riy/qiqnJMBLCPI9E3VRFKpzldoOqyOrfTDdr88SKSMwFEsLZ8LSrDzvjUjeoTam57dJVmwwVHTUTYp8+fXy8L7oBrnsl/P3fybOjtgO9ruS9nTyqdpv0mKnhfFLBl/Xkb+FVfTFak9Hob96SPnAhDUR1bpApfbmeVELWke4V1P+me7b0e6vZ/lQjUX8vulGnH6PWje6RcqoXU1tTDTEA4vbbb3fdh8UrNIxiCPCHmzzQzp0eI0JiA9tefBLoezHNV3B7nlSGpIJhMnhUpdluhNlJeXG6v6qwhVMOAw3JiQciVTZ/8zxM9zR5jFRvkFRu1Xctr29nJPvL+OXkmTEZJaqirYZFma6ve4rsQvL0wQCnMtnNb1KNA3kd3fujZl2T5bbTm0z3kZ4ZPQzMzrCya68NGjTw2Zeamurz3nWPonpvdVPn6ajvSj6T7iHz12/IdZLktXr37m08TqYul2VX72ln/LgNcZT9YSJDwyiGhOIxSmbcfKBIaLBe7YlHby3fV+CeH6mYBPoe3SgourJlN5/DbnMqk66IhvOdx0KOwnnPYK/l5HnwN4/D371VBdokA1lZWY6GhL/nUT0EUsmV/8q/6wZSly5dHLPNmcL+pHKul8Vk6DnJtn6O9JgUFhb63E/9vqtGnNt5OXZGjcfj8VHi1SQp+vsybY0aNTIanKp86O3TNHBi8viomz8D0un8wsJC2/16/6XPNTLJsdp/qX2mP4MpUaFhFEOozASH2gHFo5JaH2C92hOP7TaeQ52iRSDvRf3QB/IeA5k3oU+0dnuevzLpClMo837c1lkkZT7UvkYtW7CRF16v1yfds7/7BPIcujIq56MBEBkZGX5lwV+96/NMVMNA9VpIA0oq/urEftXbobYNOyPBKdTONL9HbqmpqcZzdG+arnyroWtSkVfXRpLPJ8uths6Z2pbuEdLDJ+227Oxs63j9WDXpg/p3aYC6DSlUZVHNWGeqs0D6FJNRpNe/P6+knphDTeKgz72Mp+9joNAwiiGhfBDiWfgiXbZAPEbxXE/xDOstsYjnyfHxSLDynZaW5qh06QqUfk83you/MBR1Pok0ht0kYjD1m/rkbDVDXigGRyCEmkRClk1dnNTte1XPBewnwbtB91jIf00j904Kr1SmVUPDqd6lweBk7KjhVVIRlkaGuhaNurmR82A2t0aClEXVUDBtdokAnI4VQvh4n1Qjz42xqteXk5FiSsfudA95rCn0LpTN6Z52SVdU2barbyk/esKHRB5UpWEUQ0JRPuN5RF9+uGVO/lgp16GEy5DEJFkNumR97nBgmhei7lPrNRDPT6NGjercy43C5eT1M41oy7/ZLeSoG0Gqp91OSTItihnPHiO7bGp2qIkzpCFhUp6DQR1Ftwtlc7v5S+WuyqNUztV01tJbIb1H+vXlYIo/5T5cm9t7yPpy8j7JZ9LnwtjJtHy3UpalwZqZmWm1CdWbF8gz+DtPNT6l4WYyVKTs+TOc7TZ/4XlOm9OAgvQ2qvVt8ixKAzaQwYl4hIZRjJCduRy58efKNJ0fr4qQOhKjZmqJdgx5sOEy8U48v/tYE+sBA33UOFrvKBiZoBwdQ800pYeCmDxxbpQMO8+Dv1AdVcExoV5bVWic0jqrXhHdm6B7jOy2SMiK2lZMCQsCuZ/e1zudL5VnqcyZ5tkEWgb1ex6MIeBGYfU390PNWib1CfX92mUqU9fZcdoi5UEKRsmXHlngD2PH33lqwhRd5k1hgU7GrJ5yPZDncnu8zO7m1CbDLWtykV1VxlTZNrVRr9frE6JZHwaiaRjFCGmBmzrBREd2suooSrjCfAIJGaqvil+slf94Jlrv3O4+6kh8uN+R07MFIxOJJkdu322gMqAO5MiwJOnxVvtmIdwbNlJhM5XN7fkm7O4vFRQ5Gi0n9KelpdVZfwWAz3FpaXUny+vPHs5BJpOHR58j51Y29frQ5z6YvG/6PAs1IYKebtpNdi41pMj0bkIxKFJSUqxRervy6EalHlqoDrzayY9bZT1Yz1ewm7+EJVIXcJsQxS6cTKL+zSm1vr+0+6b3b/IambZwh9C53WQGO3WunSrbpiyb48aNq5OaXK2jRNXBaBhFET3G267jifS9I31tp8YXzPXU3/oHL9yLHIabSI60JlpnEyjx/Jx2ipseohPOsvu7Z6DtIJ7r14Td85v6VX2SsAndY6GeZ9d3uVUyTGUMVKl0KrPpuyH/37BhQ+v/UrE03VuOupuup26NGjXyG5asy5JTO9DnmBQWFtYZ8NLfqUlOTd9QidMClqa6UD1GuvEgvXD6/dTUz/4UeCdF2N/m5DGSz6N+E3XDU/1X9YKq93Ba90bWQbDld7qe07yXQAYhQt1UmbVrL6E+q9t3rcuXurkZWDENTLvddLnQ60Ea2Hblln2J9HYl2uCbhIZRFFGFRHfny0YQqRCcSAqo/hFyanj6c5kMG72s6m/9g6s2Xp14UPwiUe/x8FzRIJ47Vbt3oH9sw/munO6prk0Sr4Q6iOHPSycVbn/zTdQQEDVG3p8CE4i3Rw+F040Bt8qaTlqau+QPgWxurufP+Dat2yLnPOl9tynUSf/uqZ4Qu/V2TAadXXkkpox+48b9kWFM/Zbp1zb9TY6K273bzp07h/x+/PUd6rOqsi3fq744qbooqd37iOTm1rsTSLlCNdxMbdtUzkDbcCgGsV059XKpfVcgSwIE+66c6kC9Pz1GcU68GEZ6Z68Llfw4REIZDLeAqtdTn0EIZ8NI/7gBdScNO3mM1FAlvYHqqHUciwXHTGEOwV4nUjISzx1XtL2B4agLXR6jYdzpo3xuyh+L9x5oGVWcZEF9FjeJX9RQI7fKRDCj13oZQzlfEkgITyQ3f+9WllOu+SL7YtnXmzZ1fpQMNZcLZ+rzgOyMEZNM2JXT9F7tFGx1TkugW7AT6dUyOWH6Psq60etIDeeSAwfRlBtdftUws2jc38lD5c8Qi7QB6c+4MxlZMuzX5NkJpSyq3uK23G689IkCDaMoo3dceqdhmuRmRyyVWlXhU5/B6/U6frzVtRX0Rqc+hzqCp44mOo3O6Zj2R7POwqUURzJDVDx7ZaJdNjehWk6YBgnCZRzblceksPurL1lOPW480gRSRokaOmaXdU0ep6cjVutbfQ+yP3FSWMOhqNk9ixtlwyn5QjxsgZbLrZdLenhMySTUaAHAbNCqoZHqvBp1XpDc7EKNnL5fdoNyTls4ZMmu7zANIqpr/vgzyuLVWxTo5q+O3RoK/uTUTaISp80ks2rKcH2fnF8mZd/0nFLO5UBCuOpUT2qh7zMZ1abkL4kIDaMoo3b4pg7YaXVhpxCScJdPzVVvurf+kVIbcjgapKle1A5CHz1RFQlZZtPkaT2cMRwGhlMMfDiur39Mwkm4ym5SPCNVr+FG3seu3OocAqdnc5qsr6+uHgpu5jQ4GWKqPPkbjdYxGXpujD99EUrdaFEVWtNEcn8fXrtQNXWxS6lYyPcAQOTl5YVNkdA3J7l1c75aN6qHPdabyeh0c55sA26VY3VtHv3d2hk1+lyZQLyC8bzZ9R0mHUAdfE1L+yMttSkES76XWD+f3bsPZPMnV8F4ZFQZkpvUNYLxAsuBXjfH6ovumtqBeozqmQ3XO3DyKOr3Mq2FFM8RKU7QMIoyumCbGqPs9GX4l7ougipk6krL/nAjoLKhy/lO8rceqy//HusOVe9gZNiG3fG6AufPsLRT+PS6dKs4B4PJQxaJNQLUZ9LfuxvZ0RVTfwa7fj/9+vp8hHDPS1H328m4/jGSbbVLly51QkJNHzsh6iqysq26bY/6MU4eU32zq3+7EFSnMsl9Xq/X53z5dzUVrto+JWq/onvEnQw9fZ0Z/T0FMxBj52kL5+YvZa2bkDhZN2pbimV/q26qpzHa4VjJtjmNwpvarK4XpKU5r89jN6E+Uls076VuhYWFxnbvxgiTx6hJBYIxjPTNX0pw1fsH+KaiVxfojYRhFOhminSK54gUJ2gYRRlVyP1lglEbjvxXn+AqG4ITqpfKLnRGHqM3LKkM6qsdJ+KmT+LVVyjXlVx9kqrJfawbEkBdd7ueVcjufnYZnNRFJSKdwQAAIABJREFUGdX1ApwMJCcviOljqt9HTdmpe9l0Y1HP5iNHmtQPtN399InyptApuTlhZ8TafVTkMfr7Vcsn36c6eKG3WaeEI3Z/l8aW+lEzvSuTwe1vFXj9PnbGlS6fbgwFk3Ll71wZHqz/XZ07onrV7J5D3a8+i3qc23oJZ7iJaXNjzLtRYExtNZLlDmRTkxpEuj5jsQWT0Ssade0GeY6UITcho/Eydy1SmzqYZTKEnPoO0/pUoQwKt27d2u8xck6b6b2onnr1mxntVOp29etP10gEYmIYzZ8/X3Tp0kV07txZTJo0qc7+2tpaUVxcLDp37iwKCgpEeXm563NNxJNhFKmRStVal8JoalRSaVU7S3/GTixHISK1qQYf4DvC63SebvRIxdjtew0krWu4UqSqCqVaTrfXl4aNSZ5kZ6wbDHYfY6dr2BnngZTHzXM5KT2qsex0fzf15k+W7N673fsLdgt3ZiS3WygT1u02GcYS74qcP4+RW/kJ5rxobHbGKbfQN1N7d0JXPgPJEgtE1ijSvwMyU1ks6tXj8Rjr1u23W54brrBMf9cxvRdZBnUdLSkDsahT9TlMiyQnKlE3jKqrq0Vubq5Yt26dqKqqEr169RKrVq3yOWbevHni3HPPFbW1tWLJkiWiX79+rs81ES+GUTSFN9zrDnCz36I9gTWYLZ5GQJNlE4JKYyS2ROnbnEIS/T2DnWc/1s8kt3gsUyy2aBnoquyY5kA6efRjXUdyk3PM4in0slGjRj4REv42fc5iMJvJA2U6Rq8nu/lyQoiw6iCB6gr+1tlKRK9R1A2jL774QgwZMsT6PXHiRDFx4kSfY2655Rbx+uuvW7+7dOkiKisrXZ1rIl4Mo3jqELhx4xbZTYj4UUq4mbdwZ3HS37/pG+Ams5XthzkO6gzwzYwZ67Ikw+Y0J1ZNmGQKv9WVcDch/JHYVAU61vUJ/DHn0c3iznKTMq8axMF45OXaQ07hb4WFhT4GjzoNQJcNILyDn02bNvX57fSM/jxEiTjPKOqG0VtvvSVuuukm6/fLL78sbr/9dp9jLrjgAvHpp59avwcOHCiWLVvm6lwT8WIYxUuHwI0bt8hvQsSPIsstsptpCQa7b4CcK+YUShMLj5E6n9DpOK/X66PsxLru69tmWlspJSXFdi0v3WME/OFVkGt66deLhtfVpKg7JauJ1eZUF/o+tZ1Lw0qtb1Nd69dUDTA1k6v6buT8U709yncvDbPMzMw6UyMCeT6nTfem2YW4+/ME0WPkgtmzZ9cxbu644w6fY84///w6htHy5ctdnSt59tlnrSwkHTp0CPNTBE+sOwFu3LhFZxMivhQAbuHf1PT0gaTV95dd0G7CfSSfRZVXf8qUmigk1u8gETa3yqnMQKcbnKoxqk9y1zN8qpnMTN4htwt3RmJTM8JG617AHwaIut80MCET/qiJceT56nFqciJ/fbw/j5L+HtV9at+ie1305EbhTEmvh4jayW9qaqqxb0s0Q0iHoXRRJhadUay3QHL3c+NWHzb1Axate9bH7GDxvqnYZUa0w0lZtjs/knNahHBegFzdYiHfibqZFpiV715XZqXyqx+rzuFwmuSue49076Rc1DlQw8RfQho3m3pvkxIfbi+WrEM1+VCg701Pgy0NJ1OyC7tndiqDHAAxGWCqAWzqV+T58rqm5RSC2QJN9uNvzlsiGklRN4yOHj0qTjjhBLF+/XorgcL333/vc8y7777rk3yhb9++rs81EU+GkalziWWaxUhuer3TOOKWLFu0s3ZJ4j1bWyJvJsMzWHQFSZ887W99q2CfQU+vrqfSVxUZpwyWavnsFt0GIBo2bOhYnmh5D2K5uVmvS26qF05mddS/o7qyqXuM9P+rmTb10C2796bLpjwvlEn+8r6mNc8i1W/J+g4mEYU05PSU5voyCE7XkKGMdsuBmJaNkPeWhpCeZVC9hvp3aUC5CQ1UjwvV06R7t1W5k9dPpPlFQsQoXfe8efPEiSeeKHJzc8WECROEEEJMmzZNTJs2TQhxLF336NGjRW5ursjPzxfLli1zPNcf8WQY1YfQA1NHqo5q+BvRUieH2hmF6ihZIGXTR+dkfv1wGZ+hxmcLEbiiHErZTaM/MkxAiOgYq3r57Twb6kfX7do6wZQjEEXM7bGyntW6lfhTJkLNKKT2b3Z15k+GnFKkm+bDeDweVx4qr9drtUXZdlRlIdEGS9R5BXbrRLlBfR+63Eu5s8OfTNqlSbdbt0z9NrkJ0XLKtqcq5PK5srKyjH2mG7kNRpZVGY21vMh3qT6f3Xp3qgLtJpGH+i51A8DuHXu9XmstRLUePR6PUUGWIXymOvc398TuWiaPZKQGh/W2qe/X+15VZtSFb6X3Rsq6fI9OGerUxamdDNdx48Y5huGqdaz+TV2fT3/fTn2EnDukrokXrGEqdatA5DIR4AKvUSaSawKpi1nKTkEIUaexy2PtUiw6LRamjjAGk57RtE91I5samduPpr9yyI8CAB9Xr77IaFpamu1Cnv7mBqj1ZDpOCHeGkexQ/CnN6sKochTItKCsOlKovj95DVn36qieqgT5G+10o4TIOlGNHvnbhMklr3fg6rM6ld9O6ZPHynrzer2W4S9XjtevrT6Pfi2T3OsTaPUwGSeZNSk3Th8c0wdTnbRt6jP0kT11MEKuOyI/tDKURAeoq5Cb+l3TBHK131IXCFbfqy5fbhWpSIQX6gsBB5NxSVVQ5XvwV3dCOPc/6gR9fX6JP6+D6Vn0326NQLs+3s6wMg1+NWrUyCgrejiZqR50z4S+iLPbwTY3AxZ270K2BbV9SaPazoujDyiaPAV29a22X4mp/oKRdT3MS5Ur1fh1U1eqbKht2l8yEn/Xle1cfx+6nqLe02lgVta7SW7VgVb922MaFNaNQTWkUE+aIetUxd97NLUxfzqTHDBW5UYvg6k+5MLceh2bsEsWkgjQMIoy/tIzOgmz2qhVhTIlJcVnsqBJKdSVZidMH1+5BkGkcPro+lOC9Bhsu+dzUmBNo6NOjd/0rlTXuIqqBKnnqiNKaielxhZLIy01NdVn5DGU0WpJICM6untc/xgWFhZav1Vl1MmDaDd6ZtqvYqcwBKuk6mVSlZhwEMizhXqukzKlK4QmhdOpPE4fOr092L0DO0+Fv/7JbrDGKXTETlkKdjMZhf5k2A61n9b7EJPXUSL7elP57AZD3PYPdoMHTka6iUDboT7oovaBdoMNdoNedmE9pj5Ibecmz7o81s3AnEnGZDtJS0vzmaxvV++qseT07CbcGLpqPblpF/p+PdGImgjC34CvSQ51vUX3YrjZpKLuFMnh8XiMWft0b498N2q5TNjJmiordsa8rFM1CYWsF7XO/WFXRrUt+TP8/PVZ6rGqd9DNoKYsS6Kl6ZbQMIoyJqVbpuMU4o9Gp+a4N414OX3AQlWYJaaGFAv0EWW9Y7Qb4QzHfe06QImbug7nSGs4nzFYOVHLYOfRkKO5pixKpnKoHyt/+FOCQpXTRB7pchqZl14fufm7TqB1aYqdd3tdf3JtGpUF6maPUvsqXSFNS0sLas0Rk7JrIhiDQB1JVlMs2/U5qiKnGhPq6LaTUhcKofRjbo83eYTc1KnbQTH1HnpWL1l/qjdAPd+kZKoDWLr32e332J+nLpB6dTPY5M+QMBk1ds8q5dBJAXcTfinRjVB1oMdJH9GNa/18fdDTaf6Ov3q2+z7YfcfUgVCnASC3fbNTGeW95EC23i+66cck6iCMei+1jTgNvITrexwLaBjFEHWURc1OEqqREy7lWe2MYincpg+P3rGHw4Nid994atjxUCY3H99ARjwD8Vzp9w9mv9tjEhUnxcqtFyHao32BvA9VuXMTLmn6bRpkUQ2OYDzm4TAI/LUZVVkxhe+4Db+KdwLtQwJB7W/0tuDPI68r7qZyBdp27ORU9/q5ua7TMbrRqYdPOxlHahpxdVDM1G/ryniggwV629QXQ5VhiqZ7mrxNaltQPWb6JufeBNs3qgMxqrEN/OFt1vuecA9o6+8kLc3X0+5v8FEtgwyn00P71Dr0990Ohz4bC2gYxRDdEyREeAQnXMIXr0KsKzh28x7qC/H6HsJFuJ8vVAUi0YmnPiRShPL+7EZR1SQR0lhRla1I1UWgRqFUmtWySeU5kh6j+oKqnKoRB0L4eqPtvOGBKoPBkpbmG0Yd6oCPWnZpQLuZfK9f32kup36fQA1btW3q8xBVg0Z6fk1GqzRK1CxyEnVgQeoNAESrVq2MhpQJ1dOo14tqLEjZUY1tafgFYzS6rT/1nejvy9+7UGXenwHnJuuc3l7UiKh4/r7QMIohJtdrfVbYwoVaR4kc9uQWykRgJLvHKBkI5f3Zhd0Af8w5ieYIZyDtWx31l8+gfj/qk1xH8lnsvKemeS8mb3g0+pZAPLxur6cq7uraRrqHxSnkSve4+ctuGGgZTc+seonVNNPqfVWjxy4MUzVg5LOawmudyq/WmV6v0gtnFy6nGuVqZs5IDgwGM/Di1uvnT//SdVw1DC+e03jTMIoxJgu/vnzYIkWy1VGyPS8hkURtT+qHXY/PjxaBDO5IpcdfcoL6QKQHhII1fPwRzwNZ0tshw8bsvET+DD/VEAilvvT61r1kKnpYmB5p48bg0BV19f6qQeg0z0/3ppnC1/zNl9Sf2c3x/tq63XGm6wdi4DslcnEj6+q9dC98vPZfNIxijMlrFC9QIY8/+E4IiRyxaF+BKNLJ4CGXRPpdhBqOaVc2u33x0Hfroed6UoOMjAxXYaPheBaT7uOUwEX3YuiRNm7mqdol11CfyV9mSP3Z/f12Uw9Ox0s59Rd6ZyfPpj7DSfb18tgNGAXr0VTPc1rzMpbQMIoxbht0LIjnka9khe+EkMgQK8U1kPvqHqN4+2YkEm48RnbvJph+OB76bl1JVsPP1JAyu1DCcGLSfZz0ITUMTp33F4iCbspoamcIOIWHRbPtBeoxMsmqXqdOz6jLqSnE2HRcIOjhmPHWj9EwigF2oRzxRjyMcBFf+E4ICS+yTcXr6KVKIpU1EdGVPTvlL5h+ONyehHCgKt2qN8Yu+UQg1wzG6+QUQWM3cT9Uj6vbd253fLDPHy1MderWY6R66UzpyYOV53iOlBKChlFMUIUynj1G8dbATSRCGQkh8UsiemHY70WGUMOiwkk0PUx6mJpUXoN59lDLbXdfu4yLgRhipnlUgb5zt6Fv8aLwmzxqbussXEagyRBy6wmLBTSMYoDJYxTtCb9uiLcGbiIRykgIiV9oZBATsZaLaHqM1CyHod4v2Lknbq4biMGi/05L+yN5g/Q6RYJIvLdQ5vOE4pkJlxHo5ACIRx2OhlGMoccoNBKhjCT+oNwQQpyIR4Ut3ERqXZlo151TuJiawdFfKu14/Sakpdln63PC6/U6pl4PlVDCJgO9RjShYRRj4nmOESH1lWRQegghwROPClu4UcPTwqmLRLvuTAPMgc7Hi+dvQrAeo2h4yOojdm0gBSQqrFy5EosWLcLKlStjXZSEpKSkBOnp6SgpKYl1UUgCUVxcjKFDh6K4uDjWRSGExCGlpaWoqqpCaWlprIsSMcaOHQuPx4OxY8eivLwckydPRnl5ecjXjXbdFRcXY86cORg3bhwAID09HQBQVVWFMWPGuOrr4/mbUFpaipqaGlRXV/utU1UnGjt2LFJSUpCSkhKXz5VwRNlAiyjx7I1JhlGpSBLPozzhgPJBSP2GbZzEA/UlesVtxrX6Sn3XiaKBXRvwCCFErI2zcFFUVITly5fHuhgkApSUlGDq1KkoLi6ulyN76enpWLhwIYYOHYqqqqpYF4cQEmbYxgkJH046QTK0tfquE0UDO5uBoXQkIajv4Q7x7N4nhIQO2zgh4cNJJ0iGtqY/P6cbhA96jAghhCQtHHklyQTlvX6SDF6ycEOPESGEEKIxdepULFy4EFOnTo11UQiJCKo3gfJeP0kGL1m0SI11AQghhJBYQYWC1HekMSTlnPJe/ygtLaUHMEzQMCKEEJK0UKEg9R3VGJKyLj1GlH1CfOEcI0IIIYSQJIHzUQjhHCNCCCGEkKSH4XSE2EOPESGEEEIIISRpoMeIEEIIIYQQQmygYUQIIYQQQghJemgYEUIIIYQQQpIeGkaEEEIIIYSQpIeGESGEEEIIISTpqVdZ6Vq1aoVOnTrFuhgAgJ07d6J169axLgaJMygXRIcyQUxQLogJygXRoUwEx8aNG7Fr1646f69XhlE8wdThxATlguhQJogJygUxQbkgOpSJ8MJQOkIIIYQQQkjSQ8OIEEIIIYQQkvR4x48fPz7WhaivFBYWxroIJA6hXBAdygQxQbkgJigXRIcyET44x4gQQgghhBCS9DCUjhBCCCGEEJL00DAKkQULFqBr167Iy8vDI488Ume/EAJ33nkn8vLy0KtXL3z99dcxKCWJNDfeeCPatGmD/Px862979uzB4MGDceKJJ2Lw4MHYu3ev8Vx/MkQSk82bN+Pss89G9+7d0bNnTzz55JMAKBfJzu+//45+/fqhd+/e6NmzJx566CEAlAsC1NTU4KSTTsKFF14IgDJBgE6dOqGgoAB9+vRBUVERAMpFxBEkaKqrq0Vubq5Yt26dqKqqEr169RKrVq3yOWbevHni3HPPFbW1tWLJkiWiX79+MSotiSQff/yxKC8vFz179rT+dvfdd4tJkyYJIYSYNGmSuOeee+qc50aGSGJSWVkpysvLhRBCHDhwQJx44oli1apVlIskp7a2Vhw8eFAIIcSRI0dEv379xJIlSygXRJSVlYmrrrpKXHDBBUIIfkOIEB07dhQ7d+70+RvlIrLQYxQCS5cuRV5eHnJzc5GWloaRI0di7ty5PsfMnTsX1113HTweD0455RTs27cPW7dujVGJSaQ488wz0aJFC5+/zZ07F6NGjQIAjBo1CnPmzKlznhsZIolJ27Zt8ac//QkA0KRJE3Tv3h0VFRWUiyTH4/EgIyMDAHD06FEcPXoUHo+HcpHkbNmyBfPmzcPNN99s/Y0yQUxQLiILDaMQqKioQE5OjvU7OzsbFRUVAR9D6ifbt29H27ZtARxTknfs2FHnGMpHcrBx40asWLECJ598MuWCoKamBn369EGbNm0wePBgygXBmDFj8NhjjyEl5Q+1jDJBPB4PhgwZgsLCQsyYMQMA5SLSpMa6AImMMCT083g8AR9DkhfKR/3n0KFDuOyyyzBlyhQ0bdrU1TmUi/qN1+vFN998g3379mH48OH4/vvvXZ1HuaifvPvuu2jTpg0KCwuxePHigM6lTNRvPv/8c7Rr1w47duzA4MGD0a1bN1fnUS6Chx6jEMjOzsbmzZut31u2bEG7du0CPobUT7Kysqywya1bt6JNmzZ1jqF81G+OHj2Kyy67DNdccw0uvfRSAJQL8gfNmzfHgAEDsGDBAspFEvP555/j7bffRqdOnTBy5Eh89NFHuPbaaykTxHqXbdq0wfDhw7F06VLKRYShYRQCffv2xZo1a7BhwwYcOXIEs2bNwrBhw3yOGTZsGF5++WUIIfDll1+iWbNmlguU1G+GDRuGl156CQDw0ksv4eKLL65zjBsZIomJEAI33XQTunfvjrvuusv6O+Uiudm5cyf27dsHAPjtt9/wwQcfoFu3bpSLJGbSpEnYsmULNm7ciFmzZmHgwIF49dVXKRNJzuHDh3Hw4EHr/++99x7y8/MpF5EmZmkf6gnz5s0TJ554osjNzRUTJkwQQggxbdo0MW3aNCHEsQxEo0ePFrm5uSI/P18sW7YslsUlEWLkyJHi+OOPF6mpqaJ9+/bi+eefF7t27RIDBw4UeXl5YuDAgWL37t1CCCEqKirEeeedZ51rkiGS+Hz66acCgCgoKBC9e/cWvXv3FvPmzaNcJDnffvut6NOnjygoKBA9e/YUDz/8sBBCUC6IEEKIRYsWWVnpKBPJzbp160SvXr1Er169RI8ePax3S7mILB4hDIGIhBBCCCGEEJJEMJSOEEIIIYQQkvTQMCKEEEIIIYQkPTSMCCGEEEIIIUkPDSNCCCGEEEJI0kPDiBBCCCGEEJL00DAihJAEYdu2bRg5ciQ6d+6MHj164Pzzz8fq1auxePFiXHjhhbEuXtTZtGkTCgsL0adPH/Ts2RPTp0+39g0YMABdu3ZF79690bdvX3zzzTcRK8fy5ctx5513AgAWL16ML774wto3ffp0vPzyy2G5z9atW633vHjxYng8HrzwwgvW/hUrVsDj8aC0tBTAsTpYvny56+vv3LkT5557bljKSgghiUhqrAtACCHEP0IIDB8+HKNGjcKsWbMAAN988w22b9/u+ho1NTXwer2RKmLUadu2Lb744gukp6fj0KFD1uKHcoX31157DUVFRZg5cybuvvtuvP/++xEpR1FREYqKigAcM1gyMjJw6qmnAgBuu+22sN3niSeewF/+8hfrd0FBAd58803cdNNNAIBZs2ahd+/eQV+/devWaNu2LT7//HOcdtppIZeXEEISDXqMCCEkAVi0aBEaNGjgo2j36dMHZ5xxBgDg0KFDuPzyy9GtWzdcc801kEvUderUCX//+99x+umn46233sIbb7yBgoIC5Ofn495777WulZGRgXvvvReFhYUYNGgQli5digEDBiA3Nxdvv/02AOD333/HDTfcgIKCApx00klYtGgRAGDVqlXo168f+vTpg169emHNmjUAjiny+fn5yM/Px5QpUwAAG/9fe/cWElXXBnD8P9a8NpiJoWFhNVlmoWNaNKbGeJgLKemgSJFiiXghHe2ENxUmZkpBkhSCkEmaiMRoFHSisqY0DzRNUVFWph2IiqywRjDnvZDZNOn0+X7vBx/m87sZZ62912nfzMOz9rKzk3nz5pGVlUVISAhpaWlcuXKF6OhoAgMDaWlpAQb/03tmZiaLFy8mPDychoaGIWvy119/4e7uDkBfXx8DAwPDrl1kZCSvX78etk6r1ZKbm4ter0ev19PR0QEMZqOMRiOhoaEYjUa6uroAqKurIyQkhAULFmAwGACUjF1nZydlZWUcOXKEsLAwbt68SV5enpLBsVgsLFmyhNDQUJKSkvj06RMwmNlxjGHu3LncvHlz2LGeOXPGKaMzY8YMbDYb7969w263c+HCBZYtW+Z0T1VVFVFRUYSEhChrm5eXR3p6OvHx8QQGBlJeXq5cv3r1aqqrq4ftXwgh/nQSGAkhxCjw4MEDFi1a5LL+7t27lJSU8PDhQ54/f86tW7eUugkTJmA2mzEYDOTm5nL16lUsFgutra3U19cDg4FIbGws7e3teHp6smfPHi5fvozJZGLfvn0AHDt2DID79+9TU1PDhg0bsNlslJWVsW3bNiwWC21tbfj7+9Pe3k5FRQV37tyhubmZ8vJy7t69C0BHRwfbtm3DarXy+PFjTp8+jdls5vDhwxQWFgJw4MAB4uPjaW1t5dq1a+zevZve3t4h8+7u7iY0NJTp06eTm5urZIt+duHCBVavXu1y7SZNmkRLSwubN28mJycHgM2bN7N+/XqsVitpaWnKVrn8/HwuXrzIvXv3lIDRQavVkp2dzfbt27FYLErQ6rB+/XqKi4uxWq3odDr279+v1PX399PS0kJJSYlTucOLFy/w9vZWAkGHlJQU6urquH37NgsXLhxS39vby+3btzl+/DiZmZlKudVq5fz58zQ1NZGfn8+bN2+AweyXq8BMCCH+dBIYCSHEH0Cv1+Pv74+bmxthYWF0dnYqdWvXrgWgtbWV2NhYfH19GT9+PGlpady4cQMYzL44shE6nY6YmBjUajU6nU5py2w2k56eDsC8efOYOXMmT548ITIyksLCQoqLi3n58iUajQaz2UxSUhIeHh5MnDiR5ORk5Qf3rFmz0Ol0uLm5ERwcjNFoRKVSOfV16dIlioqKCAsLIzY2FpvNpmRtfjZ9+nSsVisdHR1UVlY6bS1MS0vD39+f4uJitmzZ4nLt1q1bp3w2NTUB0NTURGpqKgDp6emYzWYAoqOjycjIoLy8nB8/foz4+Xz+/Jmenh5iYmIA2LBhg7L2AMnJyQAsWrTI6dk5vH37Fl9f3yHla9asUTKBjnkMNzeDwcCXL1/o6ekBYNWqVWg0Gnx8fIiLi1OySVOmTFGCJCGEGGskMBJCiFEgODiY9vZ2l/U/ZwrGjRtHf3+/8t3DwwNA2V43HLVajUqlAsDNzU1pz83NTWnL1f2pqamcPXsWjUZDQkICV69e/W1fP4/1d32dOXMGi8WCxWKhq6uL+fPnu2xz2rRpBAcHO2U7qqurefHiBampqWzatMnlvY55//r3cNeUlZVRUFBAd3c3YWFhfPz40WW7/4RjDX59dg4ajQabzTak3M/PD7VazeXLlzEajS7H/et3V+U2mw2NRvPfTUIIIUY5CYyEEGIUiI+Pp6+vz+l9kNbWVhobG0fcRkREBI2NjXz48IEfP35QU1OjZDBGwmAwKO+fPHnyhK6uLoKCgnj+/DkBAQFs3bqVlStXYrVaMRgM1NfX8+3bN3p7ezGZTEO2lv1OQkICpaWlSoDl2Ib3s1evXvH9+3cAPn36xK1btwgKCnK6Rq1WU1BQQHNzM48ePRq2r9raWuUzMjISgKioKOWQi+rqapYuXQrAs2fPiIiIID8/Hx8fH7q7u53a8vT05OvXr0P68PLywtvbWwncTp069Y/Wfu7cucNmkmBwe19xcfGwB2s45mY2m/Hy8sLLywuAhoYGbDYbHz9+5Pr16yxevBgYfK4hISEjHpcQQvxJ5FQ6IYQYBVQqFSaTiZycHIqKipgwYQJarZaSkhKXBwv8aurUqRw8eJC4uDjsdjvLly9n1apVIx7Dxo0byc7ORqfTMX78eE6ePIm7uzu1tbVUVVWhVqvx8/Nj3759TJ48mYyMDPR6PQBZWVmEh4e7/HH/q71795KTk0NoaCh2ux2tVsu5c+ecrnn06BE7d+5EpVJht9sEkwrYAAABXElEQVTZtWsXOp1uSFsajYadO3dy+PBhp+OtHfr6+oiIiGBgYICamhoAjh49SmZmJocOHcLX15eKigoAdu/ezdOnT7Hb7RiNRhYsWOAUnK5YsYKUlBQaGhooLS116qeyspLs7Gy+fftGQECA0uZIeHh4MHv2bDo6OpgzZ45TneMEvOF4e3sTFRXFly9fOHHihFKu1+tJTEykq6uLvXv3Ku9mXbt2jcTExBGPSwgh/iQq++/2OwghhBB/MK1WS1tbGz4+Pv/vofxHJpOJ9vZ2CgoK/lU7eXl5TJw4kV27dg2pMxgMNDQ04O3t/a/6EEKI0UgyRkIIIcQokJSU9D97p2k479+/Z8eOHRIUCSHGLMkYCSGEEEIIIcY8OXxBCCGEEEIIMeZJYCSEEEIIIYQY8yQwEkIIIYQQQox5EhgJIYQQQgghxjwJjIQQQgghhBBjngRGQgghhBBCiDHvb2LbV6iUw0irAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here, for example, we perform the scan over chromosome 3R, using \n",
    "# populations from Burkina Faso (BFgam), Uganda (UGgam) and Guinea-Bissau.\n",
    "# This is a good step to watch the Dask dashboard:\n",
    "pbs('3R', 'BFgam', 'UGgam', 'GW')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BFgam no. samples: 92\n",
      "UGgam no. samples: 112\n",
      "GW no. samples: 91\n",
      "finished allele counting\n",
      "found 1818277 segregating variants\n",
      "finished locating segregating variants\n",
      "finished computing PBS\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA08AAADgCAYAAADbnF28AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOyde3hUxf3/32d3s4GQAIGECCQQkhgQSUA3KKBiRAUvlULFiJeKVYs+YPoDE+ulWmjLE7Qmkhr9RrFW693YVqjSEC8Fr1RIAAkX5Y4kgRDuN81ld35/wBxnZ885e/aS7G7yeT3PPMnuuc2ZPWdm3vP5zGcUxhgDQRAEQRAEQRAEYYgl1BkgCIIgCIIgCIKIBEg8EQRBEARBEARBmIDEE0EQBEEQBEEQhAlIPBEEQRAEQRAEQZiAxBNBEARBEARBEIQJSDwRBEEQBEEQBEGYgMQTQXRh3njjDUycOLHdr7Ny5UokJycH/byPPfYYEhIScM455wTtnPPnz8ftt98etPNxioqKcM899wT9vB31G4YbL7zwAubMmeP38V9++SXOPfdcxMbGYsmSJUHMWecnNzcXf/3rX/069vvvv0dsbCycTmeQc0UEk8bGRpx33nlobm4OdVYIIuwg8UQQnZwvvvgC48aNQ69evdCnTx9ccsklWLNmDQDgtttuw4cffhjiHPrH3r17UVJSgs2bN2P//v2hzo4bWmLx0Ucf9bvDaUQgv+H8+fMRFRWF2NhYxMbG4rzzzsM///lPdfvKlSthsVjU7bGxsbjhhhuClXW/aWlpwYIFC/Dggw/6fY7f//73uP/++3Hy5ElMmTIliLkjRFJTU/Hxxx+rnwcNGoSTJ0/CarWGLE9a4k9+ZxljePbZZ5GdnY2YmBicc845yM3Nxdtvv93R2TXFiRMn8MADDyA1NRU9evTAoEGDMG3aNKxevRoAcO+992LWrFnq/q2trejRo4fmd//73/+QlJSEK664AosXL+7weyGIcIfEE0F0Yo4fP46f/exnyM/Px+HDh1FfX4958+YhOjo61FkLmD179qBv377o169fqLMS0dx88804efIkTp48idLSUtx+++1obGxUtw8YMEDdfvLkSbz//vshzO0Zli5dimHDhmHgwIGa29va2ryeY8+ePTj//PP9ur6Z8xORzW9+8xuUlpaipKQEhw4dQn19PRYsWIDly5eHOmseNDc3Y8KECaitrcUHH3yA48ePY8uWLZg+fTr+85//AADGjx+PTz/9VD2muroagwYNwmeffeb2HQA4HA4AZwZmXnjhhQ68E4KIDEg8EUQnZuvWrQCAW265BVarFd27d8fEiRORnZ0NAHjllVdw6aWXqvt/+OGHGDp0KHr16oVZs2bh8ssvV0do+b6FhYWIj4/HkCFDUFlZqR778ssv47zzzkNcXBzS0tIMG90nn3wSAwcORFxcHIYOHYpPPvlEc79jx47hjjvuQGJiIgYPHowFCxbA5XLh448/xtVXX42GhgbExsbizjvv9Dj24MGD+NnPfobevXujT58+uOyyy+ByuQAADQ0NuPHGG5GYmIghQ4bgmWee0c3r//73P4wbNw69e/fGyJEjsXLlSnXb4cOH8atf/QoDBgxAfHw8pkyZglOnTuHaa69V8xYbG4uGhgYPd8B///vfOP/889G7d2/k5uZiy5Yt6rbU1FQUFxcjOzsbvXr1ws0334wff/xRM3/yb6goCp5//nmce+65iI+Px+zZs8EY070/kUmTJiEuLg47duzwuu8PP/yAGTNmID4+Hueddx7+/Oc/u43cP/HEE0hPT0dcXByGDx+O9957zy3Pl1xyCebOnYvevXsjLS0NX331FV555RWkpKSgX79++Pvf/6577crKSlx++eXq5927d0NRFLz00ksYNGgQJkyYAAD429/+hvPOOw/x8fGYNGkS9uzZAwBIT0/Hzp07ccMNNyA2NhbNzc04duwY7r77bvTv3x8DBw7EY489prqWifnt06cP5s+fj+bmZhQWFmLQoEFISkrCfffdhx9++AHAT1aMkpIS9OvXD/3798fLL7/sVnYFBQUYPHgwevXqhUsvvVQ91uh5k0lNTcXChQsxfPhwxMfH41e/+pXbc/Liiy8iIyMDffr0weTJk9HQ0KBuUxQFzzzzDNLS0pCQkIAHH3xQfT/kZ5WXr5Zo3LFjByZMmIC+ffsiISEBt912G44ePQoA+OUvf4nvv/9eLec///nPHudqaGjA5MmT0adPH2RkZODFF19Uzz1//nzk5eXhjjvuQFxcHM4//3y1gw+Yr0d8ZevWrfi///s/vP3227j66qvRvXt3WK1WXHrppXjllVfU/Xbt2oXx48cjLi4OV111FWbPnu1WbjfddBPOOecc9OrVC+PHj8emTZvUbXfeeSdmzZqFa6+9FrGxsbjkkkuwf/9+zJkzB/Hx8Rg2bBjWrVtnKr+vvfYa6urqsGTJEowYMQJWqxU9evTAtGnTMH/+fADA5Zdfji1btuDgwYMAgM8//xzTp0/HqVOn3L4bO3YsoqKiAAAXX3wxdu7cqb43BEGcgcQTQXRiMjMzYbVaMWPGDFRWVuLIkSO6+x48eBDTpk3DwoULcejQIQwdOhRfffWV2z5ff/01hg4dioMHD+K3v/0t7r77brVj3q9fP3XU8+WXX8bcuXOxdu1aj+t89913ePbZZ7FmzRqcOHECVVVVSE1N1cxTfn4+jh07hp07d+LTTz/Fq6++ipdffhlXXXUVKisrVauI2KHhlJSUIDk5GU1NTWhsbERRUREURYHL5cINN9yAkSNHor6+Hp988glKS0tRVVXlcY76+npcf/31eOyxx3D48GEUFxfjxhtvRFNTE4AzncPTp09j06ZNOHDgAObOnYsePXq45e3kyZMYMGCA23m3bt2KW265BaWlpWhqasJ1112HG264AS0tLeo+FRUVWL58OXbt2oUNGzZo3qMeH3zwAdasWYNvvvkGFRUVmvcmwxjDsmXL0NLSguHDh3vd/w9/+AN2796NnTt34qOPPsLrr7/utj09PR2ff/45jh07hnnz5uH222/Hvn371O1ff/01srOzcejQIdx6662YPn061qxZg+3bt+P1119XXeq0qK2txdChQz2+//TTT7FlyxZUVVVhyZIlKCoqwr/+9S80NTXhsssuwy233ALgTId/0KBBeP/993Hy5ElER0djxowZsNls2L59O9atW4cPP/zQzbXr66+/RlpaGg4cOIDf/e53eOihh7B161asX78e27dvR319Pf74xz+q++/fvx/Hjh1DfX09XnrpJcyePVt9/woLC1FTU4OvvvoKhw8fxp///GdYLBavz5sWb7zxBqqqqrBjxw5s3boVCxYsAAD897//xSOPPIKKigrs27cPgwcPxvTp092Ofe+991BdXY21a9di6dKl+Nvf/qZ7HT0YY3jkkUfQ0NCALVu2YO/evWqH/bXXXnMr59/+9rcex99yyy1ITk5GQ0MD/vGPf+DRRx91E0H//ve/MX36dBw9ehSTJ0/G/fffD8C4Hvniiy/Qu3dvn++F89///hcpKSnIyckx3O/WW2/FRRddhEOHDmH+/Pl47bXX3LZfe+212LZtGw4cOIALL7wQt912m9v2iooKLFiwAAcPHkR0dDTGjh2LCy+8UK2LH3jgAVP5/fjjjzFp0iT06NFDd5/k5GQMHjwYn3/+OQDgs88+w2WXXYZx48a5fTd+/Hj1GJvNhoyMDHzzzTem8kEQXQZGEESnZvPmzWzGjBls4MCBzGq1shtuuIHt37+fMcbYyy+/zC655BLGGGN///vf2ZgxY9TjXC4XS05OZi+++KK6b3p6urr91KlTDADbt2+f5nV//vOfs9LSUsYYYytWrGADBw5kjDG2bds2lpiYyD766CPW0tKim++2tjZmt9vZpk2b1O+ef/55dvnll3ucU4vHH3+cTZ48mW3bts3t+//9738sJSXF7buioiJ25513MsYYmzdvHrvtttsYY4w98cQT7Pbbb3fbd+LEieyVV15hDQ0NTFEUdvjwYY9ra+VNPO8f//hHdtNNN6nbnE4nGzBgAFuxYgVjjLHBgwez1157Td3+4IMPsnvvvVfzPsXfkDHGALDPP/9c/XzTTTexhQsXah47b948FhUVxXr16sW6d+/OLBYLe/LJJ93uQ1EU1qtXLzW98847jDHGhgwZwpYvX67u++KLLxr+HiNHjmRLlixR85yRkaFu27BhAwOgPpeMMdanTx+2bt06zXNlZGSwyspK9fOuXbsYALZjxw71u2uuuYb99a9/VT87nU7WvXt3tnv3bsbYmTL+6KOPGGOM7d+/n9ntdnb69Gl1/zfffJPl5uaq+RWfGZfLxWJiYtj27dvV77766iuWmpqqllu3bt1Ya2uruj0xMZGtWrWKOZ1O1q1bN7Z+/XqP+zJ63rQYPHgwKy8vVz8vW7aMpaWlMcYYu+uuu9iDDz6objtx4gSz2Wxs165djLEzz4lYhs899xybMGECY8z9WWXsp/Ll93P55Zer9YLMe++9x0aNGuWWR17O8rm+//57ZrFY2PHjx9XtDz/8MJsxY4aajyuvvFLdtmnTJtatWzfGmPl6RAut/Ivv7J/+9Cd28cUXu20fOHAg69WrF4uOjma7d+9me/bsYVarlZ06dUrd57bbbnMrN5EjR44wAOzo0aOMMcZmzJjB7rnnHnX7M888w4YNG6Z+3rBhA+vVq5ep+7nyyivZQw89pH5et24d69WrF4uLi2OZmZnq9zNmzGBz5sxhTqeTJSYmslOnTrHy8nL1u969e7OVK1e6nXvcuHHs73//u6l8EERXgSxPBNHJOe+88/DKK6+grq4OGzduRENDg2aUsoaGBqSkpKifFUXxCHogRrWLiYkBANU6UFlZiTFjxqBPnz7o3bs3/vOf/6juICIZGRkoLS3F/Pnz0a9fP0yfPt3NnYhz8OBBtLS0YPDgwep3gwcPRn19van7fvDBB5GRkYGJEyciLS0NTzzxBIAzc10aGhrQu3dvNRUVFbnN8+Hs2bMH7777rtu+X3zxBfbt24e9e/eiT58+iI+PN5UfkYaGBrf7slgsSElJcbs3uaz1rDBa+HJsXl4ejh49itOnT2PHjh149dVX3VwuBwwYgKNHj6opLy9PvQfxeRH/B4BXX30Vo0aNUstt48aNbs9DUlKS+n/37t01v9PLd3x8PE6cOOHxvZiHPXv24P/9v/+nXr9Pnz5gjGk+P3v27EFrayv69++v7n/vvffiwIEDmuduamrC6dOn4XA41P2vueYaNwtR3759YbPZ1M/8dzh48CB+/PFHpKena+ZD73nTQ8zX4MGD1XdJfsZiY2PRt29ft/vXO9YXDhw4gOnTp2PgwIHo2bMnbr/9ds33XouGhgb06dMHcXFxbvkweg9+/PFHtLW1ma5HtLDZbGhtbXX7rrW1VXVX69u3r0eZ19XV4eDBg2hubgZjTM07rwcB9/J0Op14+OGHkZ6ejp49e6pWMaN3wOzzLyPnd9SoUTh69Cj+9a9/uUXLGz9+PD777DPU1tYiLS0NMTExuPTSS9XvfvjhB1x88cVu5z5x4kRAVjyC6IyQeCKILsSwYcNw5513YuPGjR7b+vfvj7q6OvUzY8ztsxHNzc248cYbUVhYiMbGRhw9ehTXXXed7lybW2+9FV988QX27NkDRVHw0EMPeeyTkJCAqKgoN3/777//XjdIgExcXBxKSkqwc+dOvP/++3j66afxySefICUlBUOGDHETBCdOnFAnVoukpKTgl7/8pdu+p06dwsMPP4yUlBQcPnxYnd8hoiiKYd4GDBjgdl+MMezdu9f0vbUXqampuPbaa00FhZCfl71796r/79mzB7/+9a/x7LPP4tChQzh69ChGjBhheu6VN7Kzs9X5fCJiuaekpOCFF15w++1++OEHjBs3zuO4lJQUREdH4+DBg+q+x48fd5ujIp47ISEB3bt3x6ZNm9T9jx07Zqqzm5CQgG7dumnOKzN63vQQy/37779XXUTlZ+zUqVM4dOiQ2zOmd2yPHj1w+vRpdZtRNMtHHnkEiqJgw4YNOH78OF5//XW339noXRgwYAAOHz7sJoR9ecfN1CNaDBo0CLt373b7bteuXarYnDBhAurq6tzmV8n0798fhw8fdisnsTzffPNNLF26FB9//DGOHTumXi9Y74DIlVdeiQ8//BCnTp0y3G/8+PH45ptvsGzZMlx22WUAgPPPPx979+7FsmXLMHr0aHTr1k3dv62tDdu3b8fIkSODnmeCiGRIPBFEJ+bbb79FSUmJ2sndu3cv3nrrLYwZM8Zj3+uvvx61tbVYsmQJ2tra8Nxzz5kOAd7S0oLm5mYkJibCZrOhsrJSN3z2d999h//+979obm5Gt27d1MnYMlarFXl5efjd736HEydOYM+ePXj66adNr8H0wQcfYPv27WCMoWfPnrBarbBarbjooovQs2dPPPnkk/jhhx/gdDqxceNGNXy7yO233473338fVVVVcDqd+PHHH7Fy5UrU1dWhf//+uPbaazFr1iwcOXIEra2tauSqpKQkHDp0CMeOHdPMW15eHpYtW4ZPPvkEra2tKCkpQXR0tGbHviOpq6vD8uXLTUWhy8vLw8KFC3HkyBHU19fj2WefVbedOnUKiqIgMTERwJlgIlqC3V+uu+46t8hhWtx3331YuHChKoCOHTuGd999V3Pf/v37Y+LEiSgoKMDx48fhcrmwY8cO3WtYLBb8+te/xty5c1XrVH19vam5ZRaLBXfddRceeOABNDQ0wOl0YtWqVWhubjZ83vR47rnnUFdXh8OHD6OoqAg333wzgDPC4uWXX8b69evR3NyMRx99FBdffLHb/MKnnnoKR44cwd69e/GXv/xFPXbUqFH47LPP8P333+PYsWNYuHCh7vVPnDiB2NhY9O7dG/X19XjqqafcticlJWHnzp2ax6akpGDcuHF45JFH8OOPP2LDhg146aWXPOYGaWG2HtHi5ptvxssvv4zVq1eDMYatW7di0aJF6pywoUOH4t5778X06dPx0UcfqfWEOAd08ODByMnJwfz589HS0oJVq1a5DTqcOHEC0dHR6Nu3L06fPo1HH33UVN70uPPOOzUD4wDAHXfcgf79+2Pq1KnYuHGj+uzI4i8jIwNJSUn4y1/+ooonRVFw8cUX4y9/+YvbfCcAWL16NVJTU90smARBkHgiiE5NXFwcvv76a1x88cXo0aMHxowZgxEjRqCkpMRj34SEBLz77rv47W9/i759+2Lz5s3IyckxFdY8Li4OzzzzDPLy8hAfH48333wTkydP1ty3ubkZDz/8sLq47YEDB1BUVKS5b1lZGXr06IG0tDRceumluPXWW3HXXXeZuvdt27bhqquuQmxsLMaOHYtZs2YhNzcXVqsV77//PtavX48hQ4YgISEB99xzj6bQSUlJwdKlS1FUVITExESkpKTgqaeeUqOSvfbaa4iKisKwYcPQr18/lJaWAjhj4bvllluQlpaG3r17e7gTDR06FK+//jry8/ORkJCA999/H++//z7sdrupewsm77zzjhoVcPTo0bjkkkswb948r8f9/ve/R3JyMoYMGYKrrroK06ZNU5+V4cOHo6CgAGPHjkVSUhJqa2txySWXBC3PN9xwA7799ltDN62pU6fioYcewvTp09GzZ0+MGDHCLTqkzKuvvqoGy4iPj8e0adMM3eWefPJJZGRkYMyYMejZsyeuuuoqfPfdd6byX1xcjKysLIwePRp9+vTBQw89BJfL5fV50+LWW29VXVPT0tLw2GOPAThjjfjTn/6EG2+8Ef3798eOHTs81ij6+c9/DofDgVGjRuH666/H3XffDQC4+uqrcfPNNyM7OxsOhwM/+9nPdK8/b948rF27Fr169cL111+PX/ziF27bH3nkESxYsAC9e/dGcXGxx/FvvfUWdu/ejQEDBmDq1Kn4wx/+gKuvvtprGRrVI59//jliY2N1j500aRKeeOIJ/OpXv0KvXr1w3XXXYcaMGZg5c6a6z3PPPYff/OY3eOCBB9CnTx8kJyfj8ccfxzvvvINBgwYBOBOsY9WqVejbty8ee+wx3Hzzzeo7cMcdd2Dw4MEYOHAghg8frjlg5Qt79+7VfYe6deuGFStWYPjw4bj++uvRs2dPDB06FGvWrEFFRYXbvuPHj0dTU5PbuS677DIcOHDAQzy98cYbuO+++wLKN0F0RhTWHjZkgiAiHpfLheTkZLzxxhu44oorQp0dIswpLy/H22+/7dUiFCwWL16MzZs3q4K1K5Kamoq//vWvuOqqq3w+VlEUbNu2DRkZGe2Qs67JzTffjGHDhuEPf/hDUM/b0tKCkSNHYsOGDeq8rPbmwIEDuPzyy7Fu3To3Vz6CIMjyRBCEQFVVFY4ePYrm5mYUFRWBMRbwiCnROdm3bx++/PJLuFwufPfddygpKcHUqVM77PozZ87s0sKJCD1r1qzBjh074HK5sHz5cixduhRTpkwJ+nXsdju2bNnSYcIJOLP0xJYtW0g4EYQGNu+7EATRVVi1ahVuvfVW1X1pyZIlaiQ0ghBpaWnBvffei127dqF3796YPn06Zs2aFepsEUSHsX//fvziF7/AoUOHkJycjPLyclxwwQWhzhZBEO0Mue0RBEEQBEEQBEGYgNz2CIIgCIIgCIIgTEDiiSAIgiAIgiAIwgRdas5TQkKC2xoXBEEQBEEQBEEQIrt378bBgwc1t3Up8ZSammq4YjhBEARBEARBEF2bnJwc3W3ktkcQBEEQBEEQBGECEk8EQRAEQRAEQRAmIPFEEARBEARBEARhAhJPBEEQBEEQBEEQJiDxRBBEp6ewsBDR0dEoLCwMdVYIgiAIgohgSDwRBNHpKSsrQ1VVFcrKykKdFYIgCIIgIhgSTwRBdHry8/MxadIk5OfnhzorBEEQBEFEMApjjIU6Ex1FTk4OrfNEEARBEARBEIQuRpqBLE8EQRAEQRAEQRAmIPFEEARBEARBEARhgpCKp+XLl2Po0KHIyMjAE0884bF96dKlyM7OxqhRo5CTk4MvvvjC9LEEQRAEQRAEQRDBJGRznpxOJzIzM/HRRx8hOTkZo0ePxltvvYXhw4er+5w8eRI9evSAoijYsGED8vLy8O2335o6Vgua80QQBEEQBEEQhBFhOedp9erVyMjIQFpaGux2O6ZPn46lS5e67RMbGwtFUQAAp06dUv83cyxBEARBEARBEEQwCZl4qq+vR0pKivo5OTkZ9fX1Hvu99957GDZsGK6//nr87W9/8+lYgiAIgiAIgiCIYBEy8aTlLcgtSyJTp07Ft99+iyVLluDxxx/36VgAWLx4MXJycpCTk4OmpqYAc00QBEEQBEEQRFclZOIpOTkZe/fuVT/X1dVhwIABuvuPHz8eO3bswMGDB306dubMmaiurkZ1dTUSExODdwMEQRAEQRAEQXQpQiaeRo8ejW3btmHXrl1oaWnB22+/jcmTJ7vts337dtXKtHbtWrS0tKBv376mjiUIgiAIgiAIgggmtpBd2GbDs88+i0mTJsHpdOKuu+7C+eefj+effx4AcN999+Gf//wnXn31VURFRaF79+545513oCiK7rEEQRAEQRAEQRDtRchClYcCClVOEARBEARBEIQRYRmqnCAIgiAIgiAIIpIg8UQQBEEQBEEQBGECEk8EQRAEQRAEQRAmIPFEEARBEARBEARhAhJPBEEQBEEQBEEQJiDxRBAEEcYUFhYiOjoahYWFoc4KQRAEQXR5SDwRBEGEMWVlZaiqqkJZWVmos0IQBEEQXR4STwRBEGFMfn4+Jk2ahPz8fI9tZJUiCIIgiI6FFsklCIKIUKKjo1FVVYVJkyahubk51NkhCIIgiE4BLZJLEATRCTGySnU0ZAUjCIIgfCFS2w2yPBEEQRABQ1YwgiAIwheio6MxZcoUVFRUoKCgAMXFxaHOkgpZngiCIIh2JZysYARBEET4k5+fj4qKCuTl5UVUUKSQiqfly5dj6NChyMjIwBNPPOGx/Y033kB2djays7Mxbtw4fPPNN+q21NRUZGVlYdSoUcjJyenIbBMEQRASxcXFaG5uDquRQ4IgCCJ8KS4uRkFBAZYsWRJRA2+2UF3Y6XRi9uzZ+Oijj5CcnIzRo0dj8uTJGD58uLrPkCFD8OmnnyI+Ph6VlZWYOXMmvv76a3X7ihUrkJCQEIrsEwRBEARBEAQRAHzAjVueImEALmSWp9WrVyMjIwNpaWmw2+2YPn06li5d6rbPuHHjEB8fDwAYM2YM6urqQpFVgiAIgiAIgiDagUhbzzBk4qm+vh4pKSnq5+TkZNTX1+vu/9JLL+Haa69VPyuKgokTJ8LhcGDx4sW6xy1evBg5OTnIyclBU1NTcDJPEARBEARBEETARNqc2ZC57WkF+VMURXPfFStW4KWXXsIXX3yhfvfll19iwIABOHDgAK6++moMGzYM48eP9zh25syZmDlzJgDQ3CiCIAiCIAiCCBMKCwtRVlaG/Pz8iHDZA0JoeUpOTsbevXvVz3V1dRgwYIDHfhs2bMA999yDpUuXom/fvur3fN9+/fph6tSpWL16dftnmiCIiCdS15UgCIIgiM5GpLnsASEUT6NHj8a2bduwa9cutLS04O2338bkyZPd9vn+++/xi1/8Aq+99hoyMzPV70+dOoUTJ06o/3/44YcYMWJEh+afIIjIJBIraoIgCILojESayx4QQrc9m82GZ599FpMmTYLT6cRdd92F888/H88//zwA4L777sMf//hHHDp0CLNmzVKPqa6uRmNjI6ZOnQoAaGtrw6233oprrrkmVLdCEESYI7oFRGJFTRAEQRCdkeLi4ohx1+MoTGvyUSfFaLVggiA6L9HR0aiqqsKkSZPQ3Nzsdf9I9MEmCIIgCCI4GGmGkC6SSxAE0RH4am0i1z6CIAiCILQImdseQRBER+GrWwC59hEEQRAEoQVZngiC6PLIEfiKi4vR3NwcVi57FCWQIAiCIEIPiSeCILo8Rm564SJawtmVMFzKiCAIgogsIrH9IPFEEESXx8hNL1xES35+Pq688ko4nU6fGpmOaJjKysowZcoUlJSURFQDSKaBQ28AACAASURBVBAEQYSWcGljfYHEE0EQXR4jN71wmf9UXFwMm82Gjz/+2KdGJtCGyYz4ys/PR0VFBfLy8iKqASQIgiBCS7i0sb5AocoJgiDCGDFsOgCfQ6gHGnbdbJh3Cu9OEARBdBaMNAOJJ4IgCIFwEwG+rlEVbMKtPAiCIIjIJ9zbFhJPZyHxRBCEN0ItVmQCtTwRBEEQRLgRbm2tDC2SSxAEcRZvc3i4/3VWVlZYRAAS52NF4sRagiAIgpCJxLlOHBJPBEF0ekTB5E2AcLFSW1sbMqGiJ/D8aWwCjbYXiWFkCYIgiPAl3F32vMJCSGVlJcvMzGTp6els4cKFHttff/11lpWVxbKystjYsWPZ+vXrTR+rhcPhCFreCYKIHOx2O1uxYgWz2+2soKBA/Ssif6+3X0fnN9TnslgsLCoqilksloDzQhAEQRDBbOPaCyPNEDLLk9PpxOzZs1FZWYnNmzfjrbfewubNm932GTJkCD799FNs2LABjz/+OGbOnGn6WIIgCE5WVhauuOIKZGVleWzjlpXS0lI3S5NR+PL2JpjuDIGeS1EUMMagKErAeSEIgiAIf9ctDBdCFjBi1apVmD9/PqqqqgAACxcuBAA88sgjmvsfOXIEI0aMQH19vc/HcihgBEF0TcSJqQDcJqnybVdeeSVsNhvi4+PR2NgIh8NB9QU6gXsFQRAEEXZQwAg/qK+vR0pKivo5OTkZ9fX1uvu/9NJLuPbaa/06liCIro1oeZItMfzz3LlzkZ+fj8bGRuTm5qKmpkZ3rk9XmkcUSgscQRAE0TmhgBF+oGXw0nMLWbFiBV566SU8+eSTPh+7ePFi5OTkICcnB01NTQHkmCCISKW2thYrVqxAbW2thxiQo9llZmZi5cqVUBRFdeOTxU5paSlcLhdKS0v9yk9HRs0LplCLJNFHEARBhC+RPDAXMvGUnJyMvXv3qp/r6uowYMAAj/02bNiAe+65B0uXLkXfvn19OhYAZs6cierqalRXVyMxMTHId0EQRCSgNcKlJQTy8/Oxe/duFBQU4IEHHlCPkcWOy+VCW1sbXC5X0PLjK2aFjL9CTev8FCqdIAiC6PJ0TMwKT1pbW9mQIUPYzp07WXNzM8vOzmYbN25022fPnj0sPT2dffnllz4fqwVF2yMIgsOj/VgsFq/R9+TIe94i0HVEpD5v0YoKCgqYxWJhiqIwq9Xqc160zh/KCIQEQRAE0VEYaYaQhipftmwZO/fcc1laWhpbsGABY4yx8vJyVl5ezhhj7O6772a9e/dmI0eOZCNHjnS7Ea1jvUHiiSC6JlxIiCKCCwGr1aopQozEiTcREYwwrN6u4W27oigMAAPgVz58EUokqgiCIAhfCPd2I2zFU0dD4okguiZ2u51FRUUxm82mWpocDof6VxZWjP1UsfP9tLb5K2zM5tlXASbmmQsnAExRFE0LVLAar0hYs4MgCIIIH8K93SDxdBYSTwTRNREtT4qisKioKAZAdb0zclGzWq0sLy+PAVBFhq+Vvq8iRctSZgaeLwAsMzNTFU78fm02m1ueg9V4hfsIIkEQBBFehHu7EZaL5BIEQXQkNpsNCQkJYIyhtbUVwJkFt10uF1paWjBhwgQ1gENhYSFKSkrQ1tYGp9OJiooK5OXlqYESjAI+eAu0wLfn5OToBnwoKyvDJ598AqvV6lMkIp4vh8PhFvjC6XSqi92KefY3cIV4j4WFhVi0aBGcTqdP5yAIgiC6LnK0vYiK5tphEi4MIMsTQXRNrFYrs9lsDADLzc11s8zgrHUGgmXJarWq27j1x+wImbdAC6J1yGhOFbc8abkNBoI/VjB5f/EeRZfIcHW/IAiCIMKbcHPjI7e9s5B4IoiuCY+OxwVRUlKSKgggzA3S+s5X0WJ2PpQ3USSKLJvNxhRF8Tr3SvzM5z3J9Z6vDZQ3MeiviyFBEARBcMLNjY/E01lIPBFE10SvUtYSShaLhVksFkPh1BGVPL+GmL+8vDzDOUviZwDqfYjn89WSpSXCwq2RIwiCIIhgQuLpLCSeCKJrohWFjosI7r7H6wfu4me1WnXP197uBaI4EV0IRUFnFH5dFoXe8mwkhux2e8ABMwiCIAgikiDxdBYSTwTRNREFgGiR4YkLKKvVypKSktR5UKGyPIniRM8tzpdFcsW5XXoCyUhYyVYvsjwRBEEQnRkST2ch8UQQXRPZEsMTDx4hJ1Fcyesj+ev+5mt+5blMRlYmPXgwB/HetOYmeZu3pDd/iiAIgiB8JRIG4Eg8nYUafoLouohzmaxWq2ptiomJcRNNoqVGFlSMuQdy8Oa6FmhgBXGtKV8j2mkJRkVRdM8hW5+0IgSSmx5BEAQRKJHQptA6TwRBdHnmzp0Ll8uFvLw8MMawfft2WK1WnD59GosWLVL3s1gsSEpK8jheURQAZ9ZGuvLKK6EoCq666irDNZIWLVoExhgWLVqEsrIyWK1WKIqirhflDb4+FGMMTqfTY50mb8dGRUUBOLPGlcVi0VzriSOv+SSuTSVui6i1OAiCIIiww981BsOGjtFv4QFZngiiayNaU2QLE3fD41Yl7u4mWqv4PlarVddKI8KvwecambU8ia6B/oYBFy1PiqKYdjPkLnpiOHcROaBGJLhfEARBEIQvkOWJIAgCP61oDgCMMdWaZLfbUVdXh5aWFkyYMAGKoqCtrQ3AGYuT1WrFxx9/jJqaGlRVVcHpdOKKK65AVlYWgDNWmtTUVJSUlOCcc85RLTMWi0W1+hQXF8PpdKrnla03okWnrKwMU6ZMQU1NDVwuF9LT01FSUgKbzeazxScvLw9RUVGora1VLUlG1NTUYNGiRWhsbHRb/Z3Dy40xpt67mfMSBEEQRKfASHWdOnWKtbS0qJ+//fZb9vTTT7N//vOfQVF1lZWVLDMzk6Wnp7OFCxd6bN+yZQsbM2YMs9vt7KmnnnLbNnjwYDZixAg2cuRI0xYlsjwRRNdEDvIgznmCYGVatGiRGpIbZ61OFotF83jZ8gQhAAU/t2y9EecwGYX/5hH/cNZqxPNis9mYxWIxZemRIwzyuV3e6kE5OITRQrxan4MBWbMIgiCIUOJ3wIjLLruMbd26lTHG2LZt21h8fDy7//772YQJE9jDDz8cUKba2tpYWloa27FjB2tubmbZ2dls06ZNbvs0Njay1atXs0cffVRTPDU1Nfl0TRJPBNE1kYM8iAJEy32PMXf3NNHlTs+VjosOLlKSk5OZOD4lhvzmroB64b9FISYmq9XqIdy0kMOU83PJx3mLosfzwsWgHP2PRyIMduTBSJhMTBAEQXRe/BZPI0aMUP9/7LHH2KxZsxhjjDU3N7tt84evvvqKTZw4Uf1cVFTEioqKNPedN28eiSeCIPyGCxNu0eEWITHSnqIo6r5iqHKbzeZm+dESIXqL2jocDnWboiiqmOHR9/QEhyhqRPEjrkOlVZ+JeedWJzmqoGhREs+thbiIcGZmpio6efQ/sXy8LSzsC2R5IgiCIEKJ3+IpKytL/X/cuHHsvffeUz9nZ2cHlKl3332X3X333ernV199lc2ePVtzXy3xlJqayi644AJ24YUXshdeeMHUNUk8EUTXRMtdTivJwomLHX6MnpVFdJFTFIVFRUWpoc0tFovbWktcpJmFCygxiRYrEb6uE883zys/h3ic3W5n8fHxbvfu7brc6iWXkXzPBEEQBBHJ+B0wIjs7G4WFhVi0aBG2b9+OiRMnAgCOHj1qdJgp2NnJxiJ88rYZvvzyS6xduxaVlZV47rnn8Nlnn2nut3jxYuTk5CAnJwdNTU1+55cgiMiFB2BwOp2oqKjQ3W/RokVwuVzqZ15PVVRUICkpCVarFbm5uR6BFPLz81FRUYG8vDwoigJFUXDBBRcgOjoa7MwgFYAzYdAZYz4FfcjNzXX7nJmZiYqKCjidTo/z5Ofnw+l0quHWa2pq0NraivXr16vH8SAX+fn5OHLkiBrUQg74UFNT41YnW61WjBo1CldccQUuuOAC2O12Nez7hRdeCEVRMHfuXNP3RRAEQRCRiKF4evHFF5GQkIDdu3fjww8/RExMDABg8+bNAa/xkZycjL1796qf6+rqMGDAANPH83379euHqVOnYvXq1Zr7zZw5E9XV1aiurkZiYmJAeSYIIjLJyspCRUWFKmz0UBQFeXl5Ht9bLBY0NjbqRpUrLi5GQUEBlixZgrlz56K5uVmNbqcoClwul/o3Ly/PMDKdvI5SWVmZmieHw4Hdu3er0f/k8/CIfnPmzEFjYyMAqOJt69atAIC1a9eqEf0cDoe6j7zehsPhcBvkmjNnDmpra7FixQrU1ta6CcZ169bp3o8/0FpSBEEQRNjiqxmrqamJuVwuPwxg7rS2trIhQ4awnTt3qgEjNm7cqLmv7LZ38uRJdvz4cfX/sWPHssrKSq/XJLc9guiayEEheOLuZzabjSmKogaCEOdCyUmcxyS6uulFoRNd5vicIaMAC2KwBHGtJj63CJI7nhy8QjxGTrwcRPc7o7lFPO98nSpxLpYYRMNMEAtfoIARBEEQRCjxe87TqlWr2OWXX86mTp3K1q5dy84//3yWlJTEEhMTTYkVbyxbtoyde+65LC0tjS1YsIAxxlh5eTkrLy9njDG2b98+NnDgQBYXF8d69erFBg4cyI4dO8Z27NjBsrOzWXZ2Nhs+fLh6rDdIPBFE14QLBj1R4Uuy2+2anXtx3pPefCj+PQ+BLu/LmKfo0ko8MIPdbnebf8SPNZrXJZ5DK+S6KMb4uRVFUeeM8cAQ8jwvfxby1YMCRhAEQRChxG/x5HA4WFVVFauoqGC9e/dmq1atYoydWX9p1KhRwc1lB0DiiSC6JkZCxJfELTB6licI0e24dUgrOp5oCZOtK6JoEdeiMpv01rLSux8xLzzwA19PStxXtmiJa2GJUQQJgiAIItLxWzyNHDlS/X/YsGFu20g8EQQRCXChEwzxJFqFxHWOuGjgIo1HsePiRRQispCLiYlxy6t8TS6++Dm5UNFzRRSPMUqyODK7j7jelfh9ZmYmudkRBEEQnQK/o+1ZLD9t7t69u9s2XyLjEQRBhIqSkhK0tLR4fB8VFeXzucrKylBVVYWamhpYrVYwxqAoihq4oba2FvHx8Thy5Aiio6OxdetWNVAEcCZyHj+Wc/r0afX/0tJSj2s2Njaq50xKSlKDOPC/DofDra7mx+ihKIoaKAI4E7zHF1wuFxISEtyi8eXl5WHr1q0eQScIgiCIyIAC9ZjHUDx988036NmzJ+Li4rBhwwb07NlT/VxbW9tReSQIggg6ra2tPh+Tn5+PSZMmISkpST2+ra1NDRuelZWFI0eOAACam5uhKIoqcqKjo7Fy5UokJSXB6XSq50xKSlL/F78XOXLkCKxWq6YoqqurMx0iPDMzE4wx5Obmwul0wuFwoK6uTnd/MWy7iBjJD4AayZAgCIKITPjgoFE0WOIMhuLJ6XTi+PHjOHHiBNra2nD8+HH1sz8dD4IgiI5GtLJo4Uunv7i4GPn5+R4ixuVyoaSkBGvXrnX7nouLvLw8NDc3A3C3CimKgv3796ufRYuUjJ6wamxsRGlpqZv1SbQIifBw5dzCVV1djYKCAt1rekMsO9ECRxAEQUQWfHAwEA+CrmK9MhRPP/74I0pLS3H//fdj8eLFaGtr66h8EQRBBIXc3FwPtzYRprFgtxZc2JSUlOieQ+9cfGFeecFbdnbBXN7gJCQkADhjIbLb7brCRhZ8jDE3K1FUVBQcDofugsBiPsV1pHyhoKDAzfXR5XKpC/ASBEEQkUVxcbHHAvC+0lWsV4biacaMGaiurkZWVhb+85//BDRCSRAEEQrKysowbdo0zW3R0dEAzogVjsPhQEFBgYfgcjqdsNlsHucQ3e7i4+MN88Ld9kQWLVqEsrIyTJkyBY2NjepcKZvN5iHU+ELlskibO3curFarKqqcTifWr1+vK4oSExMRHR2NnJwctLW14d133zXMt0xMTAwWLVrkMZcs2IvlEgRBEMGnvSxEwbBeRQKG4mnz5s14/fXXce+99+If//gHPv/8847KF0EQRFDIysrStcBwVzruzgYANTU1KCkp0ZzvI7rOKYoCm82GpqYmAGcEGA/qoEdBQQH279/vto/L5UJ+fr6aRz5XSgwkwa93+vRpTTfDlStXYs6cOaqoSk9Ph9Pp1L3vxsZGpKamoqamBtOmTQNjzCfr0+nTpzXLJ1jznrqK6wdBEEQo6CoWovZCYQY+KxdeeKGbD7/8OdLIyclBdXV1qLNBEEQHYrPZdOcL+UNmZqYaRY9bp3r27KkKJ71Id6KwkKtdxhisViusVqvP80lzc3OxcuVKOBwO1NTUuF3PrEtiUlISmpqaVEEUFRVlmA+bzabpxl1QUBCQywcnOjoaVVVVmDRpkipwCYIgiOBQWFiIsrIy5OfnB6XO5nSmuttIM5iKtidH3IuLi0PPnj3bJbMEQRDBxEhAcLc9s8TExKhWKsYY5syZgzlz5qgR9rSEE58rxcWMVn6io6NxwQUXwOl0QlEUNXFBpZUPzsqVKxEfH+8mnHj+zMLzz2ltbTW0ImkJp2BG2+sqrh8EQRAdhWjRD8b8Ji26St1taHnqbJDliSC6HoWFhZpBHjjckuQNPWuL3W7XXEcKAGJjY3Hy5El1nSY9LBaLKpamTJmCiooK5OXl4R//+IeHexw/J3BmjSa9UON8fSk9USNapqxWq8/WOdmylZmZid27d0f8aCNBEERnpDNZhToCvy1PBEEQnR0zwgk4Y20RLT4cPeEEQBU5RsIJODPvyel0oq2tDRUVFYiPj0dFRYXmvCJ+TgCGazQxxtCjRw/1s2hl4+s9AWdE0I033miYP73zc2JjY7F161avATMIgiCI0NBVrEIdAYkngiA6NWVlZW7R9AJBDuIQbLgg8Sa2zCLmt7m5GQUFBSgoKMD27dsB/CSi9AJLiBitl8UFnd58L1+hgBEEQRDBxZurXjDq3ZycHCiKgpycHL/PEQmEVDwtX74cQ4cORUZGBp544gmP7d9++y3Gjh2L6Ohojx/b27EEQRDAmdE2s9Ylb8TGxgZ0vNVqhd1u151rxcWT0bpUgfD000+jtLQUjDFYLJaglUuwoUhQBEEQHUsw6t2amhosWrTIYw5uZyNk4snpdGL27NmorKzE5s2b8dZbb2Hz5s1u+/Tp0wfPPPOMhwo2cyxBEAQAdeAlkIAGXMycPHnSQ9j4EuLb5XIhPj4ezc3NhgLJYrG0i4BijMHpdHosqusNrYAUgKeY1HJr9AdyLyEIguhYglHvOhwOzJ07F4qidGrPgZCJp9WrVyMjIwNpaWmw2+2YPn06li5d6rZPv379MHr0aLdV7M0eSxAEwUlKSjKMPme32w2P50LDZrO5BWGw2Wy6Lm9aQoIxhsbGRmRmZsLlciE5ORmAp6Wpra1NvaY30RcswWKEnhvhyZMnVStabm5u0Nway8vL0dLSgvLy8qCcjyAIgjAmGBH4qqurYbfbcdNNN6GkpERTQHUGt+yQiaf6+nqkpKSon5OTk1FfXx/0YxcvXoycnBzk5OSoi1kSBNG1OHjwoO42RVHQ3Nysu7itzWZT/+fR9rgQ04q+xzl9+rSumx+fc8QDPhhZgbwFRD19+jRGjRpluE97wqM2rVy5MmhCjouw9p5jRhAEQQQXvjB9ZmampgtgZ3DLDpl40uoQmHWr8eXYmTNnorq6GtXV1UhMTPQtkwRBdArMhOHWC3bgz2oOvD4SI+PJ21esWAEgOPOb1q9f7/ex3PoVDEjsEAThjc5geSD0WbdunTqnVssFsDO4ZYdMPCUnJ2Pv3r3q57q6OgwYMKDdjyUIomsRSAPtcDgMrUt6MMZgtVrdrFYiTqcTV1xxBQBjq1NHYBTuPBTIv1dnj9pEEOGKvyLH23GdwfJA6KMoCiwWC6xWq6YLYHst0NuRhEw8jR49Gtu2bcOuXbvQ0tKCt99+G5MnT273YwmC6FqUlZUZBnUwWkS2pqZGFUC+Bpzg6zZ1JYzCmZulrKxMnedqsVg6fdQmgghX/BU53o4LZ8uDkfAL1GLWVSxuc+bMgcViwZw5c0KdlfaDhZBly5axc889l6WlpbEFCxYwxhgrLy9n5eXljDHG9u3bxwYOHMji4uJYr1692MCBA9mxY8d0j/WGw+FonxshCCJoFBQUMLvdzgoKCoJ2PovFwgB4JEVRmKIozGq1qp+19pOPsdvtXvcLZrLZbIbbef593WY2mSmXYDUn8u+lKErQngWCIMzjb10c7Dq8I7Hb7WzFihXMbrf7tC3QcxPhh5FmCKl46mhIPBFE+NMeDYyegOCiQBZXslgwKx5CkcIpb8HuLFFngyAIxs703wC0ez/OSPgFKgrDTVRq5acjyjncykEPEk9nIfFEEOFPe1SsojiKj49nAFheXp5bx9/hcLCCggIPq01MTIzPIiI6OpoBYJmZmbr7xMTEMIvF0mHixygvviRvVrBgwS1QVqs17BtZgmhPIqWzyWmP/AJgixYtCmod09WQhZHW4FR7lLP8PETKoBiJp7OQeCKIrklSUpKbIHI4HMxut6vfZ2ZmMrvdzux2uyquuACyWq26bn9GiYsjrW38ujwfspALdvImeIKZggEXsXl5eWHfwBJEexMpnU1Oe+S3oyxPkYRZkcr3k4VRR1me5OchUgYDSDydhV46guiaAO7ubeIIGLfIiAJLtNZwAeSPgDJKvsxF4kLOn+SLMAuGFSxQuHDiIpbqbaKrEymdTU6wO+AdaYWOpLI2K1L5fryONvpdvN2/P9sjqUxFSDydhRphguia8MZc7uTzjrqRMPLXda+rpkCx2+2qWIyOjjY1eh2pjTPRftAzEToCtTxpuXlFRUWpgzvt+ZtGkpXPV8uTmXLzdv9mtvMBu0h3uSbxdBYSTwQR/rRHp0eeyySLJz3RZLVa3UbseAp0/lCwrVhGKVhzndpLPIm/txhpLzEx0XQDHEkdHqL9EJ+lzvRMRJoQDDS/Wm5evF5ob1feSCtrf+FlyqPN8vsN1PLEByoVRWE2my2i3z8ST2ch8UQQ4U97dHrsdrubYBEnzMpubVoNSXuIjI4UUB2V/Kljxd/bYrGwqKgozXMbdWa6SoeHMMZqtarWCXFeYaTTmYSgGeT3Wa+jHwmEa93ErXkAgipyyPLUCekMlShBdHbay/Kk5TPPryXO9bHZbMxisbjlob0i4mnNs4rk5M9vJs6PsFqtmsEt+NpaXZVw7YCFG/KARFRUFLNYLKHOVsB05O8fjs8a7+h3hCUj2PcvCt9Azt1e88iCKUg7W4RUEk9nIfFEEF0Tb5W63JBYrVa30bNARI6e8NIKjR7uyVvgCn86qlqdi2AJs85CV7M8+Iv8PtlsNma1WkOdrYgi2M+ar4LBKAIcr49DNefJH3EgHsPdwP0pW8BcCHHZDdpXsRaIwOtMVifGSDypkHgiiK6J6Lan1cGXR/XEQBKBhPnmFhPRDzyYAizYiVvd/D3en46qlouOeM7c3FzDDkM4jpQHG36PPLR9Z77XQJADw3SGDlxHE4z3STwHtyabrRu0xIsYLc7sufy9DyMLj78WMJ5/2aMhWPnSuhZfesPoulpikLtO+zMQpvX+RTIkns5C4okguiZyh1z8Xl7/Qhy1C1SMiOcJNHADb9DEyH96okrr+46YY+VvHStP9JctXElJSbrHdiWrTFe6V3/QE/QkOjsWrXmMZjvjRqGuFUUxfS5/35VgW54Y69g1qrQsT6LFS65rZTHI2w5FUUxdQ0Sccyi2gZEKiaezdFbx1BVGXgkiUERRwZj7Qqyi6BAbGT0hoiiKV1c+3sACP61ZFEgS8+RrCtbcKlmAyef1p1Ov9TtolbfR8V2l/utK9+oPRs8uiU7f0bIK+xoaO5iWLLMiWBY5vuShPd4xX4Sc3vUDyZdYfmJ7lJSU5FFOvJ43Eqha0RBFcSu2F5E85zBsxVNlZSXLzMxk6enpbOHChR7bXS4Xy8/PZ+np6SwrK4vV1NSo2wYPHsxGjBjBRo4caVoUdVbxRKORBOEdbrGJiYlhjLn7Z2t1Dux2u5u1iIsom81mODeHn1NuqAK1PPkinhRFURtFng9fFss1ew35O1/rWF+scgThDaNBAqP5joGEZu4ogpkPs+eS+xb+9DWCkW9frxuMfPtCMJ8hPVe7QO5BtEDxdiA5OdntvZDrYqO8yvcjRknVc9cMl/fIF8JSPLW1tbG0tDS2Y8cO1tzczLKzs9mmTZvc9lm2bBm75pprmMvlYqtWrWIXXXSRum3w4MGsqanJp2t2VvEUiQ8lQXQ0eqNlvOGQ3THk7bJg0HOD0wowwd3Q+Pf+zGHi65twMeRt4V6OL9fSspCJ9Sbw0xwkvWN9aXRFAUviqfPSHm2Ut7lyWu+OFt46peEyOOlPPvTcxczOwfHX8iTWf0ZBEsy6sxm58mlZo3zJt7/3KO7rizuh2XPK5ebLPcjfyb+3+PvwgTazwkkLrXltcpsaLu+RL4SlePrqq6/YxIkT1c9FRUWsqKjIbZ+ZM2eyN998U/2cmZnJGhoaGGMkngiC8A2jxkf2/falYyYn3hAFGuRBFBVcfDkcDo9JuXqJ593X6xq54uldmwsnbwtY6glYq9WqKUblYB1Uh0cm7dFxki3HdrvdQ/yLc+f0nh1vneWOnK9ihC/5EC3jWhHa9DrowUJ+j/XKl4srb4Mjcr3DBRMPIqF1D2Ln3cjdzx8rlSia5Doq0AECfq+yS50RWnkW3w+t++flLrdb/uRfK+CS/Hz56noZDoSleHr33XfZ3XffrX5+9dVX2ezZs932uf7669nnn3+ufp4wYQJbs2YNY4yx1NRUdsEFF7ALL7yQvfDCC6auGerKjyCI0GLkT240SqZnZYqNjfX4TnSVEzsHgSbRtY13iMwIKDEvvsx94uJPq6zE8/nSKGqVv5g/M4KTiDzay/Iki5LKFgAAIABJREFUCnajAQ69Z9kMvnSm27ND6IsAFYWFVqfWyMoiW9zN3Je8n2x50jvG7Hst1qNiJ9+o3pHrEr21luSOvxmRKpav3Db4KkTl/Ij1u9lgG3q/o5ao422dLJoCCeuvJ0C1LJuRZIEKS/FUUVHhIZ7uv/9+t32uu+46D/FUXV3NGGOsvr6eMcZYY2Mjy87OZp9++qnmdV544QV1tHbQoEHBvg2CIMIcLfcFb64Q8ve+CCAxmASvfM0ElzAz70fcx9s5xU4l/6snAkVXPHEfrcZD7rQGipgPb66IPJS61mhsZ3Zf1nIrNXNMZy0PjlwuRs+Pv8+rL53pYHYIRXHg62Kmer+9mXyK+5h179M6L8+DXBdqWZu9iQMty5M3RNEs5kG+J9n65otY1lqCwpd6SaxLeX74M5yUlOTmDqcnFo1EpCwgeTloBf7xVr+YHXgUv1MUxeO5jaR6KSzFU6BueyLz5s1jTz31lNdrkuWJILoeYmOoJYq0fLW9jeLpJd5giIKEMfcJtbIQ04vEJ+7HLTxyQyvvpzXCLN4/v19R6HFBxhttRVEMGzZv7iC+oiWQjMpYb52VSBrR9BXZrdTsMZ21PERkK4vee+nvc6pXjv5YaGS8uRLz+4qKivLJeuZPR5sju43x6xpFThMFplhX8u9F10FRvIj7eStDXwcQRGEjDk4lJSW5iSb52TAzoCb/b5Qvo/dQDO3N8wNoW8nENkRL7PG6k68tyI+R3ciNrPtGv7GekNZrD3idxfMVifVQWIqn1tZWNmTIELZz5041YMTGjRvd9vnggw/cAkaMHj2aMcbYyZMn2fHjx9X/x44dyyorK71ek8QT0ZmJpBEdbwTzXozOJbtFaLkC8b9iQ+fNSsQFEa9zZNHCEz+PnkUoMzPTrdEROwRyh1Gv4dPyoZdDyvIOhjfXFnmE0tcOupbbkBlRKt8nWZ7MHRPO5RFI/nx5bhRFCUjkewu64G+n0JsFV7ZuaNUHeogdaq08Glng5WPlAR15f/l8fFBHHgzha7XpDewYDXLxa4gdfzO/pyzU+ICWGYuemAc5cqocIMLIOqk3GMfvh59HS9iKx8rPtTxQJtbpZoLwiKJYFldaiNeRo+iJ1xR/Q18tpuFGWIonxs5E0zv33HNZWloaW7BgAWOMsfLyclZeXs4YOxOqfNasWSwtLY2NGDFCne+0Y8cOlp2dzbKzs9nw4cPVY71B4okIBuHaKelMI80ddS9anQB5ErrYGBgJJ7kxEhs42e2ERzgS3Vr0klhv6Y3+abl3yB0bsWETxRy3aIjnlMsf+Gn0WCwfM++CVmdJdJmRQ7t7S5Ho/hEsfLnnUJePt+sH8o7LHXRfkq/lYcby5A9aFmOjPMjLKshoWUZklzm9vMtWblFcyHWRaLXhZaJneRIFlDyfRrZQGdUTWnWDmWdHPKdcz3ILjZ7roPi7i8JLrv9lTwK9PMhWGVGY6ln7tRadFetgfk/ioJZRPSovPs6vLV/DCHk/USRFsoueFmErnjoaEk9EMAhXkRLpFZVIqO5FL2qQVmOg5e/OOwqi2GDMXcDI7jdyR1DutIgjgXqjk6JVQsuaJh/jLaKTkeWJN95mn3+tEWXRZUYWnN4S74QZjZ53Rnhnk/9e3tAT2h2FUT3pjyVNPt5f8eRrOOn2irYnjvhzq7G/a1Ex9tMAilgfiYt+651Pqyw5srVHFg+8TLQGW/SEq5g3sQy06hvRwi5aysX3wOzzoyVE5HvXKyN5eQhRhGidU2/wSrYyiQNuYj3G783oOea/A29nxHddSxBplbuW5wMXcnrPm+zqLt6fLPDDta9kFhJPZyHxRASDzt5J60oYuZ9421ceOWXMfW6T+L0cfU8WKmLjpuWmobU/FzBynrXmccn34UsnXOu+zYQA1jtePoevId25mJTPK3eCOst7Kne0ZEHt7Thv4ajbq5y0RDsnGJ0qf6NY6pWdr+55gZab3juglT/ZsqN1Xb371Xr3xGvZ7XZNa4oZSzCvP0SxYyRsjd51uQ4T3ci0Ah94q+eMylgsG7meNarbRAugmPiAWUxMjMfzolWO3EIjf19QUGBK9IjlwK9r9vmXB8LExO/d13ZQtA6KLu+iwI3E+pjE01lIPHknEh/wjobKyBNfyyRcytAX64W8r91u95g3pOcmw/fVmzgrdnS1GmfZpc7Idc5b2YqNtplReC2XIdHtz1/rj+y6wle8N5OMhCzvFLXXGjYdjVxOZn83jrffpL1Hh7XeMbEj50u+xY6aL9ZKnozCPnMxALh3i7yJKtFqxPPlq0VNvhf5vQZ+sjDo/V56AkHOm15nWrRCiC693sqUv2ui+6+vvws/j2jlEZNcJ3IhwfNnt7u7KPKOu9H1tOaNymJKLGdRbPA86pWnbNXT2k/vt9LLr1wussA1W+7eBqvk+zSLLJa02stQWsL9gcTTWbqKeAqkYyq+RIGeq7MS6abo9sDXMgmXMvT2fMudNa1wrHKDpCUmzIYbFq1KckOpNZlay/3J2z3JnTQzZcTzxO9L7riIHQVxP6O8mJ3jZJR4R9hud18gVYxIGOl1Fy9DX343s/W21vPjy/Fm8y93qsxEDtSqI0Qh6W8nXa8foDdhXmvwgN+XLBwA3yLjyaP2wE9BIbREs5GbLd9Hjt4plpM3McSPtVqtXudi8vPKc3eMOudG77xZa5dWJ5//PlxgmrFKcqubWH/L1xfLWrxn8bnWq5MYY27R7swkI0u8/L2vlldf9td6P/XqBLkc9DwDzFrCwwkST2fpKuIpkI6pXNGFSyc3nIhkQdmeLjq+jLiaFS2hLmOxA6M3YVt0WZFdPeTRNqOGw5tQ0yoPrY6ot3dWbIjNulOYsW5prXUiz7/SK1t/k9xQiyGJO1t9L7tGGWG23tb7fYJd74vnM3sfWs9YMAS32FHXcz/SstZoDQqI+eHvv9b+WvcmzzsRrS6ie54oCvR+F7OWDTMusuI9yPvzOpALRV+tJ0a/obclCoKd+DMvzmOSByn0oqSKg0VazxbH1zzZbLagPeNiEoWmmYXStdpyo7pC/s2NCJe23Qwkns7S2RpTPYJleeKVRKQ86IHSme9VbKz9XSyyI91/wkW0y50ksaPjy/F6a4po7SOOOnuz4PhjeRLvReta/sI7DrzxlDsXssXMWwNudi6UfF+88fc1MEAgdETdIY5+e5urZnYwQy/im7/3Ix8nT7QXB+TMWgbFc/prbdJL8to/smWMlzd/tnle5WiVsqXJyDosikejpQ/4/rIbFa+PxLLTsyqIFmy99eR8Tfw8PM9mXPt4Et3rtO63PYSDUZKD/hQUFLhZII2eGZ5nrd9MfPbDKXGLvLf9tOoPPcusXJd3pn42iaezhPuPGurOu/wSxMTEMMaYWzjMjshDqMogXDrs7QHvEIjuCb6Ws7fyMfucmOnYhfpd0MoHfz/8EZ9G9yNGyRKvpWXBCtazaWYCeiDnZMwzimAwLE3eEg8nrBe2XQutfXwtEyMLW6DoCU1R2Gs9N2aeGS2BazZPepZQOeqaHA3MqDOn1UaLZetrcBG9ZHQub5PvuWWGiwat+Uryuytah33JIx/ckK3LcrAEo8EIRVFUERVI+WmFIdcLhqB3Dj0Bx6127Vk3mL1Hrd/UbBJ/91Dfi9HzYOY+5BDuevfJn79AomeGKySezhLu4inUnXetip0xzzkSeh08ow6H2c6I3e59PYv2Ilw67FoEmje5cjMTqUi+tjfBBbiH59ZDqzOgd81AR74DRS/sdzDOL7rl8EbN277Buq/2qGt4HkURze8tJibG437bK/HnSo7+5EtZ+Fo+Zt8ns7+juJ+eVUEUo/zaWu5x3oIyyB0r+XfUWiOIl6tsPRLzrdex1uuYyvNJ2vs58SVCGU9cnGu5QYn3lZmZqZaN7Aboy/WSkpIMrQBiR9bMuQIts2CJ11Cd30zytgC6UTKzOG2ok1zGWhH8tKzccpAg8Tj+Lnmra8O5j6UFiaezhLN4Eiu2jn6wjCp0rUpZL0SyXofDjLuJ2EHn+3ak240R4fDC+zKyrSVwxd9YHDk1G23NTNhSrc6a3qi+tzVNxGfAl/IPtijQsgpp5def50Ps/MrBF9qb9nimxToiNjZWbahll6ZQdBjkwR5/nlstZLc0fs96z7bZd8nsPclilOcjKSnJrZOkZw3WEmaiVYWLBfFexAEusT2QywzwnMDvbUS/o+e9BCuJ7TcvM1ko6oWn1kv8N/AWvIE/T6Eug1CkYLtxdvZkNHjBy1MMkAF4rnnIB1H5s8zrGTOD3qE2EPgKiaezhLN4Eh/gUF7bW/InRLJoAtZ6uWR3KD3f2mBhZD3TQnzhtSqVYOTD23YtNw3xHsRIboqiqMJIdpXiv4MvCz8a5VPMF88HvxafRyDPb/B231pRhcxWuHpizV+hwI81yoO/DYJWYxRJDYtMMDs73hp5X5PD4XCrV0SBALjPZTPrgsLrAr3Rcq2BCaNnURYl3spMr9y5YNEqA7n+MjPSr2dFEcPq8/sSn99Qu2EFa46P2ec1VAJGHhyjREkv+bJ+l+wdIteLelZrI8zWreECiaezRIp4MtPZM3Kl8lUc+PoC+jpPQu5U83sUv+ejGDExMaasImbuUf5enBwqdlK0LBx6I9Na96F3PTPloTdqy5h2hCouSOTgArwjyEduuaDRck2RO0z+IP+mvPOg5QLA88vL3FvEOb3nTqty1quMxfIMhhg3Epu+CFEjwsHCGQhmOuK+CCKteqMjkijcjOogb1YEX5830eJqpnPDmKdA8dVy44uY5YMwcp3E8yMv2BlKFyx/Fg+N1GTmWaRESXS/M7O/bEHXcreXjwmFd0h7QuLpLOEsnuTK3ltUL9GiIHZKeQdMy62CI3Y4AxmxEl3rjMLParmXiJ1+3siOHDmSATAVeEC0aPC8yC6BstUD+MmlQq9y4S+1mTCvPLSpeD1xYq4cAU3vuvJvxgWxNx91LcsTT7LpXe/3A36qAvTEpuzOFEiHlo9i60WI06vYZf9qsczlQBjiMyGfx1+MKvxIagzak2A39uLARSiTuPixXl2gl8QobvwZ1+vY+zonJRw7zIHMFwl26irubO3x7lHqfEn0TPG2r958Rrndlo/TmrIRTE+QjobE01nCVTzJnVzubiALHN4QGzWyYtQa8SURRY3c4Qx2kl8grRdWK3IPoD1yKrvcmJkTIHdM2msEWxYW8nZ/r2umEyJbxuRjzXY8GfMMX83YmYrRSGz6m/gq9KJ1URbDgSSj5zqQdzTYc546G8F+TjrLfAb+ToU6H10piR4Zoc5Le6f2ePcodb5ksVjUttdoP/7uKIriNQiOOICjNXAe6QOLYSueKisrWWZmJktPT2cLFy702O5yuVh+fj5LT09nWVlZrKamxvSxWoSreNJ7iP1d5M4o2Wy2LuHK0FFJrlhCnZ9gpo6cvN0RLj5mIgsS/hOJzwSlzpfkDlxnf47a492jFLnJSByJ3kla63Npze2WkzhtQy8iJyfSBxbDUjy1tbWxtLQ0tmPHDtbc3Myys7PZpk2b3PZZtmwZu+aaa5jL5WKrVq1iF110keljtYg08UQpMhInHFyMKOkn0d2PCD7B7qSGYtFMSpGfRMjyRCnYqSPEuD/XEIWR1nbZdU5rrTMxGIzRtcRpI5FsWfKGkWawIESsXr0aGRkZSEtLg91ux/Tp07F06VK3fZYuXYo77rgDiqJgzJgxOHr0KPbt22fq2EigsLAQiqKEOhtEAFitVvX/mpqaEOaE8Mbp06dDdu3CwkJER0ejsLAwZHlobxhjQT3fnDlzUF1dDavVSvUkYYqCggK3z2VlZbBYQtbN6TDk+44koqKiQnp9X5+PYNdzWlx44YU+HzNq1CgUFxejubnZY5vD4UBxcbH6ubi4GE6nEy6XS61bFUWBzWbDxx9/DKfTqXsdh8OB/Px8TJo0CQ6HA5MmTUJ+fr7P+Y10Qlar1NfXIyUlRf2cnJyM+vp6U/uYOZazePFi5OTkICcnB01NTUG+i8AoKSkJdRYMyczMDHUWwp45c+aEOguESWJiYkJ27bKyMlRVVaGsrCxkeYgkxMZ+1KhRHdJhISKfsrIytwEK3qlTFKVTC3CxYxxptLa2enzXnr+VLNYuuOAC08dq5ctiscBut5vOs9wOyccpioJ169bpijqHw6G5rba2VvOcmZmZbtsA98E8RVEQFRUFRVFUUSS+Lw6HAwCQl5cHu92O6upqVaTl5uaipaUFJSUlOOecczr9AKFIyMSTVmMoP0R6+5g5ljNz5kxUV1ejuroaiYmJfua2a9GtWzcAwO7du0ObkQihK1UYkUxbW1vIrs0bpa44QucrFovFrbGvra2lgRzCFFOmTHEboBBH2F0uFwoKCmC322G1WpGXlxfCnAaXnJycUGchaOTl5fk0WOJwONQOvhlEq0pmZibWrVtn+lg5X1arFXPnzkV+fj6ioqJQUFDg1Qp43nnnuX2WxRxjDIqiYNq0aep38fHxan5lIcRxOp1qP+SBBx5Qv9+6datHu1NaWgqXy4XS0lLMnTsXiqJg7ty5AM60kzwPFosF1dXVcDgcqKioQFZWlsd5gDP978bGxq41QNjuToM6fPXVV2zixInq56KiIlZUVOS2z8yZM9mbb76pfs7MzGQNDQ2mjtUinOY8haMvtsPh8AgTHo75DKckhjcPdV4oGadwev87I76G2paTuKiptyiWlCjxxOdtiAv2eoPPDdFakJvX6UZLHgQjv8FKjPk356kj5hOK10hOTjZ9nNZi0XLi5zP6fXhZ8yBZ4vIn4vIYZq4nPiNiYAWtgAhyXSiu6yjnWf5OXIBWDBkuRpAUIxjz/+W1MeVnRERvLU05ErMYoVlrbpN4H+KyDJ2FsAwY0drayoYMGcJ27typBn3YuHGj2z4ffPCBW8CI0aNHmz5Wi3DqPIVjBCA9xM5LdHR0yPMZTklcCyHUeaFknDrrpNZwQZxEHIwkNsKBCjMzSRRvwUp6E6/9jWQZbgE02qPMfEm8PdJaLN4MWmvXaHWqg1XugQgxef094Ccx4E/+xHv19Xi+mL2ZfsyiRYvcjjNzfnHRZa3Et3HRYHQeM4vW663taJT4M+dNVPDfTixrvl1cVkRvPSSjbfw7q9XqEVZcTyBpHa/3PvDngv9uch860qPpeSMsxRNjZ6LpnXvuuSwtLY0tWLCAMcZYeXk5Ky8vZ4ydCVU+a9YslpaWxkaMGMHWrFljeKw3wkk8BbsREZMvlYD48Ht7EcTtWpUtr1DNVo7tce/eRreCLVq1KjRvHb1wCBVvVA5GHTt+XDA7yfy8vnYofe1Qd9YKPlzQ+j303nMz76G8ALU/zxzv4Hmrl8SRVm/PFX9OvT2vvGMkWvPF+lUuJ7FTJSZRIDD202CWOJBlVJ8qitJudY5We+DtWjzUMf9tRGuPWQu+GFmsvd5r+fxymRYUFHjcq97vIAserUXq+fOkVX58rTBv92y2XZUXoOcdbfG59fabi9eU32dx0XLx2TW7SCt/Z0URxf/m5eWpdYO8aLtYdr709+TnUbS4yFYpM5Yn8XeQy1rv+dJDT6B567uZPb/RfsBP4rczR9XTI2zFU0cTTuJJS3wwpm2R4pWstxEifn96wkbsFGiNyPi6oJlYcfEKhndwxBEPrQaGMU8ByTsPPG9anQSetEZFeOMjuhxqrUUgV7ji4m7iNm+VvJnnSRzFlMtWXPwYgNu92+120x0eXuGL+4udEj2XBHG9Bv5XzKNcqcojYXqL6fH74Q0IL1+xIRUXcua/u6IoapmLI9p6nV+xjMV74mUpPpcknNof4KdFlfXWqLNYLB4j/XrPhVFHwNtIubfBIPF4ucMkIj9bWsh5EetrrfdKLC+exG38mvwZ1hsh9mWkWuwcyx0u8f3zVs+IdYlep1BcQ8Zms7m1W+Ii3Hr3I7dbeu1DR6NV95ipn8V2x6zYE38L8fc0siYw5mnx0FvkXH6WueVCURQPASs/3/KzLD5b4m8mty//v717D4rquuMA/l0EDFVjsEIEMaEICrK7LOGlJiUIoSaaQjA0JVIgIaRDTDTkYe1Mm5Q4RLGlLQ1NY2NbxWjAOhlcJySYGJ8odgFZabStiBIeUkNQiyEFFU//0L3dXXZhQWAffD8zjux9nnt/9/U75+xd4+0ybu3R/a1blouLi3T86Ji7J+l/1l33h/KgrztmdeszbpHTJZSWVDILYfpH54fDVMuofplHIqEZaDn6zyiO3MJkDpOnW2wpeTK+SehuRMYXOXM3GuMLjv5FxtXV1aA5eChlGm7XB1MP1frL0b9I6i78xhflgR5edDeBgco3kr9mrX/R0q9htuRXt80tz1zZ9Wu3zDWjm7vp6W+zqeTHuPZO93A31DIOddtM3cQGi49uvHGXBuMHPUtvXPo1vOPtom8Nuv2uX+uqO64HSlCGU2mjq5DQv84NpeV8pAH/r6G1dL3mEkVLyjrUbTHVrcfU8kwlCENdp6lafHPXOFPLNG6pM3d9G2um9o1xEi6E+QqzoZTfuMVGx9I46p8bpn7Lxzi++nEwvi/rl0f/Idr4mmxq/+jH3dR3c0wdd6au3eaOPeN9OtC9dCD690rdturuXcO5bgy3HENZzu1czwZrtaKbmDzdYkvJk36tvLmL8UAnnn6CpLswmauNsQXmalAsqSkydyMxtY7R2G79C/RorGOgZQ7UcmU8r/FyBmtJGm0D1RIOVHM+khd0az9wjTf657P+9wEHqkEdaqz1k+yBjndjg9Xa3w5Lr1GDsXRfDDfhHKziyXifjhRz8R9sO4yvCbbW8mRcRn3Gx8RQjnNz+2UoCawusTHVam/uODWu1NWdK6bKYzzM1HOLfgWBs7OzyZcfCCEMKinNvcjAkv0xEsfISCQ+I3X/suT5aDgGS8LpJiZPt9hS8mSq1UlHl1DparJMGa2TyhpsMdnTZ83y3c5DjS3v17FKamx5HzgqU/vc1dV1WN1wLF2+JfPawwODpRUeI31cj8V5YuqcH27ibI17nvGbaAdzO2UdqURRV4bBkib96fVbrAY6L42H6faPu7u7yZYg4+8SmTrOjVvxzK3Lkm0e7jFia/eM0SjPaFYkORImT7fYSvKkX7ujeyg2Hj+U2iVbOclp5DlqjIezXY66Lxyd7kHQVNc9U91uRuvh2B6On8Faj+3ZSOx/a8ZwqMn37ZZ1JGKvO/csrWi1JEkaiKWtbeaGj0QFhz2c59Y2VpUy9o7J0y22kjyZq40hIvMc6UFyPBmoG4zx93548zbE/WE7xrq2fiS7Lg+1G5r+uody3R3ONdq4iyGfiW6yxrnPe6whJk+32EryxBsi0dDxvLFPA3UxNtedj8jW2Ov1Zzjl1n+IHsr8t7MuYHy+DtscayQy9nqMjxYmT7fYSvJERDSeDNRNhLXNRLZluA/Rt9Mdmy1PhpjIWN9AOYNMCCEwToSHh6OmpsbaxSAiIiJyKBMnTsSePXuwePFi9Pb2jvn8RCNpoJzBaYzLQkREREQOZuXKlVi8eDFWrlxplfmJxgpbnoiIiIiIzHj11VdRVFSElStXoqCgwNrFoTHAliciIiIiomEoKirCnj17UFRUZO2ikA1g8kREREREZAa7FJI+qyRPFy9eRHx8PAICAhAfH49Lly6ZnK6iogJz586Fv78/8vPzpeG5ubmYOXMmVCoVVCoVPvroo7EqOhERERGNIwUFBejt7WWXPQJgpeQpPz8fcXFxaGhoQFxcnEFipNPX14fnn38eH3/8MU6dOoWSkhKcOnVKGv/SSy9Bq9VCq9ViyZIlY1l8IiIiIiIah6ySPKnVamRkZAAAMjIysGvXrn7TaDQa+Pv7w8/PD66urkhJSYFarR7rohIREREREQGwUvJ04cIFeHl5AQC8vLzw5Zdf9pumra0Ns2bNkj77+Pigra1N+vz73/8eSqUSmZmZZrv9AcC7776L8PBwhIeHo6OjYwS3goiIiIiIxpNRS54eeughyOXyfv8sbT0y9QZ1mUwGAHjuuefQ2NgIrVYLLy8vvPLKK2aX8+Mf/xg1NTWoqamBh4fH8DaGiIiIiIjGPefRWvDevXvNjrv77rvR3t4OLy8vtLe3w9PTs980Pj4+aGlpkT63trbC29tbml/n2WefxaOPPjqCJSciIiIiIupv1JKngSQkJKC4uBg//elPUVxcjMTExH7TREREoKGhAefOncPMmTNRWlqK999/HwCkxAsAysrKIJfLLVpvU1MTwsPDR25DbkNHRwdbwhwI4+k4GEvHwng6FsbTcTCWjsXR4tnU1GR2nEyY6h83yjo7O/HEE0+gubkZ99xzD3bu3Ilp06bh/PnzyMrKkl49/tFHHyEnJwd9fX3IzMzEz372MwBAWloatFotZDIZfH198cc//lFKpuzFQL9cTPaH8XQcjKVjYTwdC+PpOBhLxzKe4mmVlqdvf/vb+Oyzz/oN9/b2NvjNpiVLlph8Dfl77703quUjIiIiIiIyZpW37REREREREdmbCbm5ubnWLsR4FRYWZu0i0AhiPB0HY+lYGE/Hwng6DsbSsYyXeFrlO09ERERERET2ht32iIiIiIiILMDkaZRVVFRg7ty58Pf3R35+fr/xQgisWrUK/v7+UCqVOH78uBVKSeZkZmbC09PT4HX4Fy9eRHx8PAICAhAfH49Lly6ZnHew2NPYamlpwaJFixAUFITg4GD87ne/A8B42quenh5ERkYiJCQEwcHB+MUvfgGA8bRnfX19CA0NlX67kbG0X76+vlAoFFCpVNJPxDCe9uny5ctITk5GYGAggoKCUFVVxVgKGjXXr18Xfn5+orGxUfT29gqlUilOnjxpME15ebl4+OGHxY0bN0RVVZWIjIy0UmnJlIMHD4ra2loRHBwsDVu9erVYv369EEKI9evXi5/85Cf95rMk9jS2zp8/L2pra4UQQnR1dYmAgABx8uRJxtNO3bh/KXshAAAN10lEQVRxQ1y5ckUIIcTVq1dFZGSkqKqqYjzt2K9//Wvx5JNPiqVLlwoheK21Z/fee6/o6OgwGMZ42qf09HSxadMmIYQQvb294tKlS+M+lmx5GkUajQb+/v7w8/ODq6srUlJSoFarDaZRq9VIT0+HTCbD/PnzcfnyZbS3t1upxGQsOjoa06ZNMximVquRkZEBAMjIyMCuXbv6zWdJ7GlseXl54b777gMATJkyBUFBQWhra2M87ZRMJsPkyZMBANeuXcO1a9cgk8kYTzvV2tqK8vJyZGVlScMYS8fCeNqfrq4uHDp0CM888wwAwNXVFXfddde4jyWTp1HU1taGWbNmSZ99fHzQ1tY25GnItly4cEH6UWYvLy98+eWX/aZhXG1bU1MT6urqEBUVxXjasb6+PqhUKnh6eiI+Pp7xtGM5OTn45S9/CSen/z+WMJb2SyaT4Xvf+x7CwsLw7rvvAmA87dHZs2fh4eGBp59+GqGhocjKykJ3d/e4jyWTp1EkTLzIUCaTDXkasj+Mq+36+uuv8fjjj6OwsBB33nmnRfMwnrZpwoQJ0Gq1aG1thUajweeff27RfIynbfnwww/h6ek5rNccM5a26ciRIzh+/Dg+/vhjvP322zh06JBF8zGetuX69es4fvw4nnvuOdTV1WHSpEkWf3fJkWPJ5GkU+fj4oKWlRfrc2toKb2/vIU9DtuXuu++Wula2t7fD09Oz3zSMq226du0aHn/8caSmpmLZsmUAGE9HcNdddyEmJgYVFRWMpx06cuQIdu/eDV9fX6SkpGDfvn340Y9+xFjaMV0MPD09kZSUBI1Gw3jaIR8fH/j4+CAqKgoAkJycjOPHj4/7WDJ5GkURERFoaGjAuXPncPXqVZSWliIhIcFgmoSEBGzduhVCCBw7dgxTp06VmkLJNiUkJKC4uBgAUFxcjMTExH7TWBJ7GltCCDzzzDMICgrCyy+/LA1nPO1TR0cHLl++DAD473//i7179yIwMJDxtEPr169Ha2srmpqaUFpaitjYWGzbto2xtFPd3d24cuWK9Pcnn3wCuVzOeNqhGTNmYNasWfjXv/4FAPjss88wb948xtJab6oYL8rLy0VAQIDw8/MTeXl5Qggh3nnnHfHOO+8IIW6+MWrFihXCz89PyOVyUV1dbc3ikpGUlBQxY8YM4ezsLGbOnCn+9Kc/ia+++krExsYKf39/ERsbKzo7O4UQQrS1tYlHHnlEmtdU7Ml6Dh8+LAAIhUIhQkJCREhIiCgvL2c87dSJEyeESqUSCoVCBAcHizfeeEMIIRhPO7d//37pbXuMpX1qbGwUSqVSKJVKMW/ePCkmjKd9qqurE2FhYUKhUIjExERx8eLFcR9LmRAmOiUSERERERGRAXbbIyIiIiIisgCTJyIiIiIiIgsweSIiIiIiIrIAkyciIiIiIiILMHkiIiIiIiKyAJMnIiIH8e9//xspKSmYPXs25s2bhyVLluD06dM4cOAAHn30UWsXb8xptVosWLAAwcHBUCqV2LFjhzQuJiYGc+fORUhICCIiIqDVaketHDU1NVi1ahUA4MCBAzh69Kg0buPGjdi6deuIrKe9vV2K84EDByCTyfDnP/9ZGl9XVweZTIaCggIAN/dBTU2Nxcvv6OjAww8/PCJlJSKyV0yeiIgcgBACSUlJiImJQWNjI06dOoV169bhwoULFi+jr69vFEs49r71rW9h69atOHnyJCoqKpCTkyP9sC4AbN++HSdOnMCKFSuwevXqUStHeHg43nrrLQD9k6fs7Gykp6ePyHp+85vf4Nlnn5U+KxQKg4SxtLQUISEhw16+h4cHvLy8cOTIkdsqJxGRPWPyRETkAPbv3w8XFxdkZ2dLw1QqFb773e8CAL7++mskJycjMDAQqamp0P3En6+vL9auXYsHHngAO3fuRElJCRQKBeRyOdasWSMta/LkyVizZg3CwsLw0EMPQaPRICYmBn5+fti9ezcAoKenB08//TQUCgVCQ0Oxf/9+AMDJkycRGRkJlUoFpVKJhoYGADcf9uVyOeRyOQoLCwEATU1NCAwMRFZWFuRyOVJTU7F3717cf//9CAgIgEajAQB0d3cjMzMTERERCA0NhVqt7rdP5syZg4CAAACAt7c3PD090dHR0W+6BQsWoK2tzeR+9fX1xZo1axAZGYnIyEicOXMGAPDFF18gLi4OSqUScXFxaG5uBgDs3LkTcrkcISEhiI6OBgCp5a+pqQkbN27Eb3/7W6hUKhw+fBi5ublSS5BWq8X8+fOhVCqRlJSES5cuAbjZQqQrw5w5c3D48GGTZf3ggw8MWobuuece9PT04MKFCxBCoKKiAo888ojBPNu2bcPChQshl8ulfZubm4u0tDTExsYiICAAmzZtkqZ/7LHHsH37dpPrJyIaD5g8ERE5gM8//xxhYWFmx9fV1aGwsBCnTp3C2bNnDVoP7rjjDlRWViI6Ohpr1qzBvn37oNVqUV1djV27dgG4mazExMSgtrYWU6ZMwc9//nN8+umnKCsrw+uvvw4AePvttwEAf//731FSUoKMjAz09PRg48aNePHFF6HValFTUwMfHx/U1tZi8+bN+Nvf/oZjx45h06ZNqKurAwCcOXMGL774Iurr6/HPf/4T77//PiorK1FQUIB169YBAN58803Exsaiuroa+/fvx+rVq9Hd3W12+zUaDa5evYrZs2f3G1dRUYHHHnvM7Lx33nknNBoNXnjhBeTk5AAAXnjhBaSnp6O+vh6pqalSt7y1a9diz549OHHihJRU6vj6+iI7OxsvvfQStFqtlNjqpKenY8OGDaivr4dCocAbb7whjbt+/To0Gg0KCwsNhuucO3cO7u7umDhxosHw5ORk7Ny5E0ePHsV9993Xb3x3dzeOHj2KP/zhD8jMzJSG19fXo7y8HFVVVVi7di3Onz8P4GYrmrnkjYhoPGDyREQ0DkRGRsLHxwdOTk5QqVRoamqSxv3whz8EAFRXVyMmJgYeHh5wdnZGamoqDh06BABwdXWVWjUUCgUefPBBuLi4QKFQSMuqrKxEWloaACAwMBD33nsvTp8+jQULFmDdunXYsGEDvvjiC7i5uaGyshJJSUmYNGkSJk+ejGXLlkkP5d/5znegUCjg5OSE4OBgxMXFQSaTGazrk08+QX5+PlQqFWJiYtDT0yO1/hhrb29HWloaNm/eDCen/9/2UlNT4ePjgw0bNmDlypVm992TTz4p/V9VVQUAqKqqwvLlywEAaWlpqKysBADcf//9eOqpp7Bp06YhdYP8z3/+g8uXL+PBBx8EAGRkZEj7HgCWLVsGAAgLCzOInf42enh49Bv+xBNPSC2Kuu0wtW3R0dHo6uqSujUmJibCzc0N06dPx6JFi6RWKU9PTymRIiIaj5g8ERE5gODgYNTW1podr9/iMGHCBFy/fl36PGnSJACQuvKZ4uLiAplMBgBwcnKSlufk5CQty9z8y5cvx+7du+Hm5obFixdj3759A65Lv6wDreuDDz6AVquFVqtFc3MzgoKC+i2rq6sLS5cuRV5eHubPn28wbvv27Th37hyWL1+O559/3mx5dNtt/LepaTZu3Ii8vDy0tLRApVKhs7PT7HKHQrcPjGOn4+bmhp6enn7DZ8yYARcXF3z66aeIi4szW27jz+aG9/T0wM3NbXgbQUTkAJg8ERE5gNjYWPT29hp8P6W6uhoHDx60eBlRUVE4ePAgvvrqK/T19aGkpERqCbFEdHS09H2Y06dPo7m5GXPnzsXZs2fh5+eHVatWISEhAfX19YiOjsauXbvwzTffoLu7G2VlZf26sQ1k8eLFKCoqkpIwXZc/fVevXkVSUhLS09Pxgx/8wORyXFxckJeXh2PHjuEf//iHyWl0L13YsWMHFixYAABYuHAhSktLAdxMwh544AEAQGNjI6KiorB27VpMnz4dLS0tBsuaMmUKrly50m8dU6dOhbu7u9T69t577w1p38+ZM8dkixRwsyvhhg0bMGHCBLPbVllZialTp2Lq1KkAALVajZ6eHnR2duLAgQOIiIgAcDOucrnc4nIRETkaZ2sXgIiIbp9MJkNZWRlycnKQn5+PO+64A76+vigsLDT7MgRjXl5eWL9+PRYtWgQhBJYsWYLExESLy7BixQpkZ2dDoVDA2dkZW7ZswcSJE7Fjxw5s27YNLi4umDFjBl5//XVMmzYNTz31FCIjIwEAWVlZCA0NNZsAGHvttdeQk5MDpVIJIQR8fX3x4YcfGkzz17/+FYcOHUJnZye2bNkCANiyZQtUKpXBdG5ubnjllVdQUFBg8Gpvnd7eXkRFReHGjRsoKSkBALz11lvIzMzEr371K3h4eGDz5s0AgNWrV6OhoQFCCMTFxSEkJMQggf3+97+P5ORkqNVqFBUVGaynuLgY2dnZ+Oabb+Dn5yct0xKTJk3C7NmzcebMGfj7+xuMW7hwodn53N3dsXDhQnR1deEvf/mLNDwyMhJLly5Fc3MzXnvtNXh7ewO4+WKSpUuXWlwuIiJHIxMD9Z0gIiIax3x9fVFTU4Pp06dbuyiDKisrQ21tLfLy8m5rObm5uZg8eTJeffXVfuOio6OhVqvh7u5+W+sgIrJXbHkiIiJyAElJSSP2HStTOjo68PLLLzNxIqJxjS1PREREREREFuALI4iIiIiIiCzA5ImIiIiIiMgCTJ6IiIiIiIgswOSJiIiIiIjIAkyeiIiIiIiILMDkiYiIiIiIyAL/AwZxtSIOGN9rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here's a scan using the same populations but on a different chromosome:\n",
    "pbs('2R', 'BFgam', 'UGgam', 'GW')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BFcol no. samples: 75\n",
      "AOcol no. samples: 78\n",
      "GW no. samples: 91\n",
      "finished allele counting\n",
      "found 436137 segregating variants\n",
      "finished locating segregating variants\n"
     ]
    }
   ],
   "source": [
    "# Here's a scan using different populations:\n",
    "pbs('X', 'BFcol', 'AOcol', 'GW')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}