ipynb/tutorials/RPCA-moving-pictures-standalone-cli-and-api.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial \n",
    "\n",
    "**Note**: This guide assumes you have installed [gemelli](https://github.com/biocore/gemelli).\n",
    "\n",
    "\n",
    "## Introduction \n",
    "\n",
    "In this tutorial you will learn how to interpret and perform Robust Aitchison PCA. The focus of this tutorial is compositional beta diversity. There are many beta diversity metrics that have been proposed, all with varying benefits on varying data structures. However, presence/absence metric often prove to give better results than those that rely on abundances (i.e. unweighted vs. weighted UniFrac). One component of this phenomenon is the interpretation of relative abundances can provide spurious results (see [the differential abundance analysis introduction](https://docs.qiime2.org/2019.1/tutorials/gneiss/). One solution to this problem is to use a compositional distance metric such as Aitchison distance. \n",
    "\n",
    "\n",
    "As a toy example let’s build three taxa. These three taxa represent common distributions we see in microbiome datasets. Where the first taxon is increasing exponentially across samples, this is a trend that we would be interested in. However, taxon 2 and 3 have much higher counts and taxon 3 is randomly fluctuating across samples.  \n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/7/72ebdf6a3303ce0a5850ce52a46befac564cc26d.png)\n",
    "\n",
    "In our distances below we have Euclidean, Bray-Curtis, Jaccard, and Aitchison distances (from left to right). We can see that the abundance based metrics Euclidean and Bray-Curtis are heavily influenced by the abundance of taxon 3 and seem to randomly fluctuate. In the presence/absence metric, Jaccard, we see that the distance saturates to one very quickly. However, in the Aitchison distance we see a linear curve representing taxon 1. The reason the distance is linear is because Aitchison distance relies on log transforms (the log of an exponential taxon 1 is linear). \n",
    "\n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/b/bc002a51edcd3e34cba1874a6aa97d7d08b6c0b5.png)\n",
    "\n",
    "From this toy example, it is clear that Aitchison distance better accounts for the proportions. However, we made the unrealistic assumption in our toy example that there were no zero counts. In real microbiome datasets there are a large number of zeros (i.e. sparsity). Sparsity complicates log ratio transformations because the log-ratio of zero is undefined. To solve this pseudo counts, that can skew results, are commonly used (see [Naught all zeros in sequence count data are the same](https://www.biorxiv.org/content/10.1101/477794v1)). \n",
    "\n",
    "Robust Aitchison PCA solves this problem in two steps:\n",
    "\n",
    "**1.** Compostional preprocessing using the centered log ratio transform on only the non-zero values of the data (no pseudo count)\n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/4/43fe1323791b5cea419e0973b8983621dbf31a20.gif)\n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/1/13b8c6f415d6ab10c81dec1a27f1f24079be398f.gif)\n",
    "\n",
    "**2.** Dimensionality reduction through Robust PCA on only the non-zero values of the data ( [matrix completion]( https://arxiv.org/pdf/0906.2027.pdf)). \n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/a/a327d5600f68b96457c227c660f533e94ee68341.gif)\n",
    "\n",
    "\n",
    "To demonstrate this in action we will run an example dataset below, where the output can be viewed as a compositional biplot through emperor. \n",
    "\n",
    "## Example \n",
    "\n",
    "\n",
    "In this example we will use Robust Aitchison PCA via gemelli on the “Moving Pictures” dataset. The dataset consists of human microbiome samples from two individuals at four body sites at five timepoints, the first of which immediately followed antibiotic usage ([Caporaso et al. 2011](https://www.ncbi.nlm.nih.gov/pubmed/21624126)). The tables needed for this tutorial are provided below:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "##### Table [download](https://github.com/biocore/DEICODE/blob/master/ipynb/tutorials/qiime2-moving-pictures-tutorial/table.biom?raw=true)\n",
    "**save as:** table.biom \n",
    "\n",
    "##### Sample Metadata [download](https://github.com/biocore/DEICODE/blob/master/ipynb/tutorials/qiime2-moving-pictures-tutorial/sample-metadata.tsv?raw=true)\n",
    "**save as:** sample-metadata.tsv\n",
    "\n",
    "##### Feature Metadata  [download](https://github.com/biocore/DEICODE/blob/master/ipynb/tutorials/qiime2-moving-pictures-tutorial/taxonomy.tsv?raw=true)\n",
    "**save as:** taxonomy.tsv\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cd qiime2-moving-pictures-tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using table.biom in [biom format](http://biom-format.org/), a raw count table, we will generate our beta diversity ordination file. There are a few parameters to gemelli that we may want to consider. The first is filtering cutoffs, these are `min-feature-count`, `min-sample-count`, and `min-feature-frequency`. Both min-feature-count and min-sample-count accept integer values and remove feature or samples, respectively, with sums below this cutoff. The feature cut-off is useful in the case that features with very low total counts among all samples represent contamination or chimeric sequences. The sample cut off is useful for the case that some sample received very few reads relative to other samples. The min-feature-frequency can be useful to remove features that only appear in a small portion of samples, which may be difficult to further asses using tools like [Qurro](https://github.com/biocore/qurro).\n",
    "\n",
    "**Note:** it is _not_ recommended to bin your features by taxonomic assignment (i.e. by genus level). \n",
    "**Note:** it is _not_ recommended to rarefy your data before using gemelli. \n",
    "\n",
    "The other main parameter of the gemelli is the number of components to use through `n-components` (i.e. the rank). gemelli relies on a low-rank assumption and therefore it is recommended to choose a value between 2 and 10. \n",
    "\n",
    "Now that we understand the acceptable parameters, we are ready to run gemelli.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "!gemelli rpca \\\n",
    "    --in-biom qiime2-moving-pictures-tutorial/table.biom \\\n",
    "    --output-dir qiime2-moving-pictures-tutorial/standalone-cli"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same command can be run using the python API through the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cameronmartino/Dropbox/bin/gemelli/gemelli/preprocessing.py:426: RuntimeWarning: divide by zero encountered in log\n",
      "  mat = np.log(matrix_closure(mat))\n"
     ]
    }
   ],
   "source": [
    "from biom import load_table\n",
    "from gemelli.rpca import rpca\n",
    "\n",
    "# import the data table\n",
    "table = load_table('qiime2-moving-pictures-tutorial/table.biom')\n",
    "# perform RPCA\n",
    "ordination, distance = rpca(table, min_sample_count=500)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To ensure the impact of differences in sequencing depth or sample sum differences are not influencing the sample distances an additional QC step can be run. To do this we will first rarefy the data at the same depth as the sample sum cut-off used above. Then we will run mantel tests between the rarefy and unrarefy distance vs. the absolute difference in sample sums.  This can be done easily through the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5571106600082708, 0.5814522872892138)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from gemelli.utils import qc_rarefaction\n",
    "\n",
    "rare_ordination, rare_distance = rpca(table.copy().subsample(500))\n",
    "t, p = qc_rarefaction(table, rare_distance, distance)\n",
    "(t, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see here that the difference between the rarefy distance vs. unrarefy distance in the correlation to the absolute difference in sample sums is non-significant. If it was significant you would need to continue on using _only_ the rarefy results but since it is not here we will continue with the unrarefy. If you are unsure about your results, RPCA works with rarefy data and you continue with that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the standalone tool is two-fold an [OrdinationResults](http://scikit-bio.org/docs/0.5.1/generated/generated/skbio.stats.ordination.OrdinationResults.html) file and a [DistanceMatrix](http://scikit-bio.org/docs/0.5.1/generated/generated/skbio.stats.distance.DistanceMatrix.html) both in [scikit-bio](http://scikit-bio.org/docs/0.5.1/index.html) format. Both the output files can both be imported and visualized from the standalone command or directly used from the python API. First we will plot the ordination biplot to explore sample separation and the features driving the difference.\n",
    "\n",
    "First we need to import the metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "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>BarcodeSequence</th>\n",
       "      <th>LinkerPrimerSequence</th>\n",
       "      <th>BodySite</th>\n",
       "      <th>Year</th>\n",
       "      <th>Month</th>\n",
       "      <th>Day</th>\n",
       "      <th>Subject</th>\n",
       "      <th>ReportedAntibioticUsage</th>\n",
       "      <th>DaysSinceExperimentStart</th>\n",
       "      <th>Description</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>#SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>L1S8</th>\n",
       "      <td>AGCTGACTAGTC</td>\n",
       "      <td>GTGCCAGCMGCCGCGGTAA</td>\n",
       "      <td>gut</td>\n",
       "      <td>2008</td>\n",
       "      <td>10</td>\n",
       "      <td>28</td>\n",
       "      <td>subject-1</td>\n",
       "      <td>Yes</td>\n",
       "      <td>0</td>\n",
       "      <td>subject-1.gut.2008-10-28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>L1S57</th>\n",
       "      <td>ACACACTATGGC</td>\n",
       "      <td>GTGCCAGCMGCCGCGGTAA</td>\n",
       "      <td>gut</td>\n",
       "      <td>2009</td>\n",
       "      <td>1</td>\n",
       "      <td>20</td>\n",
       "      <td>subject-1</td>\n",
       "      <td>No</td>\n",
       "      <td>84</td>\n",
       "      <td>subject-1.gut.2009-1-20</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          BarcodeSequence LinkerPrimerSequence BodySite  Year Month Day  \\\n",
       "#SampleID                                                                 \n",
       "L1S8         AGCTGACTAGTC  GTGCCAGCMGCCGCGGTAA      gut  2008    10  28   \n",
       "L1S57        ACACACTATGGC  GTGCCAGCMGCCGCGGTAA      gut  2009     1  20   \n",
       "\n",
       "             Subject ReportedAntibioticUsage DaysSinceExperimentStart  \\\n",
       "#SampleID                                                               \n",
       "L1S8       subject-1                     Yes                        0   \n",
       "L1S57      subject-1                      No                       84   \n",
       "\n",
       "                        Description  \n",
       "#SampleID                            \n",
       "L1S8       subject-1.gut.2008-10-28  \n",
       "L1S57       subject-1.gut.2009-1-20  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "from assets.plotting_helper import biplot\n",
    "%matplotlib inline\n",
    "\n",
    "# import the sample metadata\n",
    "mf = pd.read_csv('qiime2-moving-pictures-tutorial/sample-metadata.tsv',\n",
    "                 sep='\\t', index_col=0).drop(['#q2:types'], axis=0)\n",
    "# import the taxonomy metadata\n",
    "tf = pd.read_csv('qiime2-moving-pictures-tutorial/taxonomy.tsv',\n",
    "                 sep='\\t', index_col=0)\n",
    "# view\n",
    "mf.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can view the ordination loadings through the `ordination` result of type OrdinationResults."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "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>PC1</th>\n",
       "      <th>PC2</th>\n",
       "      <th>PC3</th>\n",
       "      <th>PC4</th>\n",
       "      <th>PC5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>L1S105</th>\n",
       "      <td>-0.090263</td>\n",
       "      <td>0.202402</td>\n",
       "      <td>0.293835</td>\n",
       "      <td>0.273304</td>\n",
       "      <td>-0.118611</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>L1S140</th>\n",
       "      <td>-0.351046</td>\n",
       "      <td>0.041167</td>\n",
       "      <td>-0.048862</td>\n",
       "      <td>-0.086615</td>\n",
       "      <td>-0.107419</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             PC1       PC2       PC3       PC4       PC5\n",
       "L1S105 -0.090263  0.202402  0.293835  0.273304 -0.118611\n",
       "L1S140 -0.351046  0.041167 -0.048862 -0.086615 -0.107419"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the sample loadings\n",
    "spca_df = ordination.samples\n",
    "spca_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "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>PC1</th>\n",
       "      <th>PC2</th>\n",
       "      <th>PC3</th>\n",
       "      <th>PC4</th>\n",
       "      <th>PC5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4b5eeb300368260019c1fbc7a3c718fc</th>\n",
       "      <td>-0.076566</td>\n",
       "      <td>0.021455</td>\n",
       "      <td>0.272938</td>\n",
       "      <td>0.357056</td>\n",
       "      <td>-0.256331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fe30ff0f71a38a39cf1717ec2be3a2fc</th>\n",
       "      <td>0.108060</td>\n",
       "      <td>-0.240542</td>\n",
       "      <td>-0.170898</td>\n",
       "      <td>-0.053301</td>\n",
       "      <td>-0.221769</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                       PC1       PC2       PC3       PC4  \\\n",
       "4b5eeb300368260019c1fbc7a3c718fc -0.076566  0.021455  0.272938  0.357056   \n",
       "fe30ff0f71a38a39cf1717ec2be3a2fc  0.108060 -0.240542 -0.170898 -0.053301   \n",
       "\n",
       "                                       PC5  \n",
       "4b5eeb300368260019c1fbc7a3c718fc -0.256331  \n",
       "fe30ff0f71a38a39cf1717ec2be3a2fc -0.221769  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the feature loadings\n",
    "fpca_df = ordination.features\n",
    "fpca_df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can directly add the metadata to the loadings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "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>PC1</th>\n",
       "      <th>PC2</th>\n",
       "      <th>PC3</th>\n",
       "      <th>PC4</th>\n",
       "      <th>PC5</th>\n",
       "      <th>BarcodeSequence</th>\n",
       "      <th>LinkerPrimerSequence</th>\n",
       "      <th>BodySite</th>\n",
       "      <th>Year</th>\n",
       "      <th>Month</th>\n",
       "      <th>Day</th>\n",
       "      <th>Subject</th>\n",
       "      <th>ReportedAntibioticUsage</th>\n",
       "      <th>DaysSinceExperimentStart</th>\n",
       "      <th>Description</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>L1S105</th>\n",
       "      <td>-0.090263</td>\n",
       "      <td>0.202402</td>\n",
       "      <td>0.293835</td>\n",
       "      <td>0.273304</td>\n",
       "      <td>-0.118611</td>\n",
       "      <td>AGTGCGATGCGT</td>\n",
       "      <td>GTGCCAGCMGCCGCGGTAA</td>\n",
       "      <td>gut</td>\n",
       "      <td>2009</td>\n",
       "      <td>3</td>\n",
       "      <td>17</td>\n",
       "      <td>subject-1</td>\n",
       "      <td>No</td>\n",
       "      <td>140</td>\n",
       "      <td>subject-1.gut.2009-3-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>L1S140</th>\n",
       "      <td>-0.351046</td>\n",
       "      <td>0.041167</td>\n",
       "      <td>-0.048862</td>\n",
       "      <td>-0.086615</td>\n",
       "      <td>-0.107419</td>\n",
       "      <td>ATGGCAGCTCTA</td>\n",
       "      <td>GTGCCAGCMGCCGCGGTAA</td>\n",
       "      <td>gut</td>\n",
       "      <td>2008</td>\n",
       "      <td>10</td>\n",
       "      <td>28</td>\n",
       "      <td>subject-2</td>\n",
       "      <td>Yes</td>\n",
       "      <td>0</td>\n",
       "      <td>subject-2.gut.2008-10-28</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             PC1       PC2       PC3       PC4       PC5 BarcodeSequence  \\\n",
       "L1S105 -0.090263  0.202402  0.293835  0.273304 -0.118611    AGTGCGATGCGT   \n",
       "L1S140 -0.351046  0.041167 -0.048862 -0.086615 -0.107419    ATGGCAGCTCTA   \n",
       "\n",
       "       LinkerPrimerSequence BodySite  Year Month Day    Subject  \\\n",
       "L1S105  GTGCCAGCMGCCGCGGTAA      gut  2009     3  17  subject-1   \n",
       "L1S140  GTGCCAGCMGCCGCGGTAA      gut  2008    10  28  subject-2   \n",
       "\n",
       "       ReportedAntibioticUsage DaysSinceExperimentStart  \\\n",
       "L1S105                      No                      140   \n",
       "L1S140                     Yes                        0   \n",
       "\n",
       "                     Description  \n",
       "L1S105   subject-1.gut.2009-3-17  \n",
       "L1S140  subject-2.gut.2008-10-28  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# merge the sample metadata and loadings\n",
    "spca_df = pd.concat([spca_df, mf.reindex(spca_df.index)],\n",
    "                    axis=1, sort=True)\n",
    "spca_df.head(2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "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>PC1</th>\n",
       "      <th>PC2</th>\n",
       "      <th>PC3</th>\n",
       "      <th>PC4</th>\n",
       "      <th>PC5</th>\n",
       "      <th>Taxon</th>\n",
       "      <th>Confidence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0160e14a78b18b903618f11bc732746e</th>\n",
       "      <td>0.043439</td>\n",
       "      <td>-0.028137</td>\n",
       "      <td>0.065841</td>\n",
       "      <td>0.186956</td>\n",
       "      <td>0.022391</td>\n",
       "      <td>k__Bacteria; p__Verrucomicrobia; c__Verrucomic...</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>01b99cb344ed2530f7d80897ffe257a9</th>\n",
       "      <td>0.008712</td>\n",
       "      <td>0.005172</td>\n",
       "      <td>-0.003146</td>\n",
       "      <td>-0.001097</td>\n",
       "      <td>-0.005254</td>\n",
       "      <td>k__Bacteria; p__Proteobacteria; c__Betaproteob...</td>\n",
       "      <td>0.777604</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                       PC1       PC2       PC3       PC4  \\\n",
       "0160e14a78b18b903618f11bc732746e  0.043439 -0.028137  0.065841  0.186956   \n",
       "01b99cb344ed2530f7d80897ffe257a9  0.008712  0.005172 -0.003146 -0.001097   \n",
       "\n",
       "                                       PC5  \\\n",
       "0160e14a78b18b903618f11bc732746e  0.022391   \n",
       "01b99cb344ed2530f7d80897ffe257a9 -0.005254   \n",
       "\n",
       "                                                                              Taxon  \\\n",
       "0160e14a78b18b903618f11bc732746e  k__Bacteria; p__Verrucomicrobia; c__Verrucomic...   \n",
       "01b99cb344ed2530f7d80897ffe257a9  k__Bacteria; p__Proteobacteria; c__Betaproteob...   \n",
       "\n",
       "                                  Confidence  \n",
       "0160e14a78b18b903618f11bc732746e    1.000000  \n",
       "01b99cb344ed2530f7d80897ffe257a9    0.777604  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# merge the feature metadata and loadings\n",
    "fpca_df = pd.concat([fpca_df, tf.reindex(fpca_df.index)],\n",
    "                    axis=1, sort=True)\n",
    "fpca_df.head(2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will plot both loadings as a biplot by using seaborn and matplotlib. To do this we will first make a plotting helper function called `biplot` which is located in the assets directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAFzCAYAAAAqiZ1iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd1hT1xsH8O9NQsIIBMIWkSEzgMgQRVFwi7ta96pd1g7rbLVDf7W7jlpbtRZXHQjWusBVF0OtA1xsRGWD7DASQsb9/RGlIEOmYZzP8/gk3HvuvW8AyZtzzzkvRdM0CIIgCIIgiNoYqg6AIAiCIAiiIyJJEkEQBEEQRD1IkkQQBEEQBFEPkiQRBEEQBEHUgyRJBEEQBEEQ9SBJEkEQBEEQRD1Yqg6gPRgYGNCWlpaqDoMgCKLTiI6OLqBp2lDVcRBER9IlkyRLS0tERUWpOgyCIIhOg6KoNFXHQBAdDbndRhAEQRAEUQ+SJBEEQRAEQdSDJEkEQRAEQRD1IEkSQRAEQRBEPUiSRBAEQRAEUQ+SJBEEQRAEQdSDJEkEQRAEQRD1IEkSQRAEQRBEPUiSRBAEQRAEUQ+SJBEEQRAEQdSDJEkEQRAEQRD1IEkSQRAEQRBEPUiSRBAEQRAEUQ+SJBEEQRAEQdSDpeoACKKjkeflQZ6RCahzwDQ2AdNAX9UhEQRBECpAkiSCqEGem4v8iZMhz8oCAKh5uEN/924wDQ1UHBlBEATxqpEkiSCeoWUylO/7szpBAgBp9B1U3b0DjVGjVBgZQahGdHS0EYvF2gXAGWR4RmekABArk8ne9vDwyFN1MJ0RSZII4hlaJoPs0aM622VPUl99MATRAbBYrF0mJiaOhoaGxQwGg1Z1PETzKBQKKj8/X5Cbm7sLwERVx9MZkU8GBPEMQ10dWnPnvrCRAfVRI1UTEEGonrOhoWEpSZA6JwaDQRsaGgqh7AkkWoAkSQRRg5prH+hu/QUsgSPYHu4w+CsYTCMjVYdFEKrCIAlS5/bs50fe61uI3G4jiBqYurrQnPIa1H2HAEwmmHp6qg6JIAiCUBGSXRLECyiKAtPAgCRIBKFiXl5e9hEREZpNbT916lTLvXv3quQ/bmhoqPbQoUNt6tt37do1jRkzZljUt8/MzMwlJyenwQ6LgoIC5g8//GD4susPHDjQLj8/n9n0iImmIEkSQRAEQbSjb775xnTp0qUtml1WWFjI3L1790vv+c+aNatw48aNL02miOYhSRJBEN1GlVSOgrJKJOeUIq+0EiKJTNUhdSqlpaUMPz8/G3t7e4Gtra1TQECAHgCsXLnS1NnZ2dHW1tZp1qxZFgqFAoCyJ+itt94y9/T0tLe2tnYKDw/XHDVqVG8LCwvnJUuW9ACApKQktpWVldOUKVMs7ezsBGPGjLEuKyur89507Ngxnb59+zoIBAJHf39/a6FQ2Oj7V2RkpGa/fv3snZycHH18fGzT0tLUACA8PFzTzs5O0LdvX4dFixb1tLW1dQIAmUyGRYsW9XR2dna0s7MTbNiwwQBQ9hB5eXnZjxkzxtrKyspp4sSJVs9f39GjR3WsrKycPDw87I8ePapbXxzFxcWMhIQETW9vbzEA5ObmMgcNGmTr6OgomD17tgVN/zfk63//+5+xra2tk62trdP69euNAGDFihU9MzIyOA4ODoJFixb1TEtLU/P09LR3cHAQ2NraOp07d44LADNnziw5duwYWfm2jZEkiSCIbkGhoBGfXYqpv0Ri/u//4rWfI3Al/ikqpSRRaqpjx47pmJiYSJOSkuIfPnwYN2XKlFIAWLVqVV5sbGzCw4cP48RiMSMoKIj3/Bg2m62IiopKWrhwYf60adNsAgIC0hMTE+OCg4MNcnNzmQCQmpqq/t577+UnJyfHa2trKzZs2FCrRyQnJ4f13XffmUZERCTHx8cnuLu7i77++mvjhuKUSCTUkiVLep08efJRXFxcwoIFCwpWrlxpBgBvv/221bZt29Lu3buXyGQyqzOULVu2GPB4PHlsbGzC/fv3E/7880/DxMRENgAkJCRobNu2LSMlJSUuPT2dc+HCBa5IJKI+/PBDy1OnTqXcvn07KS8vT62+WK5evaplb28vfv716tWre3h7e5cnJCTET5w4sSQnJ4cNKJO6wMBA/ejo6ISoqKiE/fv3G167dk1j06ZNmebm5pLExMT4nTt3Zu7Zs4c/fPhwYWJiYnxCQkJc//79RQBgaGgor6qqop5/T4m2QZIkgiC6hWJRFb469gASqbIXQK6g8dPpeJSJX12SJJFJkF2ejUMJB3D6cSiKKote2bXbgru7uzgyMlJn8eLFZufOnePq6+vLAeDs2bPaffr0cbCzsxNcv35dOzY2VuP5Ma+99loJALi6uoptbGzEFhYWUg0NDdrc3Fzy+PFjNgCYmJhUjRo1qgIA5s2bV3j9+nVuzeuGhYVpPXr0SN3Ly8vBwcFBEBQUpJ+ens5uKM4HDx5wHj58qDFs2DA7BwcHwYYNG0yzs7PVCgoKmBUVFYyRI0dWAMCCBQuqfwAXL17UOXLkiL6Dg4PAzc3Nsbi4mBUfH68OAC4uLhW9e/eWMplMODk5iR49esS+d++ees+ePSUuLi4SBoOBOXPmFNYXS1ZWlpq+vr70+dc3btzQfvPNNwsBYObMmUIdHR35s9fIHTt2bImOjo6Cx+Mpxo0bV3zlyhXtF883YMCAisOHDxssX768x61btzT09PQUz/fp6+vLGvu+EM1HZrcRBNEt0DSNnJLKWtskUgUkMkUDR7S97IosLA9bCjktBwD8/fAoNvn+DD31zjFJoE+fPpI7d+7E//3337zPP//c7OLFi6Xr16/PXbFihcXNmzfjbWxspMuXL+9RWVlZ/QFcXV2dBgAGgwEOh1Pdc8NgMCCTyShAOVmiphe/pmkaPj4+pSEhIU+aEidN05SNjY343r17iTW3NzawmaZpatOmTelTp04trbk9NDRUu2bcTCazwbjro6mpqZBIJLU6JBiMuv0TNW+7Ncbf3788IiIi6e+//+a98cYbVkuWLHn64YcfFgLKHjRNTc1X9wvdDZCeJIIgugU2i4kBNrWHbJjpaUCD/WruTohlYgQlHq5OkACgQJyPpOLERo7qWFJTU9W0tbUV77//ftHSpUuf3rt3T1MkEjEAwMTERCYUChkhISHNzvhycnLYFy9e1AKAwMBA/sCBA8tr7vfz86uIiorixsbGcgCgrKyM8eDBA05D5+vTp09lUVER6/k5JRIJFRUVpW5oaCjX0tJSXLp0SQsADhw4wH9+zMiRI4U7duwwlEgkFKDsjSotLW3wPbJv376VmZmZ7Li4OA4ABAUF8etr5+LiUpmamlod64ABA8r27NmjDwBHjhzRKS0tZQLAsGHDys+cOaNbVlbGKC0tZZw5c0Zv6NChZTweT15RUVEdR3JyMtvMzEy6YsWKgrlz5xbcuXNHEwAUCgXy8/PV7O3tJQ3FTDQf6UkiCKJb0NFQwxeTnfHLuSTcfFQIhx7a+GS8E/har+juBA3I6Lq39hR05/ngHx0drbFmzZqeDAYDLBaL3r59e5qBgYF8zpw5+QKBwKlnz55Vrq6uFc09r7W1deWePXv033//fQsrKyvJypUr82vu79Gjh2znzp2pM2fOtK6qqqIAYN26dVl9+vSpNyFQV1eng4KCHi1ZsqRXWVkZUy6XU4sXL37q6elZuXPnztT33nvPQlNTUzFo0KAybW1tOQAsW7asIDU1lePi4uJI0zTF5/OlZ86cqVun6BlNTU36119/TRs/frwNn8+X9e/fvzwhIUHjxXZubm6VZWVlzOLiYoaenp7ihx9+yJ46daq1QCBw9Pb2Ljc1Na0CAB8fH9Hs2bML3d3dHQFg3rx5+YMGDRIDgIeHR7mtra3TsGHDhM7OzuKtW7easFgsWlNTU37o0KEnAHD16lVNNze3CjW1eodGES1ENbWLrzPx9PSko6KiVB0GQRAdUIVEBnGVDGwWAzoar3b4RnJxElaFrwAN5d9dXY4utgz9FXz1ejshXimKoqJpmvasue3+/fuprq6uBe153aSkJPb48eNtHz58GNee13lOKBQyeDyeAgA+++wzk5ycHLW9e/dmtOc1v/rqKyNtbW3F8uXL2+17uXDhQvPJkyeXTJo0qezFfffv3zdwdXW1bK9rd2WkJ4kgiG5Fi8OCFkc1f/rMub2wddhvOP04FLocPYy2HA09TucYj9RVHDlyhLdp0yZTuVxOmZmZSQIDA1Pb+5qrVq3Kb+9FLp2dncX1JUhE65CeJIIgiFdMQSvAoDrWkFBV9SQR7Y/0JLVcx/pfShDtTSYB7u0HZFWqjoToxjpagkQQRP3I/1SiexEXAScWAHuHAOKSlzaXl5ZClpmJqrt3Ic99ClpCJo4QBEF0FypNkiiKGkNRVBJFUSkURa2uZ/8ciqIePPt3naIoV1XESXQh2qYAhwfkRAM73YCS9AabKsrKIDpwEE/7eyN//EQ89RkMaXzCKwyWIAiCUCWVJUkURTEBbAPgD0AAYBZFUYIXmj0B4EvTdB8AXwP449VGSXRJxi6AQgaUpAIJfzfYTFFWjtIffqz+mhaLUbxyFeQFZIgGQRBEd6DKniQvACk0TT+maboKQBCASTUb0DR9nabp4mdf3gDQ8xXHSHRFzrMB9We1KA0cGmxGi0WAovYaNvK0tDrbCILoeJhMpoeDg4PA3t5eIBAIHC9cuKDVnOOnTp1q+bIZaVu2bNG3s7MT2NnZCWxtbZ0OHjyoCwBLly7tceLECW0AWL9+vVF9BXuJzkGVPzgzADXXpsh8tq0hbwE429BOiqLepSgqiqKoqPz8/IaaEQTgtRhY+VT5/NDYBgdxM7S1wTCuXUNTfexYUFrN+ltLEIQKcDgcRWJiYnxSUlL8119/nfXZZ5+16YfsR48eqW3atMn033//TUpOTo6PiopK8PT0FAHAli1bsidPnlwGADt37jQuLy8nSVInpcofXH1Fb+pdj4CiqKFQJkmfNnQymqb/oGnak6ZpT0NDw4aaEYQSiw3MeZZzX/2h3iYMAwMY/v0XOEOGgGFiAq3588D78nMwSJJEEJ2KUChk8ng8GaAs37Fo0aKetra2TnZ2doKAgAC959vnz5/fq3fv3k5+fn42BQUFLAA4efKk9siRI3s/P9fx48d1Ro0a1TsnJ0dNS0tLwePx5ADA4/EUDg4OVcB/vVDffPONUV5enpqvr69d//797QDg2LFjOn379nUQCASO/v7+1kKhkCRQHZgqfziZAMxrfN0TQPaLjSiK6gNgF4BJNE3XW2WZIFrEdgygawWErQPKcuvsphgMsKyswP99O4zOnobO2i/BJAk4QXQKEomE4eDgILCysnL6+OOPLdatW5cDAPv379eNiYnRSEhIiLt06VLy2rVre6alpakdOHBANyUlhZOUlBS3b9++tDt37nABYMKECWUpKSnq2dnZLADYs2eP/htvvFEwYMAAkYGBgdTc3Nzl9ddftwwMDOS9GMMXX3yRZ2RkJA0PD0++efNmck5ODuu7774zjYiISI6Pj09wd3cXff3118YvHkd0HKpMkm4DsKUoyoqiKDaAmQBO1WxAUVQvAMcAzKNpOlkFMRJd3bx/lI9/z2qwCYPHA9PICAyNOmWZCILooJ7fbnvy5Enc8ePHHy5cuNBKoVAgMjJSe/r06UUsFgvm5uay/v37l1+9elUzPDy8erulpaXU29u7DAAYDAamT59eGBAQwC8oKGDeuXOHO23aNCGLxUJERMTDwMDAR7a2tpWrV682X758eY/GYgoLC9N69OiRupeXl4ODg4MgKChIPz09/dXWxiGaRWVlSWiallEU9SGA8wCYAPbQNB1HUdR7z/b/DmAtAH0A2ymKAgDZiyvCEkSr6NsA7m8Dd3YBmbeAnl6qjoggiDY2YsSIiuLiYlZOTg6rsSoTz95n6li8eHHhuHHjbNTV1ekJEyYUPy8iy2AwMHToUNHQoUNF/v7+pW+//bbl5s2b69wReY6mafj4+JSGhIQ8aeVLIl4Rld4LpWn6DE3TdjRN96Zp+ttn235/liCBpum3aZrWo2m677N/JEEi2t7on5WPu/qTmWtEiwlFVcgrrURhmQRyOfk96kju3r2rrlAoYGxsLPP19S07evQoXyaTITs7m3Xr1i3u4MGDK3x9fcv++usvvkwmQ1pamtqNGze0nx9vaWkpNTY2lm7atMn0nXfeKQCA1NRUtatXr2o+bxMVFaVpZmZWZxaIlpaW/Pm4Iz8/v4qoqChubGwsBwDKysoYDx484LT/d4BoKVLgto0oqpT/Nxhs0nPa6XC4wOR9wIk3gHt7Afe3VB0R0cnkl1Zi3d8PcCe1GPpcNr6c7AI3Cz1w2ExVh9ZtPR+TBCh7cHbs2JHKYrEwb968kuvXr3MdHR2dKIqiv/rqq8xevXrJ5s2bV3Lp0iUde3t7Jysrq0ovL69axWJnzpxZuG3bNpaHh0clAFRVVVErV67s+fTpUzUOh0Pz+XxpQEBAndVpFyxYUODv729rZGQkvXnzZvLOnTtTZ86caV1VVUUBwLp167L69OlDlvLvoEiB21ZSVFZCnpGBsm07AADaHywG09wcDHX1V3J9oo3QNLBeDaDlwGohoK6j6oiITkIkkeH7U3G4EPvf4H81JoVjS4fAUKfz/B0gBW4bN3/+/F5ubm6iZcuWdbrvBylw23Jk6mErKXJzkTdyNMR//QXxX38hb+RoKHLqzpQiOjiKAt65qXx+7mPVxkJ0KqIqGaJTi2ptk8pp5JeRzoGuwsnJyTE+Pl7jvffeIzOsuxmSJLVSxaFAQCr9b4NUiooDB1QXUHd3dy+Qc69lx/bwAKxHAvf2AQVkMiUBVMmrUF5VjsZ63NXVmBCY1Z79zWRQ0OeSoSZdRVxcXEJUVFSShoZG17v1QjSKJEmtxNCpe1uGqmcb8Yrc2Ar84QFc/KzBlbQbNeWg8nH/8LaNi+h0CsQF2B27C9/d+gaX0y+iVFJabzuuuhpWjnWEtREXAKDJZuLL15yhrU6GfBJEZ0eSpFbSfH0qGHr/lfdh6OlBa8Z0FUbUzXm8DTDUgBtbgG0CQC59+TE1cY2AYd8CpZlAUmj7xEh0eMWVxfg0YiXOPjmN2IIY/HJ3Cy6k/wOZQlZvexNdDfy2wBPHlg5B8Ec+GOpoDE0OSZIIorMjSVIrMYyNYXTxAnQ3boDuhp9gdPGfOvW+iFfI6VmCKhMDJamov/rNSwxapXw8PAGQkXEl3ZFQUoJ8ce0akKcfh6KsqqyBIwA+l4Meehow1FEHR43MaiOIroAkSa1EMRhgmhhDa9ZMaM2eBaaJCSgG+baqjJYh0P9jwMxLOVPt7p7mn4Op9t9K3JHftW18RKfAZtYdT6TD1gaDIv+3CaI7If/jia5n1I/AW9cBrikQugjIvtP8c/QeCfBtgPD1QFlO28dIdGhcNS68TQdWf82gGHjL5V3wOHXKcxFdwPXr1zWCg4PJD5eogyRJRNfEYALvxyif/+EBiIoab1+fueeVj3/NaLu4iE5Bh6OD9/t+iJ+GbMJS9+X4Y+Qu2Oraqjosop1ERUVpnj59miRJRB0kSSK6Lk194J3byue/uza/5AjfGvB4D0iPBDJvtn18RIfG4/DgwHfAsF7DYaRpDHVW51kYsrMp33+An+Pm4ZLVs5dHjpuHS/n+A/zWnnPVqlWmVlZWTgMHDrSdMGGC1dq1a429vLzsIyIiNAEgJyeHZWZm5lJZWUl9//33PUJCQvQcHBwEAQEBei87N9F9kCSJ6NrMPIGx25Sz1c4tbf7xozYoH3cNIHXdWoimaZxPPYdSiVDVoRAdUPn+A3zh/76yUOTlsUHTUOTlsYX/+8qiNYlSRESEZkhIiF5MTEz86dOnHz148ECrobbq6ur0mjVrsidMmFCcmJgY/8477xS39LpE10OSJKLr83of6D0KuPUrkBTSvGM5XOC1Z4uD3t3d9rF1cTRNY1fMH9h271c8Ka1b+Ly8qhzFlUWQyCpVEB3REZT9vMUMEknt9yKJhFH28xazlp4zLCyM6+/vX8Llcmk9PT3FyJEjS1odKNEtkSSJ6B5mnlQ+Hp4IFD1q3rF95ijXXgp5F6gkvSHNcTjxEP5JOw81Bht6nP/uYtA0jZzybGy4/SOWhS3FrthdKK4kH+C7I0V+fr1VwRva3hQNrZDOYrFouVwOABCJRC1YH4TobkiSRHQPaurA0lTl8602QJWo6cfWrOt25qM2D62rKpGUICjpMCRyCWhaAT11fq19n11djbv5d1BUWYjzqWexJ2YXRNJm/FyILoFhaFjv0vgNbW8KPz+/8vPnz/NEIhElFAoZFy9e1AUAc3Nzya1bt7QA4NChQ9VZu46Ojry8vJy8HxJ1kF8KovvQtfhvxtrB0UAj9bjqMHUDeo8GHhwA8hPbJ74uhsfmQZOlCQAw1jIBV41bvU8sE6Gwsnat0Os51yCRk8U7uxvtZUuzwOHUHvDH4Si0ly3Nauk5fX19RWPGjBEKBAKnsWPH9u7Tp08Fj8eTr169+unu3bsN3dzcHAoKCqqXRPf39y9LTk7WIAO3iReRdfOJ7sVmFDDoE+DaT8C/PwMDlzf92CkHgA1GyrpuyzOVPUxEg9LL0iCSiTDVdhoWOL1Rax+byQGDYkBB//feaKJlAop8T7sd7vx5RYBybJIiP5/NMDSs0l62NOv59pZat25d7ubNm7PLysoY3t7e9p988slTNze3yuTk5PjnbbZu3ZoNAMbGxvLY2NiE1r0SoisiPUlE9zP8e4BvC/yzAsj4t+nHaRkCI34AyrKbPwC8G9pxbxsAYIrt1Dr7NFmaWCBYCOpZ2Rg2k4OP3JZCl6P7SmNsrRJJCU6mnMDii+/i7OMzqg6n0+LOn1dkejc6xiwzPdr0bnRMaxMkAJg7d66Fg4ODoE+fPo4TJkwo9vHxIfdyiWYjPUlE98NgAO/eBn7QBXYPBFbmAtwm1tvzXg5cXA0ETQK+qARYdctXEECBuADxRfEYYuYLbbZ2nf2aapoYZTEag8x8UFxZDAMNA+jU064j235vGy6lXwBN05DTcggMBKoOiaghJCSk7nRKgmgm0pNEdE/qPGDxA+XzX+0Bef3V3etgqgHzLiifh3/TPrF1AQfj9wNAndtsNWmxtWCkaQR7vj30NfShxmzxZKZXTqqQIurpbUgVUjAZTMywnwULHUtVh0UQRBsjSRLRfRm7AJP/BCRC4NRbTT+u9whA3x6I/AYozW6/+Dqp8qpyXM64BFtdOxhqGqk6nDZXKC7EjJDXUSDOhy5HF6ZapphuT0rXEERXRJIkonvrOx8QvA7c3w/EBjf9uLlnlY9/TWufuDqxk49OAAA+dOt6yyXEFsRg4fn5kNEyfNJvNfb7H8LWYdvAYpCRCwTRFZEkiSCmHAJY6sDRmUB+Eye46FkBnu8DGdebN/i7i5PKpQhOOgwdNg9WPGtVh9NmaJrG4YRAfHZ1NQBg+/Cd8DEbrOKoCIJobyRJIggWG1jybBXubQJAUta040b9pHzcPZDUdXvmUvpFAMCqfp+qOJK2I5KKsOTKhzicdAi9dW0QPP4oemr3VHVYxEtoamq6vazNN998Y2Rtbe00ceJEq9DQUO0LFy40WOOtPa5PdHwkSSIIANDpASyMUD7f49O0hSbZWspeKAC4E9B+sXUSClqB7fd/AwD0Meij4mjaRnppOmaenoa00lTMc5yPzb5boMHSUHVYRBvZvXu34ZkzZx6eOnXqyeXLl7UjIyO5Lz+K6E5IkkQQz1kMBoZ+DTx9AIR91bRjXGYpb9WFvgeIu3cNzVu5ytItKz0/6RKLQoZlXMGHlxcDAH4YvAHT7Gd0idfVER27nc4ftyHMZcC68x7jNoS5HLudzn/5UU335ZdfGjs7Ozva2dkJli1b1gMAZs+e3SszM5MzceJEm6+++spo//79hr///ruxg4OD4Ny5c7WSpeXLl/eYPHmy1YABA+wsLCycN23aZAAAQqGQ4e3tbScQCBzt7OwEBw8erLPQV2hoqHa/fv3sx44da21paen8/vvvm+3YsYPv4uLiaGdnJ4iLiyPriHRgZLQhQdQ05HMg8QQQ/hVg6QdY+TXenqKAt/4FdroBpz8AXj/0KqLskL6/+S0AYFAPHxVH0joyhQw/R29CZFYEOEx1/DEyoFbdOaJtHbudzt9yLsmiSqZgAEBhuYS95VySBQBM6der1YtKHjt2TCclJUX9wYMHCTRNY8SIETZnz57lBgYGpoeHh/PCw8OTTU1NZUKhkMXlcuXr169/Wt95EhISNKKjoxPKysqYbm5ugqlTpwrNzMykp0+fTuHz+YqcnBxW//79HWbPnl3CYNTuf0hMTNQ4evToYyMjI5mFhYULh8MpiImJSfj666+NNm3aZLRnz56M1r5Oon2QniSCqImigDfClM//HAoIM19+jGlfwHYsEBvY9IHfXUx8YTxo0HjT+W0wGUxVh9NiRZVFmBk6DZFZEfDrORSHxwWTBKmd7Q57bPY8QXquSqZg7A57bNYW5z937pxORESEjkAgEDg5OQkePXqknpiYqN7c8/j7+5dwuVza1NRU5u3tXRoZGamlUCiopUuX9rSzsxMMHTrULi8vj52ZmVmn88HFxaXCwsJCqqGhQffq1Uvi7+8vBABXV1dxenp651kgrBsiPUkE8SIOF/gwEfjNAdhiCXwuUg7ubszkP4ENhsC+ocDKnG5X1+3n6I0AAH9LfxVH0nJxBbFYc1U54Hyl5ycY0tNXxRF1D4Xlknr/czW0vblomsbSpUtzVq1aVdCa87x4q5WiKOzcuZNfWFjIiomJSeBwOLSZmZmLWCyu0/nA4XCqBzkyGAyoq6vTz5/L5fLu9ceikyE9SQRRHwN7YPpRgJYDf01/eXstA2DkT0DFUyDxZPvH14FklmXgqegpJvWeDA6r2R/QVY6maQQlBlYnSNuH/04SpFdIn8upas725vL39y89cOCAgVAoZADAkydP1LKysup0EGhra8vLysoa7AY9e/asrkgkonJzc5k3btzQ9vHxqRAKhUwDAwMph8OhQ0JCtLOzs0mvUBdDkiSCaIhgKuD2FqruM74AACAASURBVJB0Eriz5+XtByxTPga/Bkgr2ze2DmTngx0AgGl2TUgm25i4SobMIhEOXn2CS7G5KCqXNOt4kVSEpVc+QmDiIVjpWD+b3m/eTtES9XnLzzqLzWLUWkODzWIo3vKzzmqL80+ZMqV02rRpRf369XOws7MTvPbaa71LSkrqJENTp04tOX36tG59A7cBwM3NrWL48OG2/fv3d1y5cmWOpaWl9O233y66f/++lrOzs+PBgwf5VlZW3ec/fjdB0U2Z6tzJeHp60lFRUaoOg+gK5DJgkykgKgDeuweYuDbe/vFlYP9wwGcNMOK7VxOjChWKC7Hw/Hx4mw7Emv6fv/Lrx2WW4N3dtyBXKP+O2RhrY+t8D/C5L58wlFGWjg8uKWevzXGYi+n2M7v17DWKoqJpmvasue3+/fuprq6urbpN1RTHbqfzd4c9Nissl7D1uZyqt/yss9pi0HZbWb58eY/GBnV3dPfv3zdwdXW1VHUcnREZk0QQjWGygA/igQ1GwO99gU+LAA29httbDwMMBcDV7wGvDwCdNhl72mEdTlTO5lvo3IzadzXQNA05LW9RWY9SsRTbLz6sTpAAIOVpGbKLxS9NkiIyw7ExSrkY6Pc+P8LJwBkloio8zC3D1aQ8eFrrw6WnLnS1yN2TV2FKv15FHSkpIojnSJJEdFmlpaXQ0dFp/Ym0DIG3bwK7+gM73ZWrczMauVM9+zTwixUQPBV450brr99BiaQi/JN2HlY61jDRMmn28UXiIlxMv4Cs8kyMsfSHuXYvcNlNX8uPpmlIpPI62ytldbc9J1PI8MudnxGeGQY2k40/Ru4GX52Pyio5gv9Nw96IxwCA4BvpGO/WAx+PdoC2hlqzXxvRtWzevJlUsu6myJgkoksqKChAUFAQjh49itjEWORX5KO8qonlRurT0wsY8wtQkgr8s7LxtnqWgNdHQNZNIP1ay6/ZwYU8PgUA+Mh9SbOPLa4sxurIVTiYsB9XMi7j08hVeJB/v1nn4GmysWBw7fpwhtocWBnUn2gVVxZh1ukZCM8MwxAzXwSN+wv8Z9P7yyqlCLyeWqv9mXvZEFc1nHARBNH1kSSJ6JL09fVhZGSEoqIiXIu8hqDDQfjx5g/IqchBi8fh9f8IsBwK3PgZeHi28bYjflA+7vEBFF3vjVaqkOJQwgFosjRho2vb7OPzRE+RK8qtte1IcjCEEmGzztPXUg87FvbDcCdjzB9shd3vDoC+dt1bbfGFcVhwbh4k8kqs8FiFlf0+qXWLj6IAxQu/F1qV5ZCfP4uij5dCcudOs+IiCKJrIEkS0SVRFAXPgZ5QUAqABjgKDpLyErHu2hcokbSwfAhFAXNOK58fGgsUP2m4LVsTmHpY+Tz6j5ZdrwMLzwgDAHzqtaZFxzOpujOt1RgsUGjewGltdTW4WfKxbooLFg2zgZFO7SUIaJrGkaRgrI78BADw27Ad8DX3q3MeLQ4L0/r3en4Qvj31PXYHroTi05UQHzsOyEkBY4LojkiSRHRZGtoaSNZMRC47BwAwonA0REIxJPJWzNJV0wA+Vo5bwS/WgFTccFvnGYCaFnD6/S5V101BK7D17hYAQF/DlhU6N9AwhKWOVfXXFCjMFcyHDqdlY8jYLCaYL4wTE8vEWB6+FAcT9sNCxxJB4/5CL51e9R6vwWZh/mBrbJjthjkmcjjkPQJTLgPFZIL33Tfg9POs97jWSk1NRXh4OE6cOIGDBw8iPDy8Xa5DEETLqDRJoihqDEVRSRRFpVAUtbqe/Q4URf1LUZSEoqiXDAQhiNoYFAOlhkJE695GGP8yAGBIsR/yMvNbd2I9K+XgbAA4NA5o6PYdRQFvPRuTdHpx667ZgUQ/VS6vscx9RYunzOuq62L9wK+x3GMlZtjPwq/DtsNW167NYswsy8CM0NfxqCQFsx3mYOvQ36CpptnoMTwW4HJ4B6Z8+x4AgOJqQX3CeHDnzWuzuGqSyWS4du0akpKSkJeXB4lEAnt7+3a5FgH4+vraFBQUNFozx8vLyz4iIqLOL8r169c1goODeW0Vi5mZmUtOTg6ZONUJqCxJoiiKCWAbAH8AAgCzKIoSvNCsCMASABtfcXhEF8Dj8LDMfQU4TA7KWWW4ZPgPAOBG+A3cuXOn5WOTAMBurHLxyNQrwM1fG25n4grYTQBig4C8+JZfrwP54ZZy/afBPYe06jy66nrwMx+KOY5z0Uun10uTmKa6mhmJ9y8pE53vfH7ETIfZL03mqu7cRbZVb1QE7ALb0xMmMQ9gGh8HvU1t/6dHIpHg1q1b2LNnDyoqKgAAbDYbkyZNgolJ82cJEo1TKBSQy+UIDw9PMTAwaNEAwaioKM3Tp0+3WZJEdB6q7EnyApBC0/RjmqarAAQBmFSzAU3TeTRN3wYgVUWAXYVQIkRCYTzOp55Fdnk2xLJGbhF1IQyKATs9O+wcuQu/DP0N28fsxKz5s6Cvr4+oqChcuHABCkUrxpqM2gjoWgHnPgYybzXcbvJe5eM+34Z7nTqJpKJESBVSLBAsbNHaRm1NWmOskFwhx8/Rm/BT1A9QY6hh7+j9cDZwbvR4RUUFit5bjPwJEwEA+oEHYXjyOJh8PVBMJqjGlnpoJpFIhLCwMPz555+4d+8euFwuJk2ahIEDB2Ly5MkwMDBos2t1NmefnOEvODvXZeKJcR4Lzs51OfvkTKuqCiclJbGtra2d5s6d2+tZUVt2zd6bVatWmVpZWTkNHDjQdsKECVZr1641fn7s4cOH9VxcXBwtLS2dz507x62srKS+//77HiEhIXoODg6CgICAWgulbd26VX/48OG9Bw8ebGtpaem8YsUK0+f7RowY0dvJycnRxsbGaePGjXV+wElJSWwrKyunGTNmWNja2jpNnDjR6sSJE9ru7u4OFhYWzleuXGmbTw5Ei6nyr5wZgIwaX2cC6K+iWLqssqpSBDzYiYgs5VgHChS+GLAWHsaeYFBdf0iaGpMNPpNfPdUbAKZMmYLIyEgkJibi8OHDmDZtGtjsFiwayGAAi6KBH/nKNZRW5SnXVHqRpj4wahPwzwog8QTg+ForXpFq/XJHORZprPU4lcZRUlGF208KER6fh369+XC1YuPT6x9ALBPDp4cPlnuuemkSVxkWhsI5yltpGhMnQHfjBjC0tNo81tLSUly/fh3p6ekAACMjIwwZMgR8vvJ30tjYuLHDu7yzT87wd8X8YSFVSBkAUCwpZu+K+cMCAPytxrZ4gcnU1FT1gICA1IMHD6bX3B4REaEZEhKiFxMTEy+VSqm+ffsK3NzcRM/3y2QyKiYmJiE4OJi3fv36HmPGjEles2ZNdlRUlNb+/fvT614JePDggVZMTEwcl8tVuLm5CSZNmiQcMmSI6NChQ6nGxsby8vJyys3NTTB37txiExOTWr1ZGRkZ6sHBwY89PDzS+vTp43jo0CH9qKioxMDAQN1vv/3WdOjQoY9a+j0gWk+V75L19X+3+GM2RVHvUhQVRVFUVH5+K8ecdCFimbg6QQIAGjR2xfyB0mZOte5KKIrCkCFD4OPjg4qKCuzbtw/l5eUtO5mGnrJcCQD85qgsY1Kf/s/WEgqe0mnrumWXZyGzPAPjrMZDg6WhsjjEVTLsiXiEL/96gItxudhw6RKWRLwJsUyMZR4r8InXmkYTJHlRMfInTq5OkAxDToG/YzsYWlqgaRqyzEyIToVA9riR2YtNUFRUhGPHjiEoKAjp6emwsLDArFmzMHny5OoEiQCCEgPNnidIz0kVUkZQYmCrlqs3NTWtGj58eMWL28PCwrj+/v4lXC6X1tPTU4wcObLWrIpp06YVA8DAgQMrMjMzm/TpycfHp9TExETO5XLpcePGFYeFhXEB4McffzS2t7cXeHh4OObm5qrFxcXVqQBtZmYm8fLyEjOZTNjZ2YmHDRtWymAw4O7uLsrMzHx5fR2iXakyScoEULOSZE8ALV7VlKbpP2ia9qRp2tPQsJ5P892UVF73TmWFtAJ0y/PRLkMgEGDcOGWPSGBgIFqcXJu4AhN3A+JCIPS9+tswWcCCK8rnYetadh0VC4hRLmUww2GWSuOokMhx/HYGABqaRrfBt/sLAPA/r58x1HxYg8fRNI2KY8eR69IHVdHR4L77Dno8eQS2u3KGXtHSZchxEODpoMEoXvw+Ki9caFF8ubm5OHToEI4ePYqCggI4ODhg7ty5GD16NLS1tVt0zq6sWFJcbyLS0Pam0tTUrPde+svGIqqrq9MAwGKxIJfLmzQz4cUxbxRFITQ0VDs8PFw7KioqMSkpKd7R0VEsFovrvOey2ezqgBgMRvX1mUxmk69PtB9VJkm3AdhSFGVFURQbwEwAp1QYT5ekpaaFntzaVc3HWo0HV438sQYAMzMzTJ+urF5//PhxPH78uGUncn8TsJ8E3N0NxP9dfxsrP8CoD3DtJ0CY2bLrqEhxZTGin0bB07gfdDm6Ko6GBktNBj27YHDN/oVMzEfhg0Uw0zJv8AhZVjZy+g9AyUdLQDMo0KGBUHz6IWJLE3Ey5QTS429A/NdR0OXlgBoLGq9Nhta77zQ9IppGeno6AgICcOrUKVRUVMDNzQ1vvPEGhgwZAk1NMrSkIXocvarmbG8tPz+/8vPnz/NEIhElFAoZFy9efOkvtI6Ojry8vLzB98urV6/qPH36lFleXk6dOXNG19fXt7ykpITJ4/Hk2trairt376rfv3+/7e/lEu1OZUkSTdMyAB8COA8gAcARmqbjKIp6j6Ko9wCAoigTiqIyASwH8AVFUZkURbVBMa7uQ1ddD18P+hZTbKair2FfLHFbivHWE6DGJPWontPV1cW8efNAURQuXryI6Ojols18m3YEoJjAkdeBguT628x+9jkguHONSwpOCgIAvO3S9MShvZTJ8qDtuA1qWnkoz+mHosQ5mNC3NzQ5dW+x0QoFynb+gade/UFnZePoa6b44Bdn3OOVIOBBANZf/ASiz9aBOXJa9TFsNzfobfm5ScsbKBQKpKSkICAgAOfOnQNN0/D29sabb76Jfv36tWysWzcz02F2lhpDrVavjxpDTTHTYXZWe1zP19dXNGbMGKFAIHAaO3Zs7z59+lTweLxGZ735+/uXJScna9Q3cBsAPD09y2fMmGHl7OzsNGHChOIhQ4aIpk6dKpTJZJSdnZ3gs88+6+Hq6lrn1h/R8VGtmgbdQXl6etJRUVGqDqNDkSqkqJJVQVNNs8Vr23R1MpkMp06dQkFBAXr16oVRo0aB0dzZTcJM4OdnPRprygBOPXXEzi0DbmwBFkYCFj6tD7ydiWVizAh9HT255tg+4neVxnI96xp+uK1cguBN28/w9KkR3C34sDXRhq5W7YREmpKCvBEjAakMtIkRPv+Aj2K+ss0Gn4049tNCzD+kfB9OstXCiff7Yp3nl9AztQZDo/ExV3K5HImJibh27b/afH5+frCxsWn+70wHQVFUNE3TtVbNvH//fqqrq2tBe1/77JMz/KDEQLNiSTFbj6NXNdNhdlZrBm2/jFAoZPB4PEVZWRnD29vb/vfff0/z8fERvfzIurZu3arf2KDujuD+/fsGrq6ulqqOozNS/Rxe4pVQY6hBjU16jxrDYrHw2muv4erVq0hISEBgYCCmT5/evN4AXk/l2KM/hwL7/IB3bysXlaxp2LfKJGnvYGCtDGA0ur6dyp15rFw482P3pSqLQa6Q47d7W3Ep/SKYFBO7Ru2BvoYB4FS3LV1VhdIffkT5TuUYKt2tW5A3vC+Kw5SD560eV4D7wSjMB6CggB8/sUF6L01oqdFg9DAFQ73hBKmqqgqxsbGo+SFs9OjR6NWrF/nw0Qr+VmOL2jMpetHcuXMtHj58qCGRSKiZM2cWtjRBIro+kiQRRA0URWHw4MEwMDBAZGQk9u3bh9mzZ4PLrb+yfL2s/ADftUD4eiDiW8D3i9r72ZrA68HA0RnA7R1A/w/b9DW0JZlChj/j94LNYMOe76CSGIori7H44iKIZBUY2GMQVniughqj/oS/6u5d5I9Xrnmk5u4G/T/3gcnnQ6+yGPYSPib9chNWacp1wkLe98QZgaQ6iZ1gPQFaavX/nCsrKxEdHY24uDgAgLq6OkaOHAkTExOSHHVCISEhrZu6WMOSJUsKARS21fmIjoUkSQRRD0dHR/B4PISGhiIwMBCTJ0+GkZFR00/g9z/lmkhXvgQsfQGLwbX3O00DTr0NnP0I6DNHuZRABxSZGQEAWO31mUqun1iUgE8ilBWJlrovx7Bew+ttpxCJULLqE4hPnAQA8A/sh8awoQAAWiwG9dOvWPpHGAAgfrwLTD7/H6YZ9oZu2nkkFyfDt6cf3IzcwWbW7jUsLy/HjRs3qgf06+npYejQod164UeC6E7ImCSCaERJSQmOHDkCABgxYgSsra2bfrCkDPj+2TyDFdmAtmnt/U9jgB19AMHrwPS/2ijitkPTNCadHA8AODkp9JX2mNA0jWMpf+PPOOVq5b8O2wYLHct621ZGRKBw1hwAgPq4sdD7eXP1mkfik6dQ/IGyp07NxQV6u/8AbWJUnQzJFXJUKarqrPtUUlKCyMhI5OQoiyP36NEDgwcPBo/XdStTqHJMEtG+yJikliM9SQTRCF1dXcyfPx+HDh3CxYsX4e7uDg8Pj6YlDBxt4IN4YJsA+MUaWFMK1JxVaOwCOEwG4o8CT2MB48ZLaLxqd/PuAACWuC19pQlSpawSX177DEnFSejJ7YkNvpuhpVZ39rSiuBiFb76NqlvKkjCGp06C7eEOAKiKiUH+hEmAVLlOmMGJ4+D086xzDiaDCQ3GfwlSfn4+Ll++DKFQudiqjY0N+vfvD612WImbIIiOr3NOwyCIV0hdXR1vvPEGjIyMcOfOHZw9e7bpNd8MHYGphwFZJfD3nLr7J+5WPu4d0uHquv10+wcAgK+53yu7ZnZ5FqaHTkVScRKm283EtuG/10mQaJqG6MRJ5Dj3QdWtW9B6+y3lopAe7pDn5yN/+kzkjxkLSKXQ3fATemSk1Zsg1TxfVlYW9u7di+PHj0MoFMLFxQXz58/HsGHDSIJEEN0Y6UkiiCZgMpmYNGkSrl+/jri4OBw6dAgzZsxo2sw3l5lAylng/n7g3n6g7/z/9mnygdE/A+eXAQnHAMHU9nsRzZBS8hAimQhzHOc1OEi6rdWc3v/1oO/gauhap408Owf5U6ZCnqEs+2h0+SLU7O1BSyQo3bARZVt+AQBozp4F3rq1YDQy4J6maaSlpeGff/6p3ubl5QUnJyeoqZGZoJ1ZQUEBc9euXfzVq1eTGlVEq5CeJIJoIoqiMGjQIAwZMgRisRj79u1DWVlZ0w6euBvg8IATC5RjkWrq/5Hy8cjrgFTctkG30K93lMnGBOuJ7X4tuUKOrXe24Ifb34EChT2j/6yTINEKBcoCdiG3nxfkGRnQ+WwNeqSngmVnB9G588i2tkHZll/AsraG8fWr0NvwU4MJkkKhQFJSEgICAqoTpMGDB+Ott95C3759SYLUBRQWFjJ3797djJkWBFE/0pNEEM3k4OAAHo+HkJAQHD58uGkz35gs4KMkYKOJcrD26hJA/dkgYAYTeCNMua7SlbXAqA3t/RIalVuRgyelTzDaYgw01dq3nEaJpATvX1yEcmk5Bph6Y1W/T+v0XElTHiFv5CigqgoMY2MYnjoBVs+ekCYmIn/K66CfjR8yOBIMzqCBDV5LJpMhLi4ON2/erN42fPhwWFlZddoFILuKuLPJ/DvBMWaiYjFbU0+jyn2GS5aTv12L101asWJFz4yMDI6Dg4PA19e3FAAuX77MoyiKXrVqVc4777xTHBoaqr1+/foefD5fmpSUpOHi4iI6ceLEEwaDgeDgYN7q1at78vl8mYuLiygtLY1z5cqVlOXLl/fgcrny9evXPwUAW1tbp9DQ0If29vZV27dv5+/YscNYKpVS7u7uFfv3709jschbbGdH/jIQRAuYmppixowZAIATJ04gJSXl5QdxjYG3riuf/9EPqDmuydIXMHYFrm8EhBntEHHT7YlVjpOa5VjPGKo2lFSUiPln56BcWo6P3Zbhs/5f1EqQaKkUwm++RZ6vH1BVBd2fN8Mk+jYoTS0ULnwLecNHghYKwVv/FXqkpzaYIEkkEty8eRN79uzBzZs3wWQyMXbsWLzzzjvo3bs3SZBULO5sMv/f3dEWomIxGwBExWL2v7ujLeLOJvNbes5NmzZlmpubSxITE+O9vb3LY2JiNBISEuIuXbqUvHbt2p5paWlqAJCQkKCxbdu2jJSUlLj09HTOhQsXuCKRiPr4448tzp49+zA6OjqpsLDwpZnOnTt31I8ePcqPiopKTExMjGcwGPTvv/+u39L4iY6DpLkE0UI8Hg/z589HYGAgLl++jKKiIvTr16/xmWDm3sCojcA/K4FLa4CRP/63b9YpYIsFEDQZWBTd/i+gHkKJEDdy/oWroSv46i1+j2oUTdM4kXIce+OUydjWodtgybOs1abq/n3kj1UuP6DW1xX6+/8EQ0cH5dt3oPS77wEAGhMnQPf778DQrb8+qUgkwq1bt5CcrKyjx+VyMWLEiOatd0W0uzvBMWZyqbxWpiqXyhl3gmPMWtOb9FxkZKT29OnTi1gsFszNzWX9+/cvv3r1qiaPx1O4uLhU9O7dWwoATk5OokePHrG1tbXl5ubmEgcHhyoAmDlzZtGuXbsMG7vGuXPntGNjYzVdXV0dAaCyspJhZGQka23shOqRJIkgWkFdXR0LFixAaGgo7t27h/z8fPj7+zfeO+G9XLnQ5LWfAKthgM1o5XbdXsp9/24GUiMAyyGv5kXUcDRZuSbUoj6L2+X8Elkl1l7/AglFCeih1QOb/LbUmr2mEIlQ8ulqiI8dBwDw/9wHjRHDURkWhsI58wAADBNjGAQHQc3Gpt5rlJaW4tq1a8h4PrjbyAhDhgwBn98+SR/ROs97kJq6vbkaWwuQw+FU72QymZDJZI2uHchiseiaM1slEgn17BrUtGnTCrdt29YuRXkJ1SH9zATRSkwmExMnToSzszOysrJw4MABSCSShg+gKGDueeXzg2OAkrT/9g37Rvm4zxdQNFqYvM1JZJU4+egETDRN0FPbvM3Pn12ejWmhU5FQlICpttOwfcTOWglSZUQkcmztIT52HOr+Y2CalACWlRVy+vWvTpD0D+yHaXRUvQlSYWEh/v77bwQFBSEjIwM9zXthzKSpGDN2AkmQOjBNPY2q5mxvCh6PJ6+oqGAAgK+vb9nRo0f5MpkM2dnZrFu3bnEHDx5c0dCxrq6ulRkZGZykpCQ2AAQHB1f/8lhaWkru3bunBQBXr17VzMrK4gDAmDFjSkNDQ/WysrJYAPD06VNmcnJymyR5hGqRJIkg2gBFURg4cCD8/PwgkUjw559/orS0tOED2JrAkmfjmLZYAtJK5XM1DWDas9W3b/3WrjG/6GzqWQDAMo8VbX7uG9n/4r2L7wAAvh70LRY4vQEGpfzzoygpQf7r01A4azYA5cKPeps2omTNZ8gb4gtFdja0V3+KHqmPof6s1EhNubm5OHjwIP7++28UFhbC0dERYyZNw40KUywPiscv55NQWNZI0kqolPsMlyymGrPWwmNMNabCfYZLi3tlTExM5B4eHuW2trZO169f5zo5OYkdHR2d/Pz87L766qvMXr16NXgrjMvl0ps3b04bM2aMrYeHh72RkZFUW1tbDgDz588vLi4uZjo4OAh+++03QwsLi0oA8PDwqPziiy+yhg8fbmdnZycYNmyYXUZGBpkm2QWQsiQE0cZyc3Nx6tQpAMCkSZNgbGzccOOkEODwRKD3KGDes94lmgZ+5AOVJcAnhcq1lNqZXCHHa6cmggKFk5ND2+68tBy/39uO82nnAAC7R+2DoaZyeAdN0xCHhKB48QcAAK03F0JnzWqIgo9A+MWXAADOiOHQ+3kzmC/0BNE0jYyMDJw/f776doqbmxtcXV1RUUXjo/3RSHn63/IMA+0M8L8pLtDRIB/uG6LKsiRtPbuttYRCIYPH4ykUCgXmz5/fy9bWtnLdunV5qoqntUhZkpYjY5IIoo2ZmJhgxowZCA4OxsmTJzFs2DDYNDB+BvYTgH4fALe3Abe2A17vK2/HLYwEdrgAIe8AM/5u95ivZkUCANb0/7zNzimUCPHBpfdQWlUKL5P++LTfGqg9K8siz8lB/rTpkD9JBQAYXboARXExcmztAQAUlwvD48egJnCsdU6FQoHHjx/j8uXL1du8vb3h6OiI59OtCypEtRIkAPj3YQEkUgVQu0Rbs1TKKlEuLYcaQw0sBgtsBrv69RCt4+RvV6TKpOhFW7ZsMTh8+LCBVCqlnJycRMuXLyf167opkiQRRDt4PvPt8OHDuHz5MgoLC+Hl5VX/zDf/X5QDuc98AJh5AWaeyjpujq8DCUeVi08au7RbrDRNY1O0cm0mL5P+bXLO5OIkrAxfDgD4yO1jjLQYpbyWQoGKvfsgXLsOAKC9+lNoTBiPooVvQpb8EACgF/AHNPzH1PpeyeVyJCYm4tq1a9Xbhg4dWu8UfjUmAxwWAxLZf3dw9LTYYLSy/txPt3/Avby7oCgKUoUU1rze2DJ0a6vOSXRM69aty+vMPUdE2yFJEkG0E3V1dcyfPx9nzpzB/fv3q2e+MZnM2g0ZTGDxA+AnfSCgH7AqH9AyACb+oUyS9gwGVhcre5jawYOC+wCA910/rB4n1FI0TePUo5PYHRsAAPhl6G+w4lkBAKSPHiN/1GjQlZVgGBpAPzgIFXv2Im/QYAAAd8lH0Fn6MSgOp/p8VVVViImJQXT0f0sijB49Gr169WpwqQWuuhqWjLHHhtAEAACTQWHNRCfoajav10eukCOpOBFhGVdw7tl4LQoMsCk1mGqZ4mP3pc06XyelUCgUFIPB6HrjMroJhUJBAWhisUniRSRJIoh2xGQyMX78eNy8eRMPHjzAgQMHMGvWcRyx9QAAIABJREFULHBqJAIAlOOO3o0G/vBQ3mZbnglo6AH+W4GzS4D4o4DTtHaJccNt5VpNw3uNaNV5JLJK/O/ftYgrjIOJlgk2+/4CLpsLWipF6cZNKP9tGwCAt3EDQAH5w5TXY3sPAH/bb2DWGLtVWVmJ6OhoxMXFAVAmnCNHjoSJiUnj61AB0GAzMdrFFN42BsguEcOcrwltdTUwmS9PAIUSIaKfRuFo8l/ILP9vUU82g41JNpNx6tEpjLcej9mOc19ZTTsVi83PzxcYGhoKSaLU+SgUCio/P58HIFbVsXRWZOB2IxSlZaArxQBFgWnY6FpiBPFSDx8+xJUrVwAAM2fOhI6OTt1GUX8AoYsAz8XA+O3KZQDWP/ss87lIOfutDT0RPsbHVz7CDPtZmOM4t8XnyanIwaILbwMApti+jvmCBWBQDFQ9eIB8/3EAADUXF2h/shJF8xYoD2IwYBh6CmzX/+q0lZeX499//8WTJ08AAHw+H35+fjAwMGhxbI1R0Ao8KnmEyKwInEg5VmtfX0M3jLMejz6GrtBgaVS3b21vW0dV38Dt6OhoIxaLtQuAM8hs6M5IASBWJpO97eHhQW4ftgBJkhogz8tDyedfovL8ebAsLaD782awnZ1r3QogiOaqOfNt4sSJMDExqdvo0Djg4RlgxnHAcTKQFgnsHQIMWAaM2dym8awIW4aHJckIHBsMLrv+grAvczPnBr69+TUA4Cvvr+Fm7A6FWIySNZ9B/NdRAIDu5o2oOHAQ0rv3lF//uhWar02u7hUqKSlBREQEcnNzAQBmZmbw8fEBj8dr7Uuso7yqHPfy7+JEynEkFyfV2jfLfja8ewyChY7FS3usupr6kiSC6O5IklQPhUiEks+/hPjIkeptlIYGjK9F1rolQBAtUVpaiqCgIADKwce2tra1G8gkwDfqyucfJQP6tsBOTyAnGliaplyZuw3kif7P3pmHNXGtf/w7SUhCSNj3RUAghLCJbK5VRK3UjYpUQUWvWreftS6ovcXaaqu3VbHW26q01brVpV6tC4qIioK7oIIsYRNQ2RGBsAWSzO+PEAQFBIuiOJ/n4YGcOefMe4Yk88573qUYs8/9C8PMvLHYdWmnx8tIGcIStjf56+wc+Qf0OPqQXLmK0kmTAQAsb28wepmh+o/dAAC12bOg/sVK0FQVlpmSkhJcvHgRFY1Fam1sbODh4QE1NbUXT/iKkCSJ3MpcXM+/ioNpB1oc42vx4Ws9AX30XF5ZSewpUEoSBcWLUD5JrUBWVUHSLMQYAMjaWsiKiiglieIfo66ujunTp+PgwYOIjo5GaWkp+vXr98xywWAplKEt5sB/+cCXVcDkv4EfeylyKs2/1yVy7En+AwAwVRjU6bGVkgosvLgA5ZJyuBm4498eIaCLq1EyYzLqGyPQuMuWoip0MyRQbLVp7/wdDBNjkCSJvLw8nDt3Dg0NDQAAR0dHuLi4gM1md8naaqW1SCxJwOkH4bhXcrfFMV/rCRhs8gGsNK167NZZV0Ftt1H0cF66HUkpSa1AMJlQsbODJDb2WSONBpoOVdSZomtgsViYNm0aIiIicP/+fZSUlGD06NHPIt80ewHTooB9I4A9w4HZ14ABy4FrG4HsS4Dl0H90/sr6SsTmxUCoYw9d1c75+2Q8Tceyy0sAAP/XZyE+tPBBzalwFM9T1HtjjxmNushzqApVbA3qHj8Glrs7SJJEdnY2oqKimuby8PCAvb09VFT+uRN0flUebhTcwIHU/aiXP6toYcYzg5+NP/oauEKT1XoxXIrWYTAYvxsaGtrp6ek9pRy3KXoajY7twsLCwt8BjGutD6UktQJNUxOa369Hqf8nkOUXACoq0Ph6NWitOdpSULwidDodo0ePxq1bt5CQkIC9e/ciICDgmTXFajgw+Esgdr2iGK7XWoWStMcL+KoBoL/6x/fvDEWCyvnOCzo8hiRJhD84hd/uhwEAtgzdCvM6NRQNHgLpgwcAABWXPqgLPw0A0Nj4A9QmTwYJQCQSISYmpmmuwYMHg8/nv5gOoRPUy+qR/CQJkTlncS3/aotjH5qPglevYbDVEoBOe/VzUMCBUpAoeio0Go3U09OrKCwsdGirD6UktQHd3Bx6p8NBVleDYKuCUOeB1oV+EhQUgKLmm6enJ3R0dHDx4kXs3bsXkyZNeuawPOw7IPVv4PwXQK9BwCdHgb/8gFv/BfoveaVzSmQSHM34H3TYOjBXt+jwmLXXv8b90vsw4Bhg8wc/gjj4Nwoby4cwB/RH/bXraLh7D5yAydD45mvI2Wwk3r+PmzdvNs0zfPhwWFhYvJAAsqOU1BTjduFtHBT9iYr6iqZ2HbYOJtkGwM3QvdOWMYp2oVEKEkVPpvH93eYXEqUktQFBEKDr63e3GBTvCdbW1uDxeDhx4gQOHz78LPKNIIDZ14HvNYFdg4Bl+YCqDhC5FHAOAjid3wKOyj0HAFjmtqJD/QurCzCnMbx/vJUvpql5o9TRA2RNDaCmBlRXo/7addAtLaH75z7IDA1x6+5dJCYmAgAYDAZGjhwJExOTTkeMSeVSiMpScfHhBZx/GNXi2AemQzDCfCSE2vZUeRAKCorXAuWIR0HxlmBgYIDJkxVRYSdPnkRaWmN4OlsDWNCYC26rDRB0QfH3iVmdPoeMlOHXxB0AAHsd+5f2v11wq0lB+tp9NSaeeYqSwUNA1tSAUFUFqqsBALqHD0H93Flcyc7Gnj17kJiYCB6PB19fX8ycOROmpqYdVpCe1pXhfG4U5kTNxoST4/HllS9w/mEUOAwOZjl8irARv+Ok72kEu62As14fSkF6T/juu+/0e/fubT9u3DjLzo5NS0tj7tix47VVil66dKnx6tWruyWqJy0tjWljY9Pqhzk3N1fFy8urjcKR3cP69ev1fvrpp3fGwZeyJFFQvEUoI98OHz6My5cvo7S0FAMGDAChbw9M2A8cm6rwT7KfBCQfBgoTAUOnDs9/Pf8aAGCl+7/bVVpkpAy/J/6K09nhAIDfTFZC6jERVQAIHhekuApkbS3U13wD+YSPceHGDTzevx+AQtn74IMPoKWl1SGZZKQMmU8zEPM4BqcenGhxzN3AHaMsP4KTrhNYjK6JfKN4N9m5c6deREREhkAgqH9575ZkZGSwDh8+rD1v3rxOFdGVSqVNhZPfRdavX28wa9asDhfnfX69crkcJEn+I9/B5/nss8+eeHh4CD7//PMnXTbpa4SyJFFQvGWwWCxMnToVJiYmSE5OxsmTJyGTyQCnKYBDAJB0AOjdWEJk10Cgg7nOSJLEhtvfAwD6Gfdvs19lfSVmRc7A6exwuGs4Y9clG0h9pzybR1wF1XFjwboSg3Namjj8v//h8ePHsLCwQEBAAMaPH/9SBUlcL0bM48tYHL0IH58Yh+Uxy3DqwQkQIDDVbhp+HrYdJ8aH46v+38Dd0INSkN5zAgMDez1+/Jg1btw46zVr1uhXVlbS/P39LRwcHOzs7OyE+/fv1wQUVhVXV1dboVBoJxQK7aKiotQAICQkxCQuLo4rEAiEa9as0d+6datOUFBQU8IxLy8v6/DwcB4AcDgcl8WLFxs7OTkJLly4wI2NjeW4u7vb2tvb2w0aNMgmNze3XdNlcnIya/DgwTb29vZ2rq6utnfv3mUr252dnQUODg52ixcvNuZwOC7KMV999ZWBg4ODHZ/PFy5ZssRYuZbevXvbT5482dza2tp+4MCBNlVVVQQAxMbGcmxtbYV9+vQRbN68uU2/kNOnT2v5+flVtHdtwsPDeZ6envyxY8da2tra2ivPO3Xq1F729vbCrKwsZlhYmDafzxfa2NjYz58/3wQAfv/9d63Zs2ebAsC3336rb2pq6qhcp6urqy0ALFiwwMTKysqez+cL58yZYwoAPB5PbmpqKomOjuZ05j3QXby7KjIFRQ+GRqPho48+QlxcHO7evYs9e/YgMDAQ7I/3AGkngVOfAoNDgNh1QPJfgMOkl86Z/ESxZTfHaR7oROtPhhlPM7DssqJwa7BsOKymKvIcNcllYAD88jNOZ2agtjGU387ODm5ublBVbbtkCkmSyK54gGv5V/FX+uEWx4Q69hjXezz66LuAo/JOfG9SvGEOHDjw8PLlyxqXL19ONzIyki5cuNDEy8ur8siRIzmlpaV0Nzc3u3HjxlUaGxtLY2Nj0zkcDnn//n1WQEBA76SkpNR169blhYaGGkRHR2cCwNatW9vc7qmtraU5ODjUbtmyJV8ikRD9+vWzPX36dKaxsbH0t99+0woODjY5cuRITlvjZ8+ebf7rr7/mOjo6Si5evKg2f/78Xjdu3EhfuHCh2YIFC4rnzp1btmHDhqY6V8eOHVPPzMxkJyYmppIkieHDh1tHRERwe/fuXf/w4UP2/v37HwwYMCD3o48+6r13716tBQsWlM2aNcvixx9/fDh69OiquXPnmrYmh0gkYmpoaEhVVVVJAGjr2gBAYmKi2t27d5MFAkF9WloaMycnh/3bb7/l7N+//2FOTo7KN998YxIfH5+qp6cnHTx4MH/fvn2aI0eOFG/ZssUQAK5evcrV1NSUZmdnq1y8eJHbr1+/qqKiIvqZM2e0Hjx4kESj0VBaWtr0pdO3b9/qS5cu8by8vGo6+VZ441BKEgXFWwpBEHB3d4e2tjYuXLiAvXv34pNPPoHmokwg1EihIAHA/yYD/LEAs30FIzRuIwBghPnIVo+fzjqFsPs7oFojw8ZjbNCvtyyBIgvdiNOVlcB9hUN237594eTkBCaT2ep8NQ01SCi5h/AHJ3G/9H6LYxP5/hhoPBiWGpZUQkeKTnPp0iX1yMhIza1btxoCgEQiITIzM5nm5uYNs2bNMk9JSVGl0WjIzc3tdB0pOp2OGTNmPAWAxMREVkZGhuqwYcP4gGL7SU9Pr6GtsRUVFbS7d+9y/f39rZRt9fX1BADcvXuXe+7cuUwAmD179pNvvvnGFADOnj2rHhMToy4UCoUAUFNTQxOJROzevXvXm5iYSAYMGFALAC4uLjU5OTmsJ0+e0MViMX306NFVADBz5swnFy9efKF+z6NHj1S0tbWlzeVo69o4OTlVN9/GNDIyqvf29q4GgCtXrqj169dPbGxsLAWASZMmlV2+fJk7bdq08pqaGtrTp09p+fn5TH9//yfnzp3jXblyhTthwoRybW1tGYvFkk+ePNl89OjRFZMmTWoKR9XX15eKRKJ3wjxMKUkUFG85VlZW4PF4OH78OP766y+MGTMGxv+KBf4YDGiYAxW5wIUvAZ8tbc6RW5mDJ3VP4GfjDxa95X1DIpPg2+trkFiagA+S5QjYltzieP2ns3HeQB9kZSUAYMCAARAIBC/4apAkibyqx7hecB0HUvdDRsqajlmqW2KCzUS46LtAndX19dgo3i9IksT//ve/TGdn5+aGTixdutRYX1+/4ejRo9lyuRyqqqqurY1nMBikXC5vei2RSJo0dSaTKVe+t0mSJKytrWvv3bsn6ohcMpkMPB5PKhKJUjqzlsWLFxcsX768he9QWloak8lkNu2l0+l0sra2lkaSZIeCIDgcjrz5utatW2fQ1rXhcDjy58c2l68tXF1dq3/55RddKyurOi8vr6pff/1VNz4+nrtt27bHKioquHfvXurJkyfVDx06pLV9+3b9GzdupANAXV0dTVVVVd7mxG8R1CMcBcU7gL6+PgICAgAA4eHhENXqAt7/UShIXCPg5k9AeW6b47ff+wUAMMHGr0V7YXUh/E9NQE5WHDb+5zECtiU1HWvwcEdU8FJEGRuBpNPh5eWF2bNnw8HBoUlBksgkiC+Kw7ob32L8iTFYcGEe9qXsgYyUYbTlGGz4YBP+Hn8SPw37GUPMhlIKEkWX4OXlVRkaGmqgVHSuXr2qCgAVFRV0IyOjBjqdjm3btunIZApFXUNDQ1ZVVdW03WNlZVWfnJzMkclkyMzMVElMTGw1CZ6Tk1NdWVkZ4/z582qAwmIVFxfXpgVEW1tbbmpqWr9r1y4tQGF5un79uioA9OnTp2r37t1aALBr166mSDsfH5/Kffv26VZUVNAAIDs7WyUvL69NA4aurq6My+XKIiMjuQCwe/fuVqP2HB0dJXl5eU1m3rauzcv44IMPqm/evMkrKChgSKVSHDlyRHvo0KFVADB48GDxL7/8YjB48OCqAQMG1Fy7do3HZDLlOjo6soqKClpZWRl90qRJFTt27HiUmpraZOpOT09nOTg41HZIgG6GsiRRULwj8Hg8zJgxA4cPH0ZMTAxKhYMw0NgTRH5jssYDY4AF918YV1pbipSyFAw2+QA8Jq+pPa7wNtZe/xqDY58g8HB+U7uczcaVGUEQN9YpHDVqFMzMzJqeXgurC3G78BYOiPajuqG6aZw+xwCT+JPgaugObfZri7amoMD333+fP2fOnF4CgUBIkiTR6AicuXjx4mI/Pz+r48ePaw0aNEistFZ4eHjUMhgM0tbWVhgYGFj61VdfFf/yyy8SW1tbe1tb21qhUNiqbwybzSYPHTqUtWjRol5isZguk8mI+fPnF7m5udW1JdvBgwcffPrpp+Y//PCDkVQqJT7++OOy/v371/73v/99NGXKFMutW7cajhw5spzL5coAYMKECZXJyclsd3d3AaCw4vz555/ZDAajTRPOzp07c2bPnm2hqqoqHzZsWGVrfdTV1eW9evWSJCUlsRwcHCRtXZuXYW5u3rB69eq8IUOG8EmSJLy9vSumTp1aDgDe3t5Vn3/+OXP48OFiBoMBIyOjehsbmzoAKC8vp48ZM8ZaIpEQAPDdd989Us55+/Zt7vfff1/QkfN3N0R7prR3FTc3NzIuLq67xaCgeC3I5XJERkbi0aNH0NPVwdjEyWCg0fVg+kXA0qtF/y3xm3Hx0QX8PvIP6HP0ISfl+P3+b7h54wi+/E862PXPvgPi/Sei0E4AtqoqRowYASMjIzTIG5DyJAUXcqNw6XF0i7mHmXnDu9dwCHTsoEKj8hW9yxAEEU+SpFvztoSEhBxnZ+cOh5BTtI1YLKapqanJaTQafv31V63Dhw9rX7hwIet1nnPv3r2acXFxnK1bt+a/vPeb4erVq6obN240PH78eHZ3y6IkISFB19nZ2aK1Y5QlieKVkD0pg+zRQ0izs8F0cwdNRxs0DhWZ9Cag0WgYNWoU4uPjcefOHew23YYp+UuhKq8C9gxrUdetqr4KFx9dAF+LD32OPsT1YiyOmo8BfyVj7fln976MQQOROeQDaOrrY8LQoSDUCMQV3caqyH+jpLakqZ86Ux2TbQPhbugBA7VuyZ1HQfFOcvXqVc7nn3/eiyRJqKury3bv3p3zus8ZFBRUXlpa+lbd54uLi1V++OGHvO6Wo6NQliSKTiN7+hTlq75C3fHGxH90OnQPHwKrf7/uFew95MGDBzh//jwA4JPi1dCUFgEjNwEDlgEA/kzdh8Nph/CT138hJ+UI3fspQn7IbBr/1MQE8ZP8oSOwhYGjAeKe3kZE9ukW5+hn1B8fWoyCg67jC07fFD0HypJE8b5CWZIouhSyovKZggQAMhkqvv4aOgcPgK7zzmSb7xH07t0bH3/8Mf7++2/8pb8WY0pDYXwuGOgzAw0sdRxOOwR1pgbS8hMgXxKCkPtiAAAJ4MK8IFR7GCOJdh+Pqy8CdxRzMmgMTBFMQz/j/jBWM+50vTUKCgqKnkK3KkkEQYwC8BMAOoDfSZL8/rnjROPxjwDUAJhBkuSdNy4oRQvI2hd9HGVlZYBMBlIuh7ykBJJr1wE6Hax+nlSh4NeMnp4eAgMDceDAAYTrLsOg8n0QHpqAC15fAAAcc2Rw/HQR5ASQa6aKs/59cc9KDOAeUHsPAOCk64wxVmPhpOtMJXSkoKCgaKTblCSCIOgAfgEwAsBjALcJgjhJkmTz/BI+AGwafzwBbG/8TdGN0HR0QDcxgSzv2bay2tSpoGlqQl5UhOKRoyAvU5RIohkaQv9MOOgGlP/K64TL5WLGjBn46/AhXME0lFRcxrZ7PwMEYHcmGRuXWuGBlTLKWWFNmmwbgP7GA2GhbkFZiygoKChaoTstSR4AMkmSfAAABEEcAjAeQHMlaTyAvaTCceoGQRCaBEEYkST5ToQO9lTo+vrQO34MlT/9F9L0dKj6TYDqRz4gmExUHzjQpCABgLywELXhp8GdNbMbJX4/YDKZCJwyFVHHDyINQ+BeUYjbmjfx+2xzAIC1hjU+tvGDi74LuM1SAVBQUFBQtE53KkkmAB41e/0YL1qJWutjAuAFJYkgiDkA5gBAr169nj9M0cXQjY2hseZroK4OhLo6CBoNJElC9uTFItuysk4V3qb4B9BoNIz8OBA7Ny9EuR4D6k/rUanFxGyHOfA2Hw41lVZz5lFQdIg9e/Y4SySSLrtvsFgs6fTp0xPa65OWlsYcM2aMTUZGRnJ7/VojPDycFxAQYGViYlIvl8uhq6srPXLkyAMTExPpy0e3ZO3atfpLliwp5fF4ncoUvXjxYuOhQ4eKfX19xZ0956uQlpbGdHZ2drCwsKhraGggPD09xfv27XtIp7der/F5SktL6b///rv2F198UfLy3q2f/1X/X8/PEx0dzZ03b16nbyAuLi6Cu3fvdihL+svozozbrdn3nw+160gfRSNJ/kqSpBtJkm56enqtdXnrkEskkBUXQ15R8fLObyE0Nhs0TU0QNMXbiCAIcGfMAGjN3lYqKlDzn9g9Ar6nEASBkXQRrCS2MGcogpV+T/oVAac/wb9jV0BUJmq31AAFRVt0pYL0OuZrDTc3tyqRSJSSnp6e4uLiUr1p06ZXcpIMCwszqKqq6tQ9UyqVYsuWLflvSkFSYmZmJhGJRCkikSg5PT1ddf/+/ZrPy9UWT548oe/cubPbHUkzMjJYhw8f7lRWWuW6ukpBArpXSXoMwKzZa1MAzye86kifdxJZaSnEoZtR4jsBZQs/gzQ7G6T8nShl0y50E2PoRZwGe/RHUB03FvqREaAZGna3WO8d+kQRGCQD5rWW2DlwN37y+i88DD2R/CQZK2KWYfyJMfg78xiq6t/odzcFxT8iJSWFaWdnJ7x8+XKnowvkcjnEYjFdS0tLCgDR0dEcFxcXgZ2dndDFxUWQkJDAAhQ32jlz5pjy+Xwhn88Xrlu3Tv+7777TLy4uVhkyZAjf09OTDwDHjh1T79Onj0AoFNr5+Pj0VpYVMTExcQwODjZydXW13bVrl5afn5/FH3/8oQUAwcHBRg4ODnY2Njb2AQEB5sqyKhs2bNDbsGGDHqCoP+fr62vZr18/vrm5uUNoaKjuq14vFRUVeHh4VGVkZLDCw8N5np6e/LFjx1ra2traA8A333xjYGNjY29jY2O/du1afQBYtmyZ6aNHj1gCgUA4d+5cUwD46quvDBwcHOz4fL5wyZIlxsr5WxuvvIYTJkyw4PP5wlGjRvUWi8W09taflJTEGjBgAN/W1lYoFArtkpOTWSEhISZxcXFcgUAgXLNmjb5UKsXcuXNNlXJs3LhRF1BYC59fF4fDcQEUBYf79+/PFwqFdnw+X/i8stgRunO77TYAG4IgLAHkAZgMIPC5PicBLGz0V/IEUNET/JHktbUQ/7QV1bv+AADIcnNRkngf+lGR73wkGE1NDUwHB2ht+REgCNBUVbtbpPcSFakYJnXJyGPb48KFC5g8eTJW9VsNibQOMY9j8N97P+GPpJ34I2knbDT5mOkwC0Ide8qBm+KtJSEhgTV58mSrnTt3Zg8YMKDDdb+UN9ry8nKGqqqqbMuWLY8BwNnZue7WrVsiFRUVHD9+nLdixQrTyMjIrNDQUL3c3FxWcnJyioqKCoqKiugGBgay7du3G1y+fDndyMhIWlBQwFi/fr1RTExMurq6ujwkJMTw22+/Ndi0aVMBALDZbHl8fHwaAERGRjYVLFy+fHmxso+vr6/loUOHNAIDAytWrFjRYmsrNTVVNT4+PlUsFtNdXFyEfn5+FRYWFg2dvWZisZgWExOjvnr16jwASExMVLt7926yQCCoj42N5Rw4cEAnPj4+lSRJuLq62nl7e4tDQ0MfjxkzRlVZpPfYsWPqmZmZ7MTExFSSJDF8+HDriIgILpfLlbc2XldXV5aTk8MOCwvLGTlyZLW/v7/Fxo0b9dauXVvU1voDAwMtg4ODC4OCgspramoImUxGrFu3Li80NNQgOjo6EwA2bdqkq6GhIUtKSkqtra0l3N3dBWPHjq18fl3N18/hcOSnT5/O1NbWlhcUFDA8PT0FgYGB5TRax+1D3aYkkSQpJQhiIYBIKFIA7CJJMpkgiHmNx3cAOANF+H8mFCkA/tVd8nYlZKUYtc3zDAGQl5ZC/qTsnVeSlFDZt7sXWn0F7BquII9tj8rKSpSXl0NTUxMsBhsjLEZihMVIPKzMxUHRAVzNv4J/X1kJAJhmNx0fWnxIFaKleKsoKytj+Pr6Wh85ciSrvbppreHm5lalvNGGhIQYLly40PTAgQMPG4uvWubk5LAJgiAbGhoIALh48aL6vHnzSlRUFGV2DAwMXqgEe+nSJbWsrCy2h4eHAAAaGhoIV1fXKuXxoKCgp63JEhERwdu8ebNhXV0drby8nCEUCmsBvOBv4ePjU87lckkulyvt379/ZWxsrJqFhUV5R9estAQRBAEfH5/yTz75pDI8PJzn5ORUrVQkLl26xP3oo4/K1dXV5QAwevTop9HR0Tx/f/8W5zl79qx6TEyMulAoFAJATU0NTSQSsauqqmhtjTc0NKwfOXJkNQBMmzbtydatW/UBFLW2/qdPn4qLioqYQUFB5QDA4XBItOJWc/78eXWRSMQ5efKkFgCIxWJ6SkoKm8lkks3X1Ry5XE4sXrzY9MaNG1wajYbi4mLm48ePGb169eqwT1q35kkiSfIMFIpQ87Ydzf4mAfzfm5brtcOgg25q2iIKDAAIHrebBKLoaRCSChgTz7blY2NjMXbs2BZ9eqmbY6Uzbu0BAAAgAElEQVTHv1Evq8fVvCvYcmcz9qXuwb7UPbBQt8Bsxzlw1HWirEsU3Q6Px5MZGRnVX7p0idtZJak5fn5+5f7+/lYAsHLlSpMhQ4aIo6KistLS0pjDhg2zBQCSJEEQRLtOeyRJYtCgQZWnTp1qtf5Ya87dNTU1xLJly8xv3ryZYm1t3bB06VLjurq6Vk0az3/mOvsZVPokPd/O4XCa5OqoXyJJkli8eHHB8uXLW2Reb7699jytyd/W+jshBxEaGvrQz8+vRUHf8PBwXvN1NScsLEz7yZMnjPv376eyWCzSxMTEsba2tlNuRt3pk/TeQtfRgeYP/wHRbCuK+38LQONRYdkUXQfR8OxeUlBQgKqqqlb7MelMePUahhO+4djmHQYvs2HIqczBqqtfYvyJMTiQ+iee1rX6YExB8UZQUVEhz549m3Xw4EGdHTt2dMqZtznR0dFcc3NzCQBUVlbSTU1N6wEgLCysye9n+PDhlTt27NBraFDsbhUVFdEBQE1NTab0Oxo6dGh1XFwcNykpiQUotrUSExPbrdlTU1NDAwBDQ0NpRUUF7dSpU1rKY+vXr9dbv359U8RRRESEZk1NDVFYWEi/ceMGb9CgQdWvuua2GDZsWNWZM2c0xWIxrbKyknbmzBktLy8vsYaGhqy6urpJN/Dx8anct2+frnLt2dnZKnl5eYy2xgNAQUEB8/z582oAcODAAe0BAwZUtbV+bW1tuaGhYf2+ffs0AaC2tpYQi8U0DQ0NWVVVVVNI3ogRIyq2b9+uJ5FICABITExkVVZWtqvDVFRU0HV1dRtYLBZ56tQpXn5+PrOz14kqS9JNqAgEMLgSA+nDh6Dr64OmqQWaBrXFQdF1SGQMMItEqDcQAABu3rwJb2/vdseY8kyxxHUZFvZZhBsF1/HTnR9xKO0ADqUdgLGaMeY6z4ezXh/QCOr56n2ExWJJuzoFQEf7qquryyMjIzOHDh3K53K58qlTp3Zo+0npk0SSJHg8nmzXrl05ALBy5crC2bNnW27dutVw8ODBTdaJJUuWlKSnp7MEAoE9g8Egp0+fXvLll1+WTJ8+vdTHx8dGX1+/4ebNm+lhYWE5kydP7l1fX08AwNdff53n5OQkaUsOXV1d2ZQpU0qEQqG9qalpvbOzc5PiIxKJVAcOHNj0FOPi4lLt7e1tk5+fzwwODi54FX+klzFo0KCawMDAJ3379rUDgGnTppUMHDiwFgBcXV2rbGxs7IcNG1YRFhb2ODk5me3u7i4AFNaoP//8M7ut8WlpaczevXvX7dq1S2fBggXmlpaWkuDg4BIejydva/379+/P/vTTT82//fZbYxUVFfLIkSNZHh4etQwGg7S1tRUGBgaWrlq1qjgnJ4fl6OhoR5Ikoa2t3XDmzJms9tY4e/bsMh8fH2sHBwc7e3v7GktLy05bIakCtxQUPZFvCBRXq+KceBBq+k4El8tFVVUVgoKCwGazOzVVQXUBjqYfwbncyKY2Pxt/jOk9FjqqVK2+ngJV4Lb78PLyso6IiMhis9nk0qVLjblcrmzt2rVF3S3X+0J7BW6px0EKih5KPZ0LWrmidIzyYejOnc6XPjRSM8JCl0U4Nu4EvvD4ElwVLo5mHMG/IoMwM3I6bhfegox8wbeVgoKig0RHR2ey2eyeZ7HoAVDbbRQUPRSSow+iUXmprq6Guro6kpKS4ObmBiaz01vzYNAYGGA8EAOMB6K4phh/ZxzF6exwfHtjDQBgnNV4jLf6GHqcdyOZK8W7ydGjR9VDQkJMm7fl5+czjY2NW0Q3mZmZSaKiotrdjnkb2bx58wu5AG/duqUaFBRk2bytsLCQaWho2GLNTCZTnpiY2GWJFCmo7TYKip6HrAH4lokC7VG4dDUHrEnrUFpaCi8vL0RHR8Pd3R0uLi5dcyq5DPFFcfjl3s94KlFEa6ozNbDQ5TO4G3iATutYKQSK7ofabqN4X2lvu42yJFFQ9DTqFf6QdJ4ugBzYWJiitLQURUVFYDKZuH37NpycnNDRWk7tQafR4WHkCQ8jT5TWluJU1gn8nXkM629+BwDwsfgIH9v4wVCNyrpOQUHx7kH5JFFQ9DTqFUEydK7CqdpQQ1HUNj09HSNGjAAApKWldflpdVV18S+HWfh7/El83X8NDDgGiMg5gzlRszApfCKu5l1Bg7zLg3QoKCgoXhuUJYmCoqfRoLAkMXgK3yBZjSK6WSqVwrCxjt6VK1cgEAjQmfT8HYVO0OFq4IbfRu5CWV0ZTj8Ix5H0w/jh9n8AAN69hsOf/wmMuSZdfm4KCgqKroSyJFFQ9DQaLUkq6oqEuHUVT9CrVy8AQF5eXpM16cGDB69dFG22NqYJg3B8/CmsHfAdzHhmuPDwPOadn4OPT4zD5UeX0CCjrEsUFBRvJ5QliYKip9Hok6SiqbAa1VU8gZ3DUDx8+BApKSn48MMPAQAXL16ElZXVGyk7QiNo6KPvgl+8d6BCUoGI7DM4INqP0PiNCI3fiA9Mh2CS7WSY8Xq9dlkoXp2js/o711eVd9l9g8nVlPrtvJ7QXp+0tDTmmDFjbDIyMpI7O394eDgvICDAysTEpF4ul0NXV1d65MiRByYmJh1OYqlk7dq1+kuWLCltreRIeyxevNh46NChYl9fX3Fnz/kqpKWlMZ2dnR0sLCyaEifeu3cvtatSDJiYmDjGxcWlGhkZdfoaNueLL74w/P777ws7O27SpEnmK1asKHJ1dX3l8jSdgbIkUVD0NBotSYxGS1J18WOYmioiph8+fAiCIPDBBx8AAB4/fvzGxdNgaWCyIAAnxodj/aDv0VvDCjGPL+P/LszHuOOjcSH3PCSyNhMXU3QjXakgvY75WsPNza1KJBKlpKenp7i4uFRv2rTplaqIh4WFGVRVVXXqnimVSrFly5b8N6UgKVHWblP+vI05mLZu3WrU2TFSqRSHDx/OfVMKEkApSRQUPY9GnySCqagF+DRHBDqd3mQxqqurg42NDQDg3Llz3SMjFEUvHXQdscVrK/786BCm2/8LAPDT3R/hf2oC/nNzHXIqcrpNPoq3j5SUFKadnZ3w8uXLnM6OlcvlEIvFdC0tLSkAREdHc1xcXAR2dnZCFxcXQUJCAgtQ3IjnzJljyufzhXw+X7hu3Tr97777Tr+4uFhlyJAhfE9PTz4AHDt2TL1Pnz4CoVBo5+Pj01tZ28zExMQxODjYyNXV1XbXrl1afn5+Fn/88YcWAAQHBxs5ODjY2djY2AcEBJjL5Qqj1IYNG/Q2bNigBwBLly419vX1tezXrx/f3NzcITQ0VLeV5XSapUuXGq9evdpA+drGxsY+LS2NWVlZSRs6dKi1ra2t0MbGxv63337TAoATJ07w7OzshHw+X+jv729RW1vbZHJeu3atgaOjo52jo6Odsn7dgQMHNJycnAR2dnbCAQMG8B89esQAgIqKCtrEiRMtlNdz9+7dmgsWLDCRSCQ0gUAgHDdunCUAbNu2TdvR0dFOIBAIAwMDzaVShaGKw+G4LF682NjJyUlw4cIFroeHh21MTAwHAKZMmdLLwcHBztra2n7JkiXGXXGdnodSkigoehqNliQwFVFtFXmZAABbW1sAQHZ2Nuh0Ojw9PSGTyVBY2GmLd5fDY/LgZzMRJ8aH44fBGyHQssP1gmtYFP1/GHd8NM7lRKJO+sYeHineQhISElh+fn7WO3fuzB4yZEhNR8cpa7cZGxs7xcbG8hYuXFgKAM7OznW3bt0Spaampnz99dd5K1asMAWA0NBQvdzcXFZycnJKenp6yuzZs5+sWrWqWF9fv+Hy5cvpN2/eTC8oKGCsX7/eKCYmJj0lJSW1b9++Nd9++22TAsJms+Xx8fFpc+bMaVEZevny5cVJSUmpGRkZybW1tbRDhw5pAMCKFStKVqxYUaLsl5qaqnr+/PmMGzduiDZu3Gick5Oj0plr9ejRI5ZAIBAKBALhtGnT2t3DPnbsmLqhoWFDWlpaSkZGRvKECRMqa2pqiLlz51oePnw4Kz09PUUqlWLjxo1NWWLV1dVl9+/fT507d27xZ599ZgYAI0aMqLp3754oNTU1ZeLEiWVr1641BIAvvvjCSF1dXZaenp6Snp6eMnr0aPG2bdvyWCyWXCQSpZw8eTL7zp077P/973/acXFxIpFIlEKj0cgdO3boAEBtbS3NwcGhNjExUfThhx+2qNK9efPmvKSkpFSRSJR89epV3s2bN1XRxVA+SRQUPY1GnyQwuWCoqqFerKgDKhAIIBKJcP/+fdjZ2cHe3h43b97ExYsXERgY2I0CP4MgCNjpCLFhyCZU1Vfh4sPz+D3pN/x8byt+vrcVrgZumGo3DVaa1t0tKsUbpKysjOHr62t95MiRLDc3t05py25ublXR0dGZABASEmK4cOFC0wMHDjwsKyujT5o0yTInJ4dNEATZ0NBAAMDFixfV582bV6KiotBLDAwMXqi5c+nSJbWsrCy2h4eHAAAaGhoIV1fXpht4UFDQ0+fHAEBERARv8+bNhnV1dbTy8nKGUCisBVDxfD8fH59yLpdLcrlcaf/+/StjY2PVLCwsOlTQF3i23daRvn379q0NCQkxmz9/vsn48eMrRo0aVXX9+nVVU1NTibJg74wZM5788ssv+gCKAWD69OllAPDpp5+WrVq1ygwAsrOzmb6+vqYlJSUq9fX1NDMzMwkAxMTEqB86dKgpSkRPT++F63n27FleUlISx9nZ2Q4A6urqaPr6+lIAoNPpmDFjRqvXc8+ePdq7d+/WlUqlRElJiUpCQgLb09OztqPXqSNQShIFxSsik5N4Wl2P8pp6cFkMqLEZ4LE79cD3elBaklTUoG5ihbLMRACArq7Cal9eXg6SJMFgMODk5ITExESUlZVBW1u7uyRuFS6Ti3HWvhhrNR4Z5enYk7wb8UVxiC9SZNOf6zQfXmbDwFHp9M4LxTsGj8eTGRkZ1V+6dInbWSWpOX5+fuX+/v5WALBy5UqTIUOGiKOiorLS0tKYw4YNswUUdQ4JgmjXh4ckSQwaNKjy1KlT2W3I+4Jzd01NDbFs2TLzmzdvplhbWzcsXbrUuK6urtXdnOeDKboiuILBYJDK7T0AkEgkBAA4OTlJ7ty5k3L06FGNkJAQk/Pnz1dOmDChXYWseeoQ5bVauHBhr88//7xwypQpFeHh4by1a9caA03Xs13ZSJIk/P39n/zyyy95zx9jMplyBuNFVUUkEjF//vlng/j4+FQ9PT2Zn5+fRVvX85/w0gkJglAnCMKqlXanrhaGguJd4nFZDaZuu4qp267h4y0xOHAtB+LatyCcvbbxoYvBgpaFAIDii4pGo4HFYgEAnj5V9FGWJ7l8+fKbl7ODEAQBvpYt1g36Dw6NPoJ5TgsAAGGJ2zH5tD9Crvwb6U/T0BNLLFEoUFFRIc+ePZt18OBBnR07dryyNh8dHc01NzeXAEBlZSXd1NS0HgDCwsKa/H6GDx9euWPHDr2GBsVnuaioiA4AampqMqXf0dChQ6vj4uK4Sn8csVhMS0xMZLV37pqaGhoAGBoaSisqKminTp3SUh5bv3693vr165u2syIiIjRramqIwsJC+o0bN3iDBg2qftU1K7GwsJDcu3dPDQCuXLnCycvLYwFATk6OCo/Hky9YsKBs8eLFRffu3eP06dOnLi8vj6lc3969e3UGDx7c5Hy+d+9ebQDYuXOnlouLS3XjNaD36tWrAQB2796to+w7dOjQys2bNzc5y5eUlNABhdKmVNRGjRpVGR4erpWXl8cAFNc8PT293QKTT58+pauqqsq1tbVljx49Yly6dEnjn16j1mhXSSII4hMAIgBHCYJIJgjCvdnh3a9DIAqKd4HK2npsCE9BeY3ii5QkgT8uP0BV3T+Kiu0aqhrrYxIEeIYWAABpreI71t7eHgCQmanwU2KxWODz+SgpKUFlZeUbF7WzcFQ4+Kj3aJz0PY0fhyq23+6XJiL48lKMPzEGJzKPo6q+6uUTUbwSTK5ml77BOzOfurq6PDIyMvPnn3822L9/v2ZHxyl9kmxtbYUHDx7U2bJlyyMAWLlyZeE333xj2rdvX4FM9mwHaMmSJSWmpqb1AoHA3tbWVrhz505tAJg+fXqpj4+PjaenJ9/Y2FgaFhaWM3ny5N58Pl/o6uoquH//Prs9OXR1dWVTpkwpEQqF9j4+PtbOzs5Nio9IJFLV0dFpEsLFxaXa29vbxtPT0y44OLjAwsLiHz99BQUFPX369CldIBAIf/75Zz1zc/M6AIiPj1ft06ePnUAgEP7www9Gq1evLuBwOOSOHTty/P39rfh8vpBGoyE4OLjJZ0oikRBOTk6Cbdu2GWzduvURAISEhOQHBARYubq62uro6DT9X//zn/8UlJeX021sbOxtbW2FZ86c4QHAlClTSuzs7ITjxo2zdHV1rVu1alWet7c3n8/nC4cNG8Z/9OhRu2b5/v371zo4ONTY2NjYT5s2zaL5dmdX0m6BW4Ig7gHwIUmygCAIDwB7AXxJkuQxgiDukiTZNVUyuxiqwC3F66ZUXIeZv95EcWVLy//OTz1hb9rh7+/Xw57hQPYF4BsS2TEncOOXLzDmp7PgGZqjvLwcf/31FxgMBmbOnAkAqK6uxp9//glzc/OmHErvEnXSOsQ8voyf721tahNo2WGGw0zYadu9kTxQPQGqwG334eXlZR0REZHFZrPJpUuXGnO5XNnatWuLuluu94V/UuCWTpJkAQCQJHmLIAgvAOEEQZgCoGzbFO8tPLYKvIT6OHzjYVMbh0mHgUa7D5Nvhupn361sDcUuQl35E/AMzaGhobBIS6VSyGQy0Ol0qKmpwdTUFLm5uaipqQGH8275+LAZbIy0+BAjLT5ETkUODor+xPWCa/gidjkAIEg4AyMtPoQ6U72bJaWgaB2lYznF28fLlCQxQRBWJElmAUCjRWkogOMA7F+3cBQUbyssFTpmDO6NBhmJi8mFMNHiYOVYITQ47W6jvxlqnj34szUU7ht1FYo2giCgoaGBiooKFBQUNCWZHDhwIA4fPoy4uLimRJPvIhYaFvi3ZwgkMgmuPI7FT3d/xN6U3dibshu9Naww2/FT2Os4UNald5ijR4+qh4SEmDZvy8/PZxobG9c3bzMzM5NERUVlvVnp/jmbN2/Of77t1q1bqkFBQZbN2woLC5mGhoYt1sxkMuWJiYmi1y3j+8TLlKT5AFp8m5AkKSYIYhSAT16bVBQU7wBaXBYWjbTFzCFWYNAIaKq9BQoSANQ8afqTraHwn6yrLGtqc3R0xJUrV5CamtqkJGloaEBHRwcikQj9+vUDk/mWrOUVYdFZ8DYfDm/z4Xgkfoi/0g7j8uNL+PLKFwCAQMFU+Fh+BA3Wa/H1pHiN+Pn5Vfr5+XUovL2n4OHhUdvRkH6KruVl0W3VAAxaae8H4EbXi0NB8W7BZtKhy2O9PQoSAMgbACYXAMDiKQJoasuKmw4ri91mZ7eMXh4yZAgAICGh3VJa7xxmvF5Y5rYcR8cexzLX5WAQDBwQ7ce0iEAsOD8P94rvQk52qhwXBQXFe8LLlKQtAFqrOVPbeIyCguJtRFVhQaIxFAEi5blpTYe4XG7T3zU1zxIX6+rqQlVVFXfv3oWyJEBPQoWugiFmQ3Fs/AnsGP4rhvcagcdVj7D62ir4nhiLfSl7UVZX9vKJKCgo3htepiRZkCSZ+HwjSZJxACxei0QUFBT/HLWWNTzLH6W3eG1oaAjgxQK3w4cPBwCkpPRsy74x1wSL+i7G0XHHsdL93+AwODiSfhgzzk7Dp+dmIr4oDjLyhcTAFBQU7xkvU5LaC9Xp8hopPQnZ06eQ5uaiPikJsuJiKtEdxZuFa9jsBYHq4pbKkDJf0vNba0rl6caNG2ienbenokJTwUCTQTg05gh+HbETPhYfoaimCGuuf42PT4zDrqTfUVpLRcBTULyvvMxx+zZBEJ+SJPlb80aCIGYBiH99Yr3byMrKUPntOtT89RcAgKanB72Tx8Ho1W6dQQqKroP3rCA2z8gC4oKW/kcmJiYAFJm3m5cNIAgCH374ISIjI5GZmQk+n//mZO5mDNUMMb/P/2GO0zzcLrqFX+7+F8cz/8bxzL+hxdKCuboFKuor8E3/tdBia718wp7ID7rOqH3SdeWsVHWkWFnarhNcWloac8yYMTYZGRnJnZ0+PDycFxAQYGViYlIvl8uhq6srPXLkyAMTE5NO7yevXbtWf8mSJaWtlRxpj8WLFxsPHTpU7Ovr25rrSpeTlpbGdHZ2drCwsKgjSRIcDke+e/fubGdnZ0ln59q6davOuHHjKrsimSUAbNiwQY/D4cgXLlz4pHl7R/7HaWlpzOjoaO68efPe6J74yyxJiwH8iyCISwRBhDb+XAYwG8Dnr1+8dxN5UXGTggQA8pISVPznB8ir/3FmeQqK9pE3bhFxn8VbaJq/qOiw2c+MxMoSJUqUjt2XLl16Ly2gdBod/Yz6Y99HB7Bz5G6Ms/LFU8lT3Cu5i9zKXHx2cQEelL9zkeVdQ1cqSK9jvlZwc3OrEolEKenp6SkuLi7VmzZt0n/5qBcJCwszqKqq6lRtMKlUii1btuS/KQVJibLAbVpaWkpgYGDpmjVrjF5lnv379+s+fPiwUwUpleVcWmPFihUlzytIHSUjI4N1+PDhN15gst1/OEmSRSRJDgCwBkBO488akiT7kyRZ+PrFezeR5j1+oU2WlQWy7pXrMlJQdIyGRkds1rPEiRpmCiVJ1tAipQqsrBQlGZ+PciMIAl5eXgCA3Nzc1yXpO4EeRw+zHT/FF+5fgkbQwKIzUdVQheUxy/C0rtXC5BSvkZSUFKadnZ3w8uXLnc54KpfLIRaL6VpaWlIAiI6O5ri4uAjs7OyELi4ugoSEBBagUGzmzJljyufzhXw+X7hu3Tr97777Tr+4uFhlyJAhfE9PTz4AHDt2TL1Pnz4CoVBo5+Pj01tZ183ExMQxODjYyNXV1XbXrl1afn5+Fn/88YcWAAQHBxs5ODjY2djY2AcEBJgrt7Q3bNigt2HDBj0AWLp0qbGvr69lv379+Obm5g6hoaG6rSynw1RWVtI1NTVlgMIa4+rqaisUCu2EQqFdVFSUmrLfqlWrDPh8vtDW1la4YMECkz/++EMrKSmJExQU1FsgEAirqqqI2NhYjru7u629vb3doEGDbHJzc1UAwMPDw3bhwoUm7u7utt99951Beno6s3///nw+ny/s378/PyMjg6lc2+rVqw0AIDY2lmNrayvs06ePoHltN6lUirlz55o6ODjY8fl84caNG3UBICQkxERZYmbNmjX6bfXLzc1VcXNzsxUIBEIbGxv7s2fPPotUeQXa1eIJgmADmAfAGsB9ADtJkux5YS9dDNPeHmAygfpnNyVVv49B0+zmchUUPZ/6RmulStN3H1S1GrNuVzyBmu6zB0o+n4+srCzcu3cPrq6uLaaxsrJCdHQ0zp07h08//fS9T77Yz7g/fh/5B8T1YlTVi1EjrYEmi/o8v0kSEhJYkydPttq5c2f2gAEDajs6TnljLS8vZ6iqqsq2bNnyGACcnZ3rbt26JVJRUcHx48d5K1asMI2MjMwKDQ3Vy83NZSUnJ6eoqKigqKiIbmBgINu+fbvB5cuX042MjKQFBQWM9evXG8XExKSrq6vLQ0JCDL/99luDTZs2FQAAm82Wx8fHpwFAZGRkUzKu5cuXFyv7+Pr6Wh46dEgjMDCwYsWKFSXNZU5NTVWNj49PFYvFdBcXF6Gfn19FZ7a8Hj16xBIIBMLq6mpaXV0d7dq1ayIAMDY2lsbGxqZzOBzy/v37rICAgN5JSUmpf/31l/rp06e14uPjRTweT95szfqbNm169MEHH9RIJBJi0aJFvU6fPp1pbGws/e2337SCg4NNjhw5kgMA5eXl9Nu3b6cBwLBhw6wDAwOffPbZZ0+2bNmiM3/+fLPz58+3ML/OmjXL4scff3w4evToqrlz5zYlB92yZYuuhoaGLCkpKbW2tpZwd3cXjB07tnLdunV5oaGhBsrs5Js2bWq138GDB7W8vb0rfvjhh0KpVAqxWNwp69/zvMzUuQdAA4BYAD4A7KDYgqNoB0JbG3rHjqJ81VeQFRWCM2kSOBP9QdDp3S0aRU+noVFJYj57eFKWJpFUlrVQkpRO2jKZDFKpFAzGs68DGo2GAQMG4Nq1aygoKICx8TMfp/cRGkGDrqoudFX/0UM9xStSVlbG8PX1tT5y5EiWm5tbp0zybm5uVcoba0hIiOHChQtNDxw48LCsrIw+adIky5ycHDZBEGRDQwMBABcvXlSfN29eiYqKYpfJwMDghTDHS5cuqWVlZbE9PDwEANDQ0EA0L7AaFBTUqpkxIiKCt3nzZsO6ujpaeXk5QygU1gKoeL6fj49POZfLJblcrrR///6VsbGxahYWFuUdXbNyuw0AfvvtN62ZM2eax8bGZtTX1xOzZs0yT0lJUaXRaMjNzWUBQFRUlPrUqVOb/K1aW3NiYiIrIyNDddiwYXxAYZnT09NrUtwCAgKafIXu3r2rFhERkQUA8+fPL1uzZk2LDOlPnjyhi8Vi+ujRo6sAYObMmU8uXryoAQDnz59XF4lEnJMnT2oBgFgspqekpLCZTGaLvf+2+vXr16967ty5Fg0NDbSJEyc+7YxC3RovU5KEJEk6AgBBEDsB3PonJ3tfoLFYYLr0gc6+PYBUCkJTE7R3PIMxxTtCfeP3NPOZJYmt3rI0iRLlTQAACgsLm7JvKxEIBLh27RrOnz+PoKCg1yQwBcXL4fF4MiMjo/pLly5xO6skNcfPz6/c39/fCgBWrlxpMmTIEHFUVFRWWloac9iwYbYAlIEM7TrjkSSJQYMGVZ46dSq7teOtOXfX1NQQy5YtM79582aKtbV1w9KlS43r6upatXI8b7n9J5bcgICA8kWLFlkAwLp16wz09fUbjh49mi2Xy6GqquqqXM/LzkGSJHixUWwAACAASURBVGFtbV177969VsuedMahvb3zkSRJhIaGPvTz86ts3h4eHs7rSD8AiImJSTt69KjGjBkzLBctWlT0qn5QwMsdt5u0RGqbrfPQtbVB19enFCSKN0d9K5YkTeV224tBIY6OjgCAtLS0F44xGAz07dsXdXV1KC2lwuApug8VFRXy7NmzWQcPHtTZsWPHKzvvRkdHc83NzSWAwlfH1NS0HgDCwsKaTITDhw+v3LFjh57SAbmoqIgOAGpqajKl39HQoUOr4+LiuElJSSwAEIvFtMTERFZ7566pqaEBgKGhobSiooJ26tSpphDJ9evX661fv15P+ToiIkKzpqaGKCwspN+4cYM3aNCgV476iYqK4pmZmUkAoKKigm5kZNRAp9Oxbds2HZlMYTAaNWpU5b59+3SVW1PKNXO5XFlFRQUdAJycnOrKysoY58+fVwMAiURCxMXFtZomyMXFpfr333/XAoCwsDBtNze3qubHdXV1ZVwuVxYZGckFgN27dzf9T0eMGFGxfft2PYlEQgAKC1ZlZSVNQ0NDVlVVRX9Zv/T0dKaJiUnDsmXLSqdOnVp6586df1Sx+2WWJGeCIJRaGgFAtfE1AYAkSZIqq01B0QEqJZWQkVIw6SyoNfMX6nKUliSVVixJ5S8qOpaWlrh//z6ysrLg7e39wnEnJyfcuXMH0dHR8Pf3fz0yU7w7qOpIuzwFQAdRV1eXR0ZGZg4dOpTP5XLlU6dO7dD2k9IniSRJ8Hg82a5du3IAYOXKlYWzZ8+23Lp1q+HgwYObrBFLliwpSU9PZwkEAnsGg0FOnz695MsvvyyZPn16qY+Pj42+vn7DzZs308PCwnImT57cu76+ngCAr7/+Os/JyanNMHtdXV3ZlClTSoRCob2pqWm9s7Nzk+IjEolUBw4c2KRIuLi4VHt7e9vk5+czg4ODCzobgq/0SSJJEioqKuSOHTtyAWDx4sXFfn5+VsePH9caNGiQWFVVVQ4AEydOrLxz5w6nT58+dioqKuTw4cMrfv7557ygoKDSzz77zHz58uXyuLi41EOHDmUtWrSol1gspstkMmL+/PlFrVn2tm/f/nD69OkWP/30k6GOjo507969Oc/32blzZ87s2bMtVFVV5cOGDWt+/UtzcnJYjo6OdiRJEtra2g1nzpzJ8vDwqGUwGKStra0wMDCwdNWqVcWt9YuMjORt3brVkMFgkBwOR/bnn3+2au3rKERPDPF1c3Mj4+LiulsMCgrISTnyq/Lx050fkV2ZjT66fTC/zwLovC7fltS/gcMTgHkJgKFTU/PBSXawHPIx+i1Y36K7TCbDzp07AQBTpkyBmtqLClxsbCxSU1MxadIkaGhQBWF7KgRBxJMk6da8LSEhIcfZ2ZkyI75mvLy8rCMiIrLYbDa5dOlSYy6XK1u7dm1Rd8v1vpCQkKDr7Oxs0dqxf+T1TUFB0T4VkgqsvhaCtKci1MskuFV0Ez/d+RFV9VUvH/wqtOKTpORp7ouuBHQ6vck3KS8vr9Up3dwU980rV650kZAUFBTNiY6OzmSz2T3PYtEDeO2JvCgo3mfqpLUvlLW4V3IP9fL6Nkb8Q1rxSVIizn/Q6pA+ffrg9u3bSEhIaDXDtqqqKiwsLJCTk4Pq6upWrU0UFG+Ko0ePqoeEhLSIMsjPz2caGxu3+FCZmZlJoqKi3rmsn5s3b85/vu3WrVuqQUFBls3bCgsLmYaGhi3WzGQy5YmJia06VlO8GpSSREHxGmHSWWDSWaiXPXNVMFYzBg2vKe9QKz5JgMJ5uzWfJAAwMzPD7du3XyhR0pz+/fsjJycHt27dako0SUHRHfj5+VX6+fn17ArMz+Hh4VGrDOmneLN0y3YbQRDaBEFEEQSR0fi71UJIBEHsIgiimCCIpDctIwVFV8BV4WKRy+dg0BTPIxwGB0tdg6HxuhIRShr9H1VaBnRomNm0OURHR6fp7ydPWo+U5fF4MDAwQEZGBuqozPHdjqS6HjXltZBLe34RYgqK7qS7fJK+AHCBJEkbABcaX7fGbgCj3pRQFBRdDYvBgoehJ34bsQvbvHdg2/Aw9Nawen0ZrKsKFL9pLT/aWhZ2AABS/uJNlSAIqKsrAlVzcnLanHrw4MEAgLt373aBoBSvgkwmx9PHFbiwKRanVp3H/7d372FRX+e+wL9rZpgZYBhuIlcFFRAQRRRRo9neU0nTxIamMU2jez87sWl3TpsYq0nt7sU0npw05uTJ6W41OyY1bZo0qdY02VFjFC9JqsYbKFdREeUiIHId7rPOH3MRYWBmYIYB/H6eh4e5rN+a9zcKvLPW+q03+6N8tDQwaSVyF08lSQ/AtJs3zN9X2GokpTwCYEgr/hK5mlalRbB3MKL8xiFIGwSV0o2z3I0VNh/2DTbttN3WZPuq6WnTTFfCnT17ts+ug4KC4Ofnh3PnzvVbxJLcp7W+FbvX78XV0xWou1qPE++cQcFnxRxRInITTyVJoVLKCgAwfx9QVebuhBBrhBAnhRAnq6ur7R9ANBo12a47rQ0wTam11dueTouMjARgKjXQXwK0ePFiAEBubu5goqQBaqhoRHvz7f8+hQcvobWpz+15iGgQ3PaRVgjxOYAwG09tdMfrSSnfAPAGYNonyR2vQTTsGWwvztboTUlSS/0Nm+uTLNNtAFBRUYHx48fb7Cc0NBQqlQonTpzA1KlToWQ9wiGl1ffe4Fg3xgdK1dB83n3005Upje2NLvu74af263z33vez+2tTWFiovu++++IuXLjgdGb+ySef+D3yyCOTIiMj241GI8aMGdP54YcfXoqMjHS6gsSmTZvGPvPMMzXOlN8AgKeffjpi4cKFjStWrGh09jUHorCwUJ2SkpIcExNjnYc9e/Zs/s6dO/W5ubnemzdvtv1JygEPP/xw9Pr166/PnDnTqTnempoa5Ztvvhn03HPPjbgRDLf9ZEkpl0opk218fQTguhAiHADM36vcFQfRHaWPJEnr3/9IkhDCWvC2qKio35e45557HGpHruftr0HsghjrfS9vFe56PA0aXb8VMVzGlQmSO/qzJS0tramgoCCvqKgoLzU1tfmVV14Z0MzFtm3bQpuampz6m9nZ2YnXXnutfKgSJAtLgVvLl1arlY8++mi9rQTJmanzv/71r1ecTZAAU0Hb7du3D3rGyBM8Nd32DwCrzbdXA/jIQ3EQjS4dzYCyd61Arb+lyG3fdR4TE02Luy9dsr2fkoVlau7o0aMw2lgITu6j1Wsx7/E0fPd338J9LyzFw7+/H/4RfvYPHCXy8vLUiYmJSYcPH3a6HpfRaERjY6MyMDCwEwCysrJ8UlNTExITE5NSU1MTsrOzNYApsVmzZk1UfHx8Unx8fNKLL7449je/+c3YqqoqrwULFsTPnj07HgB27dqlnz59ekJSUlJiRkbGREtdt8jIyKnr1q0Lnzlz5uS33norMDMzM+btt98OBIB169aFJycnJ8bFxU155JFHoi0/Py+//HLIyy+/HAIAa9eujVixYsWEOXPmxEdHRydv2bLFJdvzv/7668GrVq0aDwCZmZkxjz/+eNTs2bPjf/SjH0WtXbs24sEHH4yZN29eXGRk5NQdO3YEPPnkk1Hx8fFJd999d5ylPlp6evrkI0eO+ADA3/72N31SUlLi5MmTk+bOnRtvif0Xv/hFqOU14+LiphQWFqqfffbZKEuplB/84AdRAPCf//mfocnJyYnx8fFJzzzzTAQANDQ0KBYuXBg7efLkpLi4uCn//d//bfPK96HkqX2SXgLwgRDi3wGUAngIAIQQEQDelFLea77/HoCFAMYIIa4B+KWUcrtnQiYaIbx71/9U+5rKiTTX9NqnzsqS/ABAU1MTdLreG1ICplGnJUuW4MCBA7h8+TImTZo0yIDJGVq9Flq9FoHj7qwSMdnZ2ZqVK1dO2r59++W77rqrxdHjLLXb6urqVN7e3l2vvfbaNQBISUlpPXHiRIGXlxd2797tt379+qh9+/Zd3LJlS8iVK1c0ubm5eV5eXrh+/boyNDS06w9/+EPo4cOHi8LDwzsrKipUmzdvDj9y5EiRXq83bty4MeyFF14IfeWVVyoAQKvVGk+dOlUIAPv27bP+Q/30pz+tsrRZsWLFhPfff9//e9/7Xv369etvm4bKz8/3PnXqVH5jY6MyNTU1KTMzs96Z+m2WhAQAZs2a1fSnP/2ptGebixcvar/88ssilUqFtWvXRly5ckXz1VdfFZ0+fVq7ePHihB07dlzcunXrtWXLlk364IMP/B977DHrVR/l5eWqp556KubQoUMFCQkJ7ZaCuH3ZsmXLtfvuu8/bstfTrl279MXFxdqcnJx8KSWWLl0au2fPHt3169dVYWFhHYcOHSoGTCNQjp6zu3gkSZJS3gDQq5qmlLIcwL3d7j8ylHERjQo+Ib0eEuYtAW5ezu/7MJ9bH86vXr1qHVmyZcIE0+a/Bw4cwMSJE923pQERgNraWtWKFStiP/zww4u2Cqr2Jy0trSkrK6sYADZu3Bj21FNPRf3lL38pra2tVT788MMTSkpKtEII2dHRIQDg4MGD+ieffLLaUq4nNDS0q2efhw4d8r148aI2PT09AQA6OjrEzJkzrbWGVq1addNWLHv27PF79dVXw1pbWxV1dXWqpKSkFgD1PdtlZGTU6XQ6qdPpOufOndtw9OhR35iYGIcK+gK3ptv6a/Pggw/eVKlupQBLly6t12g0Mj09vaWrq0t85zvfaQCAKVOmtFy+fPm24elDhw75pqenNyYkJLQDtt+j/uzdu1d/5MgRfVJSUhIAGAwGRUFBgXbJkiWNGzduHPfDH/4w8oEHHqhfvny5m+o3OY6124hGG11on0/VXyvu91BL8nPu3Ll+2ykUCuu+SX3VfCNyFT8/v67w8PD2Q4cO2R7edFBmZmbd8ePH/QBgw4YNkQsWLGi8cOFC7scff1zc3t6uAGDZdb7fi3+klJg/f36DZc3PxYsXcz/44IMr3eLtNQ9tMBjEs88+G71r166LRUVFed///vdrWltbbf4N7vmhwx0fQnQ63W0xajQaCZjqOapUKqkwf7BSKBTo7Oy8LYC+duZXqVSy+xS8ZZquJyklnn766QrL+1daWnr+mWeeqZk2bVrb6dOn86ZOndqycePGyHXr1oUP+kQHiUkS0Wijs/17RaXxRltD/9uOxcWZrnyrq6uzu97IUudt//79AwiSyHFeXl5y7969F997773grVu39p5PdlBWVpYuOjq6DQAaGhqUUVFR7QCwbds267qfpUuXNmzdujXEsqDZMpXk6+vbZVl3tHDhwuaTJ0/qzp8/rwGAxsZGRU5OTr+r5w0GgwIAwsLCOuvr6xUff/yxdb3N5s2bQzZv3mwdAt6zZ0+AwWAQlZWVymPHjvnNnz+/eaDn7A6LFi1qPn78uF9BQYEauPUexcTEtJ09e9YXAL744gufsrIyDQD4+/t3NTc3W/ONjIyMhj/96U9jLO/n5cuXvcrKylQlJSVefn5+xh/96Ee1Tz/99PWzZ886vfbM1Vi7jWi0kOYPv34RNp/2i5yIm5f6v4o6PPxWglVTU4OxY/u+IEWpVGLWrFn4+uuvUVVV1W9bGh381H6drt4CwNG2er3euG/fvuKFCxfG63Q64/e//32Hpp8sa5KklPDz8+t66623SgBgw4YNlY8//viE119/Pezuu+9usLR/5plnqouKijQJCQlTVCqVXL16dfXPfvaz6tWrV9dkZGTEjR07tuP48eNF27ZtK1m5cuXE9vZ2AQC//OUvy6ZNm9bnhlVjxozpevTRR6uTkpKmREVFtaekpFgTn4KCAu958+ZZp5ZSU1OblyxZEldeXq5et25dhTPrkYZCRERE5+uvv17y7W9/O9ZoNCI4OLjjq6++urBq1aqb7777bnBCQkLS9OnTm6Ojo1sBICwsrGvmzJlNcXFxUxYvXly/bdu2a7m5udpZs2YlAICPj4/x3XffvVxQUKB5/vnnoxQKBVQqlfz9739/pf9I3E9IOfq2FEpLS5MnT570dBhEQ6ujBXjRB1j6EjB/Q6+nj2/9OS5l7cTK9/P6Hb5/4403AACpqamYNWtW/y/Z0YG3334ber0eK1euHFz85FFCiFNSyrTuj2VnZ5ekpKTY3leCXGbRokWxe/bsuajVauXatWsjdDpd16ZNm657Oq47RXZ29piUlJQYW89xuo1otGg3fzD18rX5tD7CtN6os9XQbzfmtZQO1Wjz8vLC1KlT0dDQgJs3ba5VJSI7srKyirVa7egbsRgFON1GNFq0m0fr1bbXtlo2lGytvwEvb9uJFADExMQgL890YUx7ezvU6t77LnWXmpqKc+fO4ciRI3jggQcGEDiR43bu3KnfuHFjVPfHysvL1REREe3dHxs3blzb/v37Lw5tdIP36quv9tqn48SJE96rVq2a0P2xyspKdVhY2G3nrFarjTk5OQXujvFOwiSJaLToMI8kqW0nQJpuSZJfmO2yIwCsO28DQHl5OWJiYvp9Wa1Wi9jYWBQXF6OxsRF+fnfG5obtBtPfJ7VP/0kkuVZmZmZDZmZmv5e3jzbp6ekt9i7pJ/fgdBvRaOHgSFJbQ9+7bgOASqWyrlkqLCx06KXT09MBAP/85z8daj+SdbR2oOrCDRx45UsceOVLVBXVoL1lWK2rJSIXYZJENFrYWZNknW6r6z9JAoDp06cDAK5cceziEp1Oh8jISJSUlKClxeHNkEekphoDdq/fi9JTZSg9VYbdG/ahuab/dV5ENDIxSSIaLazTbX2MJOlN28u03LRfiHv8+FvTcQ0NDf20vGXevHkAgFOnTjnUfqQq2F8Maby1xlYaJfL2stgv0WjEJIlotLBOt9keSVKoTGUW6q7an0ILCblV2uTq1asOvXxAQAACAwORl5eH9vZ2+weMUL5Bvfe387HxGBGNfEySiEaL9v5HkizqSu2PeigUCnh7ewMAcnJyHA5h0aJFTh8z0ky6Owa+Y24lRb7BPohfNKGfI0aPiuRpKWWR42a66qsieVqKvdcsLCxUx8XFTRlIvJ988omfn5/f9ISEhKT4+Piku+66K76srGxAFyxt2rRpbGNjo9N/M59++umI3bt3D/nVDL/+9a/HajSaGY4Uie15bgsWLIitqakZUHHZ9PT0yUeOHBn0p4bh8n4zSSIaLSwjSX2sSbJoquxVENymlBTT36/Gxka7JUosxowZA61Wi9OnT6Oz0+HNlEcU3yBvPPhKBpb/fCG+sXEhHtySYXN0aTQy3rzp0iuiXd2fLWlpaU0FBQV5RUVFeampqc2vvPLKgLaG37ZtW2hTU5NTfzM7Ozvx2muvla9YsaJxIK85GH/729+Ck5OTm999990Ae217ntvhw4eLx4wZ41TRWlcbLu83kySi0aLd/Huhj+k2AND1c+l/T1FRt7aiqaqqcvi4pUuXAgAKCkbvdi0+gd6InhWFmPQo+AR6ezqcO0ZeXp46MTEx6fDhw05npUajEY2NjcrAwMBOAMjKyvJJTU1NSExMTEpNTU3Izs7WAKY/tGvWrImKj49Pio+PT3rxxRfH/uY3vxlbVVXltWDBgvjZs2fHA8CuXbv006dPT0hKSkrMyMiYaKlDFhkZOXXdunXhM2fOnPzWW28FZmZmxrz99tuBALBu3brw5OTkxLi4uCmPPPJItOXDx8svvxzy8ssvhwDA2rVrI1asWDFhzpw58dHR0clbtmwZY+N0+pWbm6sxGAyKTZs2lX3wwQfWWneOnltkZOTUiooKVWFhoXrixIlTVq5cGR0bGztl3rx5cU1NTQIAvvrqK++UlJSE+Pj4pGXLlk2qrq62jjz98Y9/DE5NTU2Ii4ubkpWV5TPc3m9nMEkiGi0M5uoRSq8+mwSMn+xwd4GB1vqbKCkpcfg4S/23r776yuERKCJ7srOzNZmZmbHbt2+/vGDBAocvJ7TUbouIiJh29OhRv6eeeqoGAFJSUlpPnDhRkJ+fn/fLX/6ybP369VEAsGXLlpArV65ocnNz84qKivIef/zxGz//+c+rxo4d23H48OGi48ePF1VUVKg2b94cfuTIkaK8vLz8GTNmGF544YVQy2tqtVrjqVOnCtesWXPbNvQ//elPq86fP59/4cKF3JaWFsX777/vDwDr16+vXr9+vfWKivz8fO/PP//8wrFjxwp++9vfRpSUlPT9Q23Djh07gh588MHa5cuXN12+fFlrmWJ05Nx69lVaWqr98Y9/XFVcXJzr7+/f9c477wQCwL/+679O2Lx587WioqK8KVOmtGzYsMFaNNJgMCjOnDlT8Prrr19Zs2bNhOH2fjuDm0kSjRYNZXabBIyPx7UT+9HV2Q6lqv9NEIUQCAoKQm1tLXJycjBnzhyHwhBC4J577sFnn32GixcvIi4uzqHjiPpSW1urWrFiReyHH354MS0trdWZY9PS0pqysrKKAWDjxo1hTz31VNRf/vKX0traWuXDDz88oaSkRCuEkB0dHQIADh48qH/yyServbxMeUloaGivaadDhw75Xrx4UZuenp4AAB0dHWLmzJnWArWrVq2yWaNnz549fq+++mpYa2uroq6uTpWUlNQCoL5nu4yMjDqdTid1Ol3n3LlzG44ePeobExPjUEFfAPj73/8etGvXrmKlUomMjIyb77zzTuDzzz9f7ci59RQZGdl21113tQBAamqqoaSkRHPjxg1lY2Oj8pvf/GYTADzxxBM3HnrooYmWY773ve/Vms+jqampSVFTU6Osq6tTDNf3uz9MkohGi6ZKu028A03LMdoabsInKNROayA5ORlHjhwxHdPWBo1G41Ao0dHRAICsrCzExsb2W1CXyB4/P7+u8PDw9kOHDumcTZK6y8zMrHvooYcmAcCGDRsiFyxY0Lh///6LhYWF6sWLF08GACklhBD91lGTUmL+/PkNH3/88eU+4u01hGowGMSzzz4bffz48bzY2NiOtWvXRrS2ttqczen58+LMz8/x48e9r1y5olm+fHk8YEooxo0b1/b8889XO3JuPanVamt7pVIpW1pa7M5A2Yp/OL/f/eF0G9Fo0Wx/3ZDW37Q8obXOscLuERHWEXSUldkfqbIQQmDBggUAHN9CgKgvXl5ecu/evRffe++94K1btwbZP8K2rKwsXXR0dBsANDQ0KKOiotoBYNu2bdZ1P0uXLm3YunVrSEeHaRf169evKwHA19e3y7IOZuHChc0nT57UnT9/XgMAjY2NipycnH4/QRgMBgUAhIWFddbX1ys+/vhj63z25s2bQzZv3mzdd2PPnj0BBoNBVFZWKo8dO+Y3f/78ZkfP8Z133gl69tlny8vKys6VlZWdq6qqyqmsrFQXFRWpHTk3RwQHB3fp9fquvXv36gBg+/btwXPnzrWO7Lz33nuBALBv3z6dn59fV3BwcNdwer+dwSSJaLQw2E98rLtuN9Q61KVer7fednYhtmWabd++fU4dR8OXwrzo2RP96fV64759+4p/97vfhf75z3+2e8WWhWVN0uTJk5Pee++94Ndee+0qAGzYsKHyV7/6VdSMGTMSurpuzfA888wz1VFRUe0JCQlTJk+enLR9+/YgAFi9enVNRkZG3OzZs+MjIiI6t23bVrJy5cqJ8fHxSTNnzkw4d+6ctr84xowZ0/Xoo49WJyUlTcnIyIhNSUmxJj4FBQXewcHB1iBSU1OblyxZEjd79uzEdevWVcTExDhc92b37t1B3/3ud2+bmsvIyLi5Y8eOIEfOzdHXefvtty9v2LAhKj4+PiknJ8f7pZdeshbmDQwM7EpNTU146qmnordt21YCDK/32xlCSqdG3kaEtLQ0efLkSU+HQTS0fmUe4v5V3z/TjZWl+OQn38DsH/1vTFywwqFu/+d//sc6ivTEE084NfSfk5ODY8eO4f7777+tcC4NP0KIU1LKtO6PZWdnl6SkpDg27EgDtmjRotg9e/Zc1Gq1cu3atRE6na5r06ZN1z0d150iOzt7TEpKSoyt5ziSRDSaePc/E6ENMBe5rbdfv81i8uRbV8TV1zu15hFJSUkAgM8//9yp44juJFlZWcVarXb0jViMAly4TTSa+PS/pYqX1rSHUkO5zfWPNkVGRlpvX716FQEBDs90QKVSITU1FWfOnMGNGzcQHBzs8LFEtuzcuVO/cePGqO6PlZeXqyMiIm6rhTNu3Li2/fv3Xxza6Abv1VdfLe/52IkTJ7xXrVp127bulZWV6rCwsNvOWa1WG3NyckbvBmUewCSJaDTxdWwz4ZtXHP89ailPAgDZ2dmYOnWqUyFNmzYNZ86cwaFDh5CZmenUseRxRqPRKBQKxbAZ5cjMzGzIzMzM83QcQyk9Pb2loKDgjjrnoWI0GgWAPjd043Qb0Wiic2zdT2PZJae6jY2NBQAYDAZ0X3TpCI1Gg4SEBNy4cQMNDQ1OHUsed766utrf/IeEaFQxGo2iurraH8D5vtpwJIloNPGLsNtEow9Cm4NXt1lMmjQJxcXFAIDr16/ftjWAI9LS0lBQUIAvv/wSGRkZTh1LntPZ2fl4ZWXlm5WVlcngh2oafYwAznd2dj7eVwMmSUSjQZf5CmGfkP7bAfCPikVV3gmnureUGgGAy5cvO50k+fj4YPz48SgtLYXBYICPz51REHakmzlzZhWA+z0dB5Gn8JMB0WjQbt4CRK2z2zRwQiIAQDpRV02tvlXCJDc317nYzO666y4AwIkTziVoRESewiSJaDRoN29260CS5Btiulqtvdm5y/mTk5Ott1tbna8ModfrERISgqKiIrS1tTl9PBHRUGOSRDQadFhGknztNtXqzaVJnNgrCQBiYmKst69du+bUsRaWUiVnz54d0PED1dHageYbBtSXN8BwswXSOGwu1iKiYYxJEtFo4MRIktbftJdSa71zi7dDQ28VxM3Pz3fqWIugoCD4+voiOzsbnZ0urXDRp/aWDlz6shR/WbMbdZxA9gAAF45JREFU7//wH9i59lPUlfMqOyKyj0kS0WhgWZPk5cBIknnX7dZ656pNKJVKqFSmaz0qKiow0JJGS5YsATDwtU3O6jB04Mh/HYOx07QGy1DbgsP/759obRhwMXkiukMwSSIaDZwZSdJbSpM4N5IEACkpKdbbdXV1/bTsW1hYGBQKBY4fPw6jE4vHB6q9uR3GrtsTutorddakiYioL0ySiEYDJ9YkqXX+AIDmml7VD+waN26c9XZpaanTx1vcc889AIALFy4MuA9HqXVqqH29bntsXGo4VFqvPo4gIjJhkkQ0GjgxkiQUph/7myXOl3gaM+ZWbbjBLL62JFuHDx8e8LSdo7R6De7btBSB4/whFALR6VG464lZUPswSSKi/nEzSaLRwIk1SRb115wfxVEoFNDpdGhqakJbWxs6Ozut65ScIYTA4sWLcfDgQZSUlGDChAn2DxogpUqJkNhgfOs3SyEloFQrofFV2z+QiO54HEkiGg2sI0mOJUlKtQatdc4t3LboXuC2vLx8wCNBEydOBADs37/f7aNJAOAd4A2fQG8mSETkMCZJRKNBy03Td5XWoeb6iIkDfqmoqCjr7b179w54OwCFQoF58+YBMCVbRETDDZMkotGgyZxkCMeKtQfEJAz4pQICAqy3FQoFurq6BtxXQoIpjv379w+4DyIid/FIkiSECBJC7BdCXDB/D7TRZpwQIksIkS+EyBVC/MQTsRKNCI0VTjXXR5jWAHW2Gpx+KSEEgoNN2whIKaHRaJzuw0KpVCItLQ3t7e2orq4ecD9ERO7gqZGk5wAckFLGAThgvt9TJ4BnpZSJAOYA+A8hRNIQxkg0cjRfd6r5rV23nStNYmGp4yalhFbr2BSfvb6ysrIG1Q8Rkat5Kkl6AMAO8+0dAFb0bCClrJBSnjbfbgSQDyByyCIkGkkMzi3C1vpbdt0eWJIUGWn6UUxISEBYWNiA+rBQq9WYMmUK6urqBrxBJRGRO3gqSQqVUlYApmQIwNj+GgshYgCkAjju9siIXKSlvhW1pXWoKqqB4WaLe1/M4FyyM9gkSacz7cdUWloKtXrwV4vNmDEDAHD06NFB90VE5Cpu2ydJCPE5AFsfMTc62Y8OwE4AT0sp+6xKKYRYA2ANAIwfP96ZlyByuZb6Vhz8v1/i2hnTWiGfIG98++Xl0IU4vo+RU4wdTu2RZE2SGgaWJAGmn7PS0lIYjUYoFIP7vOXt7Y2JEyfi0qVLaGpqsiZhRESe5LaRJCnlUillso2vjwBcF0KEA4D5e5WtPoQQXjAlSO9KKXfZeb03pJRpUsq0kJAQV58OkVPqrtVbEyTAVFT1zN/Oo6t94FeC2eUzxn4bM41/EACg5ebAF0vHx8cDgMsWXM+ZMwcAcOzYMZf0R0Q0WJ6abvsHgNXm26sBfNSzgRBCANgOIF9K+eoQxkY0aPUVTTYea0RnhxuTJF/HPxwoVaYpsvrSogG/XEREBIDB1XDrTqfTITw8HJcuXUJra6tL+iQiGgxPJUkvAVgmhLgAYJn5PoQQEUKIT81t5gF4DMBiIcRZ89e9ngmXyDmR00KhUN6+Z1HCslj37vbs6/wC6rpBJEmWq9qys7MH3EdP8+fPBwCcPn3aZX0SEQ2UR2q3SSlvAFhi4/FyAPeab38BwLGd8YiGGa2/Ft/afA+OvXUKbU3tmHLfZERND3fvi/o5339jRcmgXjI+Ph5FRUUDruHWU2BgIPz9/XH+/HmkpaW5ZFE4EdFAccdtIjfw0qgQlhCC5T9fhG9tXoakb8RB6zfwTRf7ZTRP4emcG0nyHRtlv5EdlvprlZWVg+7LYtGiRQCA3Nxcl/VJRDQQTJKI3Eir18AnwBsKpRt/1DrMu2Zr9E4dFjB+8qBfOjzcNHpVUlIy6L4sxo4dC7Vaja+//npQJU+IiAaLSRLRSNfebPqudu6y+YBoU5Jk7OwY8Et7eXkBAPLy8gbchy3Lli0DABQWFrq0XyIiZzBJIhrpOixJknN7MPkEmvZwbW24OaiXnzZtGgCgra1tUP10Z7ly7osvvoDRaHRZv0REzmCSRDTStZu3G3ByJMmyV1LbIDaUBIDo6GgAQEWFc0V2+yOEsI4mXbp0yWX9EhE5g0kS0UhnmW5zYsdtoFuR27rBJUljx5pGpIqLiwfVT08xMTEAgIMHD0JK6dK+iYgcwSSJaKQb4EiS1jySNJjSJACgVCqh0WhcPuIjhMC//Mu/AACuXbvm0r6JiBzBJIlopLOzJkkaJQx1LTDcbIGx89b6HutI0gCL3HaXkpICAGhubh50X93FxcUBAD777DOX9ktE5AgmSUQjXT8jSe2Gdlz5+hr+8fxn2Ln2U+R8nI/WBtMCa5XWBwDQWH550CFERZn2XCovLx90X90plUrMnj0bXV1duH79ukv7JiKyh0kS0UjXz5qk5toW7Nt8GPXljTDUtuD4H8+g/Lwp2TCVRwRuXhn8ZfbBwcEAgIKCgkH31VNSUhIA4MCBAy7vm4ioP0ySiEYyoxGoM48Eid5VfEpPlvV6rCjrEjpaO633G8oGv5ZICAF/f39UVFS4fJG1l5cXpk2bhqamJtTW1rq0byKi/jBJIhrJ8ncCX20x3X5rfq+ng2ICez02ZmIglF6mH32NXyA6W12zjig5ORkA0NDQ4JL+uktNTQUAHD582OV9ExH1hUkS0Ug2diqg0gBCAUTM6vX0mAmBGJ8Wab0fMM4fScvjrWVS9FGTXBaKZV1SWVnv0avB0mg0iI+PR3V1tVuSMCIiWwZftpuIPGfMZFPNtq52YM7TvZ729tdi0U/uQruhHV2dRmh0avgEeFufD4xJRHX+SUgprWuUBkqvN9WOO3/+vHUdkSvNmDEDRUVFeP/995GamopZs3onhURErsSRJKKRTAggMRMIiAEibScNWr0G+jA/BEb535YgAYBurGn0p7253gWhCISFhaGurs4tmz92L6LrjtEqIqKemCQRjXT3/BZ44mubC7ft0fibrkprq3fNgujExEQAwI0bg997qaeEhAQEBgZCCAGDweDy/omIemKSRDTSeXkD3gEDOlRrTpJa62tcEoqlMK07dshWQIml85fBR+2DlpYWtDa6rqAuEZEtTJKI7mBavbk0iYtGknx9TXs1ZWdnu6S/7m5cqsXO//Upmj42ovOUAmc+PI+2JiZKROQ+TJKI3KS1oQ2tDa2eDqNf2gDXlSaxiImJQVtbG7q6ulzWZ0tdK4784Ti6OoxAhwDKvJDzUT46WjrtH0xENEC8uo3Ixdqa23G9oBqn3j8HoQBmPTodIXHBUHt7eTq0XtQ6fwCA4UaFy/qMi4tDSUkJqqurERYW5pI+pZRovtF7HVL3TTGJiFyNI0lELlZf0Yg9m7JQVVSD6wU1+OQXn6Op2rWFX11FoVACAG6W5LusT8u6pCtXrrisT7WPF+IWTLjtMd9gH2h0ape9BhFRT0ySiFxIGiXy917o8SBQdHDwpT/cqf5qscv60mg0AFy7LkmlUWHGd6dixsqpCBjnjwlzx+H+zffAO0DrstcgIuqJ021ELiQUAgGRfr0e94/o/dhwofBSo+VmlUv7TExMRH5+Pjo6OuDl5ZppRm9/LWY8NBVTMuKhUiuh9uEoEhG5F0eSiFwsdsFE6MN01vsB4/wRPSvKgxH1Tx8xwX4jJ02YYOqzsrLSpf0qVQr4BHgzQSKiIcGRJCIX8w3yxgMvfQMNlY0QCgG/UF2vna6Hk4DoBNRdKXRpn6GhoQCA4gsXMW7cOJf2TUQ0VJgkEbmBT6A3fAKHb2LUnT5yIgCgs60FKs3gY25tbMPV0+WABC4UFyE1PhV+YTooVcpB901ENJQ43UZ0h/P2t+yVNPgNJaWUuHq6HAdf/RKy2PQZbOdzn6K1nps+EtHIwySJ6A6nDXBdaZK2pnbkflpkulNlGjnq9G1DdbHra7kREbkbkySiO5xGby5y2zD4kSSFSgGfQPNl+fXmXy+RXbxUn4hGJCZJRHc4a5HbusGPJKm9vZD+WCq8vFWAFJBtgAjrgj5s+G6BQETUFy7cJrrDWYrctrggSQIAv1BfPPxf96O6+AbKGq8i99J5dClZPoSIRh6OJBHd4eqvmXbbvnRoF+pKiwbdn1KlhG+wD2Jmj0PijAQAQHl5+aD7JSIaakySiO5wZ/70fwChQHN1Oa59fcClfQcGBgIA8vLyXNovEdFQYJJEdIebfO8qqDRaqLTeCJs616V9CyEQGBiIqqoqSCld2jcRkbsxSSK6w0XMXASV1hfGzg4ETZzSb9uuji60NLSis6PL4f6Tk5MBAPX19YOKk4hoqHHhNtEdTqFQIuV7a1FTeAYKVd/FaA11LTj/cSGuZVcgIjkU0x5IdGhX8cjISADAtWvXEBAQ4LK4iYjcjUkSEWHighWYuGBFn8+3Nrbh8O+OofTrMgBA9YUbqL54A8t+eje0+v73QNLr9QCAc+fOWUeViIhGAk63EZFdnW1d1gTJojznOjrbHJt2i4iIQGNjI4xGozvCIyJyCyZJRGSXEICX9+1TcSq1EkIpHDo+IcG0FcCNGyxPQkQjh0eSJCFEkBBivxDigvl7oI02WiHECSFEthAiVwjxa0/ESkSAxk+NOf+Wettjsx6bDrWP2qHjIyIiAABXr151eWxERO7iqTVJzwE4IKV8SQjxnPn+hh5t2gAsllI2CSG8AHwhhNgjpTw21MES3elUahUmzY9B5NQw1FyqRXBMILwDveGldexXiI+PDwDg9OnTiI+Ph06nc2e4REQu4anptgcA7DDf3gGg14pRadJkvutl/uJGK0QeovFVwz9Cj0nzYxAQ5Q+Nr2OjSMCtaTaj0Yjdu3e7K0QiIpfyVJIUKqWsAADz97G2GgkhlEKIswCqAOyXUh7vq0MhxBohxEkhxMnq6mq3BE1EA6NUKqFUKgEACgWXQhLRyOC231ZCiM+FEOdtfD3gaB9Syi4p5XQAUQDShRB9Xj8spXxDSpkmpUwLCQlxxSkQkYsEBARg5syZEEJAq+1/ywAiouHCbWuSpJRL+3pOCHFdCBEupawQQoTDNFLUX191QohDAJYDOO/aSIloKEybNg0XLlywrk8iIhruPDXu/Q8Aq823VwP4qGcDIUSIECLAfNsbwFIABUMWIRG5lEKhwLe+9S0sXdrn5yciomHFU0nSSwCWCSEuAFhmvg8hRIQQ4lNzm3AAWUKIHABfw7Qm6ROPREtELqHVaqFScaN/IhoZPPLbSkp5A8ASG4+XA7jXfDsHQGrPNkRERERDgZeZEBEREdnAJImIiIjIBiZJRERERDYwSSIiIiKygUkSERERkQ1MkoiIiIhsYJJEREREZAOTJCIiIiIbmCQRERER2cAkiYiIiMgGJklERERENjBJIiIiIrKBSRIRERGRDUySiIiIiGxgkkRERERkA5MkIiIiIhuYJBERERHZwCSJiIiIyAYmSUREREQ2MEkiIiIisoFJEhEREZENKk8HQEQjQ21zG8puGFDb3I6kSH8E+qqhUvJzFhGNXkySiMium83t+Nlfz+LslToAgI9aibd/MBfRY3w9HBkRkfvwYyAR2XWt1mBNkADA0N6FbQcuoLmt04NRERG5F5MkIrLrZnNbr8dqGtvQ2WX0QDREREODSRIR2ZUY4Q8ftfK2x74zezz8fdQeioiIyP24JomI7AryVeOPP5iLrQcuoKaxDd+ZPR6zY4M9HRYRkVsxSSIiu5RKBcaP8cXPVySjs8sIPUeQiOgOwCSJiBzmo+GvDCK6c3BNEhEREZENTJKIiIiIbGCSRERERGQDkyQiIiIiG5gkEREREdnAJImIiIjIBiZJRERERDYwSSIiIiKygUkSERERkQ0eSZKEEEFCiP1CiAvm74H9tFUKIc4IIT4ZyhiJiIjozuapkaTnAByQUsYBOGC+35efAMgfkqiIiIiIzDyVJD0AYIf59g4AK2w1EkJEAfgmgDeHKC4iIiIiAJ5LkkKllBUAYP4+to92rwFYD8Bor0MhxBohxEkhxMnq6mrXRUpERER3JLeV9BZCfA4gzMZTGx08/j4AVVLKU0KIhfbaSynfAPAGAKSlpUknQiUiIiLqxW1JkpRyaV/PCSGuCyHCpZQVQohwAFU2ms0DcL8Q4l4AWgB6IcSfpZTfd1PIRERERFZCyqEfdBFC/BbADSnlS0KI5wAESSnX99N+IYB1Usr7HOy/GsAVlwQ7dMYAqPF0EC7A8xg+RsM5ADyPoRItpQzxdBBEw4nbRpLseAnAB0KIfwdQCuAhABBCRAB4U0p572A6H4k/6EKIk1LKNE/HMVg8j+FjNJwDwPMgIs/xSJIkpbwBYImNx8sB9EqQpJSHABxye2BEREREZtxxm4iIiMgGJknDxxueDsBFeB7Dx2g4B4DnQUQe4pGF20RERETDHUeSiIiIiGxgkuQhjhT5FUJohRAnhBDZQohcIcSvPRFrfxw8j3FCiCwhRL75PH7iiVj74mjBZSHEW0KIKiHE+aGOsT9CiOVCiEIhRLF5S42ezwshxOvm53OEEDM8Eac9DpxHghDin0KINiHEOk/EaI8D5/Co+d8gRwjxlRAixRNxEpFjmCR5jiNFftsALJZSpgCYDmC5EGLOEMboCEfOoxPAs1LKRABzAPyHECJpCGO0x9GCy38EsHyognKEEEIJ4L8AZABIAvCIjfc2A0Cc+WsNgD8MaZAOcPA8agH8GMArQxyeQxw8h8sAFkgppwF4AVynRDSsMUnyHLtFfqVJk/mul/lruC0ic+Q8KqSUp823GwHkA4gcsgjtc6jgspTyCEx/qIeTdADFUspLUsp2AO/DdD7dPQDgHfP/p2MAAsw73Q8nds9DSlklpfwaQIcnAnSAI+fwlZTypvnuMQBRQxwjETmBSZLnOFTkVwihFEKchal0y34p5fEhjNERjhYrBgAIIWIApAIYTufh1DkMM5EArna7fw29E1BH2njaSIjRHmfP4d8B7HFrREQ0KJ7acfuOMNgivwAgpewCMF0IEQDg70KIZCnlkK6JccV5mPvRAdgJ4GkpZYMrYnPitV1yDsOQsPFYz9FGR9p42kiI0R6Hz0EIsQimJGm+WyMiokFhkuRGLijy272vOiHEIZjWxAxpkuSK8xBCeMGUIL0rpdzlplD75Mp/i2HmGoBx3e5HASgfQBtPGwkx2uPQOQghpgF4E0CGufoAEQ1TnG7znH8AWG2+vRrARz0bCCFCzCNIEEJ4A1gKoGDIInSMI+chAGwHkC+lfHUIY3OU3XMYxr4GECeEmCCEUANYCdP5dPcPAKvMV7nNAVBvmV4cRhw5j+HO7jkIIcYD2AXgMSllkQdiJCIncDNJDxFCBAP4AMB4mIv8Silruxf5NX/i3AFACVNC+4GUcpPHgrbBwfOYD+AogHMAjOZDfyal/NQjQffgyDmY270HYCFM1dyvA/illHK7Z6K+RQhxL4DXYPp/8paU8kUhxJMAIKXcak5SfwfTKKQBwL9JKU96LOA+OHAeYQBOAtDD9P+oCUDSUE/d9seBc3gTQCaAK+ZDOln0lmj4YpJEREREZAOn24iIiIhsYJJEREREZAOTJCIiIiIbmCQRERER2cAkiYiIiMgGJklELiKE6BJCnBVCnBdCfCiE8DE/HiaEeF8IcVEIkSeE+FQIEW9+bq8Qok4I8Ylnoyciop6YJBG5TouUcrqUMhlAO4AnzXsU/R3AISnlJCllEoCfAQg1H/NbAI95JlwiIuoPkyQi9zgKIBbAIgAdUsqtlieklGellEfNtw8AaPRMiERE1B8mSUQuJoRQAciAaYfxZACnPBsRERENBJMkItfxFkKchal0RilM9eqIiGiEUnk6AKJRpEVKOb37A0KIXADf8VA8REQ0CBxJInKvgwA0QognLA8IIWYJIRZ4MCYiInIAkyQiN5KmCtLfBrDMvAVALoBfASgHACHEUQAfAlgihLgmhPiGx4IlIqLbCNPvcCIiIiLqjiNJRERERDYwSSIiIiKygUkSERERkQ1MkoiIiIhsYJJEREREZAOTJCIiIiIbmCQRERER2cAkiYiIiMiG/w9WlIElaEl0XQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(6, 6))\n",
    "\n",
    "# plot the biplot\n",
    "ax = biplot('PC1', 'PC2', spca_df,\n",
    "            fpca_df, 'BodySite', ax)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Biplots are exploratory visualization tools that allow us to represent the features (i.e. taxonomy or OTUs)  that strongly influence the principal component axis as arrows. The interpretation of the compositional biplot differs slightly from classical biplot interpretation. The important features with regard to sample clusters are not a single arrow but by the log ratio between features represented by arrows pointing in different directions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this visualization we noticed that BodySite seems to explain the clusters well. We can run [PERMANOVA](https://docs.qiime2.org/2019.1/plugins/available/diversity/beta-group-significance/) on the distances to get a statistical significance for this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "method name               PERMANOVA\n",
       "test statistic name        pseudo-F\n",
       "sample size                      34\n",
       "number of groups                  4\n",
       "test statistic              1.67627\n",
       "p-value                       0.057\n",
       "number of permutations          999\n",
       "Name: PERMANOVA results, dtype: object"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from skbio.stats.distance import permanova\n",
    "\n",
    "permanova(distance, mf['BodySite'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed we can now see that the clusters we saw in the biplot are significant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the biplot and PERMANOVA results that gut is very different from the skin samples. Next we can use [qurro](https://github.com/biocore/qurro) to explore log-ratios of the microbes highlighted by gemelli. For more about why log-ratios are useful you may want to read [\"Establishing microbial composition measurement standards with reference frames\"](https://www.nature.com/articles/s41467-019-10656-5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "!qurro\\\n",
    "    --ranks qiime2-moving-pictures-tutorial/standalone-cli/ordination.txt\\\n",
    "    --table qiime2-moving-pictures-tutorial/table.biom\\\n",
    "    --sample-metadata qiime2-moving-pictures-tutorial/sample-metadata.tsv\\\n",
    "    --feature-metadata qiime2-moving-pictures-tutorial/taxonomy.tsv\\\n",
    "    --output-dir qiime2-moving-pictures-tutorial/qurro-standalone"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two taxa groups whose arrows seem to be directly opposed with relation to the BodySite grouping is Bacteroides (associated with gut) and Streptococcus (associated with skin and oral). We can use Qurro to explore this relationship. To make a log-ratio we can filter by taxa who contain Bacteroides in the numerator and Streptococcus in the denominator of the log-ratio. Those features will then be summed according to thier taxonomic labels and used in the log-ratio. In Qurro the axis one loadings (or another axis) from gemelli are highlighted by if they are contained in the numerator or denominator. The log-ratio plot is contained on the left and can be visualized as a scatter or box-plot. From this it is clear these taxa can separate our BodySite groupings. \n",
    "\n",
    "![](etc/img16.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A tsv file with the log-ratios can be exported and a t-test by BodySite on the log-ratios could confirm this observation by using the `Export current sample plot data` button. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "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></th>\n",
       "      <th>test-stat</th>\n",
       "      <th>p-value</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>group one vs.</th>\n",
       "      <th>group two</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">gut</th>\n",
       "      <th>left palm</th>\n",
       "      <td>9.943842</td>\n",
       "      <td>3.749890e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>right palm</th>\n",
       "      <td>3.385889</td>\n",
       "      <td>6.933597e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tongue</th>\n",
       "      <td>20.744948</td>\n",
       "      <td>6.574794e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">left palm</th>\n",
       "      <th>right palm</th>\n",
       "      <td>-1.181070</td>\n",
       "      <td>2.587322e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tongue</th>\n",
       "      <td>1.526182</td>\n",
       "      <td>1.528810e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>right palm</th>\n",
       "      <th>tongue</th>\n",
       "      <td>1.893318</td>\n",
       "      <td>8.078723e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          test-stat       p-value\n",
       "group one vs. group two                          \n",
       "gut           left palm    9.943842  3.749890e-06\n",
       "              right palm   3.385889  6.933597e-03\n",
       "              tongue      20.744948  6.574794e-09\n",
       "left palm     right palm  -1.181070  2.587322e-01\n",
       "              tongue       1.526182  1.528810e-01\n",
       "right palm    tongue       1.893318  8.078723e-02"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import itertools\n",
    "from scipy.stats import ttest_ind\n",
    "\n",
    "# import the taxonomy metadata\n",
    "lrdf = pd.read_csv('qiime2-moving-pictures-tutorial/sample_plot_data.tsv',\n",
    "                   sep='\\t', index_col=0).dropna(subset=['Current_Natural_Log_Ratio'])\n",
    "# split data by BodySite\n",
    "lrs = {type_:df_.Current_Natural_Log_Ratio\n",
    "       for type_, df_ in lrdf.groupby('BodySite')}\n",
    "# get all combos\n",
    "ids_ = list(itertools.combinations(lrs.keys(), 2))\n",
    "# take t-test\n",
    "tst = pd.DataFrame({(id1_, id2_):ttest_ind(lrs[id1_], lrs[id2_])\n",
    "                    for id1_, id2_ in ids_},\n",
    "                   ['test-stat','p-value']).T\n",
    "tst.index.names = ['group one vs.', 'group two']\n",
    "# view results\n",
    "tst\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this we can see that gut vs. skin and oral is significantly different using this log-ratio but skin vs. oral is not."
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}