validation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/skbio/util/_testing.py:15: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  import pandas.util.testing as pdt\n"
     ]
    }
   ],
   "source": [
    "from qiime2 import Artifact\n",
    "from qiime2.plugins.fragment_insertion.methods import filter_features\n",
    "from qiime2.plugins.emperor.actions import plot, biplot\n",
    "from qiime2.plugins.deicode.actions import rpca\n",
    "from qiime2.plugins.diversity.actions import core_metrics_phylogenetic\n",
    "from qiime2.plugins.fragment_insertion.methods import filter_features\n",
    "from qiime2.plugins.diversity.visualizers import beta_group_significance\n",
    "from qiime2.plugins.diversity.visualizers import alpha_group_significance\n",
    "from qiime2.plugins.metadata.visualizers import tabulate\n",
    "from qiime2.plugins.diversity.methods import alpha\n",
    "from qiime2.plugins import taxa\n",
    "from qiime2 import Artifact\n",
    "import qiime2 as q2\n",
    "import biom\n",
    "import biom\n",
    "import os\n",
    "import pandas as pd\n",
    "import glob\n",
    "from qiime2.plugins import (feature_table,\n",
    "fragment_insertion, \n",
    "feature_classifier, \n",
    "diversity, \n",
    "metadata, \n",
    "emperor, \n",
    "deicode, \n",
    "gneiss, \n",
    "taxa)\n",
    "import pandas as pd\n",
    "\n",
    "from qiime2 import Artifact\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "import skbio\n",
    "from skbio.stats.distance import MissingIDError\n",
    "from skbio import DistanceMatrix\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import csv\n",
    "\n",
    "import qiime2\n",
    "from matplotlib.legend import Legend\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_and_create_dir(dir_name):\n",
    "    CHECK_DIRECTORY = os.path.isdir(dir_name)\n",
    "    if not CHECK_DIRECTORY:\n",
    "        os.makedirs(dir_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# import data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove blanks and katharo samples from biom table and metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import metadata\n",
    "mf = q2.Metadata.load('14385_20220208-102803.txt').to_dataframe()\n",
    "\n",
    "# import table\n",
    "bt = biom.load_table('141987_reference-hit.biom')\n",
    "mf = mf.loc[bt.ids(), : ]\n",
    "\n",
    "# directory for outputs (this will include things like figures and emperor plots)\n",
    "inputs_dir_name = 'inputs'\n",
    "check_and_create_dir(inputs_dir_name)\n",
    "output_dir_name = 'plots-and-figures'\n",
    "check_and_create_dir(output_dir_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['BLANK5.12D', 'BLANK.AssayAssure.20C.BLANK.1', 'BLANK5.12G',\n",
       "       'BLANK5.12H', 'Katharo5.2speciesmock.1E.5thousand.prev',\n",
       "       'BLANK.AssayAssure.20C.BLANK.2', 'BLANK.EtOH.20C.BLANK.2',\n",
       "       'Katharo5.2speciesmock.1F.5hundred.new', 'BLANK.None.20C.BLANK.2',\n",
       "       'BLANK.None.20C.BLANK.1', 'Katharo5.2speciesmock.1B.5million.prev',\n",
       "       'Katharo5.2speciesmock.1F.5hundred.prev', 'BLANK5.12C',\n",
       "       'Katharo5.2speciesmock.1G.50.new',\n",
       "       'Katharo5.2speciesmock.1H.5.prev',\n",
       "       'Katharo5.2speciesmock.1G.50.prev', 'Katharo5.Zymo.2F.5hundred',\n",
       "       'BLANK.EtOH.20C.BLANK.1', 'Katharo5.Zymo.2G.50',\n",
       "       'Katharo5.Zymo.2H.5', 'Katharo5.2speciesmock.1D.50thousand.new',\n",
       "       'Katharo5.2speciesmock.1D.50thousand.prev', 'K',\n",
       "       'Katharo5.2speciesmock.1E.5thousand.new', 'BLANK5.12E',\n",
       "       'Katharo5.2speciesmock.1C.500thousand.new',\n",
       "       'Katharo5.2speciesmock.1C.500thousand.prev',\n",
       "       'Katharo5.Zymo.2C.500thousand', 'BLANK5.12F',\n",
       "       'Katharo5.2speciesmock.1A.50million.new', 'L',\n",
       "       'Katharo5.Zymo.2D.50thousand', 'Katharo5.Zymo.2E.5thousand', 'N',\n",
       "       'M', 'Katharo5.2speciesmock.1A.50million.prev',\n",
       "       'Katharo5.Zymo.2A.50million', 'O', 'Katharo5.Zymo.2B.5million',\n",
       "       'Katharo5.2speciesmock.1B.5million.new', 'Q', 'P', 'S', 'R', 'T'],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf.host_subject_id.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# criteria to keep smaple  \n",
    "keep = lambda sample: 'BLANK' not in sample and 'Katharo' not in sample\n",
    "\n",
    "# identify the metadata data to keep (need to remove blank and kartho)\n",
    "keep = [ind for ind in mf.index if keep(ind)]\n",
    "\n",
    "\n",
    "#filter metadata\n",
    "mf_filtered = mf.loc[keep, : ]\n",
    "\n",
    "#filter biom table\n",
    "bt_filter = bt.filter(mf_filtered.index, axis='sample')\n",
    "\n",
    "#some features may have been removed that were unique to the samples removed\n",
    "#now there may be a lot of 0's in our table, remove them and match tables\n",
    "keep_feat = bt_filter.ids('observation')[bt_filter.sum(axis='observation') > 0]\n",
    "bt_filter = bt_filter.filter(keep_feat, axis='observation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['K', 'L', 'N', 'M', 'O', 'Q', 'P', 'S', 'R', 'T'], dtype=object)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_filtered.host_subject_id.unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### import table/tree/metadata into qiime2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/skbio/util/_testing.py:15: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  import pandas.util.testing as pdt\n",
      "\u001b[32mImported 141987_insertion_tree.relabelled.tre as NewickDirectoryFormat to rooted-tree.qza\u001b[0m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'filtered-table.qza'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#import tree\n",
    "!qiime tools import \\\n",
    "  --input-path 141987_insertion_tree.relabelled.tre \\\n",
    "  --output-path rooted-tree.qza \\\n",
    "  --type 'Phylogeny[Rooted]'\n",
    "\n",
    "#turn tables to q2 artifacts\n",
    "q2mf_filtered = q2.Metadata(mf_filtered)\n",
    "q2mf_filtered.save('filtered-metadata.qza')\n",
    "q2bt = q2.Artifact.import_data('FeatureTable[Frequency]', bt_filter)\n",
    "q2bt.save('filtered-table.qza')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### filter table to match tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #import filtered table and tree\n",
    "# rooted_tree = q2.Artifact.load('rooted-tree.qza')\n",
    "\n",
    "# #filter table to match tree\n",
    "# tree_filtered_table, removed_table = filter_features(q2bt,rooted_tree)\n",
    "# tree_filtered_table.save('tree-filtered-table.qza')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TEMPORARY!!!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rooted_tree = q2.Artifact.load('rooted-tree.qza')\n",
    "\n",
    "tree_filtered_table = q2.Artifact.load('tree-filtered-table.qza')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-onv7ajfj')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 846262f3-0e13-4e20-b4ab-6c0b1bdcc612>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAW8ElEQVR4nO3dfZQldX3n8feHURARxijo0QEy6BCVXY3BVqLrJpCNCuqAskpETdSwjHrEh5PVI6z4kPXkGE80JxoVnBiDyQoEWEMYxeBDRKJBnAFFeZB1ZFUmuAGNDoJRBL/7x60prm13T3XPVFff7vfrnHu66ndvVX3unZ777apfVf1SVUiSBLDH0AEkSUuHRUGS1LIoSJJaFgVJUsuiIElqWRQkSa17DR1gV+y///61du3aoWNI0kS58sorv1tVB8z03EQXhbVr17Jly5ahY0jSREnyrdme8/CRJKllUZAktSwKkqSWRUGS1FoyRSHJo5KcmeSCJC8fOo8krUS9FoUkH0xyS5JrprUfneSGJFuTnApQVddX1cuAE4CpPnNJkmbW957CWcDR4w1JVgHvBY4BDgNOTHJY89yxwOeAT/ecS5I0g16vU6iqy5Ksndb8BGBrVd0IkORc4Djguqq6CLgoyceAs2daZ5INwAaAgw8+eMHZzr7i2wteVtKue/4RC///q/4McfHaGuCmsfltwBFJjgSOB/YCLp5t4araCGwEmJqacoQgSdqNhigKmaGtqupS4NLFjSJJGjfE2UfbgIPG5g8Ebp7PCpKsT7Jx+/btuzWYJK10QxSFzcChSQ5JsifwPOCi+aygqjZV1YbVq1f3ElCSVqq+T0k9B7gceESSbUlOqqq7gFOAS4DrgfOq6to+c0iSuun77KMTZ2m/mDk6k3cmyXpg/bp16xa6CknSDJbMFc3z4eEjSerHRBYFSVI/LAqSpNZEFgVPSZWkfkxkUbBPQZL6MZFFQZLUD4uCJKllUZAktSayKNjRLEn9mMiiYEezJPVjIouCJKkfFgVJUsuiIElqTWRRsKNZkvoxkUXBjmZJ6sdEFgVJUj8sCpKklkVBktSyKEiSWhNZFDz7SJL6MZFFwbOPJKkfE1kUJEn9sChIkloWBUlSy6IgSWpZFCRJLYuCJKllUZAktSayKHjxmiT1YyKLghevSVI/JrIoSJL6YVGQJLUsCpKklkVBktSyKEiSWhYFSVLLoiBJalkUJEkti4IkqTWRRcHbXEhSPyayKHibC0nqx0QWBUlSPywKkqSWRUGS1LIoSJJaFgVJUmunRSHJc5Ps20yfnuQjSQ7vP5okabF12VN4Y1X9MMmTgacBHwLO6DeWJGkIXYrC3c3PZwBnVNXfA3v2F0mSNJQuReFfkrwfOAG4OMleHZeTJE2YLl/uJwCXAEdX1Q+ABwCv6zWVJGkQOy0KVfUj4BbgyU3TXcDX+wwlSRpGl7OP3gy8Hjitabo38L/6DCVJGkaXw0fPBo4F7gCoqpuBffsMJUkaRpeicGdVFVAASfbpN5IkaShdisJ5zdlH909yMvAp4C/6jSVJGsK95noySYC/BR4J3AY8AnhTVX2yjzBJnsXoeogHAe+tqk/0sR1J0szmLApVVUkurKrHAQsqBEk+CDwTuKWq/uNY+9HAu4BVwAeq6o+r6kLgwiS/BLwDsChI0iLqcvjoC0kevwvbOAs4erwhySrgvcAxwGHAiUkOG3vJ6c3zkqRF1KUoHAVcnuQbSb6S5KtJvtJ1A1V1GfBv05qfAGytqhur6k7gXOC4jLwd+HhVXTXT+pJsSLIlyZZbb721awxJUgdzHj5qHNPDdtcAN43NbwOOAF4J/DawOsm6qjpz+oJVtRHYCDA1NVU9ZJOkFatLUejjizczbaeq3g28u4ftSZI66FIUPsaoMAS4D3AIcAPwH3Zhu9uAg8bmDwRu7rpwkvXA+nXr1u1CBEnSdF3uffToqnpM8/NQRv0Bn9vF7W4GDk1ySJI9gecBF3VduKo2VdWG1atX72IMSdK4ed8Cu+kA7nw2UpJzgMuBRyTZluSkqroLOIXR3VevB86rqmvnm0WStHvt9PBRkj8Ym90DOBzofNpPVZ04S/vFwMVd1zMtk4ePJKkHXfYU9h177MWoj+G4PkPtjIePJKkfXTqar6uq88cbkjwXOH+W10uSJlSXPYXTOrZJkibcrHsKSY4Bng6sSTJ+7cB+jEZfG4x9CpLUj7n2FG4GtgA/Bq4ce1wEPK3/aLOzT0GS+jHrnkJVXQ1cneTsqvrpImaSJA2kS0fz2iRvY3Q30/vsaKyqh/WWSpI0iC4dzX8FnMGoH+Eo4K+Bv+kzlCRpGF2Kwt5V9WkgVfWtqnoL8Fv9xppbkvVJNm7fvn3IGJK07HQpCj9Osgfw9SSnJHk2o+EyB2NHsyT1o0tReA1wX+BVwOOAFwIv6jOUJGkYO+1orqrNAEmqql7SfyRJ0lB2uqeQ5IlJrmN0N1OS/GqS9/WeTJK06LocPvozRherfQ/a6xd+o89QO2NHsyT1o9N4ClV107Smu3vI0pkdzZLUjy4Xr92U5ElANaOkvYrmUJIkaXnpsqfwMuAVwBpGYys/tpmXJC0zc90l9e1V9XrgqKp6wSJmkiQNZK49hacnuTeOnSBJK8ZcfQr/AHwX2CfJbUCA2vGzqvZbhHwzcjwFSerHrHsKVfW6qloNfKyq9quqfcd/LmLGmbJ59pEk9WCnHc1VddxiBJEkDa/TdQqSpJXBoiBJas1aFJJ8uvn59sWLI0ka0lxnHz0kyW8CxyY5l9FZR62quqrXZJKkRTdXUXgTcCpwIPCn054rBh59TZK0+81aFKrqAuCCJG+sqrcuYiZJ0kC6DLLz1iTHcs/tsi+tqo/2G2tuXrwmSf3oMsjO24BXA9c1j1c3bYPx4jVJ6keXW2c/A3hsVf0MIMmHgC/hPZEkadnpep3C/cem/fNckpapLnsKbwO+lOQzjE5L/Q3cS5CkZalLR/M5SS4FHs+oKLy+qv5f38EkSYuvy54CVfUd4KKes0iSBua9jyRJLYuCJKk1Z1FIskeSaxYrjCRpWHMWhebahKuTHLxIeSRJA+rS0fwQ4NokXwTu2NFYVcf2lmonvM2FJPWjS1H4w95TzFNVbQI2TU1NnTx0FklaTrpcp/DZJL8MHFpVn0pyX2BV/9EkSYutyw3xTgYuAN7fNK0BLuwzlCRpGF1OSX0F8J+A2wCq6uvAg/oMJUkaRpei8JOqunPHTJJ7MRp5TZK0zHQpCp9N8j+AvZM8BTgf2NRvLEnSELoUhVOBW4GvAi8FLgZO7zOUJGkYXc4++lkzsM4VjA4b3VBVHj6SpGVop0UhyTOAM4FvMLp19iFJXlpVH+87nCRpcXW5eO2dwFFVtRUgycOBjwEWBUlaZrr0KdyyoyA0bgRu6SmPJGlAs+4pJDm+mbw2ycXAeYz6FJ4LbF6EbJKkRTbX4aP1Y9P/CvxmM30r8Eu9JZIkDWbWolBVL1nMIJKk4XU5++gQ4JXA2vHXD3nrbElSP7qcfXQh8JeMrmL+Wb9xJElD6lIUflxV7+47SJKHAW8AVlfVc/reniTpF3U5JfVdSd6c5IlJDt/x6LLyJB9Mcsv0cZ6THJ3khiRbk5wKUFU3VtVJC3gPkqTdpMuewqOB3wV+i3sOH1UzvzNnAe8B/npHQ5JVwHuBpwDbgM1JLqqq67rHliT1oUtReDbwsPHbZ3dVVZclWTut+QnA1qq6ESDJucBxQKeikGQDsAHg4IMPnm8kSdIcuhw+uhq4/27c5hrgprH5bcCaJA9Mcibwa0lOm23hqtpYVVNVNXXAAQfsxliSpC57Cg8GvpZkM/CTHY27cEpqZmirqvoe8LIFrlOStBt0KQpv3s3b3AYcNDZ/IHDzfFaQZD2wft26dbszlySteF3GU/jsbt7mZuDQ5qK4fwGeBzx/Piuoqk3ApqmpqZN3czZJWtF22qeQ5IdJbmseP05yd5Lbuqw8yTnA5cAjkmxLclJV3QWcAlwCXA+cV1XX7sqbkCTtHl32FPYdn0/yLEZnEO1UVZ04S/vFjIb1XBAPH0lSP7qcffRzqupCul2j0Juq2lRVG1avXj1kDEladrrcEO/4sdk9gClGF69JkpaZLmcfjY+rcBfwTUYXm0mSlpkufQpLblwF+xQkqR9zDcf5pjmWq6p6aw95OvGUVEnqx1x7CnfM0LYPcBLwQGCwoiBJ6sdcw3G+c8d0kn2BVwMvAc4F3jnbcpKkyTVnn0KSBwB/ALwA+BBweFV9fzGCSZIW36zXKST5E0a3pPgh8OiqestSKQhJ1ifZuH379qGjSNKyMtfFa/8deChwOnDz2K0uftj1Nhd98eI1SerHXH0K877aWZI02fzilyS1LAqSpNZEFgU7miWpHxNZFOxolqR+TGRRkCT1w6IgSWpZFCRJLYuCJKk1kUXBs48kqR8TWRQ8+0iS+jGRRUGS1A+LgiSpZVGQJLUsCpKklkVBktSyKEiSWhYFSVJrIouCF69JUj8msih48Zok9WMii4IkqR8WBUlSy6IgSWpZFCRJLYuCJKllUZAktSwKkqSWRUGS1LIoSJJaE1kUvM2FJPVjIouCt7mQpH5MZFGQJPXDoiBJalkUJEkti4IkqWVRkCS1LAqSpJZFQZLUsihIkloWBUlSy6IgSWpZFCRJLYuCJKllUZAktSwKkqSWRUGS1LrX0AF2SLIP8D7gTuDSqvrwwJEkacXpdU8hyQeT3JLkmmntRye5IcnWJKc2zccDF1TVycCxfeaSJM2s78NHZwFHjzckWQW8FzgGOAw4MclhwIHATc3L7u45lyRpBr0ePqqqy5Ksndb8BGBrVd0IkORc4DhgG6PC8GXmKFZJNgAbAA4++ODdH1rSojj7im8vaLnnH7Gw//c729581rtjXQvNsqvOvuLbvW17iI7mNdyzRwCjYrAG+AjwX5OcAWyabeGq2lhVU1U1dcABB/SbVJJWmCE6mjNDW1XVHcBLFjuMJOkeQ+wpbAMOGps/ELh5PitIsj7Jxu3bt+/WYJK00g1RFDYDhyY5JMmewPOAi+azgqraVFUbVq9e3UtASVqp+j4l9RzgcuARSbYlOamq7gJOAS4BrgfOq6pr+8whSeqm77OPTpyl/WLg4j63LUmav4m8zYV9CpLUj4ksCvYpSFI/JrIoSJL6kaoaOsOCJbkV+AEwfhxp9dj8TNM7fu4PfHcBmx1f53xfM719rvk+ss+VbWfPm32kr9+blZC9S97Zps3e/fku2e9fVTNf/VtVE/0ANs42P9P02M8tu2N783nNXFkXI3uX/GafPXufvzcrIXuXvHO8D7P3kH2mx3I4fDT9lhibdjI96y00Fri9+bxmrqzT5/vI3mUdZv/FabMv7PmZ2neWd67phTD7L7bNud6JPny0K5JsqaqpoXMshNmHM8n5zT6MScu+HPYUFmrj0AF2gdmHM8n5zT6Micq+YvcUJEm/aCXvKUiSprEoSJJaFgVJUsui0EiyT5IPJfmLJC8YOs98JHlYkr9McsHQWeYrybOaz/zvkzx16DzzkeRRSc5MckGSlw+dZ76a3/krkzxz6CzzkeTIJP/UfPZHDp1nPpLskeSPkvx5khcNnWcmy7ooJPlgkluSXDOt/egkNyTZmuTUpvl44IKqOhk4dtHDTjOf7FV1Y1WdNEzSXzTP7Bc2n/mLgd8ZIO7PmWf266vqZcAJwOCnHM7z9x3g9cB5i5tyZvPMXsDtwH0YDdo1qHlmP47R8MM/ZQlkn9FCrrSblAfwG8DhwDVjbauAbwAPA/YErgYOA04DHtu85uxJyj72/AVD596F7O8EDp+07Iz+gPhn4PmTlB34bUYDXL0YeOaEZd+jef7BwIcnLPupwEub1yyJ/6/TH8t6T6GqLgP+bVrzE4CtNfrr+k7gXEbVexujoUFhCexBzTP7kjKf7Bl5O/DxqrpqsbNON9/PvaouqqonAYMfcpxn9qOAXweeD5ycZNDf+flkr6qfNc9/H9hrEWPOaAHfM99vXnP34qXsrtdBdpaoNcBNY/PbgCOAdwPvSfIMds+tAfowY/YkDwT+CPi1JKdV1dsGSTe32T73VzL6q3V1knVVdeYQ4XZits/9SEaHHfdi6Q4aNWP2qjoFIMmLge+OfdEuJbN97scDTwPuD7xniGAdzPb7/i7gz5P8Z+CyIYLtzEosCpmhrarqDuAlix1mnmbL/j3gZYsdZp5my/5uRgV5KZst+6XApYsbZd5mzN5OVJ21eFHmbbbP/SPARxY7zDzNlv1HwJLp/5vJ4IdJBrANOGhs/kDg5oGyzJfZh2H2YZh9ACuxKGwGDk1ySJI9GXW2XTRwpq7MPgyzD8PsQxi6p7vPB3AO8B3uOf3rpKb96cD/YXR2wBuGzmn2pfMwu9lXUvaZHt4QT5LUWomHjyRJs7AoSJJaFgVJUsuiIElqWRQkSS2LgiSpZVHQoJLcneTLY4+1Q2daipI8svl8vpTk4QtY/jVJ7ttHNi0vXqegQSW5varuN8tzYfQ7uhRv1tarJKuq6u6x+VOBvavqzQtc3zeBqar67jyWuVdV3bWQ7WlyuaegJSXJ2iTXJ3kfcBVwUJLXJdmc5CtJ/nDstW9oBjH5VJJzkry2ab80yVQzvX/zhUiSVUn+ZGxdL23aj2yWuSDJ15J8uClIJHl8kn9OcnWSLybZN6NRvx47luPzSR4z7X28OKPR5P6hyfjmsede2Kzry0nen2RV0357kv+Z5ArgiWOvfzrwGuC/JfnMTtZxRpItSa7d8VkleRXwUOAzY8vfPrb+5yQ5q5k+K8mfNq97e5KHN+/hyuZ9P3Lh/7qaCENfUu1jZT8Y3VP+y83j74C1wM+AX2+efyqwkdFdJ/cAPspoUJPHAV8F7gvsB2wFXtsscymjv4oB9ge+2UxvAE5vpvcCtgCHAEcC2xndtGwP4HLgyYwGR7kReHyzzH6M7iz8IuDPmrZfAbbM8L5ezOjWBw8E9gauYTQ626MY3Zr93s3r3gf8XjNdwAmzfE5vGXt/c63jAc3PVc3n8Jhm/pvA/mPru31s+jnAWc30Wc1nvKqZ/zRwaDN9BPCPQ//O+Oj3sRJvna2l5d+ravyv7rXAt6rqC03TU5vHl5r5+wGHAvsCf1ejWxGTpMvNxp4KPCbJc5r51c267gS+WFXbmnV9mVFx2g58p6o2A1TVbc3z5wNvTPI64PcZfZHO5JM1uq05ST7CqNDcxaigbW52RvYGbmlefzfwvzu8j/8yxzpOSLKBUfF6CKPRvr7SYZ3jzq+qu5PcD3gScH6zHVgCg9qoXxYFLUV3jE0HeFtVvX/8BUlew9i4ANPcxT2HRu8zbV2vrKpLpq3rSOAnY013M/q/kZm2UVU/SvJJRiNpzTU+8/Rlq1nnh6rqtBle/+Ma60eYw4zrSHII8FpGezbfbw4J3WeG5adnm/6aHZ//HsAPxou2lj/7FLTUXQL8fvNXK0nWJHkQo1Grnp1k7yT7AuvHlvkmo7+kYXRoZHxdL09y72Zdv5Jknzm2/TXgoUke37x+3yQ7/pD6AKPBgTZX1fShGHd4SpIHJNkbeBbweUaHY57TvAea53955x/Dz5ltHfsx+kLfnuTBwDFjy/yQ0d7VDv+a5FEZDcP57Jk20uwZ/d8kz222kyS/Os+smjDuKWhJq6pPJHkUcHlzCON24IVVdVWSv2XUF/Et4J/GFnsHcF6S3wX+caz9A4wOC13VdCTfyujLerZt35nkdxgNn7g38O+Mhg69vaquTHIb8FdzxP8c8DfAOuDsqtoCkOR04BPNF/JPgVc076GTqrpupnVU1ReSfAm4llFfyOfHFtsIfDzJd6rqKEYDyH+U0ZCR1zA6LDeTFwBnNNu7N6Oxhq/umlWTx1NStSwkeQujL+t3LNL2HsqoI/eRNcMpsxmNfTxVzVjI0qTw8JE0T0l+D7iC0cApK+4aCi1v7ilIklruKUiSWhYFSVLLoiBJalkUJEkti4IkqWVRkCS1/j/NzG2SXqf6FAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tbl = Artifact.load('tree-filtered-table.qza')\n",
    "meta = qiime2.Metadata.load('filtered-metadata.qza')\n",
    "# Singletons out\n",
    "tbl_freq2 = tbl\n",
    "feature_table.visualizers.summarize(tbl_freq2, meta).visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "VAL_BASE_COL = 'validation_base_col'\n",
    "VAL_BASE_VAL = 'base'\n",
    "class MyDM(skbio.DistanceMatrix):\n",
    "    def _stable_order(self, ids):\n",
    "\n",
    "        id_order = sorted(self._id_index[i] for i in ids)\n",
    "        return np.array(id_order, dtype=int)\n",
    "    \n",
    "    def between(self, from_, to_, allow_overlap=False):\n",
    "\n",
    "        from_ = set(from_)\n",
    "        to_ = set(to_)\n",
    "\n",
    "        all_ids = from_ | to_\n",
    "        not_present = all_ids - set(self._id_index)\n",
    "        if not_present:\n",
    "            raise MissingIDError(\"At least one ID (e.g., '%s') was not \"\n",
    "                                 \"found.\" % not_present.pop())\n",
    "\n",
    "        overlapping = from_ & to_\n",
    "        if not allow_overlap and overlapping:\n",
    "            raise KeyError(\"At least one ID overlaps in from_ and to_ \"\n",
    "                           \"(e.g., '%s'). This constraint can removed with \"\n",
    "                           \"allow_overlap=True.\" % overlapping.pop())\n",
    "\n",
    "        return self._subset_to_dataframe(from_, to_)\n",
    "    \n",
    "    def _subset_to_dataframe(self, i_ids, j_ids):\n",
    "\n",
    "        i_indices = self._stable_order(i_ids)\n",
    "        j_indices = self._stable_order(j_ids)\n",
    "\n",
    "        j_length = len(j_indices)\n",
    "        j_labels = tuple([self.ids[j] for j in j_indices])\n",
    "\n",
    "        i = []\n",
    "        j = []\n",
    "\n",
    "        values = [np.array([])]\n",
    "        for i_idx in i_indices:\n",
    "            i.extend([self.ids[i_idx]] * j_length)\n",
    "            j.extend(j_labels)\n",
    "\n",
    "            subset = self._data[i_idx, j_indices]\n",
    "            values.append(subset)\n",
    "\n",
    "        i = pd.Series(i, name='i')\n",
    "        j = pd.Series(j, name='j')\n",
    "        values = pd.Series(np.hstack(values), name='value')\n",
    "\n",
    "        return pd.concat([i, j, values], axis=1)\n",
    "    \n",
    "def gather_between_distances(a_dm, a_md, col, id_col):\n",
    "    pleasant_dm = MyDM(a_dm)\n",
    "\n",
    "    tidy_data = []\n",
    "    \n",
    "    baseline_ids = {}\n",
    "    for hid, grp in a_md.groupby(id_col):\n",
    "        base_df = grp[grp[VAL_BASE_COL] == VAL_BASE_VAL]\n",
    "        baseline_ids[hid] = base_df.index\n",
    "\n",
    "    for (pres, hid), grp in a_md.groupby([col] + ['host_subject_id']):  \n",
    "        if hid in baseline_ids:\n",
    "            base_ids = baseline_ids[hid]\n",
    "            distances = pleasant_dm.between(base_ids, grp.index, allow_overlap=True)\n",
    "\n",
    "            distances['sample_preservation_method'] = pres\n",
    "            distances['host_subject_id'] = hid\n",
    "            distances['groupby'] = 'Distance to no preservative'\n",
    "            tidy_data.append(distances)\n",
    "\n",
    "    tidy_data = pd.concat(tidy_data)\n",
    "    \n",
    "    return tidy_data[tidy_data['i'] != tidy_data['j']]\n",
    "\n",
    "def calc_inter_group(group_condition_i, group_condition_j, dist_mat, df, id_col):\n",
    "    subjects = df[id_col].unique()\n",
    "    already_proc = []\n",
    "    pleasant_dm = MyDM(dist_mat)\n",
    "    pairs = []\n",
    "    for s in subjects:\n",
    "        pairs.append(\n",
    "             pleasant_dm.between(\n",
    "#                 df.loc[df['host_subject_id'] == s].index,\n",
    "#                 df.loc[~(df['host_subject_id'] == s) & ~(df['host_subject_id'].isin(already_proc))].index))\n",
    "                df.loc[group_condition_i].index,\n",
    "                df.loc[gropu_condition_j].index))\n",
    "        already_proc.append(s)\n",
    "    pd.concat(pairs).loc[:, 'value'].median()\n",
    "\n",
    "def fancy_plotting(df, x, xlabel, y, ylabel, hue, legend_label, col, order=None, legend_text_map=None,\n",
    "                   inter_distances=None):\n",
    "    num_col = min(len(df[col].unique()), 3)\n",
    "    g = sns.catplot(data=df, x=x, y=y, hue=hue, col=col, kind='box', col_wrap=num_col, order=order,)\n",
    "#                     palette=sns.color_palette(\"RdBu_r\", 7))\n",
    "    \n",
    "    g.set_axis_labels(xlabel, ylabel)\n",
    "    g.set_xticklabels(order)\n",
    "    g.set_titles(\"{col_name}\")\n",
    "    g._legend.set_title(legend_label)\n",
    "    \n",
    "    \n",
    "    for t in g._legend.texts:\n",
    "        t.set_text(legend_text_map[t._text])\n",
    "    \n",
    "    \n",
    "    additional_legends = []\n",
    "    additional_labels = []\n",
    "    for ax in g.axes.flat:\n",
    "        ax.set_title(ax.get_title(), fontsize='xx-large')\n",
    "        ax.set_ylabel(ax.get_ylabel(), fontsize='x-large')\n",
    "        ax.set_xlabel(ax.get_xlabel(), fontsize='x-large')\n",
    "        ax.set_xticklabels(ax.get_xticklabels(), fontsize='x-large')\n",
    "        for label, ls, value in inter_distances:\n",
    "            additional_legends.append(\n",
    "                ax.axhline(value, ls=ls, color='0.1', linewidth=2, label=label))\n",
    "            additional_labels.append(label)\n",
    "        # Create the second legend and add the artist manually.\n",
    "        new_leg = Legend(g.fig, additional_legends, additional_labels,title='Reference distance')\n",
    "    g.fig.add_artist(new_leg);\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# permanova groups and other parameters\n",
    "groups = ['host_subject_id', 'sample_preservation_method']\n",
    "use_pairwise = True\n",
    "\n",
    "def create_perma_groups(q2mat, q2meta, directory, groups, metric_used='unweighted_unifrac', use_pairwise=True):\n",
    "    for group in groups:\n",
    "        group_dir = os.path.join(directory, group)\n",
    "        check_and_create_dir(group_dir)\n",
    "        out_ = os.path.join(group_dir, metric_used)\n",
    "        group_plot = beta_group_significance(q2mat, q2meta.get_column(group), pairwise=use_pairwise)\n",
    "        group_plot.visualization.save(out_)\n",
    "        \n",
    "def create_box_plot_data(base_condition, groups, col, id_col, q2meta, dist_mat, directory, xlabel, y, legend_label,\n",
    "                         groupbycol, metric_used='rpca', fname=None):\n",
    "    # calculate inter-human\n",
    "    distances = dist_mat.data\n",
    "    total_subjects = dist_mat.shape[0]\n",
    "    values = []\n",
    "    for i in range(total_subjects):\n",
    "        for j in range(i + 1, total_subjects):\n",
    "            values.append(distances[i, j])\n",
    "    inter_distances = [('inter-human', '-', np.median(values))]\n",
    "    \n",
    "    # calculate inter-replicates\n",
    "    \n",
    "    values = []\n",
    "    hids = dist_mat.ids\n",
    "    for i in range(total_subjects):\n",
    "        for j in range(i + 1, total_subjects):\n",
    "            s_i, s_j = hids[i], hids[j]\n",
    "            if s_i[:-2] == s_j[:-2]:\n",
    "                values.append(distances[i, j])\n",
    "    inter_distances.append(('inter-replicate', '--', np.median(values)))\n",
    "    \n",
    "    md = q2meta.to_dataframe()\n",
    "    md = md.loc[dist_mat.ids, :]\n",
    "    md.loc[base_condition, VAL_BASE_COL] = VAL_BASE_VAL\n",
    "    between_distances = gather_between_distances(dist_mat, md, col, id_col)\n",
    "    out_ = os.path.join(directory, metric_used) if fname is None else os.path.join(directory,fname)\n",
    "    metric_used += '_distance'\n",
    "    fancy_plotting(between_distances, \n",
    "                   x=id_col, xlabel=xlabel, \n",
    "                   y=y, ylabel=' '.join([w for w in metric_used.capitalize().split('_')]),\n",
    "                   legend_text_map={w: ' '.join(w.split('_')) for w in between_distances[base_col].unique()},\n",
    "                   legend_label=legend_label,\n",
    "                   hue=base_col, col=groupbycol, inter_distances=inter_distances).savefig(out_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This will create emperor plots, PERMANOVA group significance plots, and taxonomy plots. Results will be stored in directory defined by `output_dir_name`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/sklearn/metrics/pairwise.py:1776: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/skbio/stats/ordination/_principal_coordinate_analysis.py:152: RuntimeWarning: The result contains negative eigenvalues. Please compare their magnitude with the magnitude of some of the largest positive eigenvalues. If the negative ones are smaller, it's probably safe to ignore them, but if they are large in magnitude, the results won't be useful. See the Notes section for more details. The smallest eigenvalue is -0.07112736748682533 and the largest is 5.948348443789557.\n",
      "  RuntimeWarning\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/skbio/stats/ordination/_principal_coordinate_analysis.py:152: RuntimeWarning: The result contains negative eigenvalues. Please compare their magnitude with the magnitude of some of the largest positive eigenvalues. If the negative ones are smaller, it's probably safe to ignore them, but if they are large in magnitude, the results won't be useful. See the Notes section for more details. The smallest eigenvalue is -0.019937968102490343 and the largest is 2.2369302704010403.\n",
      "  RuntimeWarning\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/skbio/stats/ordination/_principal_coordinate_analysis.py:152: RuntimeWarning: The result contains negative eigenvalues. Please compare their magnitude with the magnitude of some of the largest positive eigenvalues. If the negative ones are smaller, it's probably safe to ignore them, but if they are large in magnitude, the results won't be useful. See the Notes section for more details. The smallest eigenvalue is -0.006804585366340911 and the largest is 0.3896059974961698.\n",
      "  RuntimeWarning\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:47: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:48: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:47: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:48: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:47: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:48: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:47: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "/home/kalen/miniconda3/envs/qiime2-2020.6/lib/python3.6/site-packages/ipykernel_launcher.py:48: DeprecationWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 515.5x360 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 515.5x360 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 515.5x360 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 515.5x360 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAFtCAYAAAAaiCMCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeVwV1f/48dcFQcQFgRJMcQdXcMUlFZM0TcGFXDD3MCv3cK9cMhXLckMpzcQFc0lQZBEtF0wzl08qmppKaoAKBZqsApf5/cGP+Xpl8aogS+/n4+HjwZyZe+Z954LvO+ecOUejKIqCEEIIIcoMg+IOQAghhBCFS5K7EEIIUcZIchdCCCHKmHLFHYAQQoiSKyMjg+joaNLS0oo7FJEPExMTatasiZGRkVqmkQF1Qggh8nPjxg0qV66MpaUlGo2muMMRj1EUhfj4eBITE6lbt65aLs3yQggh8pWWliaJvQTTaDRYWlrmalmR5P6CRUdH07BhQwICAoo7FCFKJPkbKXkksZdseX0+0uf+HE6ePMmIESPU7XLlylGpUiXq1KlDmzZtGDx4MLVq1SqUc/3+++8cOnSI/v37U7NmzUKpsySIiopiz549dOvWjcaNGxd3OOIF8vX1xczMDDc3t+IORYgyR+7cC8GAAQP44osvWLRoEePGjaNOnTp8//339OrVi23btukcW6NGDSIiIujbt+9TnePSpUusXr2amJiYwgy92EVFRbF69WouX75c3KGIF8zX15fdu3fnKn/WvxHxYjVu3Ji+ffvi4uLC+++/z4MHDwo8PiEhgYEDB9KvXz/OnDnzgqJ8dtHR0bi4uABw4cIFFi5cWOCxQUFBLyo0vcideyFo0aJFrv+Ipk2bxvvvv8+nn35KnTp16NChA5DdfFK+fPniCFOUYIqikJ6eXuJ+N9LS0jAxMXmh55S/kdLBxMSEwMBAAGbOnMnWrVv54IMP8j3+xIkT1KtXj88//1zvc2i1WgwNDZ871udlb2+Pvb19vvtjYmIIDg7G1dX1BUZVMLlzLyIvv/wyK1aswMDAgDVr1qjlefUnpqSksHTpUrp164a9vT3t27fH3d2d0NBQALy9vfnkk08AGDFiBA0bNtSp48yZM0ydOhVnZ2fs7e1p164dkyZN4saNG7niatiwIdOmTePUqVMMHDgQBwcHunbtypYtW/J8HwcPHmT48OG0atWKFi1a4Orqyrp163SOSUlJYdmyZXTr1o1mzZrRqVMn5s+fz7///lvgNQoICGD06NEAzJ49W31f3t7e6jGRkZFMnDiRtm3b4uDgQP/+/dmzZ0+B9T5af8OGDfn555/x8fGhS5cu2Nvb4+7uzpUrV3IdHxsby6xZs3j11Vdp1qwZvXr1YuPGjejzQMnJkydp2LAhO3fuZP369epn4ebmxi+//KJzbM7vwPLlywkICKB3797Y29sTEhICQGZmJuvWraNXr17q5zlt2jTu3r2rU09UVBRTp06lc+fO6nUfO3ZsrvcWFRXFjBkz1PfVo0cP1q9fT1ZWls5xOb8bR44cwc3NDXt7e7799lvGjx9P27ZtSU9Pz/W+V69eTcOGDdXftT/++IOPPvqIN954g+bNm9OmTRs8PDyIiIjIda7Y2FhOnTqlfu7Ozs461yfn9/unn36iYcOGhIWF5Tr/7du3adSoEcuWLdMp37lzJ/3798fBwYE2bdrwwQcfcP369bw/PPHcWrRoQWxsLAB//fUXHh4euLm58fbbbxMZGcnly5dZunQp4eHh9O3bl7S0NI4dO8bgwYPp378/kyZNIjk5GQBnZ2dWr17NkCFDCAsLK/C4VatW0b9/f1xdXYmMjAQgOTmZ2bNn4+rqiqurK/v37wfIt55HXbx4kT59+jB48GC2bt2qlp88eZL33nsPgFOnTtG3b1/69u1Lv379SEpK4quvvuLMmTP07duXjRs3Eh0dzdtvv03//v3p378/v/32m1rP8OHDmTRpEj179mTq1Knq/y8RERG4u7vTp08fBgwYQFJSElqtls8//5y33noLV1dXtm/frvdnInfuRcjGxgZHR0dOnz5NcnIyFStWzPO4Tz/9lJCQENzd3WnYsCFJSUlcuXKFc+fO0atXL7p3787du3fZtWsX77//PvXq1QOgVatWAOzbt4+4uDj69++PlZUVMTEx7Ny5k6FDhxIcHIyFhYXO+f744w8mT57MoEGD6NevHyEhISxcuJAGDRqoLQwA3333HV988QW2trZ4eHhgYWFBZGQkP/30E2PHjgUgPT2dUaNGce3aNQYPHkzdunW5efMmW7du5fz58+zYsQNjY+M837ejoyNjx45l3bp1DB48mNatWwPZ//ED3Lp1C3d3dxRFYejQoVhYWBAaGsrMmTOJj4/Hw8NDr89hxYoVaDQaRo8eTVpaGhs2bGDcuHEcOHCAcuWy/wTu3buHu7s7f//9N0OHDsXGxobDhw/j5eXFX3/9xdy5c/U619atW0lMTGTw4MEA7Nixg7Fjx7Jx40batGmjc+yhQ4e4d+8eQ4YMwcLCgnr16qEoCpMnTyY8PBw3NzdGjBhBbGwsW7du5cyZM+zevRtzc3MyMjJ45513SElJYciQIVhbWxMfH8/p06eJjIykUaNG6jUcPHgwpqamDBs2DAsLC06dOsXSpUuJiYlh3rx5OjH9/vvvHD16lCFDhuDu7o61tTW2trb89NNPHD16lG7duukcHxISQrNmzdRHcI4dO8aVK1d48803qVGjBv/88w+7du1i+PDhBAQEUL9+fQC1G8vS0pL3338fIN+/DycnJ6pWrUpQUBA9e/bU2RccHIyiKDp3TIsWLcLPz49evXqp/0lu3boVd3d3/P39qV27tl6fpdCPVqvlxIkTDBgwAIA5c+aoLZbnz5/n008/ZfPmzUyaNImLFy8yd+5cEhIS+Prrr/H19cXU1JR169bh6+vLhAkTAChfvjzbtm0jISGBiRMn5nucubk5u3fvZuvWrWzYsIFFixbh4+NDpUqV1Gbyf//994nnyzF79mzmzJlD27Zt821h2LBhA3PnzqV169YkJydTvnx5pk6dyoYNG1i7di0Aqamp+Pr6Ur58eW7evImnp6f6ZfXSpUuEhIRQrVo1hgwZwv/+9z8cHBz48MMPWb58OQ4ODiQlJWFiYsKuXbuoXLky/v7+pKen4+7uTseOHbGxsXnyB6OIZ/brr78qdnZ2ys6dO/M95rPPPlPs7OyUy5cvK4qiKFFRUYqdnZ3i7++vHtOmTRtl/vz5BZ5r586dip2dnfLrr7/m2pecnJyr7M8//1SaNm2qrF27Vqfczs5OadSokXLx4kW1LC0tTWnfvr0yadIktSwqKkpp0qSJMmTIECUtLU2njqysLPXndevWKU2bNtWpT1EU5aeffnritVEURTl+/Hiu65Fj0qRJSsOGDZWIiAi17OHDh8qAAQMUe3t7JT4+vsC6/f39FTs7O8XNzU3JyMhQy/fv36/Y2dkpR44cUcs+//xzxc7OTgkLC9N5n+PHj1fs7OyUK1euFHiunN+FVq1aKXFxcWp5bGys0qJFC2XAgAFqWc7vQLNmzZTo6GidekJCQhQ7Ozvl0KFDOuW///670rhxY2XZsmWKoijK5cuXFTs7OyU0NLTAuMaMGaM4OzsriYmJOuWLFy9WGjZsqNy4cUMts7OzU+zs7JTffvtN59i0tDSlVatWyuTJk3XKL168qNjZ2Sm+vr5qWV6/i/Hx8Ur79u2VOXPm6JR37txZGTZsWK7j8/obmTNnjtKsWTPl33//1TnWxcVF6dOnj7p97tw5xc7OTtmyZYvOcXfu3FFatmypTJ06Ndf5RMEuXbqUZ3mjRo2UPn36KK1bt1ZGjBihZGZmKklJSYq9vb3Sp08f9V/Pnj0VRcn+e/z0008VRVGUQ4cOKW3btlWPefPNN5XZs2criqIoXbt2Vf8unnTc3bt3FUXJ/txHjhypKIqi9O/fX+f3+kn15Hjw4IHSpUsXdfvy5ctK7969FUXJ/vseO3asoiiKsnbtWmXAgAHKpk2blDt37uTan1PXtGnT1N9PBwcH9bhRo0apx82dO1fZs2ePcuXKFWXw4MG5rvHEiROVN954Q427a9euys8//5zn5/H45yR37kUs524kryagHJUrV+b8+fPcuXOH6tWrP/U5TE1N1Z+Tk5NJT0/HzMyMunXrcuHChVzHN2/enKZNm6rb5cuXp0WLFkRFRallP/74I5mZmUycODFX/+ejj12EhobSrFkzqlevTkJCglresmVLTExMOHHiBAMHDnzq96TVajl69Cjt27fX6esyNjZm5MiRTJ06lePHj+vVxzVo0CD1Dh2gXbt2QHbzYY5Dhw5Ru3ZtevToofM+PTw8+PHHHzl8+LDaolCQ3r178/LLL6vb1apVw8XFhZ07dxIfH4+lpaW6r0uXLtSoUUPn9aGhoVhbW9O8eXOd62ltbY2NjQ2//vorAJUqVQLg559/xsnJKc+73gcPHvDzzz8zevRo0tPTderr1KkTGzdu5Ndff6VOnTpqeZMmTWjZsqVOPeXLl+eNN94gNDRUpwUqODgYQ0NDevfurR776O9iamqq+uytg4MDFy9efMLVy5+rqys7duxg//796u/TtWvXuHr1KtOnT1ePCw0NxcjIiO7du+u8X2NjY+zt7dXrJ55fTp97YmIi7733Hlu3bsXNzY0qVaqoffH5URSFjh075upOyVGhQgW9jsuZkc3AwACtVqu+5vFHw55UT36vy8vYsWPp0qUL4eHhDBo0CF9f31zHbNy4kZdeeonAwECysrJwcHBQ9z3akmloaIhWq8333Iqi8Mknn9C5c+cnxvU4Se5FLCep59fkCDB9+nRmzZpF165dady4Me3bt1f7W/URFxfHV199xeHDh3P1c5ubm+c6/pVXXslVZmZmxh9//KFu37x5E+CJCe3PP/8kLS1Npzn/UfHx8U8KP08JCQmkpKSozbiPatCgAZDdN6uPx9+vmZkZgM61iomJoWPHjs99rkdniHq8LDo6Wie55/VI459//sndu3fzvZ45zXE1a9bknXfeYcOGDQQFBdGyZUs6deqEq6ur+gXx5s2bKIrChg0b2LBhQ571Pf755PeYpaurKwEBAfz444/069cPRVEIDQ2lffv2Ol9mkpKSWLFiBWFhYfz999961a2PNm3aUL16dYKDg9XkHhQUhEajUUc0Q/b1y8jIwMnJKc96DAxkmFFhq1y5Mp988gnjxo1jyJAh1KxZk3379vHmm2+iKAp//PGH2k2Uo0WLFixYsIBbt25Ru3ZtUlNTuXv3bq6/H32Pe1THjh3x8/Pj448/BrL/zvWpp0qVKlSqVIkzZ87Qpk2bfEe///XXX+o4kXPnznHjxg2sra11buASExOxtrbGwMCA3bt3q1888lOvXj3i4uKIiIjQaZbv1KkT27Zto3379hgZGXHjxg2srKx0vkTnR5J7Ebt69SqGhoYF/sf25ptv0rZtWw4dOsSJEycICAhQ+4Me7xN6XFZWFh4eHsTGxjJy5EhsbW0xNTXFwMCAxYsX5zkYrDD/g8vKyqJly5ZMmjQpz/1VqlR5rvrz+zab37685DfaNq9r8zRxPM9xQJ4jwrOysrCxsWHBggVPfM3MmTMZMGAAhw8f5sSJE3h7e+Pj44O3tzedO3dWB8y5u7vrtEg86vG+u/xGqeck8eDgYPr168epU6e4e/cukydP1jnO09OTEydOMGrUKJo0aULlypUxMDBg7dq1Oi1DT0uj0dC7d282bNhAXFwc1apVIyQkBEdHR6ytrdXjsrKyqFChAj4+Ps98LvH0mjRpQqNGjQgJCWHp0qXMnz+fr7/+mszMTHr16pUruVtYWODl5YWnp6c6UHPKlCm5kra+xz3qgw8+YMGCBbi4uGBgYMCECRN444039KrHy8uLjz76iAoVKtCpU6c869+0aRMnT57EwMCABg0a4OTkhEajwdDQkD59+qgDCSdOnEhYWBjt2rV7YjI2NjZm+fLlLFy4UH1CxdfXl4EDBxITE4ObmxuKomBubq7377Yk9yIUFRXFmTNnaNWqldqMmh9LS0sGDhzIwIEDSU1NZezYsfj4+DBmzJgCH0W6evUqV69excvLK9dkIPfv38/zzl0fOU21f/zxR753kQC1a9fmwYMHvPrqq890nvySoYWFBaampnmOcP7zzz8BcjVpP48aNWqo9T4qZwSuvufKq46ckeT63LnWrl2b06dP07ZtW52uhPzUr1+f+vXrM2bMGO7cuUO/fv1Ys2YNnTt3platWmg0GjQazTN/PjkMDAzo3bs3fn5+JCQkEBwcrDbX53jw4AHh4eFMmDCBiRMn6rx+5cqVuep82lnPXF1dWb9+PaGhoTg4OBAdHa2OYM5Ru3Ztjh07hp2dHS+99NJT1S+eztmzZ3W2v/nmG/Xn7777Ltfxbm5uOv9HdejQAX9//1zHHTp0SGdbn+Ps7e3VJ34qVqyY52C4/Op5VLNmzdi7d6+6nfN73K5dO7U7b86cOXm+dtOmTTrbj975T506NVc9gM5AXQcHB3bu3JmrXk9PTzw9PQuMOy/SRlVE/vnnH6ZMmUJWVhbjxo3L9zitVktiYqJOWYUKFahbt67Ovpxm/ccnisi5C3/8sabdu3fnahZ9Gt27d8fQ0JA1a9bkegTq0Tve3r17ExkZmWcTllar5f79+wWeJ+cb7ePvy9DQkM6dO/Prr7/y+++/q+UZGRls3rwZY2PjPJvRn1XXrl25desWP/74o1qW06Sds18fISEhOtc9Li6O4OBg7O3tdZrk89O7d2+Sk5NZv359rn2Koqj9yElJSWRmZursr169OhYWFuo1t7CwoEOHDuzZs4dbt27lqi8pKSnPx9vy4+rqSmZmJnv37uXAgQN07dpV50urgYEBGo0m1+/iyZMnOX/+fK76TE1Nn/i45KMaNWqEnZ0dQUFBBAcHY2RklKtFIqf/f9WqVXnW8Wg/vBBlmdy5F4Jz585hbGyMoig8ePCAS5cusX//fjIyMpg7d26Bd03Jyck4OTnRrVs3GjVqhJmZGZcvX2bXrl06/ZnNmjVDo9Gwdu1aHjx4gImJCQ4ODtSrV0+dGOL27dtUq1aN8+fPc/DgQf0el8hHzZo1mTJlCl999RVvvfUWvXr1wtzcnJs3b3L27Fl27NgBwDvvvEN4eDjTp0/nyJEjtGjRAsh+BOvAgQNMmTKlwOlFGzRoQIUKFdi2bRumpqZUrFgRW1tb7OzsdJp4H30U7ty5c8yYMSPXI37PY+zYsYSFheHp6ak+CnfkyBGOHj3K0KFD9RpMB9nXbfDgwbi7uwOwfft20tPTmTlzpl6vd3V15aeffmL58uWcPXuW9u3bU758eaKjozl48CAuLi5MnDiRX3/9lfnz59OjRw/q1q2LoaEhR44c4c8//9TpIpk/fz5DhgyhX79+DBw4kPr165OYmMi1a9c4cOAAQUFBeveF5zzytmrVKpKTk3MNZqxUqRIdOnRg/fr1pKWlUbduXa5evYq/vz+2tra5BpXa29uzd+9eVq9eTZ06dTA1NVWfdS/o+nz11VfcuHGDLl26qOMncrRu3ZrRo0fj6+vL9evX6dq1K1WqVOH27dscPXqUhg0bsmTJEr3erxClmST3QrBr1y527dqlzi1fq1Yt3n77bQYNGvTEZ2pNTEwYOnQoJ06cIDw8nPT0dKytrRk7dixjxoxRj6tVqxZz585lw4YNzJkzB61WqzbFr127lsWLF+Pn50dGRgYtWrRg06ZNLF68+Lne19ixY6lVqxYbN25k7dq1aDQabGxsdP5TL1++PJs2bWLDhg2EhIRw4MABypcvzyuvvIKLiwvt27cv8Bw5TWje3t4sWLCAjIwMJkyYgJ2dHXXq1GH79u0sX76crVu3kpaWRv369fn888/p16/fc723x5mbm7N9+3aWLVvGnj17SEpKwsbGhlmzZjFq1Ci96xk6dCgPHjzg+++/Jy4uDltbWz799FMcHR31er1Go2HFihVs27YNf39/dSKk6tWr07FjR958800ge6Dja6+9xrFjx/D398fQ0JA6deqwcOFC9XljyG6m3r17Nz4+Phw4cIB//vmHKlWqUKdOHcaPH68zGE4frq6urFq1CjMzszwHrX355ZcsWbKEPXv2kJqaSqNGjfDx8SEwMJBTp07pHDt16lTu37/Phg0bSE5OpkaNGk9M7i4uLixbtozk5GSdgXSPmjVrFs2aNWPr1q34+PigKArVqlWjdevW6vwDQpR1sp67EIUgZxGhhQsXPtOjf0KUVJcvX5ZFnUqBxz8n6XMXQgghyhhJ7kIIIUQZI8ldCCFEiZYzQLUgGzduJDU19bnOM3z48Dxn9SyNZECdEIWgXbt2OjP8CSEKjz6roW3evJk+ffqoU9fqo6QsKVsUyvSdu76rhgkhhCi5ctY7yG/J1M2bNxMXF8fIkSMZPnw4kP8Sr48vKfu4sLAwBgwYQI8ePThz5gyQvXz0o7NGvvfee5w8eVKNbenSpbi5uTFq1CgiIiIYPnw4r7/+OgcPHgR4piVgn1eZvnO/d+9ecYcghBBlxvDhw3PNIPe8nJ2d1dnl9JHXkqkjRoxg48aNbNq0CQsLC72XlM2LVqtl165dhIeHs3r1ajZu3FhgPCkpKbRt25bp06czfvx4VqxYwYYNG4iMjGTmzJm8/vrrWFpaPtUSsI8vD/0synRyF0IIUbY4ODio6wk0atSImJiYXMnw/PnzXL9+nSFDhgCo83/k6NWrV771d+/eHYCmTZsSExPzxHiMjIzUOR/s7OwwNjbGyMgIOzs79fWZmZksWLCAK1euYGBgoC7Mpe/7eRaS3IUQQujlae6wi0peS6Y+Tt8lZWfPns2lS5eoVq0a3377rU79jy4ja2hoqDOt8sOHD9WfjYyM1HUSDAwM8nz90y4BWxjKdJ+7EEKI/4aKFSuq/eotWrTgt99+U9dUSE1NVRdwepSXlxeBgYFqYs9PjRo1uHLlCllZWdy5c4eIiIinii0xMZGXX34ZAwMDAgMDCy2BF0Tu3IUQQpR6gwYN4t133+Xll19my5YtT71UbEFat25NjRo1cHV1xdbWlqZNmz7V6592CdjCUKann3Vzc1MHLQghhHh6Mv1s6SDTzwohhBBlnCR3IYQQooyR5C6EEEKUMZLchRBCiDJGkrsQQghRxkhyF0IIIcqYF/ac++zZszly5AiWlpYEBwfn2r9+/XqCgoKA7Ll9IyMjOXHiBFWrVsXZ2ZmKFStiYGCAoaGh3o+3RUREUKNGjUJ9H0IIURB9piwVT8fd3f2JK8Nt3LiRwYMHP9WqcIUtICCAixcvMnfuXLZt20aFChXo16/fU9Xx4MEDgoKCGDp06HPF8sLu3N3c3Fi/fn2++8eMGUNgYCCBgYF4enri6OhI1apV1f2bNm0iMDBQnlsXQoj/GH2XfH3a9dwLmikuMzPzqep63JAhQ546sUN2cs9vUZun8cLu3B0dHYmOjtbr2JCQEFxcXJ77nA4ODvJlQAghSrmWLVty9uxZTp48yerVqzE3N+fq1as0bdqUL7/8ki1btqhLvlatWpUtW7Zw7NgxvL29SU9Px8bGBi8vLypWrIizszNubm4cP36cYcOG0bt3b/U8w4cPp2XLlvz22284OzvTtm1blixZQkpKCubm5nh5eVGtWjWGDx9Oo0aNuHDhAklJSSxevFhnvngAb29vTE1N8fDw4NatW8ybN4+EhAQMDQ1ZuXIllpaWjBs3jgcPHpCZmcnkyZPp1q0bX331FX/99Rd9+/bl1VdfZebMmaxfv559+/aRnp5O9+7dmTRp0hOvWYmbfjY1NZWff/6ZOXPm6JR7eHig0WgYPHgwgwcPLqbohBDiv62grs7PP/+cYcOGAeDn58fMmTPzPfZZuy+KesnXBw8e4OfnR0ZGBsOHD8fHxwcLCwtCQ0NZvnw5Xl5eQHau2r59O6dPn+ajjz7Ks7s5x7Rp0xg7dizdu3fn4cOHZGVlYWRkxJo1a6hUqRIJCQkMHjyY119/nalTp3Lt2jUCAwOB7HXpb926xa5du1AUhQ8++IDTp0/j6OhY4HUqccn98OHDtGrVSqdJftu2bVhZWREfH8/o0aOpV6/eE9+YEEKIsqeol3zN2Xfjxg2uXr3K6NGjAcjKyuLll19Wj8u543d0dCQpKYkHDx7kWV9SUhKxsbHqUrLly5dXY1q2bBmnT5/GwMCA2NhY/vnnn1yvP378OMePH1eb+FNSUrh582bpS+4hISE6zSQAVlZWAFhaWtK9e3ciIiIkuQshRDHQ94572LBh6l18YSrqJV9z9imKgq2tLTt27MizjpxlXvPbfpKgoCASEhIICAjAyMgIZ2dnnaVkH30vY8eOxd3d/anqL1GPwiUmJnL69Glef/11tSwlJYWkpCT15+PHj2Nra1tcIQohhCiBCnvJ17p165KQkMDZs2eB7Dvta9euqftDQ0MBOHPmDJUrV6Zy5cp5xlWpUiWsra356aefAEhPTyc1NZXExEQsLS0xMjLi119/Vb80Pfo+ADp16oS/v79aFhsbS3x8/BOvxwu7c/f09OTUqVPcu3cPJycnJk6cqI5GzGk6+fHHH+nYsaPOcnjx8fGMHz8eyB7Z6OLigpOT04sKWwghRClQ2Eu+Ghsbs2rVKhYuXEhiYiJarZaRI0eqN5dmZma4u7urA+oK8sUXXzB37lxWrlyJkZERK1euxNXVlQ8++AA3NzcaN25MvXr1ADA3N6dVq1a4uLjQuXNnZs6cSWRkpHrnbmpqytKlS7G0tCzwnLLkqxBCiHzJkq+5DR8+nBkzZmBvb1/coahkyVchhBCijCtxA+qEEEKIkmzLli3FHcITyZ27EEKIApXh3tsyIa/PR5K7EEKIfJmYmBAfHy8JvoRSFIX4+HhMTEx0yqVZXgghRL5q1qxJdHQ0f//9d3GHIvJhYmJCzZo1dcokuQshhMiXkZHRUz1CJkoGaZYXQgghyhhJ7kIIIUQZI8ldCCGEKGMkuQshhBBljCR3IYQQooyR5C6EEEKUMZLchRBCiDJGkrsQQghRxkhyF0IIIcoYSe5CCCFEGSPJXQghhChjJLkLIQZCNxwAACAASURBVIQQZYwkdyGEEKKMkeQuhBBClDGS3IUQQogyRpK7EEIIUcZIchdCCCHKGEnuQgghRBkjyV0IIYQoYyS5CyGEEGWMJHchhBCijJHkLoQQQpQxktyFEEKIMkaSuxBCCFHGSHIXQgghypgXltxnz55Nhw4dcHFxyXP/yZMnad26NX379qVv376sXr1a3Xf06FF69OhB9+7dWbdu3YsKWQghhCiVyr2oE7m5uTFs2DBmzpyZ7zFt2rRh7dq1OmVarZYFCxbg6+uLlZUVAwYMwNnZmQYNGhR1yEIIIUSp9MLu3B0dHTEzM3vq10VERFC7dm1sbGwwNjamd+/eHDx4sAgiFEIIIcqGEtXnfu7cOfr06cOYMWO4du0aALGxsVhbW6vHWFlZERsbW1whCiGEECXeC2uWf5KmTZty6NAhKlasSHh4OOPHj+fAgQMoipLrWI1GUwwRCiGEEKVDiblzr1SpEhUrVgSgS5cuZGZmkpCQgLW1NXfv3lWPi42NpVq1asUVphBCCFHilZjk/vfff6t36REREWRlZWFubo69vT03b94kKiqK9PR0QkJCcHZ2LuZohRBCiJLrhTXLe3p6curUKe7du4eTkxMTJ04kMzMTgCFDhrB//362bduGoaEhJiYmLFu2DI1GQ7ly5Zg7dy5jxoxBq9Xy1ltvYWtr+6LCFkIIIUodjZJXp3YZ4ebmRkBAQHGHIYQQQrxQJWZAnRCiZAsLC8PPzw8ACwsLevXqRc+ePYs5KiFEXiS5C1EChYWFERoaSkJCAlBykml8fLwajxCi5JLkLkQJVpKSac+ePQkNDQVg1apVxRyNEKIgktyFKIF69uxJz549mTRpEiDJVAjxdCS5C1EKldRmeyFEySDJXYhSrCQ12wshSg5J7kKUQtJsL4QoSImZoU4IIYQQhUOSuxBCCFHGSHIXQgghyhhJ7kIIIUQZI8ldCCGEKGMkuQshhBBljCR3IYQQooyR59yFKEVWrVrF9evX1e1r164BqM+7N2jQQP1ZCPHfJcldiFLk+vXrXL34G7UqaQGoomgASLt5mr+SDIszNCFECSLJXYhSplYlLZ+0ScpVvvBMpWKIRghREkmfuxBCCFHGSHIXQgghyhhJ7kIIIUQZI33uQpQgTxoNf+3aNWyMiiU0IUQpIsldiBLk+vXrnL1wiSzT7PXZNdrsP9H/Rd4FwCA5BaoWW3hCiFJCkrsQJUyWqQVpTVzy3Gd6ZhOQ/mIDEkKUOtLnLoQQQpQxktyFEEKIMkaSuxBCCFHGSHIXQgghyhhJ7kIIIUQZI6PlhRD/GWFhYYSGhpKQkACAhYUFvXr1omfPnsUcmRCFS5K7EOI/Jz4+HshO7kKURZLchRD/GT179qRnz57qjH+rVq0q5oiEKBovLLnPnj2bI0eOYGlpSXBwcK79e/fu5dtvvwWgYsWKzJ8/n0aNGgHg7OxMxYoVMTAwwNDQkICAgBcVthBCCFHqvLDk7ubmxrBhw5g5c2ae+2vWrImfnx9mZmaEh4czZ84cfvjhB3X/pk2bpAlNCCGE0MMLS+6Ojo5ER0fnu79Vq1bqzy1atODu3bsvIiwhhBCizCmRj8Lt2rULJycnnTIPDw/c3NzYsWNHMUUlhBBClA4lbkDdr7/+yq5du/j+++/Vsm3btmFlZUV8fDyjR4+mXr16ODo6FmOUQgghRMml9537vXv32Lx5M5999hn37t0D4Ny5c8TExBRaMFeuXOGTTz7Bx8cHc3NztdzKygoAS0tLunfvTkRERKGdUwghhChr9EruV69e5c0332TLli1s376dpKQkAA4fPszKlSsLJZDbt28zceJEvvjiC+rWrauWp6SkqOdLSUnh+PHj2NraFso5hRBCiLJIr2b5zz//HBcXFz7++GOdgW+dO3dmxowZep3I09OTU6dOce/ePZycnJg4cSKZmZkADBkyhDVr1nD//n0+/fRTAPWRt/j4eMaPHw+AVqvFxcUlV3+8EEIIIf6PXsn9woULzJkzB41Go1NevXp1/vnnH71OtGzZsgL3L1q0iEWLFuUqt7GxYe/evXqdQwghhBB6NssriqLeZT/qzp07VKpUqdCDEkIIIcSz0+vOvUOHDnz//ffMnTtXLUtPT+ebb76hY8eORRacEIUlLCyMlStX8vDhw1xfVMuVK0f58uWZPHmyLCAihCgT9Eru06ZN4+233+bSpUtkZGSwYMECrl+/jlarlefOhRBCiBJGr+Req1YtAgMD2bZtG2ZmZmRlZdGvXz+GDRuGpaVlUccoxHPLWTBECCH+C/SexMbS0pIJEyYUZSxCCCGEKAR6DagLCAggNDQ0V3loaCh79uwp9KCEEEII8ez0Su7ffvstVatWzVVubm7OunXrCj0oIYQQQjw7vZJ7TEwMtWrVylVuY2NTqNPPCiGEEOL56ZXcK1eunOdyrVFRUZiamhZ6UEIIIYR4dnoldycnJ7744gvi4uLUstjYWJYuXUqXLl2KLDghhBBCPD29RstPnz6doUOH0r17d+rXr49Go+H69eu88sorTJ8+vahjFEIIIcRT0Cu5W1hYsGfPHvbu3culS5eA7MVeXFxcMDExKdIAhXhaYWFh+Pn5Adm/u7169ZJn3IUQ/yl6P+devnx5Bg4cWJSxCFFo4uPjgezkXhjkC0PptmrVKq5fv65uX7t2DYBJkyapZQ0aNNDZFqI00zu5JyUlcf78ef755x8URdHZ169fv0IPTIhn1bNnT3VehlWrVhVavYX9hUG8ONevX+fshUtkmWZ/dhpt9n99/4u8C4BBSkKxxSZEUdAruR87dgxPT08ePHiAoaGhzj6NRiPJXZR5RfWFQbw4WaYWpDVxyXOfyaXgFxyNEEVLr+Tu5eVF165dmT59Oi+99FJRxySEEEKI56BXco+JieHrr7+WxC6EEEKUAnol92bNmhEVFZXnLHXivy0sLIzQ0FASErL7LF/EYLP8znn16lV10JQMmCo8jw5Gk+sqROmgV3IfN24cX3zxBampqTRq1AgjIyOd/VZWVkUSnCg9imOw2ePnvH79Olcv/katSlqqKBoA0m6eBuCvJMO8KxFPJNdViNJHr+T+zjvvADBhwgQ0Go1arigKGo2Gy5cvF010osTLWSc9587tRQw2y++coaGh1Kqk5ZM2Sbles/BMpSKPqyyT6ypE6aJXct+8eXNRxyGEKGb5dXcIIUofvZJ727ZtizoOIUQJIc/zC1H66T2JDUBcXBy3b98mIyNDp9zR0bFQgxJCvHgFdXcIIUoXvZL733//jaenJ2fOnAH+r689h/S5CyGEECWHXku+enl5kZmZib+/PyYmJvj6+vLFF19Qp04dvv3226KOUQghhBBPQa8791OnTrFmzRqaNGmCRqOhevXqdOjQgQoVKvD111/TqVOnoo5TCCGEEHrSK7knJyers9NVqVKFe/fuUadOHRo1asTvv/9epAGKkklW2RJCiJJLr2b52rVr89dffwFQv359AgMDSU9PJygoCHNz8yINUJRMOats/S/yLv+LvEuithyJ2nLq9tkLl3SSvxCidGncuDF9+/bFxcWFSZMmkZqaWtwh5SsgIIDY2Fh1++OPPy4R//88Hpezs7P6qOmzeJrX65Xc3dzc1As1duxY9uzZQ/PmzfH29ubdd9995kBF6ZazylZaExdSmw8ktflAdTtnaU0hROlkYmJCYGAgwcHBGBkZsX37dp39Wq22SM77LPXu3r2buLg4dXvRokU0aNCgMMN6Jo/H9SLp1Sw/YsQI9ed27dqxb98+Lly4QO3atWnYsGGRBSeEeExWFrcSDfOcGe5WoiEvP8ddgRD5adOmDX/88QcnT55k9erVVKtWjcuXLxMUFMSXX37JqVOnSE9PZ+jQobi7uxMXF8eHH35IUlISWq2W+fPn06ZNG44dO4a3tzfp6enY2Njg5eVFxYoVcXZ2xs3NjePHj/Paa6/x448/smvXLgCio6P54IMPCAoKYvXq1Rw+fJiHDx/SsmVLFixYwP79+7l48SLTpk3DxMSEHTt28O677zJjxgzs7e0JDg5m7dq1KIpCly5dmD59OgAtW7ZkxIgRHD58GBMTE3x8fHItjubt7U10dDR///03N2/eZNasWZw7d46ff/6ZatWq8c0332BkZMTFixdZsmQJKSkpmJub4+XlxW+//ZYrLgA/Pz8OHz5MZmYmK1asoH79+ty/f5+PPvqIqKgoKlSowIIFC2jUqBH37t1j6tSpJCQk4ODggKIoen9met2579mzh/T0dHW7evXqvPHGG9StW5c9e/bofTIhRMESEhIwSInH5FJwnv9QtGgVzZMrEqKQZGZmcvToUezs7AC4cOECU6ZMITQ0lF27dlG5cmX8/f3x9/dn586dREVFERwcTKdOnQgMDCQwMJBGjRqRkJDA119/ja+vL7t376ZZs2b4+vqq5ylfvjzbtm3jvffeIyMjg6ioKCB7noU333wTgGHDhuHv709wcDBpaWkcPnyYnj170qxZM7788ksCAwMxMTFR64yNjeXLL79k06ZN7NmzhwsXLvDTTz8BkJKSQvPmzdm7dy9t2rRh586deb7/v/76i7Vr1+Lj48P06dNp164dQUFBmJiYEB4eTkZGBgsXLmTVqlUEBATw1ltvsXz58nzjMjc3Z/fu3bi7u7NhwwYg+0tEkyZNCAoK4sMPP2TmzJkArFmzhlatWrFnzx6cnZ25ffu23p+bXnfus2fPpnPnzlhaWuqUJycnM3v2bPr166dXHUeOHMHS0pLg4OBc+xVFYdGiRYSHh2NiYsKSJUto2rQpAEePHmXRokVkZWUxcOBAxo4dq0/YQpRBGixNMvOd591EZpUThSQtLY2+ffsC2XfuAwYM4OzZs9jb22NjYwPA8ePH+eOPP9i/fz8AiYmJ3Lp1C3t7ez766CMyMzPp1q0bjRs35vDhw1y/fp0hQ4YAkJGRQYsWLdTzPTrV8Ztvvsm+ffsYO3Ys+/btY/ny5QCcPHmS9evXk5aWxv3797G1tcXZ2Tnf93DhwgXatm2rzrbo6urK6dOn6datG0ZGRnTt2hXIXvn0+PHjedbh5OSEkZERdnZ2aLVanJycALCzsyM6OpobN25w9epVRo8eDUBWVhYvv/xyvjG98cYb6jl//PFHAP73v//h7e0NQIcOHbh//z6JiYmcPn2a1atXA/Daa69hZmaWb72P0yu5Pz5pTY74+HgqVqyo14nc3NwYNmyY+o3kcUePHuXmzZscOHCA8+fPM3/+fH744Qe0Wi0LFizA19cXKysrBgwYgLOzc4noTxGisFlYWHDjXjppTVzy3G96ZhNmxul57hOiMOX0uT/O1NRU/VlRFD755BM6d+6c6zg/Pz/Cw8OZMWMGHh4eVKlShY4dO7Js2bI8z1ehQgX15169ejF58mS6d++ORqOhTp06PHz4kE8//RR/f3+qV6+Ot7c3Dx8+fOb3Z2RkpOY1AwODfPv6jY2N1WPyeo2iKNja2qrN7vqc9/Fz5tXcnlfOfRoFJvecvnaNRsP48eN1lnrNysrizz//pHXr1nqdyNHRkejo6Hz3Hzx4kH79+qHRaGjRogUPHjwgLi6OmJgYateurX5T7N27NwcPHiz05F4c65ILIURp1qlTJ7Zt20b79u0xMjLixo0bWFlZce/ePaysrBg0aBApKSn8/vvvfPDBByxYsIBbt25Ru3ZtUlNTuXv3LnXr1s1Vb61atTAwMMDHx0dtks9J5Obm5iQnJ7N//3569OgBQMWKFUlOTs5Vj4ODA4sWLSIhIQEzMzNCQkIYNmxYoV6DunXrkpCQwNmzZ2nZsiUZGRncvHkTW1vbfON6nKOjI3v37mX8+PGcPHkSc3NzKlWqhKOjI0FBQYwbN47w8HD+/fdfveMqMLnXqlULyJ7EpkaNGjp9GTlNGgMGDND7ZAWJjY3F2tpa3ba2tiY2NjZXuZWVFREREYVyzrzIohlCCKGfgQMHEhMTg5ubG4qiYG5ujo+PD6dOneK7776jXLlymJqa8vnnn2NhYYGXlxeenp7qGK4pU6bkmdwh++79iy++4ODBg0D2HCsDBw7E1dWVGjVqYG9vrx7bv39/5s2bpzNwDaBatWp4enoycuRIFEXBycmJbt26Feo1MDY2ZtWqVSxcuJDExES0Wi0jR47E1tY237geN2HCBGbPno2rqysVKlRgyZIlAIwfP56pU6fSv39/HB0deeWVV/SOq8DkvnDhQiA70b7zzjs6zTGFLb9miaJorshLcaxLLoQQJdXZs2dzlbVr14527dqp2wYGBnh6euLp6alzXP/+/enfv3+u13fo0AF/f/9c5YcOHcpV5uHhgYeHh07Zhx9+yIcffpjr2B49eqh38QBbtmxRf3Z1dcXV1TXXax59fzn//z9u4sSJ+b7m0X2NGzdm69atT4zr0fdpb2+vxlm1alW+/vrrXK83NzdXB90BfPTRR7mOyY9efe4TJkzQ2U5OTub06dPUqVOHOnXq6H2yglhbW3P37l11++7du1SrVo2MjAyd8tjYWKpVq1Yo5xRCCCHKIr2S+9SpU2nevDkjRowgIyODQYMGERkZSbly5fD29lZHHD4PZ2dn/Pz86N27N+fPn6dy5cpUq1YNCwsLbt68SVRUFFZWVoSEhPDVV1899/n0FRYWxsqVK3n48CGZmZlqebly5ShfvjyTJ0/+T/bLP/rIVl4MUuJJSDB+wVEJIYQAPZP7yZMn1WH+hw8fJjk5mWPHjrFr1y58fHz0Su6enp6cOnWKe/fu4eTkxMSJE9VkOWTIELp06UJ4eDjdu3enQoUKLF68ODvAcuWYO3cuY8aMQavV8tZbb2Fra/us71cIIYQo8/RK7v/++686c8/x48fp3r07L730Ei4uLqxbt06vE+X3+EMOjUbDvHnz8tzXpUsXunTpotd5Clt+fTH/dU96ZMvkUrAMShRCiGKi1wx1FhYW6mNsv/zyizqgIi0tDQMDvaoQQgghxAui1517z549mTZtGnXq1CEpKYmOHTsCcPnyZWrXrl2kAQohhBDi6eh12z19+nT1uT1fX191JqG4uDgGDRpUpAEKIYQo2zZt2oSLiwu9e/dm48aNQPZ86507d6Zv37707duX8PDwAuuIjo4mKCjoiecKCAhgwYIFhRF2iabXnXu5cuXUAXWPevwZRCGEEKXb+CnTiP2n8FYXtHrJgjUrvsx3/9WrV/nhhx/44YcfMDIyYsyYMbz22msAjBo1Su88ExMTQ3BwcJ7PtP8X5ZvcY2NjsbKyUn8uSM5xQjxOpvUVonSJ/SeBG9VfK7wK7xwpcHdkZCTNmzdXW4QdHR3VBVXyotVq81xm9quvviIyMpK+ffvSv39/hgwZwvz587l48SKGhobMmjWL9u3bA9mtzh4eHkRFRdGtWzdmzJhRaG+3pMg3ub/22mscO3YMS0tLunTpkuescDkLyly+fLlIgxSln0zrK4TIi52dHStWrODevXuYmJhw9OhRmjVrRtWqVdm6dSt79uyhWbNmzJo1CzMzM51lZtPT03F3d6djx45MnTqVDRs2sHbtWgB1ZregoCAiIyPx8PBQV6+7fPkye/bswdjYmJ49ezJ8+HCqV69ebNegKOSb3Ddt2qQuL7d58+YXFpAoW2Ra35Jv1apVXL9+Xd2+du0agPqZRUdH81KxRCb+C+rXr8+YMWPUKc4bNmyIoaEhQ4YMYdy4cWg0GlauXMmSJUvw8vLKd5nZRxc2g+xlVHMWialfvz6vvPIKN27cALKnwa1cubK6LyYm5r+T3Nu2bZvnz0KIsuX69eucvXCJLNPsVhWNNvu/hf9F3sUgJYFKJkZgVFANQjyfgQMHMnDgQCB7ThQrKyt1bpWc/e+//z6Q/zKzJ0+e1NnOa12SHDnLuAIYGhrmu9xraVZgn7u+pM9diNIty9QizwmJTC4FQ1ZiMUQk/kvi4+OxtLTk9u3bHDhwgB07dhAXF6euI/LTTz+pM5Pmt8zs48ur5iyX2qFDB27cuMGdO3eoV68ely5dKpb3+KLlm9zz62fPi/S5CyGEeFYTJ07k/v37lCtXjnnz5mFmZsb06dO5cuUKADVq1FAfX8tvmdmc5vw+ffrg5ubG22+/zbx583B1dcXQ0BAvLy+dO/ayLt/k/mg/++3bt1m6dCl9+vShdevWQHZ/RlBQENOmTSv6KEspGSkuhChtrF6yeOII96eu7wm+//77XGVLly7N89j8lpmF7LFij8pZF/1Rbm5uuLm5qds5A/DKGr363MeMGYOnpydvvfWWWtatWzcaNGjA3r176devX9FGWcrJSHEhRGlR0DPpovTQaxKbM2fOMGfOnFzlbdq04bPPPiv0oMqKsj5S3CAlQV3yVZORCoBiVEHdB9bFFVqheXQk+eOjyAEaNGigsy2EECWBXsndzMyMI0eOMHLkSJ3y8PBw9XG50upJjwGB/AeelwYNGuhs51w32/o5Cd061zGl0aMjyR8dRQ45X2CEEKLk0Su5jxkzhiVLlnDx4kVatmyJRqPht99+Y9++fcycObOoYyxS169f5+rF36hVKftRiCpK9iDCtJunAfgrybDYYivJHv+yU1ZbJ+AJI8mFEKIE0iu558zes2HDBo4ePQpkP/i/YsUKunXrVqQBvgi1Kmn5pE1SnvsWnqn0gqMpG6RFRAghio9eyR2yB9CVhUQuXgxpERFCiOKjd3IX4mlJi4gQQhQPSe5CCCFUsz8cz7/xdwutPjNLa7yWrym0+kqrzMxMypV7cSlXkrsQQgjVv/F3mdXgaqHVt+T6k48ZN24cd+/e5eHDh4wYMYLBgwej1Wr5+OOPuXjxIhqNhrfeeotRo0axefNmtm/fjqGhIQ0aNGD58uVERESwePFi0tLSMDExYfHixdSrV4+3336bOXPm0LhxYwDc3d2ZP38+jRo1Us8dHR3NjBkzSE3Nfpx3zpw5tGrViri4OD788EOSkpLQarXMnz+fli1b5hnT8OHDmTFjBvb29iQkJDBgwAAOHTpEQEAAR44cIT09nZSUFL755hs+++wzrl69ilarZcKECUXW3S3JXQhRoIyMDG6lGebZlXIr0ZCXE+SRQPF8Fi9eTNWqVUlLS2PAgAG88cYbxMTEEBsbS3Bw9lMpDx48AGDdunUcOnQIY2NjtaxevXr4+flRrlw5fvnlF5YvX463tzcDBw4kICCAjz/+mBs3bpCenq6T2AEsLS3x9fWlfPny3Lx5E09PTwICAggODqZTp0588MEHaLVaUlNTuXz5cp4xFeTcuXPs3buXqlWrsmzZMtq3b4+XlxcPHjxg4MCBvPrqq5iamhbm5QQkuQshhChmW7Zs4ccffwTgzp073Lp1i7p16xIVFcVnn31Gly5d6NSpEwANGzZk2rRpvP766+pdb2JiIjNnzuTWrVtoNBoyMjKA7InEfHx8mDFjBv7+/jrTzubIzMxkwYIFXLlyBQMDA27evAmAvb09H330EZmZmXTr1o3GjRtjY2OTZ0wF6dixI1WrVgXg2LFjHDp0SF1r/uHDh9y5c4f69es/3wXMg97JPTMzkwsXLnD79m31wuWQ6WfFi/Ckx+uuXbuGjSxNWuiMjIywMU3Oc3DkwjOVMJFplcVzOHnyJL/88gs7duygQoUKDB8+nIcPH2JmZkZgYCDHjh3j+++/Z9++fXh5ebFu3TpOnz7NoUOH8PHxISQkhJUrV9KuXTvWrFlDdHQ0I0aMAKBChQq8+uqrHDx4kH379uHv75/r/Bs3buSll14iMDCQrKwsHBwcgOxV5fz8/AgPD2fGjBl4eHjQr1+/PGMyNDRUl5hNT0/Xqb9ChQo626tWraJevXpFcSl16JXcb968yXvvvad+K4LstXINDAwwMDCQ5P4Yeca7aBS07jiAQXIKVC228EQJlpCQgEFKfL4TDxmkxJOQ8N9ZMawkSUxMxMzMjAoVKhAZGcm5c+eA7M/M2NiYHj16UKtWLWbNmkVWVhZ37tyhffv2tG7dmuDgYFJSUkhMTFSXHt+9e7dO/Tlrwbdp00a9g378/NbW1hgYGLB79251bfeYmBisrKwYNGgQKSkp/P777zg5OeWKCbJXrbt48SIODg6EhYXl+147deqEn58fc+bMQaPRcOnSJZo0aVIo1/FxeiV3Ly8v6tevz44dO+jatSv+/v7cv3+fxYsXq29O/J8nJiGZtvSZ5TdbHIDpmU1Aep77hBAlk5OTE9u3b8fV1ZW6devSokULAOLi4pg9ezZZWVkAeHp6otVqmT59OklJSSiKwqhRo6hSpQpjxoxh1qxZ+Pr60r59e536mzVrRqVKlfJskgd4++23mThxImFhYbRr107t/z516hTfffcd5cqVw9TUlM8//zzPmADeeecdpkyZwt69e2nXrl2+73XcuHEsXryYPn36oCgKNWrUKLJV6fRK7ufPn2fjxo1UrVpVvXNv1aoVU6dOZfHixQQEBBRJcKVZQUlIpi0VZcW/6RquRkTg7OxMZmamzr5y5cpRvnx5Jk+eXOzLHFtYWHDjXnqBf5OyamM2M0trvUa4P019BTE2Nmb9+vV57nv8Lhxg27ZtucpatmzJ/v371e0pU6aoP8fGxqIoSr7943Xq1CEoKEjdnjp1KgD9+/enf//+esVUv359nTo+/PBDIPfysiYmJuq69EVNr+SemZlJxYoVATA3N+eff/6hXr161KxZk8jIyCINUAghxItTlp5J37NnD8uXL2fWrFkYGBgUdzgvlF7JvV69evz555/Y2NjQuHFj/Pz8eOmll9iyZQvW1qV/Wc/CJv174r/CzFjBys6hTC4YJEq/fv36/WfHhOmV3EeMGEHC/3+Wdfz48Xh4eNC7d2+MjY1ZunRpkQYohBBCiKejV3J3cfm/fqrGjRtz6NAhIiMjeeWVVzA3Ny+y4Eor6d8TQghRnJ5pEhsTExOaNm361K87evQoixYtm4E+2QAAIABJREFUIisri4EDBzJ27Fid/evXr1cHJWi1WiIjIzlx4gRVq1bF2dmZihUrYmBggKGhoQziE0IIIfKhV3KfPHkyjRs35v3339cpX7duHZcuXWLFihVPrEOr1bJgwQJ8fX2xsrJiwIABODs706BBA/WYMWPGMGbMGAAOHTqkjtDPsWnTJrnjFUIIIZ5Ar+GDp0+fpkuXLrnKnZycOHPmjF4nioiIoHbt2tjY2GBsbEzv3r05ePBgvseHhITodAeUNgYpCZhcCsbkUjAVzv9AhfM/qNvynLsQQvw3bdy4UV2kBuDdd9/Va476p6XXnXtiYmKeE9ubmJjw77//6nWi2NhYnZH1VlZWRERE5HlsamoqP//8M3PmzNEp9/DwQKPRMHjwYAYPHqzXeYvDo60R8H8z1NnWz3n/1rmOEUKIkmDC1AnExscWWn1Wllas/mp1odX3LIpquVWtVouhoeFTvWbz5s306dNHnZb222+/LfS4QM/kXqtWLY4fP07t2rV1yo8fP07NmjX1OlHOvLuPypkQ53GHDx+mVatWOk3y27Ztw8rKivj4eEaPHk29evVwdHTU69wv2uPTyuZsP/64UFhYGKGhoeqTCMOGDSv2yT6EEP9tsfGx3G59u/Aq/F/Bu6Ojo3n33Xdp3bo1Z8+excrKCh8fH0xMTLh8+TLz5s0jNTWVWrVqsXjxYszMzHReP2vWLIyNjbl+/Trx8fHMmjWLrl276r3c6rVr15g9ezYZGRlkZWXh7e1NnTp1CAwMZMuWLWRkZNC8eXPmzZuHoaEhLVu2ZNSoURw7dozXXnuNq1evsnLlSiB7nnxfX1+++eYb5s2bx4ULF3j48CE9evRg0qRJbN68mbi4OEaOHEnVqlXZsmULzs7O7Nq1i++++45XXnmFoUOHAuDt7U3FihV55513WL9+Pfv27SM9PZ3u3bvrNXW5Xsn97bff5ssvvyQ9PZ1XX30VjUbDsWPH8Pb2VmfieRJra2vu3r2rbsfGxlKtWrU8jw0JCaF37946ZTnzBltaWtK9e3ciIiJKbHJ/WvHx8cUdwjPL+YLy6Pz5vXr1Kuao/rsSEhJK/GxxQjzu1q1bLFu2jIULFzJ58mT2799P3759mTFjBnPmzKFt27asXLmS1atX8/HHH+d6fUxMDH5+fvz111+MGDGCV199Ffh/7d15WNTV/gfw9wDDogg6iuCCaIKYhguaS/YDo1REr5QLSaKYmF01B59rmaahGGZaWNFiaWkmiguZKwxpmaapYdfuYGSChvBgcmXIroBsM/P7A2dk2Bx01u+8X8/j8zDf7ZwvDvOZc77nfI5+y63u3LkTM2bMwIQJE1BVVQWVSoXLly8jPT0dKSkpEIvFWLlyJQ4ePIinn34a5eXl8PPzQ2xsrHbFuPLycrRq1QppaWkYO3YsgNosdW3btoVSqcTMmTNx8eJFzJgxA1988UWj48fGjRuHN998Uxvc09PT8dlnn+HkyZO4evUqUlNToVarMXfuXGRmZt4z/ukV3KdNmwaFQoF3330Xa9euhVqthpOTE55//nlMnz5dn0sgICAAeXl5KCgogKenJw4fPozExMQGx926dQuZmZk68+fLy8uhUqng6uqK8vJynDp1CvPmzdOrXEsWGhqK0NBQ7bcwa/7Qbd++vbmrQERWqmvXrnj44YcBAH379kVhYSFu3bqFW7duYciQIQBq08HGxsY2ev7YsWNhZ2eH7t27w9vbG1euXAGg33KrAwYMwCeffILr169j9OjR6N69O06fPo0LFy5g8uTJAICKigrtZ5y9vT3GjBkDoPZL8//93//h2LFjGDNmDI4fP45XXnkFQG1w3r17N2pqanDjxg1cvny5wVrydfXp0wcKhQJFRUX466+/4Obmhs6dO2Pbtm04deqUNhlPeXk58vLyDBPcgdoW2QsvvKBtofn6+rZogXkHBwfExcVh9uzZUCqVmDRpEvz8/LR5giMjIwEAR44cwYgRI3SurVAoMH/+fAC1zzjGjx+PoKAgvcsm49F8QakvLS3NDLUhiUSC5ORkc1eDqEUcHe9m7LS3t0dlZWWLzq//iFfzWp/lVnv27In+/fvj+++/R0xMDBISEqBWq/HMM89o88zX5eTkpPOcPSwsDNu3b4e7uzsCAgLg6uqKgoICbN68GampqXB3d8eSJUv0uqcxY8YgIyMDxcXF2t5rtVqNOXPmYOrUqfr9Mu5oUbJdkUgEZ2dnODs731ee3uDgYGRkZODo0aOYO3cugNqgrgnsQG2i/XfffVfnPG9vbxw4cAAHDhzA4cOHtecSEZEwtWnTBm5ubtoZWfv372+ytSqTyaBSqZCfn4+CggL06NGjwTGa5VY147+ys7MBAAUFBfD29saMGTMQEhKC33//HcOHD0dGRob2kenNmzdRWFjYaNlDhgxBdnY2du/ere2SLysrg4uLC9q0aYPi4mKcOHFCe3zr1q1RVlbW6LXGjRuHtLQ0ZGRkaHsHHn/8cXz11Vfac4qKivR6lKtXy72qqgpvv/02du3aherqaqjVajg6OiIiIgKvvPIKnJyc9LkMEZFFqz/IVSKRICwszKofmVmztWvXagfUeXt7Y82aNY0e16NHD0RFRUGhUCA+Pr7RmNTUcqtpaWk4cOAAHBwc0KFDB8yfPx9t27bFwoULMWvWLKhUKojFYsTFxaFLly4Nrmtvb4+RI0fi66+/xtq1awEAvXv3Rp8+fTBu3Dh4e3sjMDBQe3xERAReeOEFeHh4YNu2bTrX8vPzQ1lZGTp27Kgdk/b444/j8uXL2pZ7q1at8Pbbb9/zUahewX3VqlX47rvvsHz5cgQGBkKtVuP8+fNISkpCRUUFEhIS9LkMEZFV0LSMbDFplmd7z3uOcG/x9ZrRtWtXHDp0d5GtmJgY7c8PP/wwdu/efc8yAgMD8dprr+ls03e51RdffBEvvvhig+1hYWGNDg4+f/58g21xcXGIi4vT2fbWW281Wtfp06frjFX77rvvdPbXXTpWIzo6GtHR0Y1eryl6Bff09HQkJiZi5MiR2m1+fn7o2LEjFi1axOBODZSUlODGLXsknHNtdP/VW/bwKGEyn8ZoEiABgKi6NtmFWnzn2aGypqnTyADqD3K1xdXuzD0nnQxDr+AuFovRrVu3Btu9vb2NkhiAyFbdKwFSTs4tAFWmrhaRxWqqhWzr9IrMU6ZMwRdffIH4+HjtKES1Wo1t27ZppwoQ1SWRSNDqf5exfHBpo/sTzrnC2Qa7PO/lXgmQpFIpKvIyTV4vIrIuegV3hUKBjIwMnDp1Cv369QMAZGVl4ebNmxgzZoxOmtg33njDODUluheVClebeBQgpMcA+aV37/Hvqtov2+6OauSX2qOXOStGRBZDr+BeUFCAPn36AACKi4sBAJ06dUKnTp2Qn5+vPa6pdLJEZBj1u+0L7nTbe3b3Q69G9pN+RBX/Q07OLW1PSd2Mi0Dt71WflJ9ElkKv4F5/uD6RRbKzg08bZaOPAu73MUBJSQnsyhXaAW46xZUrUFLi2MhZxqPvugXUMiJVNdSVldpHHm7q2oZKRV4m8ktbtjAIkSVoeSaaes6cOWOIegiSTCaDVCpFTk4OcnJyIJVKIZPJzF0tIqpPpULdta3cHdVwd6zdoFZDO++dyFro1XIvKyuDvb09nJ2dtdt+/fVXJCYm4vTp0/jtt9+MVkEhYN516yWRSPDHX1Wo6DO+wT7n7EM2OQ+aiCxfs8G9qKgICxcuxC+//AJ7e3tER0cjNjYWK1euxL59+zBy5Ehs377dVHW1Ok3lXSdhYVYzATDCIx0ic2o2uK9fvx6lpaVYtmwZZDIZNm/ejJ9//hmdO3fGoUOHGiTgJ/0kJSUhNzdX+5qDd4TBlrOaEZFlaTa4nzlzBomJiRg8eDBGjRqF4OBgjBgxAgsWLDBV/QQpNzcXly78G91clQA4eMfaMasZEVmaZoP7jRs3tJnpPD094ezsrF31hh5MN9emuwCJiIgeRLOj5VUqlU56WZFIxBXgiIiILNw9R8vHxsZCLBYDqF369dVXX9UZNQ8AmzdvNk7tiIiIqMWaDe7PPPOMzusJEyYYtTJEZJk0KW/rprvVbBdKytv8Unu88mMb3KzU7dBUqoE+ZqoT0f1qNrivWbPGVPUguqfmssUBAJQ12uAjdJrpd3VnWtzv9Lt7ZeETuTjB75FAALrpbgEIJuWt2k4MkaMjxO3bw+7OrAcNj/btBXGPZFu4XiuRFTNFgiSJRKKzKh0gvBkBamc3+PX0Etx9ke1icCer0Vy2OABodW4r3B1tY61zQyZIYhY+IuGxmeDeVBYxIiIiobGZ4K7BLGJERCR0D7wqnLUIDQ1FUlIS/Pz84Ofnh6SkJOb+JiIiQbKZ4E5kTbhcMBE9CJvrlieyJlwumIjuB4M7kQXicsFE9CBsPriXlJTgxi37JhdsuXrLHh53RtgTERFZA5sP7kRkG+zKS+Dynz0QVZfrbFeLW0FUfRuAl3kqRmQENh/cJRIJWv3vcqPLrwK1S7A6c9ockVXTpI8tKSmBQlGjs699e3dIJD2YYpYExSaCe1JSEnJzcwFAJxe35rW32GxVIyIT0Py9E9kKmwjuubm5OJ+VDVUrCUTK2lv++fJ1AIBdWTnQ1py1IyIiMiyTznM/ceIExowZg1GjRmHjxo0N9p89exaDBg1CeHg4wsPD8eGHH+p97r2oWklQ0Wc8bvefgtv9p6Ciz/jaXNr2NvH9hoiIbIjJIptSqcSqVauwZcsWeHp6YvLkyQgJCWnwnGvw4MH49NNP7+tcsiyaNcABCHodcCIiS2Oy4C6Xy+Hj4wNvb28AwLhx4/Dtt9/qFaAf5Fwyj/r/N0JdB5yIyBKZLLgXFRXBy+vuVBNPT0/I5fIGx/3yyy+YMGECOnbsiFdffRV+fn56n9uUkpIS2JUr4Jx9qOFOZY22VUmGU38AkxDXARdV/A85Obd0BmcCuvfu6+vLwVxEZHImC+5qtbrBNpFIN6j27dsX3333HVq3bo3jx49j/vz5+Oabb/Q6l8jURKpqqCsrUZGXCQBwU9e+JzWv80vtzVY3IrJtJgvuXl5euH79uvZ1UVEROnbsqHOMq+vdLHHBwcGIj49HSUmJXuc2RyKR4I+/qmoH0NXT6txWuDtWteRWiLS6uSqbzZFARGQOJhstHxAQgLy8PBQUFKCqqgqHDx9GSEiIzjE3btzQttLlcjlUKhXatWun17lERIZSXFyMBQsWQKFQmLsqRPfFZC13BwcHxMXFYfbs2VAqlZg0aRL8/PyQkpICAIiMjERGRgZSUlJgb28PZ2dnrF+/HiKRqMlzW8KuvATO2YfupJkE1GKX2h3KmmbOIiJbtHXrVsjlcmzduhX/+te/zF0dohYz6STv4OBgBAcH62yLjIzU/hwVFYWoqCi9z9VX3VHZmkFPfj297ry+BYDd8kRUq7i4GOnp6VCr1UhPT0d0dDSX3iWrYxMZXOqOVq4/alsqlWoHQJlKcyvRcRU6y8NeH9uydetW7eNBlUrF1jtZJZsI7kT3i70+5iOTyZCWloaSO192JRIJwsLCjL7O/ZEjR1BdXQ0AqK6uxjfffMPgTlaHwd0MmluJjqvQWRZL6/WxRZpBbRIT/V2MGjUKaWlpqK6uhlgsxujRo01SLpEhMbiTVdF0kQNgN7nAhYaGIjQ01OQJkKKjo5Geng4AsLOzQ3R0tEnKJTIkBneyGvXT1bKbnIyhQ4cOGDt2LA4cOICxY8dyMB1ZJZsJ7prnd3VThIaFhZm5VtQS90ppy25yMpTo6Gjk5eWx1U5Wy2aCuwa/hRPRvXTo0AEffPCBuatBdN9sJrhrnt/Vl5aWZobaEBERGY/J0s8SERGRaTC4E5FeZDIZcnJykJOTA6lUCplMZu4qEVETbKZbnogeHMesEFkHBnci0ktT41aIyPKwW56IiEhg2HInuk/NLQAEcBEgIjIfBncisihJSUnIzc3Vvq6beErD19e3QVIjIrqLwZ3oPjW3ABDARYDuV25uLi5d+De6uSoBAG5qEQBosw/ml9qbrW5E1oLBnYgsTjdXZbNfmoioeRxQR0Ylk8kglUo5P5qIyITYcieT4PxoIiLTYXAno+LcaCIi02NwJyLYlZfAOfsQAEBUfRsAoBa7wK68BICXGWtGRPdDpFar1eauhLH4+vri9u3bje5bu3YtoqKiIJVKcSnzGLLym56PPGnSJO2a4aGhocjKymr0uGnTpmHdunUAALlcjrFjxzZ5zcd7e+Ht4CoAwFvHy7D/YmWjxwUEBOg8o+7SpUuT19TcEwAkJyfj1VdfbfLYwsJC7c+Guqf09HT069cPALB48WJs37690eOMdU+TJk1CRV4mlg8uxcyv/sbvxUqD3tOIESOQl5fX6HH+HezxxSR37evhn5bA1dUVbm5uDY4Vi8VwdnbG4sWLUVxcbLH/Ty4uLhg2bBiGDBkCqVRq8vdeeG8nLAluDQC4eKMGz+/93wPfkyn+nogsAVvuRNQoJycnREVF8bEKkRUSdMt94sSJ2Lt3b7PHSKVSbWuvMQnnXOHc/VFty90QmivTGOUJlSaJieZ3Zezfa0vKM1SZtoi/V6IHx6lwREREAsNueTPJL72bk/zvqtoMXO6OauSX2qOXOStGZGaWmLNfJpNh3bp1qKmp0dnu4OAAJycnxMbG8vEFWRQGdzPw9fUFUPshplAocLuidtCfuENX9Oou0e4XMplMhvfffx+VlZU6H5j8sCQienAM7mageXYrk8mQlpaGkjutEGMPXqpfnkQiQVhYmNUFUc191F1QJCwszMy1IkOxxJz9zNdA1obB3YzM9YGhUCgA1H6Imosh7p1Z74iIGsfgbkM0AbX+qG9r09QXg7S0NJPXpamxE5p9HD9BRObA4G4DuD62cdQfG1Fw5/fq2d0PANCrkWOIiEzBpMH9xIkTWL16NVQqFaZMmYI5c+bo7D9w4AA2bdoEAGjdujVWrlyJ3r17AwBCQkLQunVr2NnZwd7e/p7z1+kuro9tHPW/DFl7jwgRCYfJgrtSqcSqVauwZcsWeHp6YvLkyQgJCdFp2XTt2hXJyclwd3fH8ePH8frrr2PPnj3a/Vu3bjXrc2JrxvWxiYhsh8mS2Mjlcvj4+MDb2xuOjo4YN24cvv32W51jAgMD4e5em597wIABuH79uqmqR0REJBgma7kXFRXBy+vu6lKenp6Qy+VNHp+amoqgoCCdbTExMRCJRHj22Wfx7LPPGq2uZL00A9w4uI2IbJnJgntjKexFIlGjx545cwapqanYsWOHdltKSgo8PT2hUCjw/PPP46GHHsKjjz5qtPqS9an7iMfQg9tkMlmDefWc90xElspkwd3Ly0unm72oqAgdO3ZscNzFixexfPlybNq0Ce3atdNu9/T0BFA7t3nUqFGQy+UGC+5Cms7UWGIcW1F3gJsxBrdxXr3pCOlvksgcTBbcAwICkJeXh4KCAnh6euLw4cNITEzUOebatWtYsGAB1q1bhx49emi3l5eXQ6VSwdXVFeXl5Th16hTmzZtnkHoJdTqTJlENGQYzlJmOUP8miUzJZMHdwcEBcXFxmD17NpRKJSZNmgQ/Pz+kpKQAACIjI/HRRx/h5s2biI+PBwDtlDeFQoH58+cDqB11P378+AbP4++X0KYz1U9UExoaapbkLkT3S2h/k0TmYNJ57sHBwQgODtbZFhkZqf159erVWL16dYPzvL29ceDAAaPXj4iISAiYoc4GWOISmkRA7RiR5ORkANa7kBGRJWJwJyKzsoSFjIiEhsHdBljiEppEgO6YED5TJzIcBncLI5PJ8P7776OyshI1NTU6+xwcHODk5ITY2NgGXZfNLQ6Tk5MDb7Hx606kr7rvVy5kRGR4DO4CkZubi/NZ2VC1qm2Bi5S1/7U/X74Ou7JyoK05a0ekq+77te57FQDsyjn+g+hBMbhbmAeZT61qJUFFn/ENtrc6txVA1QPWjMiwmnq/OmcfMkNtiISFwZ3wd5UIl+RyhISEtOhRgKViqljLV1JSArtyRaOB3K5cgZISRzPUikg4GNwForkPSyhrtCk8bQVTxRKRLWNwJ7g7quHZq59gRiszVazlk0gk+OOvqia75TktjujBMLgLhEQiwdXCuwvziKpvAwDUYhcAdxfdIMPTLNbDRwEtY1deAufsQw3eq7UD6ryaOZOI7oXBXSDqL6ShCTR+Pb2Qk3MLHFBnfHwUoL+679e679VaXvD19eWXJqIHwOAuEJo5wZoPRI2wsDCkpaXh0oV/cwlNI+FjgJbTZ3lemUwGgF+aiO4Hg7tA1f1A5BKaZKmam9nAL01E94/BXWAa+0Cs/5pLaJIlYcucyPAY3O/g8z0i02PrnMg4GNzrYSuCyLLda/2FxYsX8wsD2TwG9zvYgiAiIqFgcCciq8Iv4kT3ZmfuCpDpyGQy7RKwOTk5kEql2ulGREQkHGy52yCOKyAiEjYGdxvC7kwiItvAbnkiIiKBYXAnIiISGAZ3IiIigWFwJyIiEhgGdyIiIoFhcCciIhIYBnciIiKBYXAnIiISGAZ3IiIigWFwJyIiEhgGdyIiIoERqdVqtbkrYSxDhw5Fly5dzF0NIrIh7dq1w+eff27uapCNE3RwJyIiskXsliciIhIYBnciIiKBYXAnIiISGAZ3IiIigWFwJyIiEhgGdyIiIoFhcCciIhIYBvc7lixZgpkzZ+ps+/XXXzFixAjMnz8fFRUVJinTGJYsWQJ/f38sWLCgwb6jR4/C398fffr0MVq5b775ZoN9/v7+2L9/v8HLBIDi4mK88cYbCAkJwSOPPIJhw4ZBKpXit99+M0p5mvvU/B6feOIJxMXF4a+//jJKeYB57/Hhhx9GUFAQFi9ejKKiIqOUV1FRgffeew+jR49Gv379MHToUEyaNAlffvmlUcrT3FtT/0JCQoxSLpGxOJi7Apbq5MmTWLBgAcLDwxEXFwc7O+v+HtS5c2ccO3YMxcXF6NChg3b7rl270KVLF1y/ft0o5To7O2PHjh2IjIxEjx49jFJGXX/++SemTp2K9u3bY8WKFfDz80NxcTG++OILRERE4KOPPkJQUJDByx08eDDee+89KJVKXLhwAcuXL8f169exceNGg5dl7ntUqVTIz8/HqlWrEBsbi507dxq8rJUrV+Ls2bNYtmwZ/P39UVZWhuzsbFy7ds3gZQG1f+8acrkc8+bNw549e9CpUycAgL29vVHKJTIWBvdG7Nu3D8uXL8e8efMwb948c1fHIHx8fNC5c2fs3bsXc+bMAQBcu3YNP/74I+bOnYuPP/7YKOUOHDgQt2/fxrp167BhwwajlFFXfHw8ampq8OWXX8LV1RVA7Reb9evX44UXXsDSpUvx7bffwtnZ2aDlisVieHh4AAC8vLyQk5ODpKQkVFRUGLwsS7hHT09PREREICEhAaWlpdp6GMrRo0excOFCPPXUU9ptvXv3NmgZdWnuCwDc3d0BABKJRGc7kTWx7uaoEWzatAnLly/HqlWrBBPYNSIiIpCamgpNxuE9e/Zg2LBh6Ny5s1HLXbp0KY4dO4YzZ84YtZy///4bx48fx7Rp0xoNNi+++CKKi4tx6tQpo9YDqO2xUKlUqKmpMeh1LeUei4qKkJGRAXt7e6P0anl4eOCHH37AzZs3DX5tIlvA4F7HuXPn8M477yAhIQETJ040d3UMLjQ0FH///TfOnj0LpVKJr776Cs8++6zRyx0wYADCwsLw1ltvQaVSGa2cq1evQqVSwc/Pr9H9vr6+AIA//vjDaHUAgNzcXGzfvh39+/c3eIvWnPf4008/YeDAgejfvz+CgoKQmZmJ6OhotGrVyuBlJSQk4NKlSxg+fDj+8Y9/4PXXX8fRo0fBpTCI9MPgXsdDDz2Enj17YuPGjUYbKGROTk5OmDBhAvbs2YPvv/8eSqUSTzzxhEnKXrRoEa5cuYK9e/carYx7ffCLRCKjla0JfP369cP48ePh7e2NxMREg5ej7z06OBj+iVu/fv2wb98+7NmzB/PmzcOAAQMQGxtr8HIAYNCgQThy5Ai2b9+OZ555BsXFxZBKpZg7dy4DPJEeGNzrkEgkSE5OhpOTE6KiolBYWGjuKhnc1KlT8c033+Czzz7DxIkTIRaLTVJuly5dMHPmTLz33nsoKyszShk+Pj6ws7PDpUuXGt2fk5MDAEYZ2KcJfGlpaZDL5diyZQu8vb0NXk737t3Ndo/Ozs7w8fFBr169EBsbi86dOyM+Pt7g5Wg4ODggMDAQs2bNwoYNG7BmzRocO3YMmZmZRiuTSCgY3OuRSCTYunUrJBIJpk2bhry8PHNXyaB69uyJgIAAnD9/HlOmTDFp2XPmzIFarcamTZuMcv22bdsiKCgIO3bsQGlpaYP9n3zyCTp06IARI0YYvGxN4OvatSscHR0Nfn0Nd3d3BAcHN3uPHTt2xGOPPWa0Omi89NJL2LdvH7KysoxeFlD73gUAhUJhkvKIrBmDeyPc3NywefNmdOvWDVFRUdrWkDGUl5fjt99+0/l3+fJlo5UHAJ9//jnOnDmDbt26GbWc+lxdXREbG4stW7YYrQzNtMUZM2bgxIkT+PPPPyGXy7Fo0SKcPXsWa9asMfgoclNbsWIF7O3tm7zHxMREk/TI9OzZEyNHjsT69esNfu2oqCikpKQgKysLhYWFOH36NOLj4+Hm5oahQ4cavDwioeFUuCa0bt0amzZtwoIFCzB9+nR8/vnn6Nu3r8HL+c9//oOnn35aZ1uPHj0gk8kMXpaGi4sLXFxcjHb95kyePBnJycn4/fffjXL9Ll264KuvvsKGDRuwYsUK3LhxA61bt8aQIUOwa9cuoyTrMbVOnTph7969+Pjjj7FixQoUFRVBqVSiY8eOOHjwILp3726yusyePRvPPfccTp8+jeHDhxvsukFBQTh48CCSkpJQWlqK9u3bY/DgwVizZg1OCT3oAAAFC0lEQVQkEonByiESKpGao1OIrN73338PqVSKmJgYow1yIyLrwW55IgEYOXIkNm/eDDs7OxQUFJi7OkRkZmy5ExERCQxb7kRERALD4E5ERCQwDO5EREQCw+BOREQkMAzuZJM++OADjBo1qtlj9u7dK4h58URke5jEhqxKRUUFPvnkE6SlpeH69etwcXFB165dER4ejhkzZhi0rLCwMAQFBRn0msuWLUN+fj62bdtm0OsSEdXF4E5WZeXKlTh79iyWLVsGf39/lJWVITs7G9euXTN4Wc7OzlafqpaIbBO75cmqHD16FDExMXjqqafg7e2N3r17Y+LEiXjppZe0xyxZsgQzZ87UOW///v3w9/dvcL2DBw/iySefREBAAGbOnKmTAKaxbvkLFy5g1qxZGDhwIIYNG4aXXnqpweqBP/74I5577jn0798fgwYNQlRUFPLz8/HBBx8gNTUVP/30E/z9/eHv72/UJXCJyHYxuJNV8fDwwA8//ICbN28+8LVu3LiBHTt24N1338X27dtRVlaG+fPnN7leeG5uLqZPn44BAwYgNTUVW7duhZ2dHWbNmoXKykoAtYE9JiYGffv2xa5du7B7926Eh4ejuroas2bNwvjx4zFw4ECcPHkSJ0+eRFhY2APfBxFRfeyWJ6uSkJCAl19+GcOHD4evry8GDBiA4OBgPPnkkxCJRC261u3bt/HWW2/Bx8cHALBu3TqEhobi9OnTjS6Z+tlnn2HkyJGQSqXabe+88w4effRR/PDDD3jqqafw4YcfIigoCMuWLdMeo1mqFKjt6heLxfDw8GjprRMR6Y3BnazKoEGDcOTIEcjlcvzyyy/IzMyEVCpFUFAQNmzY0KIAL5FItIEdqF2Nr127dsjNzW00uGdlZeHq1asYOHCgzvbKykrk5eUBAH799VcsWrTo/m6OiMhAGNzJ6jg4OCAwMBCBgYGYNWsW9u/fj8WLFyMzMxNDhgyBSCRq0LVeU1Oj17WbW2pBpVIhPDwcc+bMabCvbdu22p9b2oNARGRofOZOVk/T7a1QKAAA7du3x3//+1+dY7KzsxucV1JSgvz8fO3rP/74Azdv3tTpRq/rkUcewe+//45u3brBx8dH55+7uzsAoG/fvjh58mSTdRWLxVAqlS27QSKiFmJwJ6sSFRWFlJQUZGVlobCwEKdPn0Z8fDzc3NwwdOhQAMBjjz2GK1euIDk5Gfn5+di9ezfS09MbXMvFxQVLly7FhQsXkJWVhSVLlqBXr16NdskDwD//+U9cvnwZL7/8MuRyOQoKCnDmzBkkJCRoR9nPmzcPJ06cwOrVq3Hx4kVcuXIFe/fuxZUrVwAAXbt2xZUrV5CTk4OSkhJUVVUZ6TdFRLaMwZ2sSlBQEA4ePIg5c+YgNDQUS5cuhY+PD1JSUiCRSADUBveFCxfi008/RXh4OM6cOYP58+c3uJaHhwciIiIglUoRGRkJJycnfPTRR012q/fs2RM7d+5EeXk5YmJiMG7cOLz++uuorKxEmzZtAACPP/44Nm7cCLlcjoiICEyZMgVff/01xGIxAGDy5MkICAjA1KlTMXz4cBw6dMhIvykismVcz52oCTt37sQ777yDc+fOmbsqREQtwpY7USMKCwtx/PjxRhPfEBFZOgZ3okaMHj0af/75J1577TVzV4WIqMXYLU9ERCQwbLkTEREJDIM7ERGRwDC4ExERCQyDOxERkcAwuBMREQnM/wPIIt4iwKu/eQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 515.5x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "base_val = 'no_preservative'\n",
    "base_col = 'sample_preservation_method'\n",
    "id_col = 'host_subject_id'\n",
    "base_condition = lambda df: (df[base_col] == base_val)\n",
    "\n",
    "#core diversity metrics @ different sequencing depths\n",
    "for rf_depth in [38000]:\n",
    "    core_results = core_metrics_phylogenetic(tree_filtered_table,\n",
    "                                             rooted_tree,\n",
    "                                             rf_depth,\n",
    "                                             q2mf_filtered)\n",
    "    \n",
    "    # create metrics, emperor, PERMANOVA, and taxonomy output directories\n",
    "    core_metric_dir = inputs_dir_name + '/core-metrics-%i'%(rf_depth)\n",
    "    emperor_dir = output_dir_name + '/emperor-plots-%i'%(rf_depth)\n",
    "    check_and_create_dir(core_metric_dir)\n",
    "    check_and_create_dir(emperor_dir)\n",
    "    perma_dir = output_dir_name + '/permanova-group-plots-%i'%(rf_depth)\n",
    "    check_and_create_dir(perma_dir)\n",
    "    box_dir = output_dir_name + '/box-plots-%i'%(rf_depth)\n",
    "    check_and_create_dir(box_dir)\n",
    "    \n",
    "    for name, q2art in zip(core_results._fields, core_results):\n",
    "#         move emperor plots to output directory\n",
    "        alpha_dir = emperor_dir if 'emperor' in name else core_metric_dir\n",
    "        out_ = os.path.join(alpha_dir,name)\n",
    "        q2art.save(out_)\n",
    "        \n",
    "        # create permanova group plot and box plots\n",
    "        if 'distance' in name:\n",
    "#             # permanova plot\n",
    "            distance_metric = name.split('_distance')[0]\n",
    "            create_perma_groups(q2art, q2mf_filtered, perma_dir, groups, metric_used=distance_metric)\n",
    "            \n",
    "            # box plots\n",
    "            dist_mat = q2art.view(DistanceMatrix)\n",
    "            test_mat = q2art\n",
    "            create_box_plot_data(base_condition, groups, base_col, id_col, q2mf_filtered, dist_mat, box_dir,\n",
    "                                 'Subject', 'value', 'Preservation method', 'groupby',  metric_used=distance_metric)\n",
    "\n",
    "#make deicode/rpca beta diversity plot\n",
    "rpcares = rpca(tree_filtered_table)\n",
    "\n",
    "#makes a directory to put the output in\n",
    "rpca_dir_name = inputs_dir_name + '/rpca'\n",
    "check_and_create_dir(rpca_dir_name)\n",
    "for name, q2art in zip(rpcares._fields, rpcares):\n",
    "    out_ = os.path.join(rpca_dir_name,name)\n",
    "    q2art.save(out_)\n",
    "    \n",
    "#import biplot made by rpca\n",
    "bi_plot = q2.Artifact.load(rpca_dir_name + '/biplot.qza')\n",
    "empbiplt = biplot(bi_plot,q2mf_filtered)\n",
    "\n",
    "#save rpca visualization\n",
    "output_rpca_dir = output_dir_name + '/rpca-plots'\n",
    "check_and_create_dir(output_rpca_dir)\n",
    "empbiplt.visualization.save(output_rpca_dir + '/empbiplt.qzv')\n",
    "\n",
    "\n",
    "#beta diversity significance rpca\n",
    "rpca_perma_dir = output_rpca_dir + '/permanova-group-plots'\n",
    "check_and_create_dir(rpca_perma_dir)\n",
    "rpca_dist_matx = q2.Artifact.load(rpca_dir_name + '/distance_matrix.qza')\n",
    "create_perma_groups(rpca_dist_matx, q2mf_filtered, rpca_perma_dir, groups, metric_used='rpca')\n",
    "\n",
    "# # box plots\n",
    "dist_mat = rpca_dist_matx.view(DistanceMatrix)\n",
    "create_box_plot_data(base_condition, groups, base_col, id_col, q2mf_filtered, dist_mat, output_rpca_dir,\n",
    "                     'Subject', 'value', 'Preservation method', 'groupby', fname='rpca-box-plot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taxonomy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'inputs/rpca/distance_matrix.qzv'"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tax = q2.Artifact.load('taxonomy.qza')\n",
    "tax_bp = taxa.actions.barplot(tree_filtered_table, tax, q2mf_filtered)\n",
    "tax_dir_name = output_dir_name + '/taxonomy_plot'\n",
    "tax_bp.visualization.save(out_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = dist_mat.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_mat.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_mat.ids[0][:-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}