DBMI Hail Demonstration PCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# i2b2 tranSMART working session - Boston, MA.\n",
    "## 260 Longwood Ave - Tosteson Medical Education Center\n",
    "## TMEC rooms 306 and TMEC L-007\n",
    "### Thursday, June 28th"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hail Demo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://hail.is/docs/devel/hail-logo-cropped.png\" width= \"150px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Documentation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://hail.is/\n",
    "\n",
    "https://hail.is/docs/devel/tutorials/01-genome-wide-association-study.html#Check-for-tutorial-data-or-download-if-necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Running on Apache Spark version 2.2.1\n",
      "SparkUI available at http://10.128.0.7:4040\n",
      "Welcome to\n",
      "     __  __     <>__\n",
      "    / /_/ /__  __/ /\n",
      "   / __  / _ `/ / /\n",
      "  /_/ /_/\\_,_/_/_/   version devel-aa83f2a1d041\n",
      "NOTE: This is a beta version. Interfaces may change\n",
      "  during the beta period. We recommend pulling\n",
      "  the latest changes weekly.\n"
     ]
    }
   ],
   "source": [
    "# Install Hail environment\n",
    "import hail as hl\n",
    "import hail.expr.aggregators as agg\n",
    "hl.init()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Matrix Table file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# about 1.2TB data\n",
    "vds = hl.read_matrix_table(\"gs://1k-genome/1000-genomes/VDS-of-all/ALL.chr.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.mt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting to know our data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Describe the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(39706715, 1092)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#size of the matrix\n",
    "vds.count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1092"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#number of colomns\n",
    "vds.count_cols()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39706715"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#number of rows \n",
    "vds.count_rows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------\n",
      "Global fields:\n",
      "    None\n",
      "----------------------------------------\n",
      "Column fields:\n",
      "    's': str \n",
      "----------------------------------------\n",
      "Row fields:\n",
      "    'locus': locus<GRCh37> \n",
      "    'alleles': array<str> \n",
      "    'rsid': str \n",
      "    'qual': float64 \n",
      "    'filters': set<str> \n",
      "    'info': struct {\n",
      "        LDAF: float64, \n",
      "        AVGPOST: float64, \n",
      "        RSQ: float64, \n",
      "        ERATE: float64, \n",
      "        THETA: float64, \n",
      "        CIEND: array<int32>, \n",
      "        CIPOS: array<int32>, \n",
      "        END: int32, \n",
      "        HOMLEN: array<int32>, \n",
      "        HOMSEQ: array<str>, \n",
      "        SVLEN: int32, \n",
      "        SVTYPE: str, \n",
      "        AC: array<int32>, \n",
      "        AN: int32, \n",
      "        AA: str, \n",
      "        AF: array<float64>, \n",
      "        AMR_AF: float64, \n",
      "        ASN_AF: float64, \n",
      "        AFR_AF: float64, \n",
      "        EUR_AF: float64, \n",
      "        VT: str, \n",
      "        SNPSOURCE: array<str>\n",
      "    } \n",
      "----------------------------------------\n",
      "Entry fields:\n",
      "    'GT': call \n",
      "    'DS': float64 \n",
      "    'GL': array<float64> \n",
      "    'PL': array<int32> \n",
      "----------------------------------------\n",
      "Column key: ['s']\n",
      "Row key: ['locus', 'alleles']\n",
      "Partition key: ['locus']\n",
      "----------------------------------------\n"
     ]
    }
   ],
   "source": [
    "vds.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['locus', 'alleles', 'rsid', 'qual', 'filters', 'info']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(vds.row)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### See the data in more detail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------+------------+-------------+-------------+----------+-------------+\n",
      "| locus         | alleles    | rsid        |        qual | filters  |   info.LDAF |\n",
      "+---------------+------------+-------------+-------------+----------+-------------+\n",
      "| locus<GRCh37> | array<str> | str         |     float64 | set<str> |     float64 |\n",
      "+---------------+------------+-------------+-------------+----------+-------------+\n",
      "| 1:10583       | [\"G\",\"A\"]  | rs58108140  | 1.00000e+02 | []       | 2.32700e-01 |\n",
      "| 1:10611       | [\"C\",\"G\"]  | rs189107123 | 1.00000e+02 | []       | 4.79000e-02 |\n",
      "| 1:13302       | [\"C\",\"T\"]  | rs180734498 | 1.00000e+02 | []       | 1.57300e-01 |\n",
      "+---------------+------------+-------------+-------------+----------+-------------+\n",
      "\n",
      "+--------------+-------------+-------------+-------------+--------------+\n",
      "| info.AVGPOST |    info.RSQ |  info.ERATE |  info.THETA | info.CIEND   |\n",
      "+--------------+-------------+-------------+-------------+--------------+\n",
      "|      float64 |     float64 |     float64 |     float64 | array<int32> |\n",
      "+--------------+-------------+-------------+-------------+--------------+\n",
      "|  7.70700e-01 | 4.31900e-01 | 1.61000e-02 | 4.60000e-03 | null         |\n",
      "|  9.33000e-01 | 3.47500e-01 | 4.80000e-03 | 7.70000e-03 | null         |\n",
      "|  8.89500e-01 | 6.28100e-01 | 5.80000e-03 | 4.80000e-03 | null         |\n",
      "+--------------+-------------+-------------+-------------+--------------+\n",
      "\n",
      "+--------------+----------+--------------+-------------+------------+-------------+\n",
      "| info.CIPOS   | info.END | info.HOMLEN  | info.HOMSEQ | info.SVLEN | info.SVTYPE |\n",
      "+--------------+----------+--------------+-------------+------------+-------------+\n",
      "| array<int32> |    int32 | array<int32> | array<str>  |      int32 | str         |\n",
      "+--------------+----------+--------------+-------------+------------+-------------+\n",
      "| null         |       NA | null         | null        |         NA | NA          |\n",
      "| null         |       NA | null         | null        |         NA | NA          |\n",
      "| null         |       NA | null         | null        |         NA | NA          |\n",
      "+--------------+----------+--------------+-------------+------------+-------------+\n",
      "\n",
      "+--------------+---------+---------+----------------+-------------+-------------+\n",
      "| info.AC      | info.AN | info.AA | info.AF        | info.AMR_AF | info.ASN_AF |\n",
      "+--------------+---------+---------+----------------+-------------+-------------+\n",
      "| array<int32> |   int32 | str     | array<float64> |     float64 |     float64 |\n",
      "+--------------+---------+---------+----------------+-------------+-------------+\n",
      "| [314]        |    2184 | NA      | [0.14]         | 1.70000e-01 | 1.30000e-01 |\n",
      "| [41]         |    2184 | NA      | [0.02]         | 3.00000e-02 | 1.00000e-02 |\n",
      "| [249]        |    2184 | NA      | [0.11]         | 8.00000e-02 | 2.00000e-02 |\n",
      "+--------------+---------+---------+----------------+-------------+-------------+\n",
      "\n",
      "+-------------+-------------+---------+----------------+\n",
      "| info.AFR_AF | info.EUR_AF | info.VT | info.SNPSOURCE |\n",
      "+-------------+-------------+---------+----------------+\n",
      "|     float64 |     float64 | str     | array<str>     |\n",
      "+-------------+-------------+---------+----------------+\n",
      "| 4.00000e-02 | 2.10000e-01 | SNP     | [\"LOWCOV\"]     |\n",
      "| 1.00000e-02 | 2.00000e-02 | SNP     | [\"LOWCOV\"]     |\n",
      "| 2.10000e-01 | 1.40000e-01 | SNP     | [\"LOWCOV\"]     |\n",
      "+-------------+-------------+---------+----------------+\n",
      "showing top 3 rows\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# see the variants\n",
    "vds.rows().show(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------+\n",
      "| s       |\n",
      "+---------+\n",
      "| str     |\n",
      "+---------+\n",
      "| HG00096 |\n",
      "| HG00097 |\n",
      "| HG00099 |\n",
      "| HG00100 |\n",
      "| HG00101 |\n",
      "+---------+\n",
      "showing top 5 rows\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Display the samples\n",
    "vds.s.show(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = vds.s.take(1091)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find additional details regarding variant calls and genotype information here: \n",
    "https://gatkforums.broadinstitute.org/gatk/discussion/1268/what-is-a-vcf-and-how-should-i-interpret-it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Struct(GT=0|0, DS=0.2, GL=[-0.18, -0.47, -2.42], PL=None),\n",
       " Struct(GT=0|0, DS=0.15, GL=[-0.24, -0.44, -1.16], PL=None),\n",
       " Struct(GT=0|0, DS=0.15, GL=[-0.15, -0.54, -3.12], PL=None)]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# see the genotypes \n",
    "vds.entry.take(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding column fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Super_Population* will be the phenotoype of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-06-27 17:16:11 Hail: INFO: Reading table with no type imputation\n",
      "  Loading column 'Sample' as type 'str' (type not specified)\n",
      "  Loading column 'Super_Population' as type 'str' (type not specified)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Import only the phenotype Super Population\n",
    "table = (hl.import_table('gs://ines-python/ines-project/tab.csv', delimiter=',').key_by('Sample'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------\n",
      "Global fields:\n",
      "    None\n",
      "----------------------------------------\n",
      "Row fields:\n",
      "    'Sample': str \n",
      "    'Super_Population': str \n",
      "----------------------------------------\n",
      "Key: ['Sample']\n",
      "----------------------------------------\n"
     ]
    }
   ],
   "source": [
    "table.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------+------------------+\n",
      "| Sample  | Super_Population |\n",
      "+---------+------------------+\n",
      "| str     | str              |\n",
      "+---------+------------------+\n",
      "| HG00096 | EUR              |\n",
      "| HG00097 | EUR              |\n",
      "| HG00098 | EUR              |\n",
      "| HG00099 | EUR              |\n",
      "| HG00100 | EUR              |\n",
      "+---------+------------------+\n",
      "showing top 5 rows\n",
      "\n"
     ]
    }
   ],
   "source": [
    "table.show(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Annotate the VDS\n",
    "vds = vds.annotate_cols(**table[vds.s])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "struct{s: str, Super_Population: str}\n"
     ]
    }
   ],
   "source": [
    "print(vds.col.dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'AFR': 1018, 'EAS': 617, 'AMR': 535, 'SAS': 661, 'EUR': 669}\n"
     ]
    }
   ],
   "source": [
    "print(table.aggregate(agg.counter(table.Super_Population)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No module named 'dask'\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sample</th>\n",
       "      <th>Family_ID</th>\n",
       "      <th>Population</th>\n",
       "      <th>Population_Description</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Relationship</th>\n",
       "      <th>Unexpected_Parent_Child</th>\n",
       "      <th>Non_Paternity</th>\n",
       "      <th>Siblings</th>\n",
       "      <th>Grandparents</th>\n",
       "      <th>...</th>\n",
       "      <th>In_Final_Phase_Variant_Calling</th>\n",
       "      <th>Has_Omni_Genotypes</th>\n",
       "      <th>Has_Axiom_Genotypes</th>\n",
       "      <th>Has_Affy_6_0_Genotypes</th>\n",
       "      <th>Has_Exome_LOF_Genotypes</th>\n",
       "      <th>EBV_Coverage</th>\n",
       "      <th>DNA_Source_from_Coriell</th>\n",
       "      <th>Has_Sequence_from_Blood_in_Index</th>\n",
       "      <th>Super_Population</th>\n",
       "      <th>Super_Population_Description</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>HG00096</td>\n",
       "      <td>HG00096</td>\n",
       "      <td>GBR</td>\n",
       "      <td>British in England and Scotland</td>\n",
       "      <td>male</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>20.31</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>EUR</td>\n",
       "      <td>European</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HG00097</td>\n",
       "      <td>HG00097</td>\n",
       "      <td>GBR</td>\n",
       "      <td>British in England and Scotland</td>\n",
       "      <td>female</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>169.49</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>EUR</td>\n",
       "      <td>European</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HG00098</td>\n",
       "      <td>HG00098</td>\n",
       "      <td>GBR</td>\n",
       "      <td>British in England and Scotland</td>\n",
       "      <td>male</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>EUR</td>\n",
       "      <td>European</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>HG00099</td>\n",
       "      <td>HG00099</td>\n",
       "      <td>GBR</td>\n",
       "      <td>British in England and Scotland</td>\n",
       "      <td>female</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>23.04</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>EUR</td>\n",
       "      <td>European</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>HG00100</td>\n",
       "      <td>HG00100</td>\n",
       "      <td>GBR</td>\n",
       "      <td>British in England and Scotland</td>\n",
       "      <td>female</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>116.22</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>EUR</td>\n",
       "      <td>European</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 62 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    Sample Family_ID Population           Population_Description  Gender  \\\n",
       "0  HG00096   HG00096        GBR  British in England and Scotland    male   \n",
       "1  HG00097   HG00097        GBR  British in England and Scotland  female   \n",
       "2  HG00098   HG00098        GBR  British in England and Scotland    male   \n",
       "3  HG00099   HG00099        GBR  British in England and Scotland  female   \n",
       "4  HG00100   HG00100        GBR  British in England and Scotland  female   \n",
       "\n",
       "  Relationship Unexpected_Parent_Child Non_Paternity Siblings Grandparents  \\\n",
       "0          NaN                     NaN           NaN      NaN          NaN   \n",
       "1          NaN                     NaN           NaN      NaN          NaN   \n",
       "2          NaN                     NaN           NaN      NaN          NaN   \n",
       "3          NaN                     NaN           NaN      NaN          NaN   \n",
       "4          NaN                     NaN           NaN      NaN          NaN   \n",
       "\n",
       "              ...              In_Final_Phase_Variant_Calling  \\\n",
       "0             ...                                         1.0   \n",
       "1             ...                                         1.0   \n",
       "2             ...                                         NaN   \n",
       "3             ...                                         1.0   \n",
       "4             ...                                         1.0   \n",
       "\n",
       "  Has_Omni_Genotypes Has_Axiom_Genotypes Has_Affy_6_0_Genotypes  \\\n",
       "0                1.0                 NaN                    NaN   \n",
       "1                1.0                 NaN                    NaN   \n",
       "2                1.0                 NaN                    NaN   \n",
       "3                1.0                 NaN                    NaN   \n",
       "4                1.0                 NaN                    NaN   \n",
       "\n",
       "   Has_Exome_LOF_Genotypes EBV_Coverage DNA_Source_from_Coriell  \\\n",
       "0                      1.0        20.31                     NaN   \n",
       "1                      1.0       169.49                     NaN   \n",
       "2                      1.0          NaN                     NaN   \n",
       "3                      1.0        23.04                     NaN   \n",
       "4                      1.0       116.22                     NaN   \n",
       "\n",
       "   Has_Sequence_from_Blood_in_Index Super_Population  \\\n",
       "0                               NaN              EUR   \n",
       "1                               NaN              EUR   \n",
       "2                               NaN              EUR   \n",
       "3                               NaN              EUR   \n",
       "4                               NaN              EUR   \n",
       "\n",
       "  Super_Population_Description  \n",
       "0                     European  \n",
       "1                     European  \n",
       "2                     European  \n",
       "3                     European  \n",
       "4                     European  \n",
       "\n",
       "[5 rows x 62 columns]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import gcsfs\n",
    "fs = gcsfs.GCSFileSystem(project='avl-hail-ines')\n",
    "import pandas as pd\n",
    "#import all the sample informations\n",
    "with fs.open('1k-genome/1000-genomes/other/sample_info/sample_info.csv',\"rt\") as f:\n",
    "    data = pd.read_csv(f)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3500, 62)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['EUR', 'EAS', 'AMR', 'SAS', 'AFR'], dtype=object)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Super_population phenotype\n",
    "data.Super_Population.unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Legend for the phenotype Super_Population : \n",
    "   \n",
    "  \n",
    "  - AMR = Americas\n",
    "  - EUR = European\n",
    "  - AFR = African\n",
    "  - EAS = East Asian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We perform a PCA on a subset of the matrix : \n",
    "\n",
    "    - 100 samples \n",
    "    - 1 000 000 variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://hail.is/docs/devel/methods/genetics.html?highlight=pca#hail.methods.hwe_normalized_pca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1092"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import random\n",
    "indices = list(range(vds.count_cols()))\n",
    "random.shuffle(indices)\n",
    "len(indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000000, 500)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vds_reduct = vds.head(1000000).choose_cols(indices[1:501])\n",
    "vds_reduct.count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-06-27 17:21:26 Hail: INFO: hwe_normalized_pca: running PCA using 832650 variants.\n",
      "2018-06-27 17:23:05 Hail: INFO: pca: running PCA with 3 components...\n"
     ]
    }
   ],
   "source": [
    "eigenvalues, scores, loadings = hl.hwe_normalized_pca(vds_reduct.GT, k=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_pca = scores.to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[20.339878787405357, 9.524022815998423, 3.1753242869345644]\n"
     ]
    }
   ],
   "source": [
    "print(eigenvalues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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>s</th>\n",
       "      <th>scores</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NA19664</td>\n",
       "      <td>[-0.06723728812215236, -0.023979997814612047, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HG01191</td>\n",
       "      <td>[-0.03942432525221255, 0.10336143590935146, -0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HG00306</td>\n",
       "      <td>[-0.1062335298196058, 0.13428074906955415, 0.0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>HG01357</td>\n",
       "      <td>[-0.1234620889816882, 0.03570673777181793, 0.0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NA20127</td>\n",
       "      <td>[0.37091644999752693, -0.012523534110084923, 0...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         s                                             scores\n",
       "0  NA19664  [-0.06723728812215236, -0.023979997814612047, ...\n",
       "1  HG01191  [-0.03942432525221255, 0.10336143590935146, -0...\n",
       "2  HG00306  [-0.1062335298196058, 0.13428074906955415, 0.0...\n",
       "3  HG01357  [-0.1234620889816882, 0.03570673777181793, 0.0...\n",
       "4  NA20127  [0.37091644999752693, -0.012523534110084923, 0..."
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_pca.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 2)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_pca.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Merge the result of the PCA and the data about sample\n",
    "tab_merge = pd.merge(result_pca, data, left_on = 's', right_on = 'Sample')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "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>s</th>\n",
       "      <th>scores</th>\n",
       "      <th>Sample</th>\n",
       "      <th>Family_ID</th>\n",
       "      <th>Population</th>\n",
       "      <th>Population_Description</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Relationship</th>\n",
       "      <th>Unexpected_Parent_Child</th>\n",
       "      <th>Non_Paternity</th>\n",
       "      <th>...</th>\n",
       "      <th>In_Final_Phase_Variant_Calling</th>\n",
       "      <th>Has_Omni_Genotypes</th>\n",
       "      <th>Has_Axiom_Genotypes</th>\n",
       "      <th>Has_Affy_6_0_Genotypes</th>\n",
       "      <th>Has_Exome_LOF_Genotypes</th>\n",
       "      <th>EBV_Coverage</th>\n",
       "      <th>DNA_Source_from_Coriell</th>\n",
       "      <th>Has_Sequence_from_Blood_in_Index</th>\n",
       "      <th>Super_Population</th>\n",
       "      <th>Super_Population_Description</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NA19664</td>\n",
       "      <td>[-0.06723728812215236, -0.023979997814612047, ...</td>\n",
       "      <td>NA19664</td>\n",
       "      <td>m012</td>\n",
       "      <td>MXL</td>\n",
       "      <td>Mexican Ancestry in Los Angeles, California</td>\n",
       "      <td>male</td>\n",
       "      <td>father</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>20.24</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>AMR</td>\n",
       "      <td>American</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HG01191</td>\n",
       "      <td>[-0.03942432525221255, 0.10336143590935146, -0...</td>\n",
       "      <td>HG01191</td>\n",
       "      <td>PR40</td>\n",
       "      <td>PUR</td>\n",
       "      <td>Puerto Rican in Puerto Rico</td>\n",
       "      <td>female</td>\n",
       "      <td>mother</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>AMR</td>\n",
       "      <td>American</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HG00306</td>\n",
       "      <td>[-0.1062335298196058, 0.13428074906955415, 0.0...</td>\n",
       "      <td>HG00306</td>\n",
       "      <td>HG00306</td>\n",
       "      <td>FIN</td>\n",
       "      <td>Finnish in Finland</td>\n",
       "      <td>female</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>668.99</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>EUR</td>\n",
       "      <td>European</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>HG01357</td>\n",
       "      <td>[-0.1234620889816882, 0.03570673777181793, 0.0...</td>\n",
       "      <td>HG01357</td>\n",
       "      <td>CLM31</td>\n",
       "      <td>CLM</td>\n",
       "      <td>Colombian in Medellin, Colombia</td>\n",
       "      <td>female</td>\n",
       "      <td>mother</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>117.15</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>AMR</td>\n",
       "      <td>American</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NA20127</td>\n",
       "      <td>[0.37091644999752693, -0.012523534110084923, 0...</td>\n",
       "      <td>NA20127</td>\n",
       "      <td>2446</td>\n",
       "      <td>ASW</td>\n",
       "      <td>African Ancestry in Southwest US</td>\n",
       "      <td>female</td>\n",
       "      <td>mother</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>23.03</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>AFR</td>\n",
       "      <td>African</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 64 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         s                                             scores   Sample  \\\n",
       "0  NA19664  [-0.06723728812215236, -0.023979997814612047, ...  NA19664   \n",
       "1  HG01191  [-0.03942432525221255, 0.10336143590935146, -0...  HG01191   \n",
       "2  HG00306  [-0.1062335298196058, 0.13428074906955415, 0.0...  HG00306   \n",
       "3  HG01357  [-0.1234620889816882, 0.03570673777181793, 0.0...  HG01357   \n",
       "4  NA20127  [0.37091644999752693, -0.012523534110084923, 0...  NA20127   \n",
       "\n",
       "  Family_ID Population                       Population_Description  Gender  \\\n",
       "0      m012        MXL  Mexican Ancestry in Los Angeles, California    male   \n",
       "1      PR40        PUR                  Puerto Rican in Puerto Rico  female   \n",
       "2   HG00306        FIN                           Finnish in Finland  female   \n",
       "3     CLM31        CLM              Colombian in Medellin, Colombia  female   \n",
       "4      2446        ASW             African Ancestry in Southwest US  female   \n",
       "\n",
       "  Relationship Unexpected_Parent_Child Non_Paternity  \\\n",
       "0       father                     NaN           NaN   \n",
       "1       mother                     NaN           NaN   \n",
       "2          NaN                     NaN           NaN   \n",
       "3       mother                     NaN           NaN   \n",
       "4       mother                     NaN           NaN   \n",
       "\n",
       "               ...              In_Final_Phase_Variant_Calling  \\\n",
       "0              ...                                         1.0   \n",
       "1              ...                                         1.0   \n",
       "2              ...                                         1.0   \n",
       "3              ...                                         1.0   \n",
       "4              ...                                         1.0   \n",
       "\n",
       "  Has_Omni_Genotypes Has_Axiom_Genotypes Has_Affy_6_0_Genotypes  \\\n",
       "0                1.0                 1.0                    NaN   \n",
       "1                1.0                 NaN                    NaN   \n",
       "2                1.0                 NaN                    NaN   \n",
       "3                1.0                 NaN                    NaN   \n",
       "4                1.0                 1.0                    NaN   \n",
       "\n",
       "  Has_Exome_LOF_Genotypes EBV_Coverage  DNA_Source_from_Coriell  \\\n",
       "0                     1.0        20.24                      NaN   \n",
       "1                     1.0         0.00                      NaN   \n",
       "2                     1.0       668.99                      NaN   \n",
       "3                     1.0       117.15                      NaN   \n",
       "4                     1.0        23.03                      NaN   \n",
       "\n",
       "  Has_Sequence_from_Blood_in_Index Super_Population  \\\n",
       "0                              NaN              AMR   \n",
       "1                              NaN              AMR   \n",
       "2                              NaN              EUR   \n",
       "3                              NaN              AMR   \n",
       "4                              NaN              AFR   \n",
       "\n",
       "   Super_Population_Description  \n",
       "0                      American  \n",
       "1                      American  \n",
       "2                      European  \n",
       "3                      American  \n",
       "4                       African  \n",
       "\n",
       "[5 rows x 64 columns]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab_merge.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 64)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab_merge.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['AMR', 'EUR', 'AFR', 'EAS'], dtype=object)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab_merge.Super_Population.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd, numpy as np, matplotlib.pyplot as plt, seaborn as sns\n",
    "%matplotlib inline\n",
    "from itertools import cycle\n",
    "cycol = cycle('bgrcmk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scatter(table, phenotype='Super_Population', dimension=2):\n",
    "    \"\"\"Scatter the result of a PCA with two components.\n",
    "    Input:\n",
    "    - table: pd.DataFrame with n rows (ie the number of points/samples)\n",
    "    - phenotype : phenotype present in the table\n",
    "    - dimension : number of dimension for the PCA (2 or 3)\n",
    "    \"\"\"\n",
    "    for pop in (table[phenotype].unique()): \n",
    "        PCA1 = []\n",
    "        PCA2 = []\n",
    "        for i in (table[table[phenotype]==pop].index):\n",
    "            PCA1.append(table.scores[i][0:dimension][0])\n",
    "            PCA2.append(table.scores[i][0:dimension][1])\n",
    "        plt.scatter(PCA1,PCA2, label=pop, color = next(cycol))\n",
    "    plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+QFPXdJ/D3Z5ddYFnU0BAUkdlFMCKiGPcx9VR5JhY+9XBcovH0hLuFh0TNRjF1nhfvgtkkmlxt0FjGUBXBh4pGHnduQa3nBCNnSlHuHgPxcbFQThEVapddBdxdAhGXX+5+74+e3p2d7e7pme6enu7v+1XVxUxPT/e32Z5Pf/v7U5RSICIivVREnQAiIio9Bn8iIg0x+BMRaYjBn4hIQwz+REQaYvAnItIQgz8RkYYY/ImINMTgT0SkoTFRJ8DJ5MmTVV1dXdTJICKKlZ07d/Yqpabk265sg39dXR3a29ujTgYRUayISKeX7VjsQ0SkIQZ/IiINMfgTEWmobMv8iYiKdebMGXR3d+PkyZNRJyU048aNw/Tp01FVVVXU9xn8iShxuru7MXHiRNTV1UFEok5O4JRS6OvrQ3d3N+rr64vaB4t9iChxTp48CcMwEhn4AUBEYBiGrycbBn8iSqSkBn6L3/Nj8HeQ3p1G3W/qUPHzCtT9pg7p3emok0REFBgG/xzp3WlM/tVkLP3npeg81gkFhc5jnWh6oYk3ACIqyPPPPw8Rwfvvvw8A6OjogIjgJz/5ydA2vb29qKqqwg9+8AMAwAMPPIDzzz8f8+fPxyWXXIK2trZQ0sbgnyW9O42mF5rQd6Jv1Gf9Z/rRvLWZTwRE5FlbWxuuvvrqEQG8vr4eL7744tD7Z599FnPnzh3xvXvuuQe7du3Cpk2b8P3vfx9nzpwJPG3aBX+34N28tRn9Z/odv2s9AfCJgChZ0mmgrg6oqDD/TQfwkz5+/Dhef/11PPHEE9iwYcPQ+pqaGsyZM2do+JqNGzfilltusd3H7NmzUVNTg7/85S/+E5RDq+Bv5eydgnfnsfxDYuTeHKwnAiKKp3QaaGoCOjsBpcx/m5r83wA2bdqEhQsX4qKLLoJhGNi5c+fQZ0uWLMGGDRvQ1dWFyspKTJs2zXYfb731FmbPno0vf/nL/hJjQ6vgb5ezzw7elVJZ1H47j3Uy908UU83NQH/OA39/v7nej7a2NixZsgSAGeyzi34WLlyIl19+GRs2bMDixYtHfffRRx/F3Llz8bWvfQ3NfhPiQKtOXgeOHXBdP6AGit73rZtuBQA0zmsseh9EVHoH7MOC43ovjhw5gldffRW7d++GiGBgYAAigrvuugsAUF1djSuvvBKPPPII3nvvPWzevHnE9++55x7ce++92Lx5M2677Tbs27cP48aNKz5BNrTK+c84e4br+tTZqaL3fXrgNIt/iGJohn1YcFzvxXPPPYdly5ahs7MTHR0d6OrqQn19Pbq6uoa2+eEPf4iHHnoIkyZNctzP9ddfj4aGBqxfv774xDjQKvi3LGhBTVXNiHU1VTVoWdCC9O40evt7fe3f6cmCiMpXSwtQMzIsoKbGXF+strY23HjjjSPW3XTTTVi1atXQ+7lz52L58uV59/Wzn/0Mv/71rzE4OFh8gmyIUirQHQaloaFBhTGZS3p3Gs1bm3Hg2AHMOHsGWhaYf+GmF5pcW/p4YYw30Pvf/d1AiMi/PXv2YM6cOZ63T6fNMv4DB8wcf0sL0BiDEly78xSRnUqphnzf1arMHzDL5HPL5Sf/arLvwA8Afb1mM7G4XDhEZGps1O83q1Wxj5307rRtp66i1BwJrJkYEVGYtA/+gVbSHjNriIJoJkZEFCbtg7+Xjl2enK4Btg7XEPlpJkZEFDatg396dxoCn8O+KgCfG8AL64Ddw4WGM2aE02WciCgIWgf/5q3NUPDZ2kkBeGn1iMBfVQUsWpS/yzhvDkQUFS2DvzVscyBFPhUAvtUEzBuO3GedBWzZ4t5lPKzxRIioPFRWVmL+/PlDy4MPPggAqKurQ2/vcJPwbdu24Zvf/CYA4KmnnsKUKVMwf/58XHzxxXj00UdDS592TT3Tu9O4ddOtOD1wOridVvcDC+8eyv0fOWIudqy6ALfxRHRrckaUROPHj8euXbsK/t7ixYvx29/+Fn19ffjKV76Cm2++GRdccEHg6dMu59+8tTnYwG+p6QP+22RgXhozZgBOPbatLuNhjCdCRMUpx3k6DMPArFmzcPDgwVD2r13wD20IBgEwoQ+4vgmz/n0an302epOqquEu42GMJwKwHoGoUPmGei/WiRMnRhT7bNy4saDvHzhwACdPnsRll13mKx1OtAv+ToO7BaaqH9sqm3Ha5uHirLOGi3TsxhMRMSuKi8V6BKLC5RvqvVhWsY+1WEM32028nr1u48aNuOyyyzBr1iysWLEi8NE8LdoF/5YFLaiurA71GAMT7J8ususBGhuB5cvNgG9RCli7FqitBSZPLjz3Hta45ERJlm+o96AZhjFiZq4jR45g8uTJQ+8XL16Md955B9u3b8fKlStx6NChUNKhXfBvnNeI2664LdRjVH7uMHR0zuotW8yAn+vzz4G+vsJz76xHICpcvqHeg/aNb3wDTz/9NABgYGAAra2tuPbaa0dt19DQgGXLlmH16tWhpEO74A8AWz7cEtq+q6UG414fPRas3RCxXoOy19x7WPUIREnmNtS7H7ll/itXrgQA/PSnP8VHH32Eyy+/HFdccQVmzZqFpUuX2u7jRz/6EX7/+9/jM7tKRJ+0G9IZACp+XuG/c1c2a1fHUuYQD7tHttU0DGD16tFNOOvqzJy9FyJAvuG8rTL/7KKfmhpg3To2HyW9FDyks81Q73GYlY9DOhdoxtkzghvTB4KKt+7A4AtrHLeorbUPvi0to4O1Ey+5d+sYcRyXnChKdkO9J52WxT4tC1pQEcCpG+MNGNuedg38gHPxTmOjmSs3DPfjFDKrUGMj0NFhPiV0dDDwE5E9LYM/AAzC/5RotdW1OPJ/8kdXt1x7YyPQ2wu0tgKplFm8YxjmImKuC6PYhv0BiPSmZbFPUGP4dx47kHdMUK+5drsiG7t6giDk1g1YLYqy00FEyaZlzj+w8v6jM2ybalrscu1OOe5SdtBifwAi0jLnXymVGFAD3r9gBfjsbH7O5C3ZDMMsysnlluMu5UBv7A9ARFrm/AsK/JZ/bgWOpgAl5r85k7dkcxrR0y3AlzIgsz8AUWk8//zzEBG8//77AICOjg6MHz9+RPv/06dPl3QoZ4uWwT91dqqwL/QbZqD/TQfw80HzX4fADzgHUbcA7yUgB1VJazeuUCEtiojIm7a2Nlx99dVoa2sbWnfhhReOGPOnutocbmbx4sXYtWsX/vSnP6GlpQVdXV2hpi2Q4C8iC0Vkr4h8JCIrbT6/RkTeEpEvROTmII7ph12PPldjTnnetLraPoim02bQtmO1x3cLyEHWCVhNTK3WRWG1KCKKjRCavx0/fhyvv/46nnjiCWzYsMHz98IeytniO/iLSCWAxwD8WwCXAPiPInJJzmYHAHwHwP/0e7wgNM5rxLpvrUPq7BQEgglVE9y/UH3c874nThy9zgrcAw6lTYsW5Q/IQVfSsj8AUUZIrS02bdqEhQsX4qKLLoJhGNi5cycAYN++fUNFPnfdddeo74U9lLMliJz/VQA+UkrtV0qdBrABwA3ZGyilOpRS7wABNK4PSOO8RnT8lw4M3j+IyTWT838hR0WFGaCBkSNz9vWNvm7sAne2LZmhhrIDckuL+T0rI+I0DAQraYl8Cqn5W1tbG5YsWQIAWLJkyVDRT3axz2OPPTa0famGcrYEEfzPB5BdONWdWRcbeYdu7R/dBdcaZ8cwRo/MmXvd5AvQuZ/bZURshgAHUPpKWnYOo8QJobXFkSNH8Oqrr+L2229HXV0dHn74YTzzzDNwG0utVEM5W8qqwldEmkSkXUTae3p6SnZc16FbFYCX7IdU7ew0c/p2sq8bpykdh46fc3i7jIhSo28AImYaShWEOVkMJVIIzd+ee+45LFu2DJ2dnejo6EBXVxfq6+s9VeKGPZSzJYjg/zGA7NmFp2fWFUwptU4p1aCUapgyZUoASfOmZUELxKmvrtXSp0BKeQvKhQz1rBRQWTnyPVC6IJzYzmF8nNFbCM3f2tracOONN45Yd9NNN2HVqlWevh/mUM5DlFK+FpgdxfYDqAdQDeBtAHMdtn0KwM1e9nvllVeqsLW+06pSj6aUPCCq9pe1Cg9g5PLjGoV5rcoMs8UtNTXun7e2jk5XKmW/rYj7vlKpcP+/nI4vEu5xQ9XaOvqPVFNj/4eh2HjvvfcK+0Jrq/kDEjH/jcnf3+48AbQrDzHWd85fKfUFgB8A+COAPQCeUUq9KyK/EJHrAUBE/kZEugH8BwD/KCLv+j2uX7mTNh8/fRzVldUwxhsQCIwxKRjbnTtyZTOM4crfXP39I3Pr2VIp56Ge7eb3zTf1QtiVv4nsHJbYxxkqiIbN3wIp81dKbVFKXaSUulAp1ZJZ9zOl1ObM6zeVUtOVUhOUUoZSam4Qx/XDbtLm0wOnUVtdi8H7B9Hb3IHVtzeOCsK5amrMAdg6OpwrZQcGRgfz6mrg+HH7kga7Zp9e5tzJV7fgVxI6h+WW8KhOjnVBeiqrCt9S8jJpc74mmrlt8Z1ywNZ22UM2K+U+T29uRsTpyaKU4t45zK7CukuS+DhDlJ+2wd/LpM1umb9Uyvy8uXk4aLvljLODeW0tcObMyO3ylTTY7TuX05hCfmXnlpubzbTE8enY7ma+UrWgX2L+OENUBG2Dv5dJm90yf9m5x1tvNQOk15yx002ls9O5wUn2vp2EkVm1yy0vXQpMnhy/RjF2/+9taMT3VIwfZ4iK5aVWOIql1K19Uo+mVOs7I2v47RqCOC2G4f24huGtlZDV4CC3IcKdd3proBJEAwanlkdxbBTjdC5ht5Ki0iu4tU9M+WntE3mQd1pKEfy9yA2gbgHb6/6qqrzdUKyAbRfoFyxQqrLSfF9Zad4Qco8TRAvGqJuXBomtOvVRDsG/oqJCXX755UPLqlWrhj7r6elRY8aMUWvXrh3xnSeeeEJdeumlat68eWru3Lnq+eefdz0Gg78PheaO/Qb/fDeQ3PbzXtv85waxoHK5+dIbtzb+MW3OTQUqh+A/YcIEx8/WrFmjrr76anXNNdcMrevq6lIzZ85UR48eVUop9dlnn6n9+/e7HiPSdv5xVsxwBcboYX5GrM/XWbSQFoQzZrj39s3mdTyhQlsw5qtojlujmECac7NHcOKkDx9G3Y4dqNi2DXU7diB9+HCox2tra8MjjzyCjz/+GN3d3QCATz/9FBMnTkRtbS0AoLa2FvX19aGlQevgX0z/ntWrgaqqkeuqqsz1Xm4mTsEyt4+A1eCkkOCaHdiD6pBlVTTb3fRi0yjGa7D2sh0HOEqc9OHDaNq7F52nTkEB6Dx1Ck179/q+AZw4cWLEjF0bN24EAHR1deHgwYO46qqrcMsttwytv/zyyzF16lTU19fju9/9Ll544QW/p+bOy+NBFEspin2KHa7AqeggX1FLa6t9ZW9NjVlmb7dPu3Jqp3RnF+mEUb4dyyITr/8RXrdjrXEsFFLsk9q+XeG110Ytqe3bfaXBqdjn4YcfVj/+8Y+VUkq9/fbbKjvWDQ4OqjfeeEP98pe/VBdeeKG6//77XY/BMv8iefkdFxLw3G4mTi2HDCN/EI2ytU/seQ3WXrdL5ABHyVNI8BebwI/XXlPy2mu+0uAU/L/61a+qadOmqVQqpVKplKqqqlIffPDBqO3efPNNdemll7oeg2X+RSp26sQVK+xLB9yKWpx6C9fW5i93zi2nXrPGW38CDYcrGc1r5YfX7RI5wJHeZowdW9B6Pz744AMcP34cH3/8MTo6OtDR0YH77rsPbW1t+OSTT/DWW28Nbbtr1y6kQuzaPya0PcdA9hSJ1iTqVm9cALj7bvs6gbVrh99bNwTA/G5T08jvWDeTZcvs01DsEDKNjZoG80LNmGE/DVpusPa6ndsfmWKpZeZMNO3di35rhiYANRUVaJk509d+rTJ/y8KFCzF+/HjboZ4XL16M5cuX495778Unn3yCcePGYcqUKXj88cd9pcGVl8eDKJao2/m3tto/3edr7+5U1OLUsauQzmFUhKDL/K1ttS9PK2+FNvVsPXRIpbZvV5Ip6289dCiklAWLZf4hKKQ9vpci33IP/uUYzwJLk9cdleN/AhWlHNr5lwKDfwjy9WwttLFHOdcTlmPPVz9pYgwnBv/8MVbrCl83hdTf5Rb52jUXL+d6wnKcz6TYNAXSDJ+duBLBjIPJ5fv8vNwholiizvl7HdQtt6mmU47Va/PMKJTjU0mxafLdDL8cH4OoYPv371c9PT1qcHAw6qSEYnBwUPX09NgO/wCPOX9RZXp3bGhoUO3t7SU5Vjpt3+Ind/2iRcCWLfYtgyx1dfaNRlIpc3unlkVRcktzR0epU2MqNk0VFWbEziViNnkN7cBUVs6cOYPu7m6cPHky6qSEZty4cZg+fTqqcoYcEJGdSqmGvDvwcoeIYinlwG5OGT23smOnz8oxF51POWZ2i02T75x/HP+ARFnACl9vnIKFYbjfFJw+i2vv/3KsJC0mTb5vZHH9AwatHC8I8oTB36NCW/UYhnt8uPNO+89yx9u38DcWPF//p+X4GFRq/D+INQZ/jwptz++2uI2/b5dx5G+sTOl+R+bTT6x5Df7aV/haTQOzmxWK2FcaWiorgYGB0eutSd3tvmtX4ehUt1hZCaxfXx6VwaQh37XmFCWvFb5at/O3WvP095sBFzADeL774cCA84BwhbTndxrXZ2CAQ8QnTpz6DpRzpxQKjLbBP7szEDAc0FtazBuAG2sUTbtRNZ1GCl20yHvHLyD6TlYUoLhNAJNvuNtixekGqAMvZUNRLGGX+bsVa7p18PJSJu91/H279WxdmEBxLEMPut6DFVwlA1b4usvXnDu72WZl5cgbQ6Hy3Wis/ZdrbNC9/tM3HfoO5PvBxPEGGFMM/nmU8lr0cqMp10xROactNpIe+Lw8KutwAywTXoO/tmX+bsWaQRdN5qs/syZJzzczVxTKcdC32AmrDL1cOE1TBwxfLKxELj9e7hBRLKVo529XnBHWxOdh5p7DLJZhhi0g5V525id9+XpKOk1izUfIUIDFPsUJ6wk9rN9+2L+ppJdYkPJ/EeXrKZlvmjsKlNfgr30nr1xx698S9iCUdp3gamrKp1iKAuD3IrK7SCy8WEqOnbyKFLeiSaeOYsVODJ+rnOsjKCB+L6LsiwQY2WOylBcL+xEUhDn/HHHL6XL4efItCRdR3H64IWLOv0iF5nSLzWwElUlJekMSKoFSX0Rh5NDZLK1wXioGoliinsbRi2LryYKupGU9GvlWqosorBYKbi2NNAO29ime199BsS1h2IKGYieom0MYF39rq3Pw1/BHxeDvkddxeOyu9WLbwLPtPMVKkLn1MC5+w3AO/k5pLHaauBg8YjP4e2B3TTtdm3YZCOb8SQtBXrBhXPxufQzsFHMzi1EntZIGfwALAewF8BGAlTafjwWwMfP5GwDq8u2zFMG/kFm87DIm5VLmTxSqIHPrQV78bpNmuwX/Ym5AMcqxlSz4A6gEsA/ATADVAN4GcEnONisAPJ55vQTAxnz7LUXwL2T+Xqe/cbFPgjF5giQKPvAFcfG7DSZnLYZh/123H75TumJUVlvK4P+3AP6Y9f4+APflbPNHAH+beT0GQC8yfQyclihz/rl/Z+bKSWtRPqo63Sjy5firqwtvqeH2w09gzj+Idv7nA+jKet+dWWe7jVLqCwDHABgBHNsXp+bNd9wxsrOi1VyYHQZJS1F1816xAli2zH4GNLfex4YBPPmkc/rsfvh2E3dn9xNIYIeasurkJSJNItIuIu09PT2hH8/pml6zZvhvbU3UXu4z7xGFqrHR7O07OGj+G3bgT6eBxx93Dshu462cOJF//+PHD782jNHHsVg3GStYGFl51v5+4O67h4NC3IaX8PJ44LYgpsU++YodY/SUR5Q8bsU6TkNEe/mh3nnn6OIdp6n0cvfT2qpUVZV9EZPbnKx2ASbESj+UsMx/DID9AOoxXOE7N2ebuzCywveZfPsNM/h7KcKMUf0OUfLkq5RVyr1zl13QdZtRzO0mY3G7IbndQHIDTMh1KCUL/uaxsAjABzBb/TRn1v0CwPWZ1+MAPAuzqee/ApiZb59hBn8vw4879Rthzp+oBNwqZb0GZK+Vtm5Ldg690O/a3SBEnG8U1uc+nwS8Bn8tR/V0GrM/W1WVWQ9w+vTwOk0HCSQqPbtROkXM1hhr1gyvW7HCvm4gmzU6qZcffrnwEWw4qqcLL2PznzkDTJzIceyJImHXGuPpp0cG/nQaWL8+f0C3Km3LdVIOOyUYkVTL4N/SAlRX59/uyJHSNnCIStwaKVACeLno8rUwcps4PpsV9L3+8MtFUDMyOdAy+Dc2mrn6fCoqkh8Iradru6bURKEI6qLzEhxz2+JbbbfLQWWl+VRjzXyWK+QnFS2DP2Dm6vMZGEh+IOQcGFRyQV10+YKjCLB8+fATQ3NzeQX/c84x6zDOOWf0ZyXoQKZt8Pd6Uw06EJZbEUvYcwATjRLURWfX6zabUsCWLcXvP2i5RU59fcDatea/2QyjJBWM2gb/fNdNtqCumXIsYonbhPWUAEFddLkTx9vJ/vFGfVFnNx10U1tbkgpGbYO/XWMCw2G0oaCumXIsYkngkCVU7oK86KxKYacbwKRJw4/ax48Xtm+ngBC2zs7SHMdLZ4AoliimcQx78MJy7TXM4aWp5IK+6Ox+vNXVo4djqKgorIOX345dxSy5HdkKBM7kVZwwAyHHC6LYiVPOIDetbtM75lusuQCiCP4+g4LX4K9lD9+o2HVaZK9hKltxv2CL7dE7Zgzwve+ZlcWlKoLJJWL2byjqq+zhG7nclj1ANMOiExWlHCupCuFUWWcYzq09DMMM/OvXRxf4gZJ0MmLwD4lTyx5Aj17DlABxbwfsVLG8evXoXFhrq/lD7e01c/xeeg6HqRSdjLyUDUWxBFXmH1WRJcv3KfaScBEXEwCCGMEzwrJ/lHAax7KVr119mB2u4p5pIkpEO+BiZiCLuj9AthADRqKDv1uRZdgdrth5imIvqrl7o+Y0xy9g1gkYhvneMMzK4WxeBo6rqBg9no+1/1whBoxEB3+33HfYdVlJyDQRlXzu3nLgNJy0VSfQ22v+f/T2Ak89NXK7J5907nCWSpn7GBgwK5Szv3fHHaUPGF7KhqJYgijzdyuyLEWHqzg1kSaigBTbWzSggAG283dvptzcbN+Sy5r0h4ioaOm0GWQOHDCLblpaSvbU5LWd/5h8G8RZ9kiudn8DuxsDi2WIyLfGxrIvIkt08Afc/wbjxw8Hf8Mwm/+W+d+LiCgQiQ/+duyKg06ciC49RESllujWPk7i3mudiMgvLYM/O2ARke60DP7sgEVEutMy+LMDFhHpTqvgb43ls2yZ2dLH6qWtS691IiKLNq19clv49PWZuf2nn2bQJyL9aJPzZwsfIqJh2gR/tvAhIhqmTfBnCx8iomHaBH+28CEiGqZN8Nd1XgoiIjvatPYBYjHQHhFRSWiT8yciomEM/kREGmLwJyLSEIM/EZGGGPyJiDTE4E9EpCEGfyIiDfkK/iIySUReFpEPM/9+yWG7l0TkqIj8wc/xiIgoGH5z/isBbFVKzQawNfPezsMAlvk8FhERBcRv8L8BwPrM6/UAvm23kVJqK4DPfB6LiIgC4jf4T1VKHcy8PgRgqp+diUiTiLSLSHtPT4/PpBERkZO8Y/uIyCsAzrX5aMQ0KEopJSLKT2KUUusArAOAhoYGX/siIiJneYO/Uuo6p89E5LCInKeUOigi5wH4NNDUERFRKPwW+2wGsDzzejmATT73R0REJeA3+D8I4O9E5EMA12XeQ0QaROR31kYi8i8AngWwQES6ReTvfR6XiIh88DWev1KqD8ACm/XtAG7Pev9v/ByHiIiCxR6+REQaYvAnItIQgz8RkYYY/ImINMTgT0SkIQZ/IiINMfgTEWmIwZ+ISEMM/kREGmLwJyLSEIM/EZGGGPyJiDTE4E9EpCEGfyIiDTH4ExFpiMGfiEhDDP5ERBpi8Cci0hCDPxGRhhj8iYg0xOBPRKQhBn8iIg0x+BMRaYjBn4hIQwz+REQaYvAnItIQgz8RkYYY/ImINMTgT0SkIQZ/IiINaRH804cPo27HDlRs24a6HTuQPnw46iQREUVqTNQJCFv68GE07d2L/sFBAEDnqVNo2rsXANA4dWqUSSMiikzic/7N+/cPBX5L/+Aglu/ZM+IJgE8HRKSTxOf8D5w6Zbt+AMDSPXvwD3v24NpzzsGOv/6VTwdEpI3E5/xnjB3r+vkggK1Hj9o+HTTv3x9iyoiIopP44N8ycyZqKoo7TaenBiKiuEt08E8fPmxb5u9VvqcGIqK4SmyZf24rn0JVi6Bl5syAU0VEVB4Sm/P3k+MHgCqwspeIkstX8BeRSSLysoh8mPn3SzbbzBeRHSLyroi8IyKL/RzTK7/l9Z8rxeaeRJRYfnP+KwFsVUrNBrA18z5XP4B/UErNBbAQwG9E5Byfx80riPJ6tvYhoqTyG/xvALA+83o9gG/nbqCU+kAp9WHm9ScAPgUwxedx8/LTysfSydY+RJRQfoP/VKXUwczrQwBcC8lF5CoA1QD2OXzeJCLtItLe09PjK2GNU6di+bnn+toHABb9EFEi5Q3+IvKKiPw/m+WG7O2UUgqActnPeQCeBvBdpZRtTaxSap1SqkEp1TBliv+Hgy19fb73waIfIkqivE09lVLXOX0mIodF5Dyl1MFMcP/UYbuzALwIoFkp9eeiU1ugIDppsaMXESWR32KfzQCWZ14vB7ApdwMRqQbwvwD8k1LqOZ/HK0iNiO99sKMXESWR3+D/IIC/E5EPAVyXeQ8RaRCR32W2uQXANQC+IyK7Mst8n8f15HPlWArl2SLDCCAlRETlxVcPX6VUH4AFNuvbAdyeed0KoNXPcaIURL0BEVG5SWwP36CwzJ+IkijRwb9ZlSMSAAAFDUlEQVS2stL3PiYFsA8ionKT6OD/+EUX+d9JAJXGRETlJtHBHwDG+AzeR774IqCUEBGVj0QH/+b9+/GFzxY/bOpJREmU6ODvt7K2CuCY/kSUSIkO/n5y7RNE8Ps5czimPxElUqKDfzEje1YCaJ0zB8e//nUGfiJKrMRO4wgMz8S1dM8ez9/5xjnnMOgTUeIlOucPFD4V47ajR0NKCRFR+Uh88C/UQNQJICIqAS2CvzHGe+kW+/MSkQ60CP6rZ89GtcfOXk3TpoWcGiKi6GkR/BunTsWTF1+M1NixEACpsWPROmcO7pw2bSinXwngzmnTsCaIISGIiMqcqADGvA9DQ0ODam9vjzoZRESxIiI7lVIN+bbTIudPREQjMfgTEWmIwZ+ISEMM/kREGmLwJyLSEIM/EZGGGPyJiDRUtu38RaQHQGfU6XAwGUBv1IkIWNLOKWnnAyTvnJJ2PkB5nFNKKTUl30ZlG/zLmYi0e+lEESdJO6eknQ+QvHNK2vkA8TonFvsQEWmIwZ+ISEMM/sVZF3UCQpC0c0ra+QDJO6eknQ8Qo3NimT8RkYaY8yci0hCDvwciMklEXhaRDzP/fslhu5dE5KiI/KHUafRCRBaKyF4R+UhEVtp8PlZENmY+f0NE6kqfysJ4OKdrROQtEflCRG6OIo2F8HA+/1VE3hORd0Rkq4ikokhnITyc0x0isltEdonI6yJySRTpLES+c8ra7iYRUSJSfi2AlFJc8iwAfgVgZeb1SgAPOWy3AMC3APwh6jTbpK0SwD4AMwFUA3gbwCU526wA8Hjm9RIAG6NOdwDnVAfgMgD/BODmqNMcwPlcC6Am8/rOhPyNzsp6fT2Al6JOt99zymw3EcD/BfBnAA1Rpzt3Yc7fmxsArM+8Xg/g23YbKaW2AvisVIkq0FUAPlJK7VdKnQawAeZ5Zcs+z+cALBDxOP9lNPKek1KqQyn1DoDBKBJYIC/n85pSqj/z9s8Appc4jYXyck5/zXo7AUC5V0R6+S0BwP8A8BCAk6VMnFcM/t5MVUodzLw+BGBqlIkp0vkAurLed2fW2W6jlPoCwDEARklSVxwv5xQnhZ7PbQD+d6gp8s/TOYnIXSKyD+ZT9n8uUdqKlfecROSrAC5QSr1YyoQVYkzUCSgXIvIKgHNtPmrOfqOUUiJS7jkTSjgRWQqgAcDXo05LEJRSjwF4TET+E4CfAFgecZKKJiIVAH4N4DsRJ8UVg3+GUuo6p89E5LCInKeUOigi5wH4tIRJC8rHAC7Iej89s85um24RGQPgbAB9pUleUbycU5x4Oh8RuQ5mpuTrSqlTJUpbsQr9G20AsDbUFPmX75wmArgUwLZMqem5ADaLyPVKqbKZmJzFPt5sxnBOZDmATRGmpVhvApgtIvUiUg2zQndzzjbZ53kzgFdVpuaqTHk5pzjJez4icgWAfwRwvVIqDpkQL+c0O+vtvwPwYQnTVwzXc1JKHVNKTVZK1Sml6mDWzZRV4AfA1j5eFpjl3lthXpSvAJiUWd8A4HdZ2/0LgB4AJ2CWA/591GnPOY9FAD6A2VKhObPuFzAvTAAYB+BZAB8B+FcAM6NOcwDn9DeZv8XnMJ9i3o06zT7P5xUAhwHsyiybo05zAOe0GsC7mfN5DcDcqNPs95xytt2GMmztwx6+REQaYrEPEZGGGPyJiDTE4E9EpCEGfyIiDTH4ExFpiMGfiEhDDP5ERBpi8Cci0tD/B/kTquKBUSHCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter(tab_merge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scatter3d(table, phenotype='Super_Population'):\n",
    "    \"\"\"Scatter the result of a PCA.\n",
    "    Input:\n",
    "    - table: pd.DataFrame with n rows (ie the number of points/samples)\n",
    "    - phenotype : phenotype present in the table\n",
    "    \"\"\"\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    for pop in (table[phenotype].unique()): \n",
    "        PCA1 = []\n",
    "        PCA2 = []\n",
    "        PCA3 = []\n",
    "        for i in (table[table[phenotype]==pop].index):\n",
    "            PCA1.append(table.scores[i][0:3][0])\n",
    "            PCA2.append(table.scores[i][0:3][1])\n",
    "            PCA3.append(table.scores[i][0:3][2])\n",
    "        ax.scatter(PCA1,PCA2, PCA3, label=pop)\n",
    "\n",
    "    plt.legend()\n",
    "   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXeYZGWZ9/85oXLnMN093dPTcSLD5GEQdo3ggAoSFpAVs6uurvuyrgK6IgZecXd90Z/4W3cXXTANwX0VlJVVEZbMzAAjaQKdc6zqyumE94+ac6aquqq6qrt6EvW5rrl6uuqpc546fc733Od+7iDouk6JEiVKlDj5iCd7AiVKlChRIkFJkEuUKFHiFKEkyCVKlChxilAS5BIlSpQ4RSgJcokSJUqcIpQEuUSJEiVOEUqCXKJEiRKnCCVBLlGiRIlThJIglyhRosQpglzg+FJaX4kSJUoUjpDPoJKFXKJEiRKnCCVBLlGiRIlThJIglyhRosQpQkmQS5QoUeIUodBFvRIlSpQomHg8zsjICJFI5GRPZVmx2+20tLRgsVgW9XmhwHrIpSiLEiVKFEx/fz/l5eXU1tYiCHkFHJx26LrO7Owsfr+f9vb29LdLURYlSpQ4NYhEIme0GAMIgkBtbe2SngJKglyiRIkTwpksxgZL/Y4lQS5RokSJU4SSIJcoUeINw69+9SsEQeDw4cMADAwMIAgC//AP/2COmZmZwWKx8JnPfAaAW265hebmZrZs2cKGDRvYu3fvss2vJMglSpQ4JZkNRPnT8ByzgWjRtrl3717OP//8FFFtb2/noYceMn+///772bhxY8rnrr/+eg4ePMgDDzzAJz7xCeLxeNHmlExJkEuUKHHK8cDBUc771h95/53Pcd63/siDB0eXvM1AIMCTTz7JD3/4Q+655x7zdafTyfr16zlw4AAA9957L1dddVXGbXR3d+N0OvF4PEueTyZKglyiqOi6jqIoKIpCgSGVJUoACcv4hv98iUhcwx9ViMQ1vvCfLy3ZUn7ggQfYs2cPa9asoba2lueff95875prruGee+5heHgYSZJYuXJlxm288MILdHd3s2LFiiXNJRslQS5RFDRNIx6PE4lECAQCBAIBvF4vXq+XQCBAOBwmFouhKAqappXEukRWRjxhLGKqNFlEkRFPeEnb3bt3L9dccw2QEOBkt8WePXv4/e9/zz333MPVV18977O33347Gzdu5JxzzuFLX/rSkuaRi1KmXokloWkaiqKgqiqQCPsRRRFBEEzRVVUVRVFSPmeMkyQJWZYRRdH83BshPKpEdlqqHcQ1LeW1uKbRUu1Y9Dbdbjd//OMfefnllxEEAVVVEQSBT3/60wBYrVa2b9/Ot7/9bV577TUefPDBlM9ff/31/P3f/z0PPvggH/3oR+nt7cVuty96PtkoWcglCkbXdTRNIxqNEo1GzZNbTLNqDHE1hDf5nyAIaJpGLBYjFAoRCATw+Xz4fD78fj+hUKhkUb9BqS2z8Y9XnI3dIlJuk7FbRP7xirOpLbMtepu/+MUvuO666xgcHGRgYIDh4WHa29sZHh42x3zuc5/jW9/6FjU1NVm3c8kll7Bjxw7uvvvuRc8lFyULuUTeGEJsiCSwaIs22+eMfaiqSiwW4/Dhw6xbtw5gnqiXLOozl0u2NHNeVx0jnjAt1Y4liTEk3BU33HBDymtXXHEF3/zmN83fN27cOC+6IhM333wz1157LR//+MfnGSFLpVTLosSC6LpOPB4nGo0iy7IpgNkENRaLFU0k9+/fz86dO9F1PeVf8vaTBdr4aYh1iVODQ4cOsX79+pM9jRNClu+a18lYspBLZEXXddP/Ozs7y9TUFOvXr88qdIqiMDQ0xPj4OA6HA5fLhcvlwul04nQ6l2RN5GNRpwt1uqukJNQlTnVKglxiHslCbIicLCdOlUxiFovFGBwcZGpqysxoisfjBINBgsEg09PThEIhIFGe0BBqQ6yXS6gNyz7dYi8JdYlTlZIglzAxYoiTrU1DLI1FuGQikQgDAwPMzs6yevVqzj33XARBIB6PY7FYcDqd1NfXp2w/HA6bQj0zM0M4HEbTNBwOB06nM0WoJUla9HfJJdSQsObj8TiDg4PU1dXhcrlKQl3ipFMS5BIpyRxAxogJURRNQQ6FQvT39+Pz+Whra2Pt2rWmaOVakxAEwXRfpAt1JBIxhdrtdhMKhdA0jUgkQm9vb1GFOvlnNBpN+b6GUKd/Jj08z4gUKVGimJQE+Q2Mpmn4/f6UULRsIiMIArFYjJdeeolwOEx7ezsbNmwoiigJgoDD4cDhcFBXV2e+rus6+/bto7KyklAohMfjIRgMomkaNpstxfXhcrmKKtTJcwDmuXB0Xc9pUZfEusRiKAnyG5DkZI6jR4/S1tZGRUVF1vFer5ejR4/i9/vZvHkzNTU1OQUnfXFtsRiWa7JIG9uPRqOmRT06OkooFEJV1XlC7XQ6Tf93trkuNIfkn+mfy5X0IstySahLFERJkN8gJC9yJccQi6KYUZR0Xcfj8dDX14coirS2tjI6Okptbe2Jnvo8BEHAbrdjt9tT5mOE3BlCPTY2RjAYRFVVrFbrPIs610JlPnPI9FnjWBqJMwBHjhwx3Tql7MSThyRJbNq0yfz9mmuu4cYbb6StrY0DBw6YN/7HHnuMf/7nf+Y3v/kNd911F5///Odpbm4mEonwiU98guuvv37Z5lgS5DOchZI5kn3DxviZmRn6+vqw2+2sXbuW8vJyIpEIIyMjBe33RIuMIAjYbDZsNltKtlW6UI+Pj5tCHYvF0HWdqqoqysrKcDqdi25Qacwh+SdAMBhEkqR5SS/JlJJeMhCcgblBqFoNrrqFxy+Aw+Hg4MGDBX/u6quv5o477mB2dpa1a9dy5ZVXsmrVqiXPJxMlQT5DSRdi4+JOv8ANQdZ1ncnJSbMZ5aZNm3A6nfPGLYSxj1Mp1TmXUL/yyitUV1ejaRoTExMEg0EURcFiscyzqJci1MY8Co2lzuSjfkNEfrx8PzzwNyBZQI3DpXfApitP6pRqa2vp6upifHy8JMgl8sOwBkOhEHa7PasQGwiCwNTUFEeOHKG6upqtW7dmLJpyqolsMTBcCJWVlbhcrpT3ki3qqakpAoEAiqIgy7JpSRtCbbValzyPUtJLEsGZhBgr4cQ/gAc+Ax1vWZKlHA6H2bJli/n7TTfdlLGyWzaGhoaIRCKcffbZi57DQpQE+QwhOZnD5/MxODiY88RRVZWRkRFGR0eprq5mx44dOYUlXwv5TMFqtWK1Wqmurk55PT3hZWBggHg8jizLJ9SiTk966e3tpbOz88wQ6rnBhGWsJJXblCyJ15cgyNlcFpmOTfJr9957L48//jiHDx/mjjvuWJYqbwYlQT7NyZTMIctyVvFMTm9uampi1apVlJeXL2jlnYkW8mKwWCxUVVVRVVWV8no8HicUCpkJL4ODg8RiMSRJIhKJMDw8nGJRL0UgMwm11+vNGUt9Wgl11eqEmyIZNZ54fRmora3F4/GYi3putzslssfwIR84cIALL7yQSy65hMbGxmWZS0mQT1NyJXNksmaT05tbWlrYvXs3kiTR39+ft2/4TBTkYi0+WiwWKisrqaysTHk9Go1y8OBBJElidnaWoaEhU6jTLepiCHXyT4P07MRkDGHOFKJ30nDVJXzGD3wm1YdchIW9TLzlLW/hJz/5CV/72tdQVZWf/vSnvPe97503bseOHVx33XV897vfTakSV0xKgnyaYfgUk4U420IdJNKb+/v78Xg8tLa2cu6556Zk4RWyWFeicCRJwmKxzGsJpCiKaVG73W6Gh4eJRqNIkpTin3a5XNhstgXjvnNxWia9bLoy4TMuYpRFug95z5493HbbbXz5y1/mU5/6FJs3b0bXdfbs2cP73//+jNu44YYb2LZtG1/84hcpLy9f8pzSKQnyaYKmaQSDwZTY2WwXhiiKxONxXnnlFfx+P+3t7axbty7j+Deab/hEk80Cl2WZioqKeQk5qqqaPmqPx8Po6CiRSARRFFOSXVwul7loq2naogo0FZr0MjMzQ3l5uZlss+xC7aorqlVsdLVJp7Kykp///OcZ3/vQhz7Ehz70IfP3lStXMjExUbQ5pVMS5FOY9GSO/fv386Y3vSnnie/3++nr6yMQCNDV1cXGjRtzjjfE+43KcrthCnWJSJKUVagNi9rr9TI2NmYKtd1uJxaLMTMzkyLUiyWbULvdbhwOR0rSS/JnSkkvS6ckyKcguZI5sp3cXq+X3t5eNE2jvb2dQCCQUsAnG8W2kCORCIODg4iiSFlZWV6P3Ceb5ZxbsXzUkiRRXl4+7zFZVVW8Xi89PT14vV7Gx8cJh8OIophSk9rlcuFwOJY0F03TTKs4nYWSXjRNM+PhDU7lc+JkURLkU4h8kzmSx7vdbvr6+pBlmc7OTnNRKd+TPVvqdKGEw2H6+/vxer00NzebqdcjIyOmb7TYi1inA8udsShJklmYqbOz03xd0zTTovb7/UxMTBCJRADm+ajtdnteLg9VVbMWcFooltqYU/pncv18I1IS5FOATAXhc53gANPT0/T39+NwOFi3bt2iFxgMH+RiCYVC9PX14ff76ejoYP369Sk3FANjESsQCKREG6TH75aVlS05frcQTjWXxWLI5EM2nlDKysrmjQ2FQubfYnJyknA4Ee+byaJO3m4uQc5G8nmczVed/jdQVTVjq7A3glCXBPkkkkmIF7JUxsfHGRwcpKKiYl5682JYrMsiGAzS19dHMBiko6Mjt686HsISmqFSkKhobABRBl1H8A6hBmYIalF8uMybTHLqciwWw+v1phQDKjbLeaGn35iWax/5LurlEurk5gFTU1PzhDoajRIOh7O6LQol23ExnqjShfqNYFGXBPkkkJzMsW/fPnbt2pXzBNc0zaxcNjc3lzW9eTEUKsiBQIC+vj7C4TAdHR3U1dXlviAic8g9/w1qFEHT0V31RCUX8dcewh4aw17biq2sicqV29HWvMn8mJG6PD09nVIMKLlqm+GjXkod5OXGCCE71feRHMWRjNEkwChzOjw8TCgUQtf1eRb1UttxJZPLTXcmC3VJkE8gmZI5soXiQOLRbXh4mNHRUVasWEFVVRUdHR1FTd3MV5D9fj+hUIjXXnuNzs7OBWsio+sIoRmk/j9CLIg4fRhpdD8e9yyPqRuJYEXAwu6JZ2lzxZEGHycemERd+y6QrGbqstVqNUtXpldtS6+DbAh0IZ1FltulcKJcFsu1D1EUzS4vAwMDbNy4ESi8HVexKFSoJycnaWxsNG8Uv/rVr7jssss4dOgQ69atY2BggPXr17N27Vrzc/v27ePnP//5CS27aVAS5BNArmQOYwU6WTzi8ThDQ0NMTEywcuVKzjnnHGRZ5sUXXyzIms1HDBYSZJ/PR29vL4qiYLPZ2LlzZ85tCjNHkTxDSENPYz36X+AfR1CP1STQ4Y/K+czp5VSKIWxaiF+GNnF2pJ9Vc5N0CveB1YXadeHx7R0TYuOYZavaZhSsDwQCKS2gltOKy4eT5UNebgShsHZcVVVVZvRHrqp1ycfLHXEzFhhjZdlKauw15CLbMfZ4PDQ1NZm/7927l/PPP5+9e/fy1a9+FYDOzs6MNS5OZNlNg5IgLyPJnTkge1adsVgSi8UYGBhgenqaVatWmenNBpIk5bSok0kWslxki7JIDqPr7OykurqaZ599dt7NI2Vbr/wn4jPfRZ7pARSS96zr8IiymZ+o7ySOiKypSGg4iCIJGkeUVqZGJ3hT5R9SBDnf75qtYH26FWd0v072i4ZCISwWy7KI2pkqyNkQhMztuIaGhrDZbGb4WywWMw2BTAL92/7fcsuzt2ARLcS1OLfsvoWL2i9a0twCgQBPPvkkjz76KO95z3tMQV6IE1F206AkyEUmW2eObBelLMuEQiF6e3vxeDysXr2arq6ujBdYIf5eY+xCF2p6lMXc3By9vb1AwnJILqKTsZ5FyI0wegDx+R8h9T8GJLY1o1cwotVjFRQ6hDF61JXcr74ZPzZUZGJYiGDBQZRV6iTt4iT74u1sDnhTTsql1NDIZsUlL2BNTEyYER+QGhJWVla25CSLN5og58IQ3XQMkTYMmNnwLLc8ewtRNUpUTSSg3PLsLexq3EWto/CONcbxf+CBB9izZw9r1qyhtraW559/ntraWnp7e82U6vPOO4/vf//7KZ8/EWU3DUqCXCSSY4gPHTpk+j1zXYzBYBCfz8err75KZ2dn1vRmg0Is5HzF2xjn8Xjo7e1FkiS6urrmFcnJtE1h/CDyA3+F4B0BLeGO0XX4nbaNh9TdlOthysUQ1bqffdpaDtNKDIkQDqzEiWAD4A/sokWbRtMk1KiTtYem6ax3sqrakdd3LZTkBayxsTHWr1+P1WrNGLtrPGY7nc4UH3W+yS5niiAvZ3p9slCrqsq0ZxqLaDHFGEAWZPpn+7FX21OyAnMVQ0o/9nv37uVv//ZvgUT7pr179/KZz3wmq8viRJbdNCgJ8hLJlMzhdrtzXiBGenM0GsXlctHV1TWvnGMmDPdGPhi+6YXm7vV68Xg8iKJotmvKtX/TWg3NYNl7BUT9SduDW2PX8Gv9TUSwEsOCRVVQkZBRENAJ4UBAJ4YFBxF0BGTiDNKAioW+YRv2iSO0Vju5YH09a9O8I7quE1d1LFLxU3KzhYQlpy2nJ7ukC3V6ssuZJMgnIppF13VWlq0krqWm8yu6QkddBy67K8WiNp5Ek8NGjXNf13Vzzm63mz/+8Y+8/PLL5mK6IAh8+tOfzjqXE1l206AkyIskVzJHtgtwbm6Ovr4+NE2jo6ODmpoaDh06lLf1kY/IGuQSb13XmZ2dpa+vD4vFgtPpTKmClQ1bdAb55XsQHeVIz/8HRP3MaS5e15uJaRL3q3/GbzkXFRHQ0RCJYkEHZGQkNEAlRqL2sgCIaARwEseKVRbRdYgq0D8b5qleN9Sp7Dx2E5gJxNh7YIxJf4Rym8z7djTTWrM0Kzofd0i2tOV8kl1UVV324k0nQpAXkxSyGHRdp9pWzS27b5nnQzYW9rK5PoxrUtM0fD4fs7OzxONxjhw5wi9/+Uve97738cMf/tAc/+Y3v5nh4eGM8wiFQszOzvLyyy/T3Nw8r+zmxMQEMzMzCEKi/nhbWxs2m23J378kyAVSaDJHrvRmKMwNsVSXRXIDU4fDwYYNG3A4HOzfvz/r3J/tc/NEzyyWyAzveu1r2OURBEFgOibzb8q1PKpuIYrEHGUEcAE6Ahr6sVNLQEdEQ0VERTr2vo6MQgwbumTDKouocY2oAoIATglEAUa9UUZsie+rajo/3TeKPxqnwm5hyBPitt/18OWLummqXNqj5GIt2GwV24xi9YFAgLm5OTN+fDn69MGZY4UbCILARe0XcU7TOXlHWRifM5KH6uvrcblczMzM0NLSwn/913/x+c9/PmX8FVdckbGusXHNVlVVsXHjRg4dOsRnP/tZzjvvPLPsptPpZP369UiSxNTUFCMjIymp64ulJMh5kqkzR7YT1Hgkmp2dNdOb169fP+9RGJbHL5w+Vtd1pqen6evro6ysLCXDL7nWQDr7Bjz8+LlhqhwW9Nf38f9F380N6r20i+P8u/I+/qR2EEFmlkqiWI4JsYjOcUtKPybJmDEXAlbi2AWw2eyEVAGbRUTTIRw/VvNAB6soIAsCspiYYyimMhuM4bBKPNXrRtF0oorGP/6+l1vetYZKx4lLt16I5GL1VquVQCBAe3t7Sgz15OSk2VA1PdnFKG+ZL2eahWzcXGrsNXkJcTaMOVssFh5//PF573/2s5/ls5/97LzXg8EgV111FWvWrEnM41iIZXLZzeSbcFlZGW63e9HzTKYkyAug6zrBYJBQKER5eXleFrGqqjz33HNUVVUtmN68nBayqqpMTEzQ399PRUUFZ5999ry55LKsnu51U+mQsVskXo5UMKa38DfxT/MZ+QEm9SqmqGaGKmJY0BEAHREVzRRknYQQG/80BAR0wYZgk2iocjDqDdNcYWdsLkwwdsyNIQqsKLdR5bDQfaylncMiIksCh8YTPmuXTUISBIJRlWf6POzZuCKv45LOiaxlkalP30LJLulCnUkUzzRBLtZ3MWpiFEosFktpaWa1WgkGg1nHz8zMzHtKWiwlQc5CcjKHz+djcnKSDRs25BxvpJZqmsZZZ52V10JdoSKbT+1iI/72lVdeoa6uji1btuBwLOxrfWnEyzN9bhxWiXesq8ciC4TjGk/3TuPRVqKiE0fkW8pVVAkhZqkghnxMjBP/DGtYTOgzhu0tCaDrCS+yDiiqxoAnTCiq8qdwwJyDRQSbLNBSbWfn6iraxCl0XccqiVy+uZGbf3OUmKqj6Rot1Q5kESZ8EUIxFad1cYKx3LUscglMPskuwWDQTFnWNA273Z6ykKiq6glZ1DsdQuuSORE3kdnZWYLBYEqm31IoCXIamZI5LBZLVtFUFMXs3tzQ0MCOHTs4dOhQ3ndmSZLmFfvONTaXeGuaxsTEBAMDAwiCQEdHBy0tLVnHq5rO/kEPz4zFGdo/zEMvT2KVRSZ9Eb7/WD8N5VamAlH8EQVNl9HQ8OEiiJMZASLICXVNwvhV0EESIXZMkTX9+Hu6DuG4jo6a/vFjLgrY3FLJVdtW8vJL0+Z7m5or+MibVvGLF8aoLbMRVzVeGvUz4Yvy6Otu/uq8Vs7rXPwj7nKwWP9urmQXIxMuEAgwMzOD2+1menqa8vLynNXalsLJcFksFcMdVChWqzWlpnMsFsvo6/f5fIyPj7N27dqiHeeSILNwMocsyymtbCB7ejMsvxsiHU3TzCpwtbW1bN++nZGRkZw3BV3X+f//p48ne2bxB1TmevuoK7OiotM/E0LTIRJXUTUNxXQxJ046HQhrxwXWKkFcTRLjYz8t0nFBNr+vABU2EX8ssV2R41Y0QESBuK5S5ZCRxPkRK+86awV2WeTxnln2DfppLLeysspBVNH41yeH6Kp30VCx+NVuXdf549FZXhnz01xl5z2bGrDJi7/Yir3glikT7siRIzQ0NJj+6vRqben1jxdTqP50FGRN0xblsjAyOKPRKBaLBbfbTUdHR8qYUCjE4OAg3d3dRS0X+4YW5FydOZJJFuRoNMrg4GDW9Ob08QtRqHgnL8AZbpKhoSHq6+vZsWOHaREstAA46A7zdK+bujIrfp9OXNMYdIcBzbRmY4qGxnGBNQTXeE0SQNUhdmz6AmC3iCiajqbpxNWE+0HTdCySiI6OXRZQFA1ZAIVUMTb2YZEERj2R468d8/MGowpf/+3rPN7jRicRidFVl6hOZpNFAlGY8EcLEuR0AfjBE4Pc/+IEmp5wxPzuaC+XnOMjqoZZXb6aLXVbkMX8L5sTFQFhxESnrxEslOyS7KPOlezyRnJZCIJAa2srR48eBaCurg6Hw8Ho6Cgul4uqqiqGh4dRVdXMarVarXR3dy95zm9IQTYW3oyYxYXihyVJIh6Pc+jQoQXTm43xy2EhG2M1TWNkZITh4WFWrFjBzp075z2apUdZ/O7QNPsG3NS5bPzF9pVM+iKAzqGJAO6Qji7oRBQNUTi2FCeALArE1IQwSSLEk9RTJ+FakAVQ9OOvOSyJOraGYJbbZUIxlTUrynhtIoAsClQ4RGYDMVwixFTN3K4sQJUNqmwag2MT9Paq5oKXxWLhm//dwyNHZrEkHNL4YxqvjPvZtqqSuKqh61DnKvwR1SASV7n/xQnssogoCqgEORJ8lt6ZVXTVV/Ga+zU0XWNXw668t3kiym/mqvaWT7JLeo++5NKmRrLL6WghL2XOlZWVbNq0KeW15uZm8//F8hmn84YS5OQY4oMHD3LWWWeldCbIRHIh9nzSmyFhIRczoy55/j6fj6effpqmpiZ27dqV9XEpebs/fW6En+0bJhhTCEZV/uXxPhoq7Ez7o6iajk0ERAFNSli7CetXQNMTLgWHVaSt1slr44EUn68GoCfHUIA3HMdplWmpsvGXO5vZ0FTGn0b9HJ0M0lHnZMQTwR2KoetQ6bAgi9BW6+TwRIAKu4VAVCEYU3jrWS1UVlqZnp5mdHSUI6/38ExP5FgyiQiigE0SCMYUZoJRQODanStprlpcTLKiKYwHZ9DkceJ2DwICgu5AEDQsghOLaKHeUU+fr69gQV5uC3kxor9QskswGExJdlEUxWxwmizUxWY5BVmSpBSRveaaa7jxxhuBRKREU1MT3/ve9/jkJz9pjvnRj37E7bffbtZ8ufXWW7n00kuLMr9MvCEEOVMyhxFFkU3QfD4ffX19xGIxOjo68Pl8NDQ05LU/SZIKclksNFZRFIaHhxkZSSRl7N69O6dvbHA2xDf+MMLBsRBWud/8znPhOOox3++MP0aty8q4N4ImgNMiEVVU0z2haToOq0RDhY1gVOHoVNAUY7sMMSUhwBrgsibiiKOKhg5saHSyvbWS8zprqCuz0r0i1TpTtITV7QnFsckiZTaJHzw5xF3PDKMDZTaZP/Z4efuGtTidTrq7u7HabFT+6Xn8c5HEPPREamytHS5fFaOxuoyWiih/6hvnu09NMuaN0b3CyZcv6mZFeXYXhq7reGNefjXwK8aCY0gNrxAKVSNrDagEsIpOmo/V1IhpMexSYYJ/uiVtZEt2MeqcCILA9PQ0AwMDxONxZFlOsaaLlewCoLjdKKOjyM3NyDWFL9amh705HI6MNSsA7r//fnbv3s3evXtNQR4ZGeHWW2/lhRdeoLKykkAgwPT0dMbPF4szWpAzFYQ3Ttxsfl6Px0NfXx8AHR0dZsyo4U/KB8PFke/YbBayoigMDQ0xNjZGc3Mz27Zt4/DhwznF+M4nB/jeo71EDD9C9Pi2bbKAekwMY6qG0yphlUXKrTphJeEZrnZKVDutTPqjbG+tZMfqan7yzAB1Np2ZCChawnXhsEoIQFRNPC7LAuiySFzV6J+NMOCO8NvXZvj+1WfNKxIkiwmBqis7bmG9OuanpdpOmU1G13UOTQZ4vMfNimPV3kRB4JPnt/LN3/Xij6oICLhsMl+9dB07WysSC1keH1/4TR+BqIoswPMDET750+f57ntWU1FeRkywMuqLU19mTbGkHxl9BH/Mj6qrrKmFbRh5AAAgAElEQVSrY8gTQA1KuOzV7Gh1MRebxK9IoMPbmt+W19/V4HQT5Gzouk5FRUVKaB4kFreNhcRsyS7pXUjywffQQ0ze/BUEWUZXFBq+/jUqLr64oG0UEg64d+9evv3tb3PttdcyMjJCS0sLU1NTlJeXm+6eTK6fYnNGCnImIc61UGfUdujv78disdDd3Z0x0Dvfi0uSJLPDbz5j090byREczc3NnHvuuabI+yMKT/bMYpVFtq6qxCIdP+H+NOLlp88NE1XSg8kSqJpuWrnHNJH6chsNljhHvDoVDguNFTYkUSCmamyql5kdG0LTobLMBVKM6UDiRtNUYWN3Rw2P98ziD8cTPmVAVaHMnjitAlGFf39qiK+9O7u/LRRT+c6j/Tw7MIcswqpqB06rhK4n3B8rklyAF25YQUu1g/95fRZJFHjn+npW1yYWsSoqKuj3g4ZEuSOxf5uu44loTPqjPNs3y3efm0scB13gyo3lvH9HE9PRaV6afYkqaxWKpmCTZVbXONnQVY1DctBe0U53VTcxNUaNvYZqW/W875CLM0WQs/ljLRYL1dXV85JdkoV6fHwcWZYJBoMIgpBSpS1TpTbF7Wby5q+gRyLm+Tr55Ztx7t69KEvZIBwOp9Rsuemmm7j66qsZHh5mfHycXbt2cdVVV3Hvvffyuc99js2bN9PQ0EB7eztvf/vbufzyy3nPe96z6P3nwxklyLk6c6QjyzLxeJzJyUn6+/txuVxZ05vh+CJZPosEi42yiMfjDA4OMjk5mTGCY9wX48tPBtDEV9F1ndU1Tv7uHV3IksCGxnKGPWFTGDN+B0HAISdC0URRIK7pvG9HC5usUxwM13LfC+PEVR1/OIoSi7OuWsTZsZ7fjx1F08EuC1Q6JFpry7jl3etY21DGI4enueOxxBPFlD+KXRbNGsYWUcAdyv2k8LXfHuWZvjkkIZE63TsToqPWgSTCuoYytNnZlGy6DU3lbGjKXJHOaZFQdRCTRFAXBFavauYLv3sRXZSwSCKCpvHLw0Ga6v7EYd8zTOlT9Gg9VFmqiGgRVEklHA0j2SU21mykxl7Dvsl9PDXxFE7ZyflN59PozK/q15kiyIXsQxAErFYrNTU1pkU9NDSEyzW/Uptx7hvirOs60eHhhGWcvE1ZTrgvChDk9Os/m8vi3nvv5aqrrgISfuWPfOQjfO5zn0OSJB5++GH279/PI488wvXXX8/zzz/PLbfckvccCuWMEGRN05ibS1g/RoxlrovACAWamJigvr6ezZs3L5jJZohsPoJcaGcPRVE4evQo09PTtLa2cu6552Y8+f/p9z14ozqypCKgc3DEy8d/+iJOq0RtmY3PX9CFIMyP7QWwSfDeLU080eNGFOHijQ1cvnUlzVUOvnrfMJOKm1qHQCAUot5p5W/ecza7O+vRdZ2rtge598AouqZT77Lwf67cxMpjj/zvWL+C3e01eMJxXhjy8H/+0IuiJvauaDrnd2S/gBRN5+neORwWAZtsQYgoROIq/qjK37+jgw1N5bwym9dhBKCr3sl5HdU82esmrmpYZIErtzYhCkLi92NPE5KYCMHrUfqps9expXULB6YPMB4Yp8pWRZ1YR020huZYMyOHRnhUe5Sh+BCN5Y2E1BC/Hvg1V3VeRaVtfs3odE7Uot5y76NYURa5itQb14xeX4+e5vLT43Gkpqa8v2shKfF79+5lYmKCn/3sZwCMjY3x+uuv093djSAI7Nq1i127dnHBBRfw4Q9/uCTI2TDqAGiaxtTUlFlNKxvJcbtWq5XVq1fT1taW174KEdl8LeRoNMrAwAChUAin05lRiMfmInjDcdrrnBwc8RJVIZo0D0XV0HSJ8bkIv311kmt3tvCvj/cTSMvI+MxbO7nrmWH8kTga8LP9IzRW2nni9Vme64ug6hEkUaBrRRl3fXSnKV6CIPDJP2/nL7Y10zM4Qq1DNMXYoMwuU2aXaamy4w0r/PS5YTRN54otjfzFtiayIQogS4loDkmESruEwypx04WdvGVNnbn/fBEEgVvetYbHjs4y4o3QWZcQaE1P+LwDUQWrJKJqOqBR4ZCQYhIO2cF5jecx5B/iglUX0FXVdfz4Kgr7X9lPnVRHNBQlGokyHZvmCf8TrK1Zu2D36+VsQJr+3ZeT5Q57M857SZJwrFxJwze+zuSXbwZZBiVOzc03o7hcxI614DKEPdn9kT7ffCz6o0ePEggEGB0dNV/7yle+wt69e/nYxz7GxMQE27ZtA+DgwYOsXr26WF85I6e1IBsnofGIlG0hLT29edeuXUxPT+edsgzFTfaIRCIMDAzgdrtpa2vD5XLNS3HWdZ1v/NcR7n9hDIskYJNF5jI8/kvmSafTNxPim+/dyLs3NfKvj/XwdO8M1eUu/vrN7fTPhgjGFERRQATiisYdj/YSiSmgg80iISAw7IlweCLApuZUH3ptmZVotT0lpTQdTYf1jeXctKeb9fUOaspyJ2iIgsBfndfKD54cJKIkhKuj1sG5aVZ1IdaOJAq8fV1d6msCfPuy9fzd/z1EXE0UOLr+7V00NIg84n6EJrWJsBKm1lHLqvLUnmmyLKPLOtPRaSpcFTTXNyOHZNatXEcNNWZBoGAwOK/ORFlZWUGP+t6ol6HAEE7ZSUdFxwkR8nw5UQuHBhUXX4xz9+6MURZGZq2RS5CcYSuKIrFYjFgsZrpAjOOY7kPes2cPDoeDyy67LGUeF154Iddddx3nnHMOX//615mZmcFut1NfX88PfvADNE2jv7+fUCiELMt0dHQUpRYynOaCDMd9uxaLxUwVNUhfHEtObzYWGfKlGIIciUTo6+tjbm6O9vZ2s83TwMDAvLGPHZ3h/x4cJ6ZoxFUIRjMLvHisPCWCwNnHRHRllYMvXNDB622qeQIefeJY+BsCOjq6phNXtGMx0wpCUkvSbAKY3n8vmZii8eG7X+C1CT+ikIhLvusDW1izIvcK+zU7VtJaY+eFYR8ryqy85+ylpSpnY1NzBb/+5A7GfVFqXVYq7DKKVseAbQBZlGl2NbOtfhs2KfXCemH6BQb8A4wFx9B0jVp7LResuoC2qjZkUZ5XEChTnYmxmTFeUF9gRpuhtbyVK7qvoK6sLkVwB3wDfOel7xBTY2hobK/fzkfXfxRRODWy405UYkjyMZFrajL6jA2XZPoNwsi8NZKJIpEIhw4dQhAEVq1aldcTrhFNcujQISwWCy0tLXR0dKS4NKemppBlmU2bNuF2u4tWCxnOAEE2MBbp4LgrYGZmJqtPthCBLXR8emJIOBymr68Pn89He3s769evX9D66ZkOElM0jGFGIUuLmFiQ045ppigI6MCmlRX8r7cdPynSU6ffvraeHzzeTySe2KbNInHNrhZeHfNzcMhDXNUQRYHmSgfrGjMvmuVqOPqLF8Z4ZcxH5HjhC7704BHu/9i2nN8T4E0dNbwpi695KU1O07FbJNprj6cWy6JMt72bnW07M47XdI2Hhx6mtayVZlczwXgQT9TDWTVnZUyfzlRn4tDhQ/zS+0uGo8PYBTsHZg7Q6+7l8vLLscpW05L+t8F/Q9MSgq/rOgemDrBrxS42120uyndfKici43CpvnAjgqOyshJRFLFYLLS2tubtaoREIphRfQ8StZDn5uZSBHlubo6VK1cCUF1dzdDQUNH8+GeMIFssFiKRCK+99hpzc3O0tbXR3d2d9SQqVJAXUwQoFArR19eH3++no6ODDRs2ZPyjGZZn8lzb65xYZZFwTEUQEqnMrdUOHEQZCQrUl9m49dL1rKpxomo6jRWpdQiS+4pNTU0x1dfHV/68hvuPxgnEVPZsbODj57cRVTRu+cU+hgIC65oq+V9v78KaxULNVR9jyB1KEWOAcV9+oX8nk1xir+oqiqZgES1YBSsO2YGiK2h6/i2Z5uJzDIeGaSxrRBAEqvQqZiOztJ7VSoOt4XjB+sAkVs3KbHAWURQJ62EGJgdos7Qh22WcVucp5cJYDopZmzrZoi/Ess+nFnI8HjfHGDcBRVGKkhBz2guyIAgEAgF6e3txu92cddZZeVmgy2khh0IhwuEwL730Ep2dnWzcuDHnfAyxTxbkt6+t5+KzGvjNyxNYxERx9juu2cxs38ucc845C1orgiAQDod59tlnqaysZOvWrdjtdi79s9RxTqvEh7dWsWLFipQyj9m2me2i2byqEscB0ez6IYuwMYulXQjFtJALxSJaWFe9jsOew9TaawkqQRyyI++QNziW5g3o6KZbSEdHFmQkSWJWmMVn9dFd102Pr4dQPERIDSHpEpOhST7yxEcIq2FqpBqubbyWdTXrTB91rvKawXiQR0YeIagE2bViF52VxXmkXgr5WJGnQh2LpbDUc/W0F+TBwUHGx8dpbW0lEonknd68GEFeaBHQuDFEIhEsFgvnnHNO3okk6WncgiBw66Ub+Pj5bXjDcbrqXbhsMvsG54t3Mka7pt7eXuLxODt27FgwpC+lm/QC47JZyHs2rODgcCIxRUBnhVPgoxtlRkdHTQFZTCnE5Wahv88lbZfgkB30entpcDRwUetFuCz5Z56ViWXsqNvB8+7nsUgW4mqc7fXbqbfX85vB3/DU+FOIgkhEiTAZniSiRJAECVmQuXv4blPIA1qAn7l/xjdWfmNeec1wOMzg4KB5nBVR4aOPfZTZyCyqpnLna3dyTfc17Fyxk821m0+KpW21WnG73dTU1GTdf7HrWCzGYs2nFrLFYjEtaWNxUZZlM8HMbl98j8dT7wopkNbWVlpaWtB13SyFlw+LEeRsi4B+v5/e3l5isRidnZ3U1NTwzDPP5H1y5XKHtNU6M45NP0nS++Zt3ryZgwcP5tUpJN9efbmsVUVRuLJTYKfDxoqmFtoaa8wFromJCQKBAKqqmlEIhng4nbkfxU/2Y7pdtvOetqVlZ/1F+1+woX4Do4FRGp2N7G7czWRkkqfHn6bB2YAoiEyGJlE1la11W+nx9jARPt6/TRREVF1lJDTC9/q+l6hn0rCbyzZchq7qPPn8k3h0DyF3iNhIjF9N/YqRwIj5eR2dHx/+MQ8OPMiO+h28e/W7ccgO1lWvM632XMe5GE8odXV1zMzM5KwFoShK1jjlQvF6vdhstoJrT+i6ztjYGNFoFEmSGB8fp76+3sxzgMT1PjU1RW1trdne7fDhwwDY7facTSEW4rQXZENMCr1w87UKDTKJps/no7e3F0VRTCFOppBU68U2OjU6Sff29lJWVpaxb16h2yxknFFvw3hKedufvQlBEIjFYtjt9pQ2VulRCIaVJwhCSjxvWVlZih9vOV0WJ6KnnizKrCpbxcNDD/P7kd/z6NijvLP1nYiCaEZRWEQLOjp9vj6CSnDeNnR04lqcQf8gTouT+3vvJ6JGqLXW8h/T/4Er6MIhO/jg2g/y1OBT6Gn5mioq0UiU3w3/jifGnkASEvHX7qgbSZC4sutKPrT2Q7ijblY4VqREmxTDcpVlmcbG3K6enp4eqqqq5l1Hi+GGG27gyiuvZPv27QV/dmBggCuuuAJVVfnIRz7Cl770JW6++WZ27NjBJZdcQiQS4brrruPFF1+kpqaGe+65Z14B+8Vy2gvyibKgki1qr9dLb28vuq6nFCBKxlhUK3ZmnzHWeDzq7e3F6XQuSogNFmMhq6rK8PAwo6OjNDc3p6R55wqbS49CMLYVCoUIBALMzs4yODhIPB7HYrGYCyiGaJ9uRdJ1XSeiRvjuy98lqkapt9czGZ7k/p77cckufDEfZZYy4lqcOnsd48FxZFFGREQ7lm9p/JQFmRpHQqxkUebX/b9G1VWsgpUKawWBeIAvP/dlgmrmJ7kAid6FOjohJYRfSTSMVXSFe47ew/2v349NsiGJEt845xtsq09EyJwof+xim5Jmwuv15tXTMhMXX3wxF6cVMvra175m/t9ut3P//fcvaX7ZOO0FOZ3lSiOVZZlQKMTzzz+PIAh0dnZSWZk9ddZYeV2OVGu3281rr72Gw+HgrLPOWlQ1rWQKsZBVVWVoaIjh4WGamppSYrsXS7bavLFYjEOHDpnlR0OhELqu43A4TLfHQp0uFuJEpDXPRGcIKkFqbYlF02prNbORWf6y+y95ZPQRpsPTtJW38dmzP8sXn/0i3piXBkcD46FxAkoAWZBpLWsloh6PWhkLjJmCChDwBvDGvSj6wm64sBqe95pK4iYv6zKxeIybnrqJr6z8ClVlVdjtdlRVJRqNmjfH5SDf6yUfvF5vRkPpVOeMEmSLxVJQ+EmmcLNMuN1uXn/9dYLBIDt27Mir5Xexu4bouo7b7WZmZoZIJFIUITbIlfCRvP+pqSmmp6dxOp05i+Mb45d64VqtVhwOBw0NDebNT9M0wuEwgUAgpdOFJEkpIn2qLCLquk65tRxN11A1FUmUUDQFHZ3Oyk621G9JOVZf3P5Fvn3w28S1OC1lLexcsZNPbPwEwXiQLzzzBWYjsynWrcFsrICiH1kQRRFJlJDERN3jlWtXUivW4na7URSFw4cPE41GkWV5XleRYoR8FdNC9vl8OQ2mU5WTf8YukZTMnmPJIfmeHIYbIlPnA0MAe3t7sVqtrF27liNHjuQlxlBcQXa73fT09GCz2airq6OxsbFoYgzZm6dC4jhMTk7S19dHRUUF1dXVefUOK6YVlewCMVoMuVyulIgaRVHMco9LWUQsNpqmUeeo4+LWi3lo8CFEIVHY6PKOyym3Jp4IkufTVdnF/979vxn0D+KSXXRVdiEIAuXWcr65+5s8PPQwDw0+hDvqNrMui4VFSFw3qq4iI1PvrDcjSkKhEOvXrwcScbhm/HRSDWSbzZbSo69QF1MxLWS/35/3tXoqcdoLcjKGhZwvmQQ52Tdrt9vZsGEDZWVlZniLgRYKo4fDCA4HonN+JEMxUq2TbwjGPHp6evJu+WR8n4UESBTFeXVAksPnKisr2bZtW6Jw/KFDee/7RCLLMlVVVQUvIiqKMi8ZoJgYx/+S9kvYWLOR6cg0jY5G2irasn6m2ladse5yla2Ka7qvIRgP0uvrLaoYX9hyIc9MPoMkJGqa3LD1BlOM09dCLBZLxmMdi8XMm6LhYtI0zWymmhw/nemcLKavWtf1kxKHvFROe0HOZCHni7FoBMejFfr6+nA6nfNcAoKQaKyp6zrK2BihJ55Ei8cQLVYc55+HXF8Puo54LAaxUAs5ed4ej4eenh4sFgvr1q1L8a0W6m/OV5CTG6IaNySXy8WWLVvM0LloNFrQzaAYLCUxJJ9FREVRePXVV80nq3S3x1IXEZOPf2dlZ1ESNC7ruIzfDP6GyfDkkrdlMBOZ4a633cVEaIKWshbqHfXme/lUThMEwUw5Tk4w0nXddDEFAoF5Xa+Tj3Wx1n9OViJRMTjtBTmZZIHNh/Qi9WVlZWzatGletIKu64SffpqKf/t3hr/zHZTpGYhEjrVdlhH37sV+1kbksnJsZ2+i7J3vLNhCjkajeDwes3dZuhAb5LsAZ2w3Hx+5sU2Px8Prr7+O3W7P6KM+mVlzxSR5EXFkZIStW7cCmBZeIBAo2iLiciwy19hr+M553+GT//NJfHFfUSxld9RNk6uJJtf8cqlLsVwFQcDpdM67ppJvim63m6GhIUKhEC+88EJK+ONi1gKSO8mfbpwRgmwIRSEiaDzOvvbaa9TW1mYtUh/rH8D/+9/j3bsXu9tNpq1rgQAhv5/K915K5MWDSLW1SI2NeVuykUiE4eFhvF4va9euzSjEBoupqbHQCR0OhxkbGyMQCLB+/fqs+y/kZlAsTuRNIL3LBSx9EXG5on50QafSVokkSMzF5hIV/BYpzHbRzobqDVnfX46wt0yRNfv372fTpk2miyl5LcDwTxvH2+l0ZjU0AoFAzmvoVOaMEGQDi8WyYHqzrutMTEzQ39+PJEm0trbOK1Kv6zqBp55i+ss3g9ttvp7zspqbI9bTg7Wri/jwMFJz84LC6fV66enpIR6PU1VVxebNC1f2yuTvzTU2l4D6/X56enqIRCJUVlam1IrNxJliISez0PdZ6iLicglyo6MRh+wgrsQpk8oIqkFU8q9qZmATbXRXdfM3m/4m65gT1SJKEISs/uloNGoKtdvtTnl6SfdPe73e03JBD84QQU62kLOlN2uaxsTEBAMDA9TU1LBt2zZmZmbmiWbouefw/PBHRJ97rrBJqCp6IIgeCiHX1iHLctZGp4YQC4JAd3c3uq4zMjKScWw6S8nqMwgGg/T09BCLxejqSnTHGB8fX/T2lpNT9bEzeRHRF/Px8uzLALS52hBjormIGAqFOHDgQEr0QXom4mKwy3Zu3nEzX3/u68TUGBYSFemiStRMJFnwOwgyn9/6eS5cdWHGcqIGxQxHy7WPbFa4IAjY7XbsdnuKf9p4egkGg/j9fiYmJvj5z3/OH/7wBwBuv/12Nm3axDnnnFOwxex2u7n66qsZGBigra2N++67b15c88GDB/nUpz6Fz+dDkiS+9KUvcfXVVxf4zVM5IwTZIJMPWdM0xsbGGBoaora2lu3bt5u1TpMLBum6judH/4H3Rz9CDwQWtX/BZkNubMSxezfBgH+e+yQ5w6+rq8uMk/T7/QW7IfIdmyyg4XCY3t5egsEgXV1d5snt9XrzztQ7GSynVb7U7zQbmeVr+7/GXCxR66DKVsXNO26mra4NSDyGb926NWsmYrLbw+l0FuQaaK9o59azbuX52ee5c/jOhKjaIBKL4I67qbHVmCnQgViAcms5k6FJNBLWaJOzibc2vzWnGAOmy2A5WYxbJPnpZcWKFQB885vf5K1vfSv33nsvNTU1PPzww2bfzEK47bbbePvb386NN97Ibbfdxm233ca3vvWtlDFOp5Mf//jHdHd3MzY2xvbt23nnO9+56AxBOEME2bioksPekvvn1dfXs2PHjnlWieFz1lWVyb/7HKHHHlv0HOQ13dR+6UtYGlYgWCxI4ZApnD6fj56ennlCbFBoREahi3pGpxKv10tXVxd1dandKk6G5Xum8ED/A3hjXmrtiZubO+Lmwf4H+fD6D5tjcmUiLmURcSQwwqG5Q0S0SMKDfMw9YpWtiPHjPRFjaoxyazkfWPsBvvfy98wElagW5ejc0QWL4J8Il4WiKEWzwmOxGN3d3Xzwgx9c9DYeeOABHjumBx/84Ad5y1veMk+Q16xZY/5/5cqVrFixgunp6ZIgG8iyTCwWM1N7jf552RJFDIs6+LvfE3rqqUXvV2hspPnnP0dK2o+Rav3iiy+iqipdXV1Z/1DFzuoz0HXd7P3V0dGRtU50sQXZ6ALudDqXbFktt996qdv2RD1YxON/d4towRP15PXZpSwi/mH6D+zt3YumaOjoOC1O/PFE9p6AwPlN59Pj7SGmxrBJNm7ecTM/PfpTqmxVOOTE4nUgHuDR0UcXFOQTUcuimPtYSh0Lg8nJSZqaEhEnjY2NTE7mDjHct2+fWe1xKZwxgqyqKuPj43g8HqqrqxdM7YXjFrIyPg4FJJSYSBIVH/sYtZ/4K4Skk8nv95vdbLds2bJgTn0x/MLJxONxBgYGmJqaorm5mc2bc9fALbTyXTaMBVOjBGg0GjXryRpCUl5ennOF/HRja91W/jTzJxyaAx2dqBpla93WRW8vn0XEV4de5a7+u7BgQUJCEzQi8QgXNF2AR/GwqW4T721/LzE1xk+O/IReby//PfTfaLqW0u1E0zWs0sK+7HyLZC2FYlrIXq93wWYLAO94xzuYmJiY9/qtt96a8vtCIXTj4+Ncd9113H333Us+r88IQZ6dneXVV1+lqakJp9OZ913KzNRbuybRbryAGGZqamj5jx9hTYrQCAQCZtTEqlWrGBsby6vASbEsZEVRGBwcZGJiwqwTXVVVVVBiyGIwkmqM8onbt29POYmTH80HBwcJHWvlnrzIlW2ha7kt5KX6kN/a/FbmonP8dui3QCJp4y3NbzHfL9bcjUXE8spypm3TOMed2CV7orYEMhE9whZhC07BieAWeMr/FP809E944h5kUabaVo1DTtw0fDEfAHbJzrtXv3vBfeeTGLJUimkhGy3TFsJY/MtEQ0MD4+PjNDU1MT4+bvqo0/H5fLzrXe/i1ltvZffu3Yues8EZIcgVFRVm1bFMd7xsGILsPO88Kq99H967f5xzvJD4EM63vZXaz30Oy7H6rslCbNRFjsfjDA8P5zWPQkQhk3gml8JsaWkxS2H29vbmvVi3WEFOTibZvHkzTqcTXddTFlezPZobYUzJC11WqzVFpE/1MDtBELi883Iu67jM/L2YBONB7u25lx5vDxE1wlR4KuESis4l0qt18KgewlqY7099n4tXX8zVnVfz1//z13jiHkREVE1lNjxLtVzNnvo9hMUwdpudd7W/i9by1gXncCJcFsW2kJfqsrjkkku4++67ufHGG7n77ru59NJL542JxWJcdtllfOADH+DKK69c0v4MzghBttlsBWXoGSSLW+3f/R1VH/4wk7d8lUjS4p5QXc2qB36FVFHBaF8fcU2j8VhxHaNlUzQapaurK0VwjPKbxSbZQtY0jZGRkaylMPNdAFyMhWy4ZQRByJlMkmuf6QtdyfUQDKH2eDxIksT09DTl5eWmUBejulgxySbESxFoXdf5zkvf4ejcUTRdYzgwjF2yU2mtRNEUpkJTyMiE9bBZ2vOXfb8krsVxx9yIgogkSCAk3BOCLLC2fi1d1i4CgQDTr08zpU8tuIh4IlwWp0otZIMbb7yRq666ih/+8IesXr2a++67D4ADBw7wgx/8gDvvvJP77ruPxx9/nNnZWe666y4A7rrrrgXj+XNxRghyMvmW1DTGJiNVV9P0ndsJPPww4SeeQKqto+rDH0I6FhUhu1xEj1W56u3tJRwOm0Kcvq1i+WXTEUURRVEYGRlhcHCQhoaGrDWJ8w2RK0SQNU3jpZdeIhqN0t3dnfXEX0xCRKZ6CAMDA9jtdlwuF36/n+npafr7+4CvcnYAACAASURBVM3qYulhY4Xu82RZ4Kqu8tDAQ+yb2odTdnJl55WsqVqTMsYddXN07igRJWL2xwtqQcJKGF3QQU/EEldZqnBanOjohKNhfnLkJ6h64u+u6IrZpqnKVsWb2t6EXU7UW4koEfp8fRAHm27Luoh4ImqYKIqypF50yRSjFnJtbS2PPPLIvNd37NjBnXfeCcD73/9+3v/+9y9pP+mcEYKcXmAoW0nNfLdVftFFlF900bz3FEVhfHycmZkZOjs7qa2tPaGxuUYpTGNxZ+fOnTm/Z76CnI+fNhqNmjehdevWndDvns2ajkajpjU9PT1NKBQyF8WShfpUs6YBHux/kIcGH6LCUsFcdI7b/3Q7/7DjH2h2NZtjZEFmOjxNRImAgJmJp+s6mq4hIhLQA8SUGC6rC0/UQ0AJzCvLaXS5DikhBvwDrKtex0RogpuevQlfzIema2yv385N229KWNSkLiJGo1FeffXVZS1neqpFWZwszghBTia57U++LGTNhUIhent78fl8OBwOtm3btjz1CbLMI7kUZlVVFU6nk7Vr1y64vXzTrHN9l3g8Tn9/PzMzM3R0dODxeE6oGGe7WSRnb6VXcjN808nWdLKQlJWVmSUgT1ayy9MTT1NprcQm2bBjZzo8zSH3oRRBtkgWREE052i0ddLQTKsXIKbFGA2OmscpU00LVVcJxoPcdfgubjv3Nu54+Q5mI7O4LIn07v1T+3lk5BEuXHUhkJqJaCQ9LFTOdCmZiMV0Wfj9/tOyOD2coYJcaE3kbCdDKBSir6+PQCBAZ2cnbW1t9PX1LctFbPiGk+dhlMLs6emhrKzMLIX59NNP57XNpURPGK2axsbGaG1tZffu3YiiSH9/f97uiOWq45ALSZKoqKhIqWVgCIlhTU9OThIOh5EkiUgkwsjIiCkkxUwRzvXUYZNs+GN+M5NOR58XgiYJEtXWaqqsVYTVMHOROcLa8fZLZq89ZKyilbgWTwlrS5kLOp6Yh35ff6KDdWDE3J8gCKi6ynAg9yL0Ynoi5puJWMzi9JqmnRLdYhbD6TnrNJZSE9kYn/wHNFKMDSHeuHEjgiAQDofzDk8z5pWvPztdkI0uIXa7fdENTAvJ6jNIznBcuXJlSvNSyM+9YViexfDPFmM7yUJSX3+8zq+iKOzfvx9BEJicnKS3tzflsdxYRLTb7Yu6seS6IV3RcQX/8uq/EFJC6Og0OBrYXp/aIdkm2dizeg+/HfotqqaioVFnq8Mf9xPVjhfRUlEJq+G8qr0NB4f52yf+lnpHPe45N7Igo6MjCRKdFYtLaihGJmKxLGT9WM3y05UzQpDh+IW72K4hkBDivr4+M47REOJMY/PBENlCBNnr9fL6668jy7LZJWSxFFL3Qtd1xsfH6e/vp76+PmtijWF1n47dGNKRZRlJkmhuPu4mSLam/X4/4+PjGTPlysrKFjwGuQR5a/1Wbth6A6+6X8UhOzi38VyzQ0cy13Zfi67r3H3kblQ9Icobajbw8uzL5oJdoWU3X5x5kWprNTWOGvwxP5qu8Y6Wd/DnK/+8oO0sRCGZiMa1V1FRsaSeiIYYn6pFqRbijBFkg8VYyMFgkKGhIXw+H52dnWzYsCHjH3SxgpzPopKmabzyyitIksSaNWsWLB9YaCeQXNuZmZkhFAoxNzeXUnwpEye6BOeJTgzJZU0b1t74+DiBQABN0+ZZe8nW9EJ/o3w6iHiiHh4dfZQ6ex0zkRmiapQB3wCiIFImlhFQCy+EpaMT1+I4JAffevO3sIk2auw1mccW+dhny0Tct28fjY2NhEKhJfVEDIVCRe03eaI5YwQ52ULOVvYynUgkwtz/Y+/L4+Oqy/Wf2ZNM9q3Z97V7mqS0LFoovYBAQUHEi1IBuVyvFqtcEH4goqItUrislR1RkUJBhCtQBIEC0jRtWmibZplp9sk6S2bf5/z+yP1++c5kljMzZ9I2zfP59CPC9MyZ7Tnved/nfZ7paWi1WtTX14ckYvY5ogEfAidLJWazGTU1NSgpKYl43FiimYKBXepITk6mIZbxHHO+IlRmH6n22GpaKpXSwaHX641LQaCxauDlvMhNyoVYJIbWoYXD60C2IhtSkRRuuxsuzhW0Sk4WJ8PNueHhZn8HXV4XrB4rksXJ+Iv6LxixjKA2oxZX115NvS6AuTEWAma+01lZWX5ytViGiNPT06fsQA+YR4RMwKdCZt3PUlNTkZub63e1Fgrh1pyJcsNms6GmpgYKhSJoYkm44/KNZgqE2WyGSqUCALrUwXdQON8q5HjAxhOxq7Vut5vekrtcLhw6dIiGfUYbBZUhz5jxoIAPOUk5SJWmwulzYvOyzbjv0H3wcT4oJAo4vc5ZpOzjfMiQZ0Dn1M06rptzY3nOctx78F6MWceQJEnCgHkAI9YR3NV8Fz2vuSLkYIh2iLhnzx50dXXBYrHg0KFDaGxsjEnbzMcLmcBkMmHx4sW4/PLL8dhjj8X8Wgnmh8MLg3C5ek6nE11dXTh48CCys7OxZs0aZGZmJmSjDghOyA6HA52dnfjiiy9QUFCA1atXIycnB1KpNCpbTT694cChns1mw+HDh9HV1YWqqiqsWrWKDmL4Et/pWiFHA5lMhqysLBQWFiItLQ0tLS1oaWlBZWUlTbTo6enB/v37cfDgQfT29mJ0dDSoL3ZZWhkuKr8IZpcZRpcRDp8D/7HkP9Cc34yHzn4Iy5KXYVHyIijECkgggfT/aiwRZgar065pP4kcQW5SLi4svRATtgnqc5Elz0LPdA/0zi9TcuZibTpakCFiYWEhamtr0dTUhNWrV+O73/0uVq5cCZ/Ph4cffhhf/epX0d7eHvXxiReySqXC+vXrsW3btpCP/fnPf46vfEW43vu8qZDJFT1Ym8DpdKK/vx96vR6VlZVoaGjw81C22+2zjhcOfOVc7Lm4XC709fVBr9cH7VPHYzwf7nFer5cudRiNRtTW1gbVERNFSKQf38lcsZ5sILFEQOjeKammA5UIbDX9jbJvoDWvFQaXAcXKYhQpiwAARcoiXJtzLSxFFjzZ+SR8Ph8mHZMQcTODPqfPCfH/1VyElMUQI0mahEfOeQQezuPno0z+WSr6khbmgpDZ9ykeFBQUoLa2FmazOSyJRgIfL2QA6OjowMTEBC688EIcOHAg5udjMW8ImYCtkFkirqioQH19/awPPtZBHZ8JMEmT7u3thVarDXkO7HGjOYdI8Pl8sFgs6OjoQGVlZUg/ZIC/eiISIft8PupHTXp9rPVmtAsDp/IFgM+Fm1TT7C2xz+ejt+QGgwHDw8PUxtSaasVY6hgdcgEzySE6hw5eeOHmvrw7FEE00+pQ5MDutaM8tRyLUhbhR8t+hEJlIbycF/WZ9egydEEiksDLeXF24dnIUGT4ncup5PQ2V17IPp8Pt9xyC/785z+HdY2LFvOGkNmK1+VyoaenBzqdLiwJAtETMnl8JEL2eDwwGAx0WEcWK0JBSE9kstSh0WggEokiPjefY0Z6HCubKygoQGtrKx3KmM1mGAwGDA0N0S1K1ihIqPXbkw2xLsaIxWL63rAI1PVarVbYbDZ8euxTiCGG3et/p0d6ymKxGHnyPDx97tN+5yMRSXBb0234x/A/MGIZQXVGNdaXrPc7xlyZ08+1sVC8Xsg7duzA1772NV5D+GgwbwgZmPnC9vf3w2KxoLy8HLW1tRGJKFZCDgV2wy0tLY36EkcCqab5IBR5By51nHHGGThw4ACvCiceQtZqtVCpVMjMzKT+Gm63Gz6fj96mE7CObsQsyGazhdT5xmMNygeJTiMR8kITTNfb3t6OnOwceMY9MyGnnP93iNhvXll1ZdBzUUgUuLTi0pDPOVfWm0JWyKWlpREfF68X8t69e/HJJ59gx44dsFgscLlcSE1NjatVAswjQjYYDDh8+DAqKiqgVCr9xP7hIBQhs1aYZMNtamoqZAp2IOKpkElSR7ClDr6Ew9edjm0hGI1G9Pb2Qi6XUy9kPn8/0NENmPlRkvRgVudLHk/0vokO2xQSc7E6LhKJcEbJGUDnTIsiSZwEp3eGlJUSJbJkWbgs+zKU6EqwX79/VgspkkZ+rloWJ5P1Jh8v5BdffJH+8x/+8AccOHAgbjIG5hEhZ2RkYO3atRCLxRgcHOT992L1viAgqdbBrDATlZVHHssmdWRkZERc6ggHvhWySCSCw+HAF198Abfbjfr6+ohLLHwglUqRkZHhpyEl/WibzebXSyUm9qTtEU/LI5GEmWhCJhfGnKQcfLXwq2ibaINELEG2Ihs+zodr66/FxRUX08y/wFCAgYGBoDamycnJlIRPtQrZZDLNiRdyojBvCJmVeBF1AZ8POVojefJ40jMdGBhAbm5uUCvMaKrvaAnZZDJhcHDQL6kjHvBpDbhcLuj1ekxNTaGxsdFPG5oIiMViKBQKSCQSv9tQYrvJtjzYvmtaWlrMq7dCYi4ImRx/y4oteOTII1AZVZCKpNhUvwlfLf6q3+PDhQKYzWbqkMcuX3i9XhoAkSgbUyFJ32Qyxb0YwscLmcX3vvc9fO9734vrOQnmDSEH85zg8yFH+4ORSCTQ6XTo6+uj+XGhqtJEVMhms5kO65YvXx51UkcohKuQvV4vBgYGMD4+juTkZBQXFyecjMMhWMvD6/XSgRe7esvKx4jKg11tTiQSTchsOyFdno6vFH4FYoiRKk9FeXo5r2OwLaRgNqZDQ0OwWCw4cuSIXzVN7k6IjWk8EDq+KV5z+hOJeUPILIj0Tch+I2kPDA0NQaFQoKmpKeIWkJCEbLPZoFar4XA4kJ+fD4VCIRgZA8EJmR0SFhcXY+3atdR+c67AV/YmkUhmtTw4jqPyMaPRiJGRESofI4Ti9XoT1iedS0J+f+R97FTvRIo0BW6rG73Tvbir5S4UphTGdGxiY0paQvn5+X6hAGazmdqYsncnsdiYnmw95BOJeUPI7Bc/2r4wEP7Ho9froVKpkJycjIqKCmp6EglCtCzYpY6amhrk5uZibGyMtyKDL1hC5jgOk5OTOH78OHJzc/2GhHwJ8mSQsbGeB+wyBrlFN5vNcLlc6OjogEgkiotUgmEuCfkDzQdIl6XTeKZJ2yQOaw+jsCw2QiZg2wnhQgFYr+lAG9PAUIBACBnfFE9a0MmAeUPIwJdkEa3jW6iliOnpaahUKshkMixZsgSpqamYmpqCwWDgddx4KmS3242BgQFMTU3NWuqI5rgAf2c4juNgMBjQ29sLpVKJVatWzfqhzPXqdCIWQ+RyOXJycpCVlQWtVouWlhZ6i04qP0IqycnJtJoObHlEwlwSslQkpYb1BGJR/FU/n7uHUHcnoUIBWJImfWohesin6gIRi3lFyASxeiKTL4XJZIJarQaAWSqCRA7qiDNYsKQOFtGQIt8NPI/HA7VaDYVCQS8+wXAqb86FQ6ikEda7V6PRwOl0+rU8yC19MNISaiU4FNjjX1JxCZ4+9jScXic8Pg/S5emzDO9jQaxkGc7GlFz4SK/farXCZDJRs694QgHIc5+qmJeEHIsnssfjgdvthlqthsfjQU1NTdBeVDSEHO0Xw+l0oq2tLWhSB4tYNMuhjuVwOKBWq6HT6VBUVITa2lpex5srnMgLQCg3N3axRafTwWaz0fYIS9QcxyVUw8se/4xFZ0ApVaJjqgPJ0mScV3xeSI/jaCC07C2YvPHo0aMoLCyE1+sNGQpAlDPhzsXhcPB2TDxZMa8ImfVE5ruQQdDd3Q2v14uamhq/TahARNsuiASSJN3X1wefzxcyqYNFLBVyIDweD/r6+miCtlKp5NV7i2Zz7kRk6kWDWM8v2MYcG65KWh5OpxMSiQQcx1Gi5mO5yReB7YSlOUuxNGepIMcm4Jt4E+9zkPeGvfCxoQAajQZWqzWsjen09LQgmvgTiXlFyAR8q1iHw4Hjx4/DYDCgrKwMVVVVEX8s0W72hQIbYJqeno7m5mZ0dHTw0npG2wphCZQsW2g0Gr+WyNDQEO/V6WjuPuLFqdIiCdbyGBsbg9VqhVKpnNXyYCu/UC2PSJiLLbq5iOsKVYWHCgVglTPkPd2zZw+OHDkCq9WKjo4OLFmyJKFeyENDQ/j+979PTbTefvttVFRURP18gZhXhMwaDIUjDafTib6+PkxPT6OqqgoKhQJKpTJqS02+CKzEyLBQoVDEFGAabcuCbPWx5j/sRiH7uEg4VQjyZAGp+oK1PCwWCwYHB/1aHixRR1J5nGpObEI8RyjlTFVVFXbt2oU333wTO3bsQGdnJx577DG0tLREdS7EC/n222/Htm3bsG3btqDWm9deey3uvPNObNiwARaLRbDPYV4RMkEo0nS5XBgYGIBWq/XzRXY4HLxJlq/nA/t4UmWwSR0NDQ1BdcRCZeWxj9Xr9Th69CgyMjKCbhSSx/F5D/g8N1FqeDwe2lclt+zRSsnmYvU4kccPdv6RWh5TU1Po6+ujKg+WpNmWx1xVyHMV4RQPcnNzUVVVhTPPPBO/+93vYj4OHy/kY8eOwePxYMOGDQAQVxBxIOYlIQdWyB6PBwMDA5iYmEB5efks5UK0Q8BoIJFIaCXkcDhQW1sbcpOIb1Ye3wrZZDJBp9PB4XBErMSj8bIIRWRWqxW9vb3gOA6NjY2QSCRUpcBKyUgPkBB1pN71XIacColoetSRVB6BmX1paWn080okaZ7scwAWc+WF3Nvbi8zMTHzjG99Af38/zj//fGzbtk2QO4l5RciBqSGsL3BZWRk1HwqEVCqNOjWED5xOJ2w2G44ePYq6ujrk5uaG/XIT46J4vYvJVp/T6URmZibKysoitkXisd90uVx+iSTZ2dlwu91wuVz0NryoaCblgvQAAz2S2ZXctLS0uGRPJwviVVmEy+wzm80YHR2F3W5HR0cHAMxSeQjhPXEqfQZz5YXs8XjwySef4NChQygrK8O3vvUt/OEPf8ANN9wQ+8n/H+YVIbNgJWRr164Ne/WKpUIOVzmwSx1JSUlobGzkNf0lfdxYVRYkJspgMNCtvt7eXt6Vb7QVMquZrqyspEEApGdNFAY+n49W9BzH0ZTr/Px82gJyOp3U4GZ8fJxWglKpFBzHwWKxxDz8CoVT1ctCJpMhOzsbdrsdmZmZKCkpocGfxHCpv7+fbsCxJB3thS7R75GQWm2TyYSCgoKIj4vXC7mkpAQrV65EVVUVAODyyy9HW1vbAiEHguM4jIyMYHBwEBzHzRpchYJMJotKykaUC4EkH2ypo7u7W/Dw0sAvsNfrxeDgIMbGxmYlpASqLEIh2py+sbEx9PX1obCwkLaACPGKRCJIJBL63pPjEqImJM3+O6lUipycHOTk5FDSdbvd0Gg00Ov1fsMvtqdKjOxjxcnSsoj1+OS9IsGfgU5uZFsusOXBvn9KpfKEJUufbPFNfLyQW1tbMT09jampKeTl5eGDDz6IengYCvOKkEmkzerVq7F//37eH3S0FTKx4CTHZz2RCUGR/xaNvWcsK9EajQaDg4Mhl0n4qif4EjJZLRaJRGhtbYVMJoPP56OvMdhtHvmxsz96tv8ZWEWT/yUTdbfbTRdWWN8EYmRPQkHZFedEWUVGg0QPxHw+X9iCI9S2XGCwKtHsB7Y8pFJpwlsWQju9zYUXskQiwfbt27F+/XpwHIfm5mbceOONQryE+UXI6enpqKurAxC6ig2GWFND5HI5XeoINOEhSIQFJ8dxcLvd2Lt3L3JycsIuk0RT+Ya7PSUDO7fbjaysLDQ2NkYk4kjPx/4vAUvQdrsdGo0GmZmZcLvdtCIkpEtSYUgoKNmcI8br7O16oEIBOHVbFgSxEn6oYFVW5dHf3w+3200louR9FLq3f7JVyHy9kDds2IDDhw/H9VzBMK8ImQWpehNByBKJBFqtFmNjY3SpY648kaenp6mkbPXq1RFXReMNL3W5XFCr1TCZTKirq4NUKkV/fz9tNURLxHzOg6hi9Ho9ampqkJWVRUk6sN1BkJycjJSUFDohZwNWyQAs0IdCKIexUDhZCTkYWPN68h46nU4cPXoUaWlp9I6EbXmwK82xnofQ1punshcyMI8JmUjf+PzootEWT09P85aSAcKZEZEK1efzobGxEUeOHOH92vj2pVlCZvvSVVVVfpptu92OtrY2eoubnp7OS74WCRzHYXR0FENDQygtLUVra6vfD529uLLkzPal2eGhXC5Hbm4ucnNz6XFY683JyUmYTCZ0dHT49VRTU1MFIbpTiZBDHV8ulyMvL493yyOavD7g5ItvOtGYV4QcLDVEKFgsFqhUKvh8PuTm5iI/P5/Xhl28adJOpxNqtRpmsxl1dXV0mSAazTKf/jipkNmNPqJQIWRNFCBkgGGz2ajWeWBgAC6Xi9pVEpLm692g1+uhVquRlZWFlpYWXkoT8voI2L50sGqatDzI7brH48GxY8ewdOlSOvjSaDSwWCwA/HuqfDbnAjEXhJzI44eSYEZqeQS2jViSDmx5CNmycLlcCb/rSTTmFSED/gZDQix72O12qNVq2O12utRBJEV8IJVKYbPZeD2WJWSPx4P+/n5MTU2hqqoKixcv9vsik8fy0SzzHeo5HA7s27ePbvSRgR07ZGPPgaywBrYJiJXi8PAwnE4nTTchJM0alVutVqhUKohEIixdujSubMBwfWm2kiavx2q10jsjcn7k7xKCITIydnOO3TwMl0ozHypkvmQZKq+PlTOynsjk/SMtkHgxX9b55x0hE0RbIZMqkHzBWb8LoukNXDzhg2hVFi6XC4ODgxgZGUFpaWlQP2Qg/t4wC4vFgu7ubtjtdpxxxhlISUmJaWDHTvVZnwGn0wmTyUQ9cG02G12C8Xg8qKqqQmFhYcLISywW+72HxOmO6LXZdg1bSaekpPj1R8nmHFlqISnYoXLmTnVCjrd6ZRNGAj2RyR2JXq+Hy+XCxMSEn1ImNTU1phbYqbTIEgzzjpBjrZAJyYrFYlqZsn4XLOJtQwQDx3EwmUwYHx9HaWlpRA013+OGI2QSD2U2m1FdXY3jx48jKSlJ8IGdQqGgfUiS0zc8PIzc3FzIZDLodDoMDQ3RKotU0kL1cgk4jsP4+DgGBgZQUlKC2tpav9cX2JcONjxUKBRISkqiFxw2tZnNmWOHyhKJJCFa30T7LSfKepN1cfN4PEhNTUVOTo6fUmZwcHDWBme4GCin03lKRzcRzDtCJpDJZFHlzkkkErpdF64yBWa+UHz9lvkQp06no1FRixYtQk1NjSDHBYITMpsiXV1dTSVsdrsdIyMj9NZT6GpDq9VCrVZTiWDgBYdUTiaTCcPDw7SXm5qaSkk6LS0tpqrNaDSit7cXaWlpIXvUofrSkYaHUqmUmgWxSy3Hjh2jy0KBJvZCLLWcTC2LWEGqcDYolSBUqCoxrmcXW4xGo5/p/amKeUvIfJc9yFKHXq9HUlJS2KQO9tjR9JBDPdZsNqO3txcSiQTLli2D0+kMamYSDHxbFuymHlExDAwM0BRpcrvu8/mwbNky2vu1Wq10K45UrbESCBmISqVSrFixIqRUL5j/LVkEIZtmRGkSqPAINQQkQ1GHw4HGxsaonbliGR4CM3dqMpkMBQUFdG2eze0jSy3EcJ0lab6V3snesuCDcIshoUJVA43rX3vtNbz22msAgPvvvx8rV67E2rVro/6s+Xoh33bbbXjrrbfg8/mwYcMGPPzww4IVL/OOkFlP5HCkySZ15ObmoqCgAAUFBbx1y/Foi+12O1QqFZxOJ2praykBEUOkWI8bDGSop9Pp0Nvbi6ysLFqdBg7sQpGhyWSCRqOB2WwGgFkkHeoHRbw1zGaz3+uMBsECNMkiiMlkogM3t9uNlJQUStBKpRITExMYGxtDdXU18vLyBPvRRFpq8Xq90Gq1MBqNKC8vp4UB63kcuNQSSp0QbiFjLgj5ZPRbDrxw33XXXVi/fj2eeOIJFBUV4R//+AdKSkrQ2NgY1XH5eCF/9tln+Ne//kWXQs4++2zs2bMH69ati+q5Qr42QY5yEiJchUxundPT02myskqlSsigjiVvt9tNE0oCB4XkuNFm5UWC3W6HTqcDx3G0OuU7sAtFhoSkx8bGYDabacXK9n7HxsYwOjo6y1tDCIS6vSU9SI1GA51OR6f5JpMJABLuIicWi2G1WtHT0wO5XE6VKqSCJn8A/+EhWWohxjiBHhTsUgvriJdowiQ65ERCqNVpq9WKyspKXHPNNbjmmmtiOgYfL2SixXe5XHRjlh1gx4t5R8jhUkNIH1Eul89a6hCqDRHsfLxeL/r7+yMSlJBbfeRW3WQyISUlBStWrKAVcTwDO7FYPMu3l0jETCYTBgcHYTAYIJFIkJmZCYfDAb1eL8jiSDiQ1zI2NgapVIozzzwTCoXCT+Gh0WjgcDggl8v92h0pKSlxkzT5jPV6Perq6mbdDYTrSwcbHsrlcuTk5PipE9hwVa1WC4vFgkOHDvlppSMFgUb7mk6mtJBwmCsv5LVr1+Lcc89FYWEhOI7Dj370o6gr8XCYd4RMwJImu9QRKqkjGlUGX0ImPVsyAIzUnxaiQiYDu4mJCVRVVaGmpgaHDx/2I2KhqypyvPHxcSQlJeGss86CXC73Wxwh3ghsvzQ9PT2sjpcviIxtenp6VgAA6UEGRigRkp6cnITNZqPVNNvy4PM+cRyHyclJ9PX1oaSkBK2trbwSX4Doh4fkIpeZmQmxWIz29nasXLmS9qVJECjHcbQ1Qsg6FrOlk7VlEQxz5YWsVqvR1dWFkZERADOeFp988gnOOeecGM/cH/OOkMmbSIZVR44cgc1mQ11dXdg9d6lUyluVEaldwHEcbYtkZ2dDqVSisrIy4nHjqZBZ57fi4mKsWbMGIpEIHo8HHo8H7e3tVK0Qz4AuEEQ6RxZn2Mo52OII0fFOT0/TxZGkpCRKhISk+VSsZKtwcHAQpaWls2RsocCuVBMQ03ez2YzBwUGak8YqPALfM9KeUCgUaG5ujusOINrhKv5rZQAAIABJREFUodvtpp8/Id7ApRbWKMjj8dA4KDapJdz7NRcqC0AY7bDRaKT+xOEQrxfy66+/jjVr1tB22UUXXYS9e/cuEHI4kAQLm82G+vr6iEkdQPRtiFAgbRGFQoGVK1ciOTkZOp2O13Gjzcojj9VqtVCpVMjOzg46sGttbfXr/Y6MjMQtKSNyromJCVRWViI/Pz/ie8wmYLA6XmICRIaHDocj7HYf8OX7nJ6ezmvVOhKI6Xtgzh0habJSTaw+XS4XnE4n6uvrkZOTE9dzh0Ko4SG52BcXF89aagG+fJ+DXQxJWvPIyAhcLhdt3xCiZt/nuWhZCIW58kIuKyvD008/jTvuuAMcx2HPnj3YsmVLXM/LYt4RstVqxf79+1FZWQm9Xs+LjIH4vS/ICrDX6w3aFuGztRVNpSCRSKgxDisnCzWwCzWgY93QzGYzOI7zq6QDSZq9RS8sLMTq1avjjiki231sRcL2fsfGxmC32yGXy5GSkkIvJqHaT0KBbREAX1bkfX19yMjIQEpKCvr6+qBSqWgrRiijpWBwOp3o7e2F1+udJR9kK+jA4SEw8z6TpZbA95msNk9OTvqtNlssFtjt9hNqYM8XQhgL8fFCvvLKK/HBBx9g2bJlEIlEuPDCC3HppZcK8RIAzENCViqVNDtvZGSEt71frIRMbtlNJhNqa2uDVkukvSCUzaDD4YBGo4HNZsOKFSuQkZER08BOLBZHVFEQ3S/pQ+r1emo5msghHbvdR86rr68P4+Pj9PV2dnbScFBChkLHPBFYLBb09PQgOTkZq1ev9nvtROEhlNFSIEgSzsjICGpqavwGfQR8+tKEpD0ej58NQOBSi8fjoevhExMTGBwc9FtqYQ3s44GQ5khCLIbw8UKWSCR48skn43qecJh3hMwOrYj0LVGErFKp6Ip1Y2NjyC8XOXa8X2DiEzw5OYm8vDykp6cjIyND0IFdMBWFzWZDT08PTCYTUlNTYbVacfDgQT+pW3p6umAXHBakH3/8+HEsWrQIZ555pt9rJL1fk8mE/v5+WK1WvxXs9PT0uCo8dmBYX18f9EdPyCpeo6VgMJlM6O7upvrxaFoI4frS4eK0yIVaJpOhoaGBRpyR4WFggjjb8ojmIi20Of2p7oUMzFNCJoi0HMKCLyH7fD6MjIzAarWioKAg7Io1QbTRTIFgB3YlJSVYs2YNTWy2Wq208hJaX0sUG1NTU6iurvZr/7BuaJOTk34/ULbdEU9vl3hAy2QyrFy5Mqi1YrDeL6nw2AGdSCTyU3dEGmqyvhfRDAwJojFaIqb5rMLD4/Hg+PHjsFgsWLx4cdRbZ6EQKU6LVNNGoxEOhwM+nw9ut9tvqYU1WwqVIM6SdCjtt9DxTQuEfJKCGAxFk5UXiTTZzb68vDxkZ2ejoKCAV+UVzWYfeS7yBZ6amqJqDXZgp1QqkZycjK6uLj+lAvkTj5yMJaPi4uKgfWLWbrGoqIj+PaJHJiTt8Xj81pzT09MjkrTb7UZfXx+MRmNQTW8kSKXSWX69fLYOSb+ctCdSUlIEGRiyCGzFAP6m+VNTUzAajXC5XMjIyEBBQQElykT1cVmSJhcCs9mMFStWQKFQhF0PJwni7FIL6UuTiw6x2GRJOiUlRdAK2eFwREzPORUwLwmZIJoKOVz1o9froVKpkJqaSuOajhw5EpVELdrUELvdTje+SHXIDuykUikN/gx1exwLSU9PT0OlUsWkXmAToQNJmiwzsHpktpKWy+V+iSFlZWWoq6sTrOrns3XY09MDu90OACgsLAzaq00EyBJIUlIS9Ho9cnJyUFlZSatpoY2WQoFc/EtLS/3e+1B96VDDQ5lMNmupJVBWyHpRazQa+r2J5fWQ45zq1pvAPCfkaNOkA0HMf8RiMZYuXQqlUkn/WyLTpDs7O6mkKj09PeLALtjtMUvSwTS/gSRN/DW8Xi8WL17s91rjAUvSbH81cGnE4XDA7XYjNTUVlZWVyM7OTvgPjPTL09LSMD4+Tr2vMzIygrZiEqWiIK0hrVaL+vp6ekeQnJwsqNFSKDidTvT09AAAVq1aFfbCHc3wMNxSCwBMTk5Cq9X6LVDFkyC+QMgnKeJJDSG3XCqVCna7PeQts1BZeQQkIcRsNqOoqAglJSUAEPPAji9Jkz6hz+dDSUkJiouLBdmei3RuZAiWlZUFlUoFsViMsrIyeDwe2lYIVCoItdnHwmw2o6enB6mpqX53BIGtmHBbh6yKIloQ69XCwsJZGYKB4GO0RNpEfM6NzCaGh4dDqjf4INbhIbnQFRUV+S21kL60VqudlSDOJrUQAhayF32iMT9eRQjIZDJ6C8oHIpEI3d3d1PwnnEOYUN4XxLB9cHAQZWVlyMvLo7Iu8gUWamAXqPkdHR3F4OAgioqKkJKSArPZjKNHj84JEbKLJcSNjYDtR7KbfUNDQzQ3jW13xCInI0ZPZrM5op45lIqCTQ8ZHBz0O7dIW4cOhwO9vb3gOC7kwJIPIhktBZ4bOS+JRIK+vj6kp6ejtbVVcEKLNDwkBUFlZSW8Xq9fyyNYgji71KLRaKjZEsl/JCk3sfbZd+3ahXvuuQddXV1ob2+nuZGB2L17N3784x/D6/Xi+9//Pm6//faYni8URFFmUZ0SwVXExlKv12NiYiKi+YfP58Pg4CDUajVqampQUVER8Qc+NDQEkUiE0tLSiOdD9NDl5eX03xE5l0qlQm5uLqqqqiCRSHDs2DFwHIfc3FxkZGQkROvLBopWVFTMuiVkK2nyh5B0INlEC47jqGVmQUEBysrKovoRsedG5G6kFRN4AQn2GbLr1uXl5YJGRwVuHZrN5llbh6mpqdBqtRgbG6OOf3MBcudnNBoxNDQEi8UCuVzu974JZbQUDh6PB2q1GlarFY2NjZRIg7U8yHmLxWJalLDfFafTiS+++AIvv/wy3nvvPWRlZUGpVOLRRx/FypUrozqvrq4uiMVi3HTTTdi+fXtQQvZ6vairq8N7771HfUteeuklLF68mM9T8HpT532FHMkTmRi2kyEOX99cvmnOwEyF7HK56P83mUzUA6GpqclvYEc2DPV6PY2xIXrfeKVkfANF+bQ7SLUaDUlbLBa6Vt7U1BQToYc6NzIAI6vhRPPLnpvb7UZPT0/Y1JB4EGnrcGpqCseOHYNYLIZSqYRer4fH45kTIhSJRLDb7ejv70dBQQGam5shFosFNVqKBNKeKS0t9XM8jKUvDcz8rlpaWuis6Pnnn4fJZIrpc+Xj2Nbe3o6amhrqmXH11VfjjTfe4EvIvDAvCZndQgpGmqz5T1ZWFlpbWyGXy3Hs2LGoSJZvO4QMAMltKhnYpaWlzRrYKRQKFBYWBh2AkR6h1+udRdLhbjlZGVmgGxpfhCLCYC0FQtLkj0gkQl9fH0wmE+rq6gSP2mGTJQgRsiQ9PT2N3t5eOjSUSCR04zCR/sgEYrEYWq0Wdrsdq1evhlKppFI3k8mEiYkJGvyaiK1Dt9sNlUoFh8Mxa+U6WqOlWNJj3G43ff/5tmeiccR7++23odFoAMBvoUloaDQavzvikpIS7Nu3T9DnmJeETBDOE1mhUGDFihVz4okMzBjCaLVa1NbW0i8/n4FdsP4lO8iZmJig1qKBHhQikQgajQYjIyMoLy8XVEZGzi2YWRAhab1ej56eHjgcDiiVSuTl5cHtdlNTm0SCXNzcbjd0Oh1qampQUFBAK0IidSMeGSwRRtqe4wu2PVJRUeEXmEukbuyqvdBbh0Q739/fj8rKSixatIjX6+JrtMR+58jALbAwINr9aJ4/FAJJenJyErfccgvEYjEefvjhiH8/nPVmMCOhE4F5TcisusFms6G3txcej4fKyQIhNCGTrb6BgQHI5XKsXr2aunPFM7BjBzlECcCaxI+NjaGzsxMOh4MOR5RK5ZzYKRKSdjqdmJ6eRn5+PioqKuB2u2EymfyGTIGVtJAkTdpCgUOrYIsZZJGBpH7b7XbIZLJZ1Wo0n5PVakV3d3dUyyXRbh2GsgUFZmSM3d3dkMvlgrRnAo2WAP/v3MTEBNRqNVVOpKSkwGAwCGJNGgiO4/Daa6/h/vvvxy9/+Ut8/etf5/XZhLPe5IPi4mIMDw/T/z8yMkKjuITCvCTkQE/krq4uersezipRKCkbGVyp1Wrk5eVhxYoV6O/vp1Ifcm5CxxqRqnhychLp6eloamryk5Gx22lsJS3kBhhpy/h8Pr8+tUwmmxVTZLfbYTKZaL/c5XL5ybViIWm32w21Wg2bzcY71FShUEChUPjdtgfrrbIrzqFImk0OCeV9EQ0ibR0GWqmmpqbCZrNR7w2W3IUGu61JQAbkw8PDyMjIgNvtRkdHh2BGSxMTE/jpT38KpVKJDz/8cM6GosCMja1KpUJ/fz+Ki4uxc+dO/OUvfxH0OealysLnm4m0HxwcxPHjx7Fs2TIUFBRE/AKQ7DI+ZvIOhwOdnZ1obm72+/ekJZKUlITa2looFAo4HA7s27ePJhBnZGTELHMKBTZQNFyflv0xm0wmWnEFmtdHS9JerxeDg4OYnJyMWT3AkjQhQ0LSbLUajKRZTW1FRQWvzztasH1fs9lMB2Dk3LxeL62aSktL53RRwefzUdMfqVRK/bJZQ/p4vUUiwW63o6urC8nJyaitraV3JaGUMdEYLfl8Prz66qt44IEH8Otf/xqXXXaZoO/v66+/js2bN2NqagqZmZlYuXIl3n33XYyOjuL73/8+3n77bQDA22+/jS1btsDr9eL666/HnXfeyfcpeJ3svCRkm82GtrY2FBcXQ6PR4Mwzz+T14U1OTtJKOhI8Hg86OjpwxhlnAPhy083lcqGurg5paWl+Awi2d8lKtTIyMuK6Zff5fBgeHsbo6GjMfTrSH2RJmm/vkvVHLioqQmlpqaAVdyBJm0ymWQnTANDf34+MjAxUVVXN6ZIA6VH39fXRBYVgZkGJ9BNm/ScaGhroXQFrAEXIkDWAYpND4gHHcfQ7WF9fz3tozBotkQsc+97JZDJkZmZCq9XiJz/5CdLT0/HQQw8lLBAgwTh9CZlUyFKpFPv378fKlSt5VQZ8dcvAzJdw7969aG1tRV9fH/R6PV0mIa2JcN7EgVpfo9HoRzQZGRlhKxpWz5ufn4/y8nJB+8Okd0nOjwyY2J6vz+eDSqVCcnIyampqEj6oIyDKE4PBQLcN5XL5rFZMos/H5/NhaGgI4+Pjfu0wdjhnNpv93jtCNkKRNOs/QRJEwoFV7ZDzi2fr0Gq1oquri14M4/0OklmD2WzGu+++i0cffRTT09NYs2YNvv71r+OSSy6hw+1TDKcvIZOKFAAOHTqE+vr6kJpbFmazGf39/Vi+fHnEx/p8Pnz88ceQyWQoLy+nzf14BnaBPxaTyTRL4paenk4HlElJSaipqUn4qjMB6UdPT09jbGwMTqcTKSkpyMrKCttXFRqsYTu5KwBAt9OCVdJ8neb4Ynp6Gj09PcjNzUVFRUVEImIvcCShg+3BRtsqYv0n6uvr43b3I3ch5BwjbR2SXvHU1BQaGhoEl5uNjY1hy5YtyM7Oxq9+9SsMDw/j0KFDWL9+vaC63znE6UvIAGhg6dGjR1FWVsbrC0N6YKtWrQr5GHKLfvz4cTidTpxzzjmQSCQJWXUG/CfZBoMBWq0WPp8PWVlZyM3NjVoTGg9YIiR9WkLS5A+57WRJUCgZGTDTo+/p6UFWVhYqKyvDtifYCxwhmnhJ2uVyQaVSwel0oqGhgdeFPhQ8Hs+sfj4brBqMpIXyn4iEcFuHCoUC09PTyM3NRU1NjaDfPZ/Ph5deegmPPvoofvvb3+Liiy+eF6ZBON0J2eVygeM4dHd3Iy8vj1ffye1249ChQ1i9enXQ/07IICUlBTU1NTh48CDOOOOMhBAxi8BA0dzcXErSRqORTtlJJZORkSF435JYkObk5KCioiIsEbLDL0LSROtL/kS7kMESYX19fcxudOFu2cORNGsNKoSmNhTYfj6ppAHQFoxWq0VmZiZqamrm3FDH6/Wit7cXBoMBWVlZcDqdVMctxPr16OgofvzjH2PRokV44IEH5oXhPIMFQuY4DsePH4dSqaRyq3AgfeEzzzzT79/bbDaoVCp4PB7U1dUhNTUVHMfhwIEDdD2X9HyFHmixgaLhfB9CDeZYklEqlVH/UMiw0ufzoa6uLuaKkB3gmEwm2O12+t6F859gq/Kqqipe6dbRIhJJy2QyjI6OIiMjA9XV1XNOhORipNVp4UxywugwIlWSitKMUvq9E9obORAGgwE9PT10cMt+BqxEkHy20WQd+nw+vPjii3j88cexbds2XHTRRfOlKmZxehOy2+2mfS6JRELtLCPhs88+o4RMVo71ej3dsGMHdmxCBqlUWdE+ueWM5ctlNBqhUqmgVCpRXV0d04AqsJ1gtVp5txNYj96ampqETLbZW2KTyQSHw+Hn5MZxHAYGBpCdnY3Kyso5jaTnOA5ms5mqF4iBfrTpJ/FCr9ejt7cXBQUFmEyaxIB5AHKJHE63ExXJFSgWF1My5LjwieGxgJgBEV0331SOYINN4pFBziszMxM6nQ4333wziouLsX379riTo09iLBCyz+eLSlsMzBDymjVrMDQ0BI1GE/XAjlSqRqORkiCpFkglHa6n6nA4oFarqXxOqCw1AjLFJhcRdnWYVFvT09M0vqmkpGTOIuCJ/wRZFCEWi2w7IVEOeIEgc4KSkhKUlJRQj+3AoSuJqGKrQSFImvg/uFwuNDQ0wC1x4x9D/0B+8swdgpfzQu/Q45KKSyCXzLwfbAIKuUvis94cCsSNsKysDEVFRXFXrezW4TvvvIPHHnsMer0ea9aswWWXXYaNGzfSzdN5iAVC9vl8UWmLOY6jyolFixZRKVm8A7tQJMhqkGUyWchA0USDdSObnJwEx3FIT09HZmamXzsh0SCr5hqNxs8jOZgVKJ9lkVhAorNIRFak1x1KGRNr2Gso/4lp5zT+OfxP5Kfkw+Qywewyw+qx4uraq6GUhe6ns0NhUq2SlBH2/WNJmrjieb1eNDQ0CP7Zj4yMYPPmzaioqMA999yDoaEhdHR04Nxzz+UlOT1FcXoTMvFENhgMGBsbiyiVITImq9WKtWvXUkvMRA3siDctIUKbzYaUlBTk5+fTVIhE3w4TuFwuHD9+HFarFfX19UhNTQ1Lgom4XTcYDOjt7aV5cpESocMti8RyfqS9NTExgbq6urhWjtkcwWhImvhPKBQK1NbW+v13j8+DPZo9GDANQGPVwOV1IU2ehrMLz8bawrUQi/jfxYRaGCEzBoPBgOrqasGrVZ/PhxdeeAFPPfUUtm/fjvPPP38+9opDYYGQyZpwX19fSG0x0fR6vV7U19ejq6sLS5cuhUwmC7vYIQTYQNHKyspZPV/ihcxW0kJLjEhFGkk5wJIguZCwt+vkT7QDLxKX5Xa7eevFQ51fuHZCuPMjFwOyYJOIFk04kk5LS4PdbofRaERDQ0NIdYHdbcfTx56GDz7kKHJQnFoMnUOHc4vPRW5yfJ4OxArA6/UiNTUVVqt11lZfPBfhoaEhbN68GdXV1bj//vvDJrTMU5zehOz1eqkHcTDPCRLhYzAYUFdXh5ycHHAch8OHD8Pn8yE7Ozsh8jHAX7lQW1sbUsIVODQkt5vsUClWZYdOp4Narea92BDp/MifYDagwY5NVr7HxsZQVVXFOxggnvNjK0Ey0BwfH6e35nMdI0/aE2q1GjKZjJphhaqkfZwPrx9/HTnJObQinrJN4eyis7EoZVHM5zA6Okp1zawHSbCLCMnr47sR6fP58Pzzz+OZZ57BAw88gPXr159OVTGLBUImVfKBAweo5wRZd9VoNKioqPCzryR/Agcj0QzlwoEEmZI161iUCz6fz+8HYjab6cYXqaTDydvsdjt6e3sBAHV1dYKTELkdJlU0cZhjLyKkRRLPxSDe8xscHIRWq4VMJoNEIhFcnRAJrP9EY2MjvSiHuogQEhzwDmDcM46clBzYPXZIRBKsL12PJEn0ZlVkEYro6vnc3YRbtmElbikpKRgcHMSPfvQjNDQ04L777hN8QH2K4fQmZJ/PB7fbTbXFa9eupW5YBQUFqKiooI5YkfrEwYZyrP440tCLXSrg6zkQDdi1XKPRGFTeJpfLMTAwAJ1Oh9ra2oTaMgaCtI50Oh00Gg08Hg+Sk5Ppe0fuROaicjKbzeju7vYzIgq8iBB1QqB8USiSjsV/gm5rGg3o1HZC69IiMyUTK/JWoDC7MKp2UaxmQOGOx66tP/XUU3j77bdhsViwceNGXHHFFTj33HPj2mokiBQy+vHHH2PLli04fPgwdu7ciSuvvJL+N4lEgmXLlgEAysrK8Oabb8Z9PlFggZBJWsjHH39Mrf6IpjfegR0ZepEfMRl6sf1eqVQaMVA0USBifaPRiKmpKVitViQnJ9NU64yMjDnzwGBNeIh6Iph5EXsnIrQvBqlITSZTxJRpcs5EohWq0o/WplRo/4lQ7Rj2HANJ2mKxoKurC5mZmYKYAQWiv78fmzdvxuLFi3HTTTehu7sbHR0d2Lx5c9xm7nxCRgcGBmAymbB9+3Zs3LjRj5BTU1Pp5uMJwELIqdVqpVrOVatWQalU8nJi44NgGW7kVm5ychK9vb3Uca6oqAi5ublzemsul8shl8thMBiQkZGBVatWwefz0R/w8PAwTe0IlN8JCbLYkJ+fj9bWVvoeBDNeZ+9EiCl8vCvX7LZjWVkZ7xgrsVhML1wErJf08PCwX7uI3YYMJGniPzEyMiJY0jTxsA5MjWG/g2z+YlpaGqxWKywWCxobGwU3A/J6vXjmmWfwxz/+EQ899BC+8pWvQCQSYdmyZfjmN78pyHPwCRmtqKgAgDnTzguNeUvITqcTR44cQW1tLRwOB5KTk3ll2MUKkn0nl8thMpkgFouxfPlyyGQyGI1GjIyM0B8w249OhDuay+WCWq2G3W6nYaoEgRcRMt3XarXo6+vzG3qRRZFYLiRsckhgsGYoyGSyWTlzbKU/OjpKDW7Y9zBUpWmz2dDT0wO5XC5IjJBEIglK0qSKZgNB2Qp1cHAQGRkZfhekRCBUtNfk5CTUajXkcjlEIhGOHTsWVUhuJPT19WHz5s1YsWIFPv3005h9RiIh3pBRh8NBU6pvv/12XH755Yk4zbgwbwk5KSmJDvKkUin6+vrookMivAh8Pl/IQFG2GmGlbUT7y1aBhGBiIWnWrJ6P7wMbUsoGqLLZfL29vXQll1TS4W7V2faEENVgYCoymyZtMs1EGDmdTr98vtTUVIyOjmJqagp1dXUJNakJljXn8XgwPT2NwcFBmM1melFWq9V+fsiJ7pl7vV6aNt7U1ESJktUhsyG50ZK01+vFU089hRdffJFWxSczBgcHUVxcjL6+Ppx33nlYtmwZqqurT/Rp+WHeEjJZX87MzERDQwMMBgOVGHEcR1UJQgyUtFotjh8/jpycHL9AzWCQSqWzgixdLhftRWs0Gj+CISQYqZWg0+mgUqmQl5eH1atXx1yJsR69pOfH3qoPDQ2FNC4ijnCLFi3C6tWrE3LbKBKJgraLHA4HraK1Wi0lSqPRSD/vuerfm0wmqNVqFBYWoqmpCWKx2M/boa+vD1arFVKpNGE9c9YMqLm52e+47GccLCSXDSwNRdJqtRqbN29Gc3MzPv30U0EGdpEQb8goeWxVVRXWrVuHQ4cOnXSEPG+Heu3t7bjllluo2L65uRmtra1YsWIF5HL5LFWCVCr1m/rzqVItFgtUKhVdsxUqJy/YEgb5cZBzJK0EstgiFotRW1s7Z1pattI3GAyYnp6GSCRCfn4+cnJyYur3xgOn0+m34JOUlET7qUTDLVQ7JhQC/ScifRZsz1woL2mPxwOVSgW73R6VGVAwsCRNfi+33XYbJBIJhoeHcccdd2DTpk2CytnCqSg8Hg/KysqQlZWF7u5ulJaW4q233sKSJUsAAC+88ALuvfdeAMCiRYuwZcsWOtQzGAxISUmBQqGAVqvF2rVrZ/WfE4zTW2VB4Ha70dnZiba2Nuzfvx+ff/45xGIxmpqasGrVKrS2tqKurg5er9ePAEnfOdjAi2+gqJAIlGaZTCZqwl9YWIjCwsKYneXiOSeyblxbW4u0tLRZ76FCofB7D4VWdrD2nJEM20P5OsS7aMNxHMbHxzEwMBC3V3Jg9iJr/hRpsEnMgMrLy1FYWCj4d6G3txc333wz6uvr0dzcjCNHjmBychK7du0S5Ph8VBTPP/887rnnHuj1elxyySV46aWXcPfdd6O+vh4///nP8dRTT+Haa6/F+Pg4MjMzUVhYiM7OTnz22We46aabqNR1y5YtuOGGGwQ5b55YIORg4DgOFosFHR0dlKR7e3uRm5uLlpYWNDc3Y/Xq1cjPz6d+E2yVKhKJ4HA4UFpaKniOHd/zJz/+oqIipKenz5KOsZV+oqpUrVYLtVpNTZhCBaCynh2BnhiRcgMjwWQyobu7G1lZWTFLuAIXbYJZqIbb1iTLFSRlPBFtEbZnzl7o2Cp6YGAAPp8vIWZAHo8HO3bswK5du/Doo4/O8gsXCnv37sU999yDd999FwCwdetWAMAdd9wx67Hf+973cMkll9AK+KWXXsJHH32EJ598EgBw0003Yd26dfj2t7+dkHONAQuyt2AgP7Z169Zh3bp1AGaIY2xsDO3t7Whra8PTTz9N4+ybm5vR3NxMI5vWrVuHwsJCmM1m7N+/n/ZSCQkmMlPOaDSit7cXaWlpaGlpoT/+YNIxo9GIsbGxWR7D8dpXkk0/kUiElStXhm3TsP1eknvHygOnpqaoNCuaTTmy9m6xWLB48eK4bplDydsClRPEy5e90A0PD2N8fFyQ5YpwUCgUyMvL86v+iQ5+bGwMOp0OMpkMqampGBkZEfRupLu7GzfffDPOOuss/Otf/xKsLRcM8agogv1djUYj+DkmGqcdIQeDSCRCUVERLr/8ciqF8Xq96Onpwa7qPLZeAAAb50lEQVRdu3DjjTfSIeCRI0fQ3NxMq2mRSEQrF2LkzW7xCeHf63Q6oVar4XA40NjYGJaAAqVjbJVKJv/EtIitUvlM1AcHBzE5OemXsBwtiDxQqVTOUnaQoVzgEgY7eCXWlOXl5aivr0/IxS+YcoIdynV3d2N6ehpyuZzeSdlsNkGzA/lgdHQUMpkM55xzDqRSqd+y0vDw8Cz1CdnY5AOPx4PHHnsMf/3rX/H4449TxdICEosFQg4BiUSCxYsXw26343//93/R1NQEm82GQ4cOob29HY888giOHTuGtLQ0StCtra1YtmyZn3Z2aGgILpfLjwD5uraxErJYDXhCVansRJ3Inlg/DFbaRirZgoKChKgn2Kk/QWCVSnrmCoUCpaWlc9K3Z0GGbVNTUwCANWvWQCaT0Yvx+Pg471iqeMCu4RNTLILk5GQkJyf7fc6EpA0GAwYHB3nZqHZ1dWHz5s346le/ik8//TShVTGLeFQUxcXF+Oijj/z+LrkDPpVw2vWQhQTHcdDpdGhvb8e+ffvQ3t6O4eFhlJWVobW1lbY7MjIy6AIGmfiz0jtCgOwPl/RoiSVkonvVbC+VjaNyu92Qy+Wora1FZmbmnDt1kSgpnU6HqqoqiMVi2pMWYpOPL4j/RKT0DNLvJefodDoFaxnZbDZ0dXVBqVTGHHIaykt6YGAAnZ2dMJvNOHDgAJ566im0trbGdJ7BEMmDwul04jvf+Q7+9re/YcmSJXjllVfwzW9+E/fffz8uu+wy1NfXA5i5ED7xxBOzesh6vR7Nzc04ePAgAGDVqlXo6OiYU8+WCFgY6p0I+Hw+9PX1Yd++fdi3bx8OHDgAq9WKxYsXo6WlBS0tLVi+fDmkUmlQ6V1ycjLMZjOSkpJQX18/55aQwJckODU1hYKCAnAc50eAgfLARIFoqwsLC1FaWhq0Mmc13GTgRQiQnGc8LSOn04nu7m6IxWLU1dVF/XrZKjXYYJPPyjrHcRgaGsLY2BgaGhoEz53jOA4ffPABHnjgATidTojFYthsNvzyl7/Exo0b4z4+H/XEjh07cPjwYWzcuBE33HADzGYz7rjjDlxzzTVYvXo1nnnmGWzcuBH79+/H17/+dRgMBiQlJaGgoACdnZ0AgOeeew6//e1vAQB33nknrrvuurjPXUAsEPLJApfLhcOHD1OSPnLkCORyOZqamihJZ2dn45133kFVVRWNWCfSu2gWROIFaU+EIsFA1YTT6RRMNUFA1q45jqOaYr4IJECj0ejXMw9luhPsOMS8Xyj/CfbYoWKfWJ25VCqlZkBZWVkJCXp1u9146KGH8NZbb2HHjh1oaWkBMEOi5LONF3zUExdccAHuuecerF27Fh6PBwUFBZiamsLg4CAuueQSHD16NO7zOMFYUFmcLJDL5ZR4f/jDH9KKc//+/di7dy+uu+46DA4OoqWlBatWraKPzc3NpQSo0+mo1wRRJBACFKKnS3wfZDIZmpqaQlaCCoUC+fn5Qf0wgqkmojlHkmAyOjoaMwmKRKKgvVSbzQaj0Uh9HVgj/YyMDD97TUKCifKfCDbYDNbXt9vt4DgOxcXFYfXVseLo0aO4+eab8W//9m/45JNP/D5ziUQi2PYdH/UE+xiypKXT6QDMOMg1NTUhPT0d9957L8455xxBzutkxAIhnwCIRCJkZGTg/PPPR39/P84991z8/Oc/h8ViQVtbG/bt24dHH30Uer0edXV1lKCJzIyQy8jICO31xiq983q96O/vh06ni8n3IZQfBlm1Zs8xnIm+0WhET08PsrOzBSdBlgDZVWFyjhqNhio7OI6Dx+NBTU0N8vPz58w1jHVvUyqVMBqNKC0tRXZ2NiwWi985xmMBCsxUxQ8++CB2796N3//+91i1alUiXpIgKCwsxNDQEHJyctDR0YHLL78cnZ2dgrvVnSxYIOQTjBtvvJH+c1ZWFkpLS6ldocfjQVdXF9ra2rBr1y78v//3/8BxHFasWEFJmkRTBUrvSK83lIE+x3G0oi0uLkZra6tg5MP6XJSUlACA3yYk8XIg2lmS3xavpjjWcwRmhkI9PT3IyspCcnIypqamMDAwkFCP5kB4vV7q2bx06VJqBsT2jFn1SThfkVDneOTIEdx888246KKL8MknnwiW1h0OfNQT5DElJSXweDwwGo3IycmBSCSi393m5mZUV1ejt7eXtlbmGxZ6yKcQyG1tR0cH2tvb0d7eTkmEld4VFRX5DbuMRqOf9E4ul2N8fJyqJ+bKqD7wtYyMjGBgYADp6enw+Xx0IMcODRPdM3e5XFCpVCH9J4L5TRBlBzlPIWRhxDe6uLgYJSUlUftXEDMtco7shUSv16O6uhoPPvgg3n//fTzxxBNYuXJl3OdMEElBYbVasWjRIuTk5CA/Px8WiwWvvvoqlixZgq1bt+LZZ5+F2WxGS0sL3nrrLezcuRN//etf8corr2BqagrZ2dmQSCTo6+vDOeecgyNHjpxM6gm+WBjqnQ4gBuxkYLh//36MjY2hsrKSGio1NTUhLS2NutK5XC6aJRe4fDEXt+hWqxXd3d1QKpWorq6mpMsO5Ai5uN3uWf1oIdoZ7Ao6H6tSFsGkbbEuYLjdbqhUKjidTkGDVsmFRKvVYvPmzVCr1VAqlbjiiiuwbt06fO1rXxPkefgqKP7+979DrVbDaDQiPz8fR44cwQ9+8AO888476OnpQX9/P1paWlBQUIDs7Gzs3LkTVVVVeO2113D33XdDJpNBLBbjl7/8JS699FJBzn2OsUDIpyt8Ph9UKhXa2trQ3t6Ojo4OTExMwOPx4Fvf+hYuv/xyLF26lErv2Cw51oRdaF0v6Vfr9XrU19fzWu4gdwWssgOYvcUXzYXEZrOhu7sbycnJqKmpibsKZ+0/WW1vpG1Iom1OlBmQ0+nE/fffjw8//BBPPvkkFi1ahI6ODuh0Onz3u98V5DniUVBs27bN77Hs4+YhFlQWpyvEYjHq6+tRX1+PTZs24dprr0VDQwOuuuoqqFQqPPHEEzh69ChSUlL8VB2LFy/26/WShA7SRohHekf61UVFRWhtbeVNPuywi/VnJn3UgYEBvzw+co7B1pjJ5uPExISgxvWssqOgoABA+G3IlJQUGAwGSCQSrFq1KiEto88//xw//vGPcfnll+Pjjz+mn5lQlTFBPAoKjUaDNWvW+P3dU9F/QkgsEPJpgIcffngW+XAcB4PBgP3796OtrQ2vvvoqBgYGUFJS4jcwJJpok8kUk/TO4XCgp6cHIpEorJwuGoTymgi1xpyRkQGRSIS+vj7k5uYKOsAMhWCZd8QPZHh4GCkpKXC5XPjiiy+CenbECqfTifvuuw+ffPIJnnvuOZqyvIBTAwuEfBogWCUoEomQnZ2NCy64ABdccAGAL/2N9+3bhz179mD79u0wm81+Bv/Lly+HQqGgbQQixyLSO9YNbWRkBGNjY3GZEfFFsDw+h8MBg8GA48ePU19hs9mMoaEh3gsiQsHhcKC7uxsymQxr166lFSu5IyFJIuxALly1HwwHDx7Eli1bcMUVV+Cjjz6ak4SUeBQU8SaAzEectj1kvV6Pb33rWxgYGEBFRQVeeeWVoMR14YUXoq2tDWeffTb+/ve/n4AzPbFwu904evQo7UcfPnwYEonEz+C/trYWAGirQ6vVwmg00iWSzMzMhK9ZBwNJXib+EwCohpuQIGtQTxZEhKye2cRpvhcm1kKVmNSHC3Z1OBzYunUr9u7diyeffJImaMwFPB4P6urq8M9//pPKJ//yl7/4ncPjjz+OI0eO4IknnvBTUHR2duLf//3f0d7ejtHRUaxfvx4qlWrOPcbnCAtDvXC47bbbkJ2djdtvvx3btm2DwWDAfffdN+tx//znP2Gz2fDkk0+eloQcCI7jYDab/Qz+SZbf4sWL0d3djXPOOQf/9V//BalUOkuNwK4HJ6pCJW0S0ksPp3ggCyLkHIVMBidmQKmpqaipqYmLaAL9MJxOJ5599lkoFAp8+umnuOaaa3DXXXclrOIPV8C8/fbb2LJlC7xeL1auXInDhw9Dp9Ph2muvxUMPPUQDHYxGI6RSKcrLy7Fnzx7k5+fjN7/5DZ577jlIpVI89NBDuOiiixJy/icBFgg5HOrr6/HRRx+hsLAQY2NjWLduHXp6eoI+9qOPPsL27dsXCDkEOI7DI488gv/5n/9Ba2srtFottFotamtrqePdqlWrkJKSMisrkE20jld6J5T/BJsXaDQao9YecxxHo63q6+sFNwMCZoIC7r77bhw5cgRVVVXo7++HxWLB3r17E7LswaeA0ev1aGlpwYEDByASidDc3IyOjg5kZWVh3bp12L59+7xd6OCBBZVFOExMTNBV34KCAkxMTJzgMzp1IRKJUFFRgcOHD9PNN6/Xi+7ubuzbtw9/+9vfcPfdd8Pr9WL58uV0iWXVqlUQi8VB0zlYg38+0jviP5GZmRn36nWwZHBWezwyMkK1x4GZi6wZUKKGhyTA99vf/jYefvhhWhX7fL6EDSvfeOMN6je8adMmrFu3bhYhv/vuu9iwYQN93zZs2IDdu3efTDFKJz3mNSGff/75GB8fn/Xvf/Ob3/j9f5FINOc+v/MNl112md//l0gkWLJkCZYsWYLrr78ewMwt/MGDB9He3o6HHnoIXV1dSE9P99syDJTeEcVEqA0+om02GAxoaGjwM7kXEoExSqyvsFarpQM5juNQVFSE3NxcRHn3GRF2ux333nsvDh48iBdffBENDQ1+/z2RyhE+BUykGKXrrrsOEokEV1xxBe66666F31wQzGtCfv/990P+t0WLFmFsbIy2LIh7mRDgY8Z97bXXoqOjAzk5OXj55ZdRUVEh2POfrEhJScHZZ5+Ns88+G8AMqWm1Wmrw/+c//xkjIyMoLy/3k96lp6fTClWv12NgYAAejwcymYyu5TY1Nc2ZYgLwN1VKTk7G9PQ0ysrKkJ2dDbPZ7BdFxfajY5W1tbW14b//+7/xne98B9u3b0/I4CuRBcyLL76I4uJimM1mXHHFFfjTn/6Ea6+9Nq7znY+Y14QcDhs3bsQLL7yA22+/HS+88MKsCi9WeL1e/PCHP/RbJd24caPfKumzzz6LrKwsqNVq7Ny5Ez/72c/w8ssvC/L8pxJEIhHy8vJw8cUX4+KLLwYwc9t9/Phx7Nu3D++99x62bt0Km83mZ/Cfl5eH119/HWeffTbKyspgs9nQ0dHh5ygXzzCOL7xeL9RqNSwWi58ZULDAVKPRiP7+fhpEENiPDnWeNpsNv/71r/H5559j586dqKurS9jribeACRejRORsaWlpVFmxQMizcdoO9XQ6Ha666ioMDQ2hvLwcr7zyCrKzs3HgwAE88cQTeOaZZwAA55xzDrq7u2GxWJCTk4Nnn32W6naDIZ5V0oVbuOAgCxRtbW3YuXMnjh07hhUrVqCxsZGSdHV1NfWZZlNY2GGckNI7YgZUUlKC4uLiqD47krlIzjMw5UQqlSIzMxOfffYZbr31VmzatAmbN29OmByMjwT01ltvxRtvvIHJyUkUFBRg48aN+N3vfgdgxq/46quvxuTkJN3IlMvlNEYpPT0d09PTyM3Nhdvtxre//W2cf/75+M///M+EvJ6TFAsqixOBV199Fbt376aE/qc//Qn79u3DY489Rh+zdOlS7N69m1pTVldXY9++fYKmUsxHvPnmm9i9ezeN6dm/fz/NMuzr60NRURHVRhODf7fbTbMMhZDesWZAjY2Ngji9Baac/OxnP8OxY8fgcDhw/fXX4+KLL8batWsT1pLho6DQ6XTYsGEDxsfH4XK50NvbSwuYq666Cr/97W9x9dVXY926dTh27BjS09NpjJLVasVXvvIVuN1ueL1enH/++XjwwQfnq944FBZUFguYX9i4caNfxtuGDRuwYcMGAF/mzu3btw979+7FI488AoPBENTgnwzjSHoIK70LZ/pOzIAqKipQUFAg2B0N8cJISkqCSqXCxMQEbrnlFqxfvx4dHR14+eWX0dramjBC5qOgyMnJwcGDB6kElCgpmpubYTKZaNjo1q1b/e4QAUCpVKKjoyMh5z7fsEDIAiOeVdIFxA6RSITy8nKUl5fjqquuAjCjJ+7s7MS+ffvw8ssv4/bbb4dIJJpl8C8SiWaZvrPSu+TkZAwMDFBtbSJ0vlarFffccw+6u7vx6quvorq6GsDM3dSmTZsEfz4W8UhAdTodMjMz6cViwSAoPsxNPs1phNbWVqhUKvT398PlcmHnzp2zknvJQBGYaXGcd955UVVbu3fvRn19PWpqaqiFIYuPP/4Yq1atglQqxauvvhrfCzqFIZVKsWLFCvzHf/wHnn76aezfvx8ffvghvvOd78BkMuG+++7DWWedhUsuuQTbt2/HwYMHkZ2djTPOOAMrV65ETk4OxsbG0NHRQYNIh4eHodVq4XK5BDlHjuPw8ccfY8OGDVi8eDHee+89SsZC4vzzz8fSpUtn/XnjjTf8HrcgAT2xWKiQBYZUKsVjjz2GCy64AF6vF9dffz2WLFmCu+++Gy0tLTTm/Lvf/S5qamqoGTdf8FFxlJWV4Q9/+AO2b9+eiJd4yoI4sK1bt45O/zmOw8TEBDX4f+655zA+Po7CwkJMTk7ivPPOw09/+lNkZmbSFBaDwUCld0ql0i8pOpq+qMViwd133w21Wo3XX38dlZWVCXrliZOA5uTkYHp6Gh6PB1KpdMEgKE4sEHIC8LWvfW2W7+yvfvUr+s9JSUnYtWtXTMdub29HTU0NqqqqAABXX3013njjDT9CJprmuQroPJUhEolQUFCAyy67jEof//jHP2Lbtm249NJLYTabceWVV8LlcmHp0qXU9W7FihWQSqXU9Y7ojiOFuQJfVsW33347fvCDH2DHjh0n9LOKRwIqEolw7rnn4tVXX8XVV18tqIT0dMQCIZ9i4GMIvoD4sGzZMhw4cAApKSn03zkcDnz++edoa2vD448/js7OTiiVSj+D/4aGBmq+RMJcbTYbZDIZ0tPTcejQISxZsgRPP/00BgYG8MYbbyR0IYivo2FHRwc+/vhj/PrXv8batWvxyiuvAAAuueQS7Nmzh7ZQxGIxhoeHYbFYUFJSQiWg9913H66++mrcddddaGpqwg033JCw1zTfsUDIC1hAAJqammb9u6SkJKxZs4YmXHAcB71eTw3+X3nlFQwODqK0tHSWwb/b7YbBYMDu3bvxq1/9Ci6XC62trXjllVdw6623Jqxnu23bNqxfv57K2bZt2xbU0fDOO+/ET37yk1mOhrm5uXj++eepgiIUqqqq0N7eLvj5n45YIORTDHNl6h1p/fvBBx/EM888A6lUiry8PDz33HMoLy8X/DxOVohEIuTk5ODCCy/EhRdeCGBmy3BgYABtbW348MMPcf/998NsNqOurg6Tk5NITk7GRx99hNLSUhw/fhzHjh1L6ACNj5wNANavX++3YbeAEwiO46L5s4ATDLfbzVVWVnJ9fX2c0+nkli9fzh09ejToYzdt2sTt2rUr6ufweDxcVVUVd/z4cfocnZ2dfo/54IMPOKvVynEcx+3YsYO76qqron8xpwFcLhd34MAB7he/+AXn9Xrn9LkzMjLoP/t8Pr//H4gPP/yQu/jii/3+3aZNm7i6ujpu2bJl3JYtWziHw5Gwcz0NwItjF6Y+pxhYFUdjYyOuuuoqquJ48803AcxssJWUlGDXrl246aabok6QYAeHcrmcDg5ZnHvuubTHumbNGoyMjAjzAucZZDIZmpubcc899yRkcJdIOdvWrVvR3d2N/fv3Q6/XB62uFyAsFloWpyAiqThaW1vjIshoB4fPPvvsfE56OKmRSEdDsiyiUChw3XXXLcgo5wALFfIC4sKf//xnHDhwALfeeuuJPpUFBIBdQIpFjjY2NgZgpq35t7/9DUuXLhX8HBfgj4UKeQGzwHdw+P777+M3v/kN9uzZM+cBpqcz+MjZPv/8c3z44Yfo6urCL37xC9TV1WHPnj0AZoZ9N954Iw0H0Gg06OnpmSVnu+aaazA1NQWO47By5Uo88cQTJ+Llnl7g22zmFoZ6pw34DA4PHjzIVVVVcb29vXE91zvvvMPV1dVx1dXV3NatW2f999///vfc0qVLuRUrVnBnnXXWrOHi6Yhbb72Vvldbt27lbrvttlmP6enpoZ+NRqPhCgoKOIPBwHEcx33zm9/kXnrpJY7jOO6mm27iduzYMUdnflqDF8cuEPICguKtt97iamtruar/3979hFRWhnEc/z4h7orAQUYo0BFdWHgkM2lhis7V4OIQJOHGP4t2d9CFAxWDm+Eqrdy5iaBMchFJXrFFUDIrGf+AwhSBobu0WmjqanTw12KOYv7Je/Xe6x17PiCe976v933P5vFw3uc859YtxeNxSVJ/f78SiYQkqampSYWFhQqCQEEQqLW1NeU5ksnm2NraOjxOJBJqaWm5xFldD+Xl5VpbW5Mkra2tqby8/Ny/qays1PLysvb391VQUKC9vT1J0szMjJqbmzO6XicpyRjrtyzcqc7bOPyvzaRkJfMY+MFLU+FZRTQvfJN6dba5uTl2d3cpLS316mw5zgOyuzLJZnMMDw8zNDTE7u4u09PT2VzilUnX++3W19fp6OhgZGTEa5s8Bzwgu5wXi8WIxWKMjY0Rj8cPMweus3Sks21vbxONRhkYGDh85Nurs+U2/5fprkyqj4G3t7czMTGR8jzn1Y8+MD4+jpmxsLCQ8hzpsLGxQSQSoaysjEgkwubm5okxS0tL7OzsUFVVRWVlJX19fYfpbN3d3ZSUlBz2RSIROjs7/1WL4mh1NrhYOpzLoGRvNss39VyaJZPNcTSLY3JyUtXV1SnNkczGoSRtb2+rrq5OtbW1mp+fv9gJXVKy2ROzs7NqbGxUcXGx8vPztbq6KkmKRqNqbGyUJI2OjiovL+9w0zUIAi0uLkqSVlZWVFNTo9LSUrW1tfkj0dnhWRYu952XzdHT06OKigoFQaCGhoYz63ac5XgWweDgoAYHB0+M6+3t1dTUlOrr668sIF8me0K6eO0SlxVJxdhU3zrt3HPFzNqAdyV9GLY7gFpJd4+MeQO4L+l9M3sI3JOU9fsWZva3pJfDYwM2D9pnjH8LGAFek7RvZl8CbwNPgJ+AjyU9yfzKXbr4pp77XzOzF4AhoDtL8/0I3Dyl6/7RhiSZ2ZlXS2ZWBIwCXZL2w48/Af4A8oHPgI+AB6d/g8tFHpDddfc78OqR9ivhZwdeBF4HHobpYzeBSTO7k4mrZEm3z+ozsz/NrEjSehhw/zpj3EvA9zy7qn905LvXw8MnZvYFcC+NS3dZ4FkW7rqbB8rMrMTM8oF2YPKgU9KWpBuSiiUVA4+AjATjJEwCXeFxF5A4PiA8h++AryR9e6yvKPxtwHvAzxldrUs7D8juWpP0FLgL/AD8Cnwj6Rcze2Bmd652dSd8CkTM7DfgdtjGzN40s8/DMR8A7wDdZrYU/lSFfV+b2WPgMXADiGd3+e6yfFPPOedyhF8hO+dcjvCA7JxzOcIDsnPO5QgPyM45lyM8IDvnXI74B3LhXlpyhh9zAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter3d(tab_merge)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Hail",
   "language": "python",
   "name": "hail"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}