{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulate the Sequences\n", "This notebooks reads a phylogenetic tree and simulates sequences down the tree." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "**Rule 2: Document the Process, Not Just the Results.** Here we describe the steps how to produce the dataset.\n", "\n", "**Rule 4: Modularize Code.** To keep the code clean, we have written helper functions in a separate Python script, [helper.py](./helper.py).\n", "\n", "**Rule 7: Build a Pipeline.** This notebook describes the entire workflow, and its modularity makes it easy to change models or model parameters.\n", "\n", "**Rule 8: Share and Explain Your Data.** To enable reproducibility we provide a `/intermediate_data` directory with files produced by the workflow.\n", "\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define Sequence Simulation Parameters\n", "The Jukes-Cantor (JC69) model [(Jukes & Cantor, 1969)](http://dx.doi.org/10.1016/B978-1-4832-3211-9.50009-7) does not have any parameters, as it assumes fully neutral sequence evolution. We will need to simulate a root sequence, so we will specify that we want a sequence of length *k* = 300 (similar to the length of *Alu* sequences).\n", "\n", "We will also include import statements here to keep the notebook clean and organized." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "INPUT_TREE_FILE = \"./intermediate_data/dualbirth.tre\"\n", "OUTPUT_SEQUENCE_FORMAT = 'fasta'\n", "OUTPUT_SEQUENCE_FILE = \"./intermediate_data/sequences.fas\"\n", "K = 300" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from warnings import simplefilter\n", "simplefilter(action='ignore', category=FutureWarning) # suppress SciPy warning\n", "import matplotlib.pyplot as plt\n", "from pyvolve import Evolver,Model,Partition\n", "from pyvolve.newick import read_tree\n", "from random import choice\n", "from seaborn import distplot,regplot\n", "from treeswift import read_tree_newick" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulate the Root Sequence\n", "We will be evolving a root sequence down the previously simulated phylogenetic tree, but we first need to create our root sequence. For the purposes of this notebook, we will simply randomly generate a sequence of length *k* = 300 with equal probabilities of each nucleotide." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "root_sequence = ''.join(choice('ACGT') for _ in range(K))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulate the Sequences\n", "We will now run [Pyvolve](https://github.com/sjspielman/pyvolve) to simulate the sequences. To do so, we need to perform the following steps:\n", "\n", "* Use Pyvolve to read the input tree\n", "* Create a `Model` object, which specifies the parameters of our model (JC69)\n", " * In Pyvolve, specifying a nucleotide model without any additional state frequency or transition rate information will default to the JC69 model\n", "* Create a `Partition` object, which specifies the details of the evolutionary unit (the model and the root sequence)\n", "* Create an `Evolver` object, which specifies the full details of the actual evolutionary process (the tree and the partitions)\n", "* Evolve the root sequence down the input tree" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "tree = read_tree(file='./intermediate_data/dualbirth.tre')\n", "model = Model('nucleotide')\n", "partition = Partition(models=model, root_sequence=root_sequence)\n", "evolve = Evolver(tree=tree, partitions=partition)\n", "evolve(seqfile=OUTPUT_SEQUENCE_FILE, seqfmt=OUTPUT_SEQUENCE_FORMAT, ratefile=None, infofile=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## (Re)Load the Datasets\n", "Pyvolve automatically writes the sequences to file, so we will need to load them to be able to analyze them. Also, Pyvolve loads the simulated tree in its own format, but because we will want to use more advanced features than Pyvolve's tree class allows, we will reload the tree using [TreeSwift](https://github.com/niemasd/TreeSwift). Note that, to reduce clutter in our notebook, we have written our code to parse a FASTA file in a separate Python script containing helper functions, [helper.py](./helper.py)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from helper import read_FASTA\n", "sequences = read_FASTA(OUTPUT_SEQUENCE_FILE)\n", "tree = read_tree_newick(INPUT_TREE_FILE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can explore the simulated sequences by printing a couple of them." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ ">0\n", "GAAGGTGGAACACGTTGACGACTACATTAGCCGTTTAGAAGCTCCTTGTGTACCACGTCATACACATATAAACGTGCTTGGCACTCACAGCGCAACGTATGCCTGGGGAGGGCAGGCTATTTGAGTAATAATGCGGTTGTCATCGAGGAAACTGCGATGTAGTCCGCACCAGTGAACTGCATCAGTAAAGAATTAACTTTGGAGCAGCCCATTGAGAATACGCACTGTGGTTCTGCTATTACCATTATGTCTGAATGACGCACACTTCTACGTCTGCCACCCGAAGCGCAGAGCTCTTCC\n", ">1\n", "GAAGGTGGAACACGTAAACGACTAGATTTGCCGATTAGAAGCTCCTGGTGAAACACTCCACTCTATTACAAACGTGCAGGTCAGTGAAAAATCAAAGAATGTCTGCGATAGACATCACTTTTGAGTAATAAAGGGGGCATGAGCGTGGAAACTGCGACGACCCCAGCATCAGTGAACGGCATCAGACATGAATTCACATTAGAGCAGCACATTGCCCATACGCGCTGTCGTGGTACCATTACCATTATGACGCAATGGCCATCATTTCCACATATGGCACCCGCAACGCCTGGTTCTTCA\n", "...\n" ] } ], "source": [ "for seq in sorted(sequences.items())[:2]:\n", " print('>%s\\n%s' % seq)\n", "print('...')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare Pairwise Distance Distributions\n", "We can compare the pairwise distance distributions between the simulated datasets. Specifically, we will show the pairwise distance distribution of leaves on the simulated tree, and we will show the pairwise Hamming and JC69-corrected distributions of the simulated sequences.\n", "\n", "For two sequences $u,v$ with true JC69 distance $d(u,v)$, the Hamming distance $h(u,v)$ is expected to be *at most* equal to $d(u,v)$, but as $d(u,v)$ increases, $h(u,v)$ will deviate further and further from $d(u,v)$. Specifically, the following relation holds in expectation:\n", "\n", "$h(u,v) = \\frac{3\\left(1+e^{-\\frac{4d(u,v)}{3}}\\right)}{4}$\n", "\n", "This is because, any time an already-mutated site mutates again, the JC69 distance increases, but this subsequent mutation is unobservable to us, so the Hamming distance does not change. An intuitive way to help internalize this phenomenon is to realize that, by definition, the theoretical maximum possible Hamming distance between two sequences is exactly 100%, whereas JC69 distance is theoretically unbounded. This relation allows for a simple correction of Hamming distances to JC69:\n", "\n", "$d(u,v) = -\\frac{3}{4}\\ln\\left(1-\\frac{4h(u,v)}{3}\\right)$\n", "\n", "Note that, to reduce clutter in our notebook, we have written our code to perform pairwise-distance-related computations in a separate Python script containing helper functions, [helper.py](./helper.py)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from helper import distance_matrix_to_list,compute_hamming_distance_matrix,compute_jc69_corrected_distance_matrix\n", "distance_matrices = dict()\n", "distance_matrices['Tree'] = tree.distance_matrix()\n", "distance_matrices['Hamming'] = compute_hamming_distance_matrix(sequences)\n", "distance_matrices['Jukes-Cantor 69 (JC69)'] = compute_jc69_corrected_distance_matrix(sequences)\n", "pairwise_distances = {k:distance_matrix_to_list(distance_matrices[k]) for k in distance_matrices}" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xt8znX/wPHXe8yZlGNMCGFmFpOQQ1JRdwelIiWdpHOpOz86SSqKznW7dVLdcivqLqkoORRSaM5FJCZnM6cN296/Pz7XNTPXtuvarmuba+/n47HHdV3f4/uaud7X5yyqijHGGOOPiKIOwBhjzMnDkoYxxhi/WdIwxhjjN0saxhhj/GZJwxhjjN8saRhjjPGbJQ3jNxHpJCK/B+laB0TkzGBcK8D7fi0iNxX2fQuLiJzh+d2WCtL1xonI457nXUUkMRjX9VwvaH9PpvCIjdMoeURkI1ALSAcOAl8D96jqgaKMq6BERIFDgAKHgQRgvKpOzue1mqjqH8GNMv9EZADwDpDi2bQTmAM8p6pr83Gt21T1vADO6Qr8R1WjArlXlvOL3e/UBM5KGiXXZapaCWgNxAOPFeRiIlI6KFEVXCvP+2oKTABeF5EnizakoFroeX+nAN1xCWSJiMQE+0bBKq2Y8GJJo4RT1S24kkYMgIjcLCJrRGS/iGwQkTu8x2avnhCRjSIyRESWAwdF5HYRmZZl/zoR+STL680iEud5riLS2PP8EhFZ7bnnFhF5OMs5/xCRBBHZKyILRCTWz/e1S1U/BO4EhopINc/15ojIbZ7njUVkrogki8guEZns2T7Pc5llnqqe60TkVBH5UkR2ikiS53nmN27PdZ8Wkfme9zFTRKpn2X+eJ/69nt/DAM/2siIyRkQ2ich2T3VQeT/eX7qqrlfVu4C5wHDP9Rp4frelPa8HeP4d94vInyLST0SaA+OA9p73t9dz7AQR+ZeIfCUiB4HzPdtGZr23iAzz/L42iki/bL+D27K8HiAiP+byO83+99Tcc429IrJKRC7Psm+CiLwhItM972WRiDTy7BMReUlEdojIPhFZEYokahxLGiWciNQDLgF+9WzaAfwDqALcDLwkIq1zuURf4FKgKjAL6CQiESJSBygDtPfc50ygErDcxzXeAe5Q1cq45PW955yzgXeBO4BqwL+BL0SkbABv8XOgNHCOj31PAzOBU4Eo4DUAVe3s2d9KVSt5qrcigPeA+sAZuG/4r2e73vW431lNz3t/2PM+6uMS82tADSAOV3UGMAo4y7OtMVAXeCKA9wfwKdAp+0YRqQi8CvT0/G47AAmqugYYhKfUoqpVs72HZ4DKwI8+7lUbqO6J8yZgvIg0zSvAHH6nWWONBKbh/j1qAvcCE7Nduw/wFO7f6w9PnAAXAZ1xv8dTgGuB3XnFZPLHkkbJ9T/PN8wfcd9UnwVQ1emeb7CqqnNx/4lP+EDK4lVV3ayqKaq6AdiP+wDsDMwA/haRZkAX4AdVzfBxjaNAtIhUUdUkVV3q2T4Q+LeqLvJ8s34f11Zxrr9vUlWPAruA03K4b32gjqqmqqqvD0nvdXar6lRVPaSq+3EfWF2yHfaeqq5V1RTgY9zvAdwH8XeqOklVj3qulSAi4nmPD6rqHs91n8V9OAbi7xzeH0AGECMi5VV1q6quyuNan6vqfFXNUNXUHI55XFUPe/4+puM+pAvqXNyXilGqekRVvwe+xH0p8fpMVX9W1TRgIsd+v0dxSa4Zrp12japuDUJMxgdLGiXXlapaVVXrq+pdng86RKSniPwkIns8SeUS3DfLnGzO9nou0BWXNObiGmq7eH7m5nCNqz33+ctTXdTes70+8JCnumKvJ556QB1/36TnG2wNYI+P3Y8AAvzsqQ65JZfrVBCRf4vIXyKyD5gHVJXj6/23ZXl+CPchiCfm9T4uWwOogGuT8L6/bzzbA1EXH+9PVQ8C1+FKFVs9VTvN8rhW9n/P7JI81/X6iwD+PXJRB9ic7UvFX7j35uXz9+tJMK8DbwA7RGS8iFQJQkzGB0saJpOn2mcqMAao5am2+Ar3wZqT7N3vvEmjk+f5XPJIGqr6i6pegauW+B/uWzq4D7BnPMnN+1NBVScF8LauANKAn33cd5uq3q6qdXBVYG+Kp53Fh4dwjevtVLUKLilC7r8br81AIx/bd+GquVpkeX+neBq6A9EL+MHXDlWdoaoXAqcDvwFveXflcK28ulOe6qn28joDV9IB1xOvQpZ9tfO4VlZ/A/VEJOtn0hnAFn9OVtVXVbUNEI2rpvpnAPc2AbCkYbIqA5TFdeVME5GeuPriQMwFzgfKq2oi7sOsB65N4tfsB4tIGU/j7CmeqqR9uCoVcB9wg0Sknaexs6KIXCoilfMKQkRO8zTSvgGMVtUT6rhF5Bo51pidhPvA9N57O5B1HEll3Af8XhE5DQikR9ZEoLuIXCsipUWkmojEeb5Vv4VrN6rpiamuiFzsx/srJSINReQ1XJJ+yscxtUTkCs+H/GHgQLb3FyUiZQJ4H15Pef7dOuHav7ydHRKAqzylssbArdnOy/47zWoRrvTwiIhEiuveexnw37yCEZG2nr+RSFziSuXY+zRBZknDZPLUqd+H+6afhKuL/yLAa6zFfTj94Hm9D9gAzFfV9BxOuxHY6Kn2GQT085y7GLgdV/WQhGv8HJBHCMtE5IDn2Ntw7QU5NSy3BRZ5jv8CuN/TLgOuN9L7nmqja4GXgfK40sFPuGokv6jqJlz120O4aqQEoJVn9xBPrD953v93uBJNTtp74t2Hq/qrArRV1RU+jo0ABuO+xe/Blfbu9Oz7HlgFbBORXf6+F1wVUZLnmhOBQar6m2ffS8ARXHJ437M/q+Ec/zvNpKpHcEmiJ+53/CbQP8u1c1MFl3yTcFVau4EXAnhPJgA2uM8YY4zfrKRhjDHGb5Y0jDHG+M2ShjHGGL9Z0jDGGOO34jLJXFBUr15dGzRoUNRhGGPMSWPJkiW7VNXvAaVhlTQaNGjA4sWLizoMY4w5aYjIX4Ecb9VTxhhj/GZJwxhjjN8saRhjjPFbWLVp+HL06FESExNJTc1plmdj/FOuXDmioqKIjIws6lCMKTJhnzQSExOpXLkyDRo0wC1fYEzgVJXdu3eTmJhIw4YNizocY4pM2FdPpaamUq1aNUsYpkBEhGrVqlmJ1ZR4YZ80AEsYJijs78iYEpI0jDHGBIcljUJQqdLxC7FNmDCBe+65p9Du//fff9O7d+9Cu58pOZ5/Htq1A1thoeSwpFEC1KlThylTphR1GCYMff01/PwzbNiQ97EmPFjSKGLTpk2jXbt2nH322XTv3p3t27cDMHz4cG666SY6depE/fr1+fTTT3nkkUdo2bIlPXr04OjRo4CbOmXo0KHExcURHx/P0qVLufjii2nUqBHjxo0DYOPGjcTExACulHPVVVfRo0cPmjRpwiOPPJIZyzvvvMNZZ53FOeecw+23316opSFz8lGF5cvd89mzizYWU3jCvsvtcR54ABISgnvNuDh4+eVcD0lJSSEuLi7z9Z49e7j88ssBOO+88/jpp58QEd5++22ef/55xo4dC8D69euZPXs2q1evpn379kydOpXnn3+eXr16MX36dK688koAzjjjDBISEnjwwQcZMGAA8+fPJzU1lZiYGAYNGnRCPAkJCfz666+ULVuWpk2bcu+991KqVCmefvppli5dSuXKlenWrRutWrU64VxjvLZsgT173PPZs+G224o2HlM4SlbSKCLly5cnIUuymjBhQubEiomJiVx33XVs3bqVI0eOHDcGoGfPnkRGRtKyZUvS09Pp0aMHAC1btmTjxo2Zx3kTUMuWLTlw4ACVK1emcuXKlC1blr17954QzwUXXMApp5wCQHR0NH/99Re7du2iS5cunHbaaQBcc801rF27Nri/CBNWli1zj/Xru6ShCtbBLPyVrKSRR4mgKNx7770MHjyYyy+/nDlz5jB8+PDMfWXLlgUgIiKCyMjIzC6fERERpKWl+TzO+9zXcdmPByhVqpTPY4zJi7dq6u674ZFHYO1aaNq0aGMyoWdtGkUsOTmZunXrAvD+++8XWRxt27Zl7ty5JCUlkZaWxtSpU4ssFnNyWLYMGjQATy2ptWuUEJY0itjw4cO55ppraNOmDdWrVy+yOOrWrcuwYcM455xz6NixIw0aNMiswjLGl+XLITYWGjeGunUtaZQUomHUwTo+Pl6zL8K0Zs0amjdvXkQRnVwOHDhApUqVSEtLo1evXtxyyy306tWrqMMqVuzvyUlJgUqV4NFHYcQI6NHDNYr//HNRR2YCJSJLVDXe3+OtpGEyDR8+nLi4OGJiYmjYsGFm7yxjslu9GjIywNvBrkYN2LWraGMyhSNkDeEi8i7wD2CHqsZ4tk0GvE1lVYG9qhrn49yNwH4gHUgLJAua/BszZkxRh2BOEt5G8NhY91i9uiWNkiKUvacmAK8DH3g3qOp13uciMhZIzuX881XV/gyNKYb+/ts9nnGGe6xeHfbvh8OHIUvnPBOGQlY9parzgD2+9onrO3otMClU9zfGhM6uXVC58rEE4e3DsXt30cVkCkdRtWl0Arar6roc9iswU0SWiMjA3C4kIgNFZLGILN65c2fQAzXGnGjXLqhW7dhrb9KwKqrwV1RJoy+5lzLOU9XWQE/gbhHpnNOBqjpeVeNVNb5GjRrBjtMY48Pu3ccSBVjSKEkKPWmISGngKmByTseo6hbP4w7gM+CcwokuNLJPjZ5dgwYN2BXk/23btm2jT58+NGrUiDZt2nDJJZfke1qQCRMm8Le3ErsAPv74Y6Kjo2nRogXXX3995vYhQ4YQExNDTEwMkyfn+GfBAw88wLx58wDo2rVr5lQsBw4c4I477sh8r127dmXRokUA7N27l969e9OsWTOaN2/OwoULAVi2bBnt27enZcuWXHbZZezbtw+AFStWMGDAgAK/13C3a5cljZKqKEoa3YHfVDXR104RqSgilb3PgYuAlYUY30lPVenVqxddu3Zl/fr1LFmyhOeeey5zBt1A5SdppKenH/d63bp1PPfcc8yfP59Vq1bxsmdKl+nTp7N06VISEhJYtGgRY8aMyfwAz2r37t389NNPdO58YqHztttu47TTTmPdunUsWbKE9957LzMJ33///fTo0YPffvuNZcuWZY6xuO222xg1ahQrVqygV69evPDCC4CbvysxMZFNmzYF9H5LGksaJVfIkoaITAIWAk1FJFFEbvXs6kO2qikRqSMiX3le1gJ+FJFlwM/AdFX9JlRxFpY5c+bwj3/8I/P1Pffcw4QJE447JiUlhZ49e/LWW28B8J///IdzzjmHuLg47rjjDtLT00lPT2fAgAHExMTQsmVLXnrppRPuNXv2bCIjI4+b4bZVq1Z06tSJAwcOcMEFF9C6dWtatmzJ559/Drjp05s3b87tt99OixYtuOiii0hJSWHKlCksXryYfv36ERcXR0pKCrNmzeLss8+mZcuW3HLLLRw+fBhwJaYhQ4bQunVrPvnkk+Nieuutt7j77rs59dRTAahZsyYAq1evpnPnzpQuXZqKFSsSGxvLN9+c+M89derUzAkbs1q/fj2LFi1i5MiRRES4P+eGDRty6aWXkpyczLx587j1VvenV6ZMGapWrQrA2rVrMxPQhRdeeNy0KZdddhn//e9/T7iXOSZ7m4ZnnkusWTH8hazLrar2zWH7AB/b/gYu8TzfAIRkTu4HvnmAhG3BnRo9rnYcL/co+ESIBw4coE+fPvTv35/+/fuzZs0aJk+ezPz584mMjOSuu+5i4sSJtGjRgi1btrBypSt8+ZrFduXKlbRp08bnfcqVK8dnn31GlSpV2LVrF+eee27mLLnr1q1j0qRJvPXWW1x77bVMnTqVG264gddff50xY8YQHx9PamoqAwYMYNasWZx11ln079+ff/3rXzzwwAMAVKtWjaVLl55wX2/VWMeOHUlPT2f48OH06NGDVq1a8dRTT/HQQw9x6NAhZs+eTXR09Annz58/3+fqg6tWrSIuLo5SpUqdsO/PP/+kRo0a3HzzzSxbtow2bdrwyiuvULFiRVq0aMHnn3/OlVdeySeffMLmzZszz4uPj2fUqFHHrTVijjlyxHWvzVrSiIyEqlWtpFES2IjwYuKKK67g5ptvpn///gDMmjWLJUuW0LZtW+Li4pg1axYbNmzgzDPPZMOGDdx777188803VKlSJaD7qCrDhg0jNjaW7t27s2XLlsxqq4YNG2au+9GmTZvjpl/3+v3332nYsCFnnXUWADfddFNmOwPAddddd8I5AGlpaaxbt445c+YwadIkbr/9dvbu3ctFF13EJZdcQocOHejbty/t27f3mQC2bt1KoB0d0tLSWLp0KXfeeSe//vorFStWZNSoUQC8++67vPnmm7Rp04b9+/dTpkyZzPNq1qwZlDaccOXtVpt9qjQb4FcylKip0YNRIsiv0qVLk5GRkfk6NTX1uP0dO3bkm2++4frrr0dEUFVuuukmnnvuuROutWzZMmbMmMG4ceP4+OOPeeqpp7jssssAGDRoEC1atMhxedeJEyeyc+dOlixZQmRkJA0aNMiMJfuU6SkpKQG/z4oVK/rcHhUVRbt27YiMjMxMOuvWraNt27Y8+uijPProowBcf/31mQkpq/Lly5/wOwNo0aIFy5YtIz09/YRkExUVlXlfgN69e2cmjWbNmjFz5kzAlYKmT5+eeV5qairly5cP+L2XFN7EYEmjZLKSRiGpX78+q1ev5vDhw+zdu5dZs2Ydt3/EiBGceuqp3H333YBbKGnKlCns2LEDcKv9eRdLysjI4Oqrr2bkyJEsXbqUevXqkZCQQEJCAoMGDaJbt24cPnyY8ePHZ15/+fLl/PDDDyQnJ1OzZk0iIyOZPXs2f/31V56xV65cmf379wPQtGlTNm7cyB9//AHAhx9+SJcuXfK8xpVXXsmcOXMA2LVrF2vXruXMM88kPT2d3Z6vrsuXL2f58uVcdNFFJ5zfvHnzzHtm1ahRI+Lj43nyySfxTr65ceNGpk+fTu3atalXrx6///474Epv3qov7+81IyODkSNHHtf+s3bt2szlcc2JvIkha5sGWNIoKSxphFhaWhply5alXr16XHvttcTExHDttddy9tlnn3DsK6+8QkpKCo888gjR0dGMHDmSiy66iNjYWC688EK2bt3Kli1b6Nq1K3Fxcdxwww0+SyIiwmeffcZ3331Ho0aNaNGiBUOHDqV27dr069ePxYsX07JlSz744AOaNWuW53sYMGAAgwYNIi4uDlXlvffe45prrqFly5ZERET4XFI2u4svvphq1aoRHR3N+eefzwsvvEC1atU4evQonTp1Ijo6moEDB/Kf//yH0qVPLABfeumlmUkn6+8V4O2332b79u00btyYmJgYBgwYkNnQ/tprr9GvXz9iY2NJSEhg2LBhAEyaNImzzjqLZs2aUadOHW6++ebMa8+ePZtLL700z/dUUln1VAmnqmHz06ZNG81u9erVJ2wrTAkJCdq2bdsijSFcdOzYUZOSkjQ1NVWjoqJ07969Qb9HamqqtmvXTo8ePepzf1H/PRUH//qXKqj+/ffx2x9+WLVcOdWMjKKJy+QPsFgD+Jy1kkYIjRs3jr59+zJy5MiiDiUsjB07luXLlxMXF8ddd90VkkWiNm3axKhRo3yWdoyTW/VUaiocOlT4MZnCY/8zQmjQoEF+Vd0Y/3gbtNesWROyezRp0oQmTZqE7PrhwDtZYZYOZ4BbU8O7P4f+ECYMWEnDGBOQ7KPBvWxUeMlgScMYE5DskxV6WdIoGSxpGGMCYiWNks2ShjEmIJY0SjZrCA+h3bt3c8EFFwBuqvJSpUplToXx888/Hzd1hTEni+yTFXpVrQoREZY0wp0ljRCqVq0aCQlugsThw4dTqVIlHn744eOOyez7HGGFPlP8HT4MBw74LmlERMCpp9qSr+HOPqmKwB9//EF0dDT9+vWjRYsWbN26la+//pr27dvTunVrrrvuOg4ePAjAL7/8QpcuXWjTpg09e/bM95oYxgRDTqPBvU45xc2Aa8JXiSppPPAAJAR3ZnTi4uDlfMyD+Ntvv/HBBx8QHx/Pjh07GDVqFLNmzaJChQo888wzvPLKKzz00EPcf//9fPHFF1SvXp2JEyfy+OOPHzenlDGFKafJCr2qVAEfa2iZMFKikkZx4p1oD2DBggWsXr2aDh06AHDkyBHOO+881qxZw6pVq+jevTvgVsOLiooqspiN8ZY0vIsuZVe5siWNcFeikkZ+SgShknUKcVWlR48efPjhh8cd8+uvvxIbG8sPP/xQ2OEZ45N3zS/PAognqFIFtm4tvHhM4bM2jWKgQ4cOzJ07lw0bNgBw8OBB1q1bR3R0NFu2bOHnn38GXAlk1apVRRmqKeGSk91jTtN+ValibRrhzpJGMVCrVi3eeecdrrvuOlq1akWHDh1Yu3YtZcuWZcqUKQwePJjY2FjOPvtsFi1aVNThmhIsr6Rh1VPhL2TVUyLyLvAPYIeqxni2DQduB7zLzw9T1a98nNsDeAUoBbytqqNCFWdhGT58eObzxo0bZ3bF9brwwgu58MILTzivdevW/Pjjj6EOzxi/eJNGTqsMW0N4+AtlSWMC0MPH9pdUNc7z4ythlALeAHoC0UBfEYkOYZzGGD8lJ0OFChAZ6Xt/lSqQkgJHjxZuXKbwhCxpqOo8YE8+Tj0H+ENVN6jqEeC/wBVBDc4Yky/JyTlXTcGxEoi1a4SvPJOGiJwlIrNEZKXndayIPFaAe94jIstF5F0ROdXH/rrA5iyvEz3bcopvoIgsFpHFO3fu9HmMetaONqYg7O8o76RRubJ7tKQRvvwpabwFDAWOAqjqcqBPPu/3L6AREAdsBcbm8zqZVHW8qsararx3XqesypUrx+7du+0/vCkQVWX37t2UK1euqEMpUv6WNKxdI3z50xBeQVV/FpGs29LyczNVzZwDQ0TeAr70cdgWoF6W11GebfkSFRVFYmIiOZVCjPFXuXLlSvzgyuTknMdogCWNksCfpLFLRBoBCiAivXGlhICJyOmq6j23F7DSx2G/AE1EpCEuWfQBrs/P/QAiIyNp2LBhfk83xmSRnAz16+e835JG+PMnadwNjAeaicgW4E+gX14nicgkoCtQXUQSgSeBriISh0tAG4E7PMfWwXWtvURV00TkHmAGrsvtu6pqI9qMKQasTcP4kzRUVbuLSEUgQlX3e0oBeZ3U18fmd3I49m/gkiyvvwJO6I5rjCla1qZh/GkInwqgqgdV1fv9YUroQjLGFEdHj7oxGJY0SrYcSxoi0gxoAZwiIldl2VUFKNldSIwpgfKaQgSgUiX3aEkjfOVWPdUUNw1IVeCyLNv346YCMcaUIP4kjVKloGJFa9MIZzkmDVX9HPhcRNqr6sJCjMkYUwz5kzTA5p8Kd/40hP8qInfjqqoyq6VU9ZaQRWWMKXYsaRjwryH8Q6A2cDEwFzfYzgqfxpQwljQM+Jc0Gqvq48BBVX0fuBRoF9qwjDHFjb9Jo3Jla9MIZ/4kDe8kx3tFJAY4BagZupCMMcWRlTQM+NemMd4zG+3jwBdAJeCJkEZljCl2LGkY8CNpqOrbnqdzgTNDG44xprhKToby5XNegMnLkkZ4yzNpiEhVoD/QIOvxqnpf6MIyxhQ3eU0h4uVt01CF4yfHNuHAn+qpr4CfgBVARmjDMcYUV/4mjSpVIC0NUlNdycSEF3+SRjlVHRzySIwxxVogSQNcFZUljfDj1zgNEbldRE4XkdO8PyGPzBhTrOQnaZjw409J4wjwAvAonoWYPI/WKG5MCZKcDPXq5X2cd00NSxrhyZ+k8RBugN+uUAdjjCm+Ai1p2AC/8ORP9dQfwKFQB2KMKd6sesqAfyWNg0CCiMwGDns3WpdbY0qOo0fh0CH/u9yCJY1w5U/S+J/nJyAi8i5uPY4dqhrj2fYCbm2OI8B64GZV3evj3I24SRHTgTRVjQ/0/saY4PEmACtpGH9GhL+fz2tPAF4HPsiy7VtgqKqmichoYCgwJIfzz7d2FGOKB3+nEIFjJY0DB0IXjyk6uS33+rGqXisiKzjWayqTqsbmdmFVnSciDbJtm5nl5U9A74CiNcYUiUCSRoUKbiS4NYSHp9xKGvd7Hv8RonvfAkzOYZ8CM0VEgX+r6vgQxWCM8UMgSUPErRVuSSM85dh7SlW3ep7epap/Zf0B7irITUXkUSANmJjDIeepamugJ3C3iHTO5VoDRWSxiCzeuXNnQcIyxuQgkKQBrorKqqfCkz9dbi/0sa1nfm8oIgNwpZd+qnpCtReAqm7xPO4APgPOyel6qjpeVeNVNb5GjRr5DcsYk4tAk4aVNMJXjklDRO70tGc0E5HlWX7+BJbn52Yi0gN4BLhcVX2O/RCRiiJS2fscuAhYmZ/7GWOCw0oaxiu3No2PgK+B54D/y7J9v6ruyevCIjIJ6ApUF5FE4Elcb6mywLfi5kz+SVUHiUgd4G1VvQSoBXzm2V8a+EhVvwn0jRljgic/ScNKGuEpx6ShqslAsog8BmxT1cMi0hWIFZEPfI2vyHZ+Xx+b38nh2L+BSzzPNwCt/IzfGFMIkpOhXDkoU8a/4ytVgsTE0MZkioY/bRpTgXQRaQyMB+rhSiHGmBLC3ylEvKx6Knz5kzQyVDUNuAp4TVX/CZwe2rCMMcVJoEnDGsLDlz9J46iI9MUt+fqlZ1seqwQbY8JJfkoaljTCkz9J42agPfCMqv4pIg2BD0MbljGmOMlP0jh0CNLTQxeTKRq5TSNSRVX3qepqIHNGW0/imFQo0RljioXkZKhb1//jK1VyjwcPHpvA0ISH3Eoac7xPRGRWtn0Bz3prjDl55aekAVZFFY5ySxqS5Xn2NcEFY0yJkd+kYT2owk9uSUNzeO7rtTEmTKWluWqmQHtPgZU0wlFuI8JrishgXKnC+xzPa5vkyZgSIpAFmLyseip85ZY03gIq+3gO8HbIIjLGFCuBTiECVj0VznKbRuSpwgzEGFM85SdpWPVU+PJnnIYxpgQrSEnDkkb4saRhjMmVVU+ZrPJMGiJSqjACMcYUT/lJGhUrukcraYQff0oa60TkBRE1aFf9AAAgAElEQVSJDnk0xphiJz9Jo1QpqFDBkkY48idptALWAm+LyE+eNbltYgBjSoj8JA2w6dHDVZ5JQ1X3q+pbqtoBGIJbgW+riLzvWWPDGBPGkpOhbFn3EwibHj08+dWmISKXi8hnwMvAWOBMYBrwVYjjM8YUsUCnEPGykkZ4ym1wn9c6YDbwgqouyLJ9ioh0Dk1YxpjiIr9Jw0oa4cmfNo3+qnpr1oQhIh0BVPW+nE8DEXlXRHaIyMos204TkW9FZJ3n8dQczr3Jc8w6EbnJz/djjAmygpQ0LGmEH3+Sxqs+tr3m5/UnAD2ybfs/YJaqNgFmeV4fR0ROw7WdtAPOAZ7MKbkYY0LLqqdMVrktwtQe6ADUyDJZIUAVwK+xG6o6T0QaZNt8BdDV8/x93LodQ7IdczHwraru8cTyLS752OJPxhSy5GQ4/fTAz7PqqfCUW5tGGaCS55iskxXuA3oX4J61VHWr5/k2oJaPY+oCm7O8TvRsO4GIDAQGApxxxhkFCMsY48u+fVY9ZY7JbcLCucBcEZmgqn+F4uaqqiJSoLU5VHU8MB4gPj7e1vkwJsiSkqBq1cDP81ZPqYLYsm1hI7fqqZdV9QHgdV8f7Kp6eT7vuV1ETlfVrSJyOrDDxzFbOFaFBRBFluVnjTGF4+hRtwDTqfloUaxUySWMQ4eOTStiTn65VU996HkcE+R7fgHcBIzyPH7u45gZwLNZGr8vAoYGOQ5jTi7p6bBokfsk7tChUL6+JyW5x/wkjawz3VrSCB859p5S1SWex7neH2A5kOR5nicRmQQsBJqKSKKI3IpLFheKyDqgu+c1IhIvIm977rkHeBr4xfMzwtsobkyJNH061KkDHTvCeefB2We7BBJiwUga1oMqvOQ5uE9E5gCXe45dAuwQkfmqOjjXEwFV7ZvDrgt8HLsYuC3L63eBd/O6hzFh74svoHdvaNECXnnF1ReNGAGXXw5Ll0Jdn31EgsKbNE47LfBzbSGm8OTPiPBTVHWfiNwGfKCqT4rI8lAHZowB1qyBa65xJYsZM461SLdvD+ec4/bNnQuRkSG5fbCqp0z48GdwX2lPg/W1wJchjscY45WRAQMHugaBL744vgtTdDS89RYsXAgTJ4YshIIkjSqeubD37QtePKbo+ZM0RuAapv9Q1V9E5EzcfFTGmFB6+2348UcYOxZq+RjO1KcPtGwJY8a4xvEQ2ONpSbSkYbz8mRr9E1WNVdW7PK83qOrVoQ/NmBLs8GEYPtw1eg8Y4PsYEXj4YVi1Cr7+OiRhWEnDZOdPQ3gN4HagQdbjVfWW0IVlTAn3/vuwdSt8+GHuXWv79IFhw1xp5JJLgh5GUpKrHctPk4l3FLkljfDiT0P458APwHdAemjDMcaQlgajR7uG7m7dcj+2TBm4/XZ46inYtg1q1w5qKElJ+StlgFvuNSLCkka48adNo4KqDlHVj1V1qvcn5JEZU1J9/DFs2ABDh/o3gO+qq1ybxue+xskWTEGShoirorKkEV78SRpfikjwy73GmBNlZMBzz7neUZf7OVNPTAw0bgyffhr0cAqSNMCSRjjyJ2ncj0scqSKyT0T2i4j9GRgTCtOnw8qVrpQR4c9/T9xX+l694PvvYe/eoIaTlJS/gX1eljTCjz+9pyqraoSqllPVKp7XVQojOGNKFFV49llo0MA1cAeiVy/XFjJ9elBDspKGyS7PpCHODSLyuOd1PRE5J/ShGVPCzJ0LP/0EjzwCpf3po5JFu3ZQs2bQu97u2WNJwxzPn/Lvm0B74HrP6wPAGyGLyJiS6tln3SC+m28O/NyICOjSxSWeIA30O3LETWtuScNk5U/SaKeqdwOpAKqahFvVzxgTLIsXw7ffwuDBUK5c/q7RpQskJsLGjUEJqSAD+7wsaYQff5LGUREpBShkDvbLCGlUxpQ0I0e6uaUGDcr/Nbp0cY9z/Vq5IE/BShrJyUEJxxQT/iSNV4HPgJoi8gzwI/BsSKMypiRZuNCNsRg8+NjcG/kRHe26OhWzpHHokGujN+Ehz9Y2VZ0oIktwa2AIcKWqrgl5ZMaUBKqu4btWLZc0CiIiAjp3hnnzghJasJIGuOnRC3IdU3z403uqJdASt5b3HEsYxgTR55+7mWyHDw/OmqidO7vR5ImJBb5UQRZg8rL5p8JPjklDRE7xrNr3P1zPqX7A5yIyW0RsnIYxBZWcDPfc40Z033prcK553nnuceHCAl8qmCUNSxrhI7eSxtPAYqCJqvZS1SuBs3Brdj9TGMEZE9YeecTNZPvuu8FbeS821l3rl18KfCnvWhpZ134KlCWN8JNb0ugO/J+qZvaUUtV0YJhnX76ISFMRScjys09EHsh2TFcRSc5yzBP5vZ8xxdL06TB+vGvHaNs2eNctWxZatXJdeAsoKcmt812QfGZJI/zk1hB+RFVP6POgqmkicji/N1TV34E4AE9X3i243lnZ/aCq/8jvfYwpthIT4aabIC4Onn46+Ndv29YtAZuR4f/8VT4UdAoRsKQRjnJLGuVE5Gxcj6msBCgbpPtfAKxX1b+CdD1jire0NOjXD1JTYfLk/A/ky03btvCvf8HatdCsWb4vY0nD+JJb0tgKvJjDvm1Bun8fYFIO+9qLyDLgb+BhVV3l6yARGQgMBDjjjDOCFJYxIfL0065L7AcfwFlnheYe8fHucfHiAiWNPXsK1nMKLGmEoxyThqqeH8obi0gZ4HJgqI/dS4H6qnrAs5bH/4Amvq6jquOB8QDx8fHBmXTHmFCYN88ljQED4MYbQ3ef5s3dsnm//AI33JDvy+zY4WrQCqJSJfdoSSN85L/Cs+B6AktVdXv2Haq6T1UPeJ5/BUSKSPXCDtCYoElLg7vuctOev/56aO9VujS0bl3gHlTbt7sxhwUREQGVK1vSCCdFmTT6kkPVlIjUFnHrXHqmYY8AdhdibMYE17//DatWwdixwRnEl5c2bSAhAdLT83X64cNuGEnNmgUPxeafCi9FkjREpCJwIfBplm2DRMQ7W1tvYKWnTeNVoI9qkOZ7NqawJSfDE09At25w5ZWFc89WrSAlBdavz9fpO3a4x4KWNMBmug03ObZpiEjr3E5U1aX5vamqHgSqZds2Lsvz14EQl+GNKSTvvedalUePdkuzFoZWrdzj8uX5anD3Jo1glTQsaYSP3HpPjc1lnwLdghyLMeEnIwPeeAPatz/Wq6kwREe7BoXly6F374BP3+5pabSkYbIrst5TxpQI33wDf/wBI0YU7n3LlYOmTWHZsnydHszqqVNOCcr8iaaY8GeW2woi8piIjPe8biIiNlLbGH+88QbUrg1XX134946NdSWNfLDqKZMTfxrC3wOOAB08r7cAI0MWkTHhIikJZs6E/v2hTBGskNyqlVv6NR9dl7Zvd0M9vOMsCsKSRnjxJ2k0UtXngaMAqnqIE6cWMcZkN22aG59RFKUMcCUNgBUrAj51x47glDLAJY39+13zjjn5+ZM0johIeY6tEd4IyPeEhcaUGJ9+ClFRhdsAnpU3aeSjiioYA/u8sq7eZ05+/iSNJ4FvgHoiMhGYBTwS0qiMOdkdOAAzZsBVVxVoptkCiYpyMw7mI2kEs6ThnfTQu6iTObn5s0b4tyKyFDgXVy11v6ruCnlkxpzMvvnGzWR71VVFF4NIvhvDt28P3jIf3kkP9+xxs6iYk5u/X4HKAUnAPiBaRDqHLiRjwsBXX7mv2N7lV4uKN2kE0KCQkQE7dwavpJE1aZiTX54lDREZDVwHrAK8f3kKzAthXMacvFRh1iw4/3woVapoY4mNhYMH4c8/oVEjv07Zs8dNWRWsNg1LGuElz6QBXAk0VVVr/DbGHxs2wKZNMGRIUUdy/HQifiaNYI7RAKjmmTDIkkZ48Kd6agMQpFXvjSkBZs1yj92KwUw7LVq4to0A2jWCOYUIHGsI323zVIcFf0oah4AEEZlFlq62qnpfyKIy5mQ2axbUqeOm8ShqFSpAkyYBTScSzClEwM1oUqGClTTChT9J4wvPjzEmLxkZMHs29OhReDPa5iU2Fn791e/D81M9dTjtMAsTFzJ341wS9yWyK2UXqkrdynWpX7U+Vao+wJ49kdi44JNfrklDREoBF6lqv0KKx5iT2+rVrutRcaia8mrVCqZMcWNH/JgXZPt2N7SkWrU8D2XF9hW8+NOLTF45mZS0FCIkgpoVa1KjQg0UZd5f80hKTYL0i/jvLztoMGchg9sPpnLZykF4Y6Yo5Jo0VDVdROqLSBlVPVJYQRlz0po/3z126lS0cWSVdTqR9u3zPPzvv13VVG5jEtfuXsuDMx7kq3VfUSGyAjfG3silZ11Kl/pdOKXcKccdu+PgDrp8mc72fZEMnzucN355gye7PMldbe9CiktpzPjNn+qpDcB8EfkCOOjdqKovhiwqY05WCxa4ep0zzyzqSI7JOp2IH0njzz+hYUPf+1KOpvD0vKcZs2AM5SPL80y3ZxgUP4jTyp+W4/VqVqxJi/ogq09nxm0/M+S7Idzz9T18v/F7JlwxwUodJxl/ksZ6z08EYP+6xuRmwQLo0KH4tGcA1K8PlSv73YPqzz99j0lc8vcSbvzsRtbsWkP/Vv0Z3X00tSvV9uua1aq5hvC2ddsyq/8sXvrpJf757T85d9e5fNn3SxqemkOWMsWOP9OIPAVuXQ3PDLfGGF927HALLg0cWNSRHC+A6USOHoXNm48vKKVlpDH6x9EMnzucmhVrMvOGmVzY6MKAQjjtNJc0VEFEGNx+MHG14+j9cW+6f9idH2/+kdMrnx7oOzNFwJ9FmNqLyGrgN8/rViLyZkFvLCIbRWSFiCSIyGIf+0VEXhWRP0RkeV5rlhtT5BYudI8dOuR+XFHwJg3VXA/btMl1APNWT/2x5w86v9eZx2Y/Ru/o3qy4c0XACQNc0jh61LXFe3Vr2I2v+33N9gPbufg/F5OUYjMangz8Gdz3MnAxsBtAVZcBwZp76nxVjVNVX3NH9wSaeH4GAv8K0j2NCY0FCyAyEtq0KepIThQb61ZC2rQp18P+/NM9nlE/nbELxtJqXCvW7FrDR1d9xKSrJ+XadpGbnEaFt4tqx//6/I/fd//OlZOvJC0jLV/XN4XHrwkLVXVztk3pIYgluyuAD9T5CagqIlZ+NcXXggUuYZQrV9SRnMjPtTU2bHCP9y+6jIe/fZhuDbux4s4V9G3Zt0C3z23+qe5nduedy99h3l/zGD5neIHuY0LPn6SxWUQ6ACoikSLyMLAmCPdWYKaILBERX5XAdYGsySrRs+04IjJQRBaLyOKdO3cGISxj8uHIEfjlF+jYsagj8S0mxj3mkDRUlUWJixj15WSIOMLu0sv55JpP+KLPF0RViSrw7fOatPCG2Bu49exbefaHZ/luw3cFvp8JHX+SxiDgbtwH9hYgDrgrCPc+T1Vb46qh7s7vdOuqOl5V41U1vkaNGkEIy5h8+PVXOHy4eLZngFs+r2HDE5LGvsP7GLd4HG3Gt+Hcd85ly19lqHb6QdY/sJbe0b2DNo7Cn5luX+35Ks1rNOeGT29gx8EdQbmvCb4ck4aI1ANQ1V2q2k9Va6lqTVW9ASjw/wxV3eJ53AF8BpyT7ZAtQL0sr6M824wpfhYscI9+jIMoMll6UC3+ezG3fXEbp489nTun34mivHnJm8SUuZzW0adSIbJCUG/tbdPIbdLCCpEVmNx7MntS9jB4xuCg3t8ET24ljW9FpEH2jSJyM/BKQW4qIhVFpLL3OXARsDLbYV8A/T29qM4FklV1a0Hua0zILFjgvsmfXnyb3TS2JVNK/U6Ht86l7Vtt+e/K/9KvZT9+uf0Xlg5cyp1t72TTxlI5DuwrCO9Mt3lNWhhTM4ah5w1l4oqJzPhjRvADMQWW2ziNwbg2h0tVdR2AiAwFrge6FPC+tYDPPEXf0sBHqvqNiAwCUNVxwFfAJcAfuJl2by7gPY0JDVWXNM4/v6gjydGybcu4p/IUfrxGabz3b17p8QoD4gZQpWyVzGP274ddu3IeDV4Qgcx0O7TTUCavmsyd0+9k5V0rg17qMQWTY9JQ1a9E5DDwtYhcCdyGq0LqrKoF6lCtqhuAVj62j8vyXHFtKcYUb5s2uQmbimF7hqoyev5oHv3+UU4rcwpvfw433/MkEe1uPeFYb3fbUM2A4h0Vnpdypcvx73/8m67vd+XpuU/zXPfnQhOQyZdcG8JVdRbuG/4c4EygW0EThjFhx9ueUcySxpH0I9z2xW0MnTWUa6Kv4fd7fufW38oTsSJ7TbDjTRqhKGmAawz3dyGmLg260L9Vf1786UX+TPozNAGZfMmtIXy/iOzDVRNVAS4AdmTZbowBN7NtpUrHurUWA0fTj3LFf6/g3YR3ebzz425gXqUaLsYcut2uXu0eGzcOTUzeqUT89Wy3ZykdUZpHvnskNAGZfMkxaahqZVWt4nkso6oVs7yuktN5xpQ4CxZAu3ZQ2p/5P0NPVblz+p1888c3jLt0HCPOH3Gs62xsrFvFz8d0IgsWuMUGvY3WwRZo0qhbpS5DOg5hyuop/PDXD6EJygTMrxHhxpgcHDjgPoSLUdXU6PmjeefXd3is02PcEX/H8TtjY10d0bZtx23OyHBJI5RjE6tXdw3tgXi4w8NEVYniwRkPkqEZoQnMBMSShjEF8fPP7hO3mCSNGX/MYOisofSN6cuI80eceEAO04n8/rsrBYTybZx+ulvU8OhR/8+pEFmBZ7s9y5KtS/h41cehC874zZKGMQXhbQQ/99yijQNISknili9uIbpGNO9c/o7v0dwtW7rHbEnD+zZCWdKoW9fVim0NcLRVv9h+tKrVike/f5Qj6baAaFGzpGFMQSxYAC1aQNWqRR0J9359LzsO7uCDKz+gfGR53wdVq+Y+vRMSjts8f77b1bRp6OKrU8c9bglwXocIiWBU91FsSNrAvxf/O/iBmYBY0jAmvzIy3BoaxaBq6tM1nzJxxUQe6/QYberkMTV769Zurqws5s8P/YKDdT3TjQaaNAAubnQx5zc4nxHzRrDvsHXeLEqWNIzJrzVrYO/eIk8aB44c4L6v7yOudhzDOg3L+4TWreG33+DgQcA1Tq9dG/oJeguSNESE0d1Hs+vQLsYuGBvcwExALGkYk1+zZ7vHLgWdVadgnvvhObbs38Ibl7xBZKnIvE9o3do1LixbBsCECW5zt26hixFc9VfZsvlLGuDWF78m+hrGLhzLtgPb8j7BhIQlDWPy6/vvoUGD0A2h9sP6PesZs3AMN8beSId6fpZ4vCsLLlnCnj3wzDPQsye0bRu6OMFVfdWp42Zcya9nuj1DaloqT899OniBmYBY0jAmP9LTYc4cuOCCIg3jwRkPUqZUGUZ3H53jMSkp8OKL0LcvPPUUzP69Dhk1asHSpYwc6VaBff75wom3Tp38lzQAmlRrwsA2Axm/dDzrdq8LXmDGb8VjCKsxJ5tlyyApKfR1OrmYtWEW09ZOY3T30Zxe2feU7IsWwdVXuw/qunVh8mRQFRqUT6DUR0dYfwRuvbXwZkCpW/eENviAPdHlCd5f9j6PzX6Myb0nBycw4zcraRiTH7Nmuccimg49QzP457f/pP4p9bmv3X0+j9m50yWMyEjX/JKY6EoVH30ELeruJfbIEsaOOsorBVodJzB167oE5mMWE7/VrlSbh9o/xMerPuaXLb8ELzjjF0saxuTH999D8+ZFtujSpBWT+HXbrzzT7RnKlS53wv6MDLjhBtczaupU6NrVba9UyVVTfTl6NZ9yFYO7LqVixcKLu25dOHQIkpMLdp2HOzxM9QrVGfLdELQgGcgEzJKGMYE6dAjmzoXu3Yvk9qlpqQz7fhitT29N35Z9fR7zwQcwcya8/LLrLHWC+Hj3uHhx6AL1oSDdbrOqUrYKj3d+nNkbZzNz/cyCB2b8ZknDmEDNmOFal6+8skhu//rPr7MpeRMvXPgCEXLif+HUVHjiCdcb6o47fFwAoF49V0r66afQBpuNd1R4QXpQed3R5g4aVm3IkO+G2GSGhciShinRVN2X7XHj4Ii/0xr9739u/vBOnUIamy97UvbwzA/P0LNxT7o19N0I/8YbsHkzjBqVywhvETdfViEnjWCVNADKli7LyG4jWbZ9GZNWTCr4BY1fLGmYEmvbNjdkoW1buPNO12h8+HAeJx09CtOmwWWXuRbmQvbMvGfYd3hfjl1s9+2DZ5+Fiy/2o2PXuefCH3/4v5xeEOR3/qmc9InpQ1ztOB6b/RiH0/L6xzPBUOhJQ0TqichsEVktIqtE5H4fx3QVkWQRSfD8PFHYcZrwN2QIrFwJb74JL70EX37pEkeu7arz5rmutkVQNfVn0p+8/svrDGg1gJa1Wvo85rXX3BTnI0f6ccF27dzjokXBCzIP5cu7xZiClTQiJILR3Uezce9Gxi0eF5yLmtypaqH+AKcDrT3PKwNrgehsx3QFvgz02m3atFFj/PHjj6qgOnTosW1jx7pt06blcuLAgarly6sePBjyGLPrM6WPlh9ZXhOTE33u37dP9bTTVC+91M8LHjigGhGh+vjjwQvSDzExqpdfHrzrZWRk6AXvX6DVRlfTvSl7g3fhEgJYrAF8zhZ6SUNVt6rqUs/z/cAaoG5hx2FKLlW4916IioJHHz22/d573YwgTz2VQ2kjORkmToQ+faBChUKLF2DeX/P478r/8nCHh6lbxfd/lzfecKWMJ/wtl1es6BZlKsSSBkCjRrAuiIO5RYRR3UexO2U3YxaMCd6FjU9F2qYhIg2AswFff7XtRWSZiHwtIi1yucZAEVksIot37twZokhNOJk3z41KHjGC48YoREbCsGGuYfzrr32c+OGHbmbYu+4qtFgB0jLSuPfreznjlDP4v/P+z+cxe/fCmDHQowecc04AFz/3XJc00tODE6wfWrRws+rm2X4UgPg68VzX4jpe/OlFtu4PcJUnE5AiSxoiUgmYCjygqtknyF8K1FfVVsBrwP9yuo6qjlfVeFWNr1GjRugCNmHjnXegShW47roT9/XvD/Xru0n8jqPqGj/atj02xqGQjFs8juXbl/PSxS9RIdJ3CWf0aNeefULceTnvPFeCyraSXyjFxLgctXZtcK87sttIjqQfYcRcH8vcmqApkqQhIpG4hDFRVT/Nvl9V96nqAc/zr4BIEaleyGGaMLR3L3zyCVx/ve8apjJl4L773IJ8K1dm2TFtmls/o5BLGdsPbOfx2Y9z4ZkX0qtZL5/HbN7sBvHdcEMOA/ly4x0qPmdOQcIMSAtPvcGqVcG9buPTGnNHmzt4a+lbrN0d5IxkMhVF7ykB3gHWqOqLORxT23McInIOLs7C6xdoTn6JiTB9Okya5ObRWLIEDh9m0iQ3+O3WW3M+9cYbXVXVO+94Nhw6BPffD9HR0K9foYQPrpPK7dNuJ+VoCq/1fM33mt/A//2fKwj51WMqu7p1oUmTQk0aTZtCqVLZknKQPN75ccqVLsewWX4sRmXypShmue0I3AisEBHvQsXDgDMAVHUc0Bu4U0TSgBSgj6eV35icbd8Ob7/tVhX6448T91eowLuRvxJbvyZtWlYAyvi8TI0arkftBx+4AXJln3sONm50s/4V4tiM9xLeY9raabx08Us0re578e6pU90EhE884arV8qVrV/j4Y1dnVKpUvuP1V9myLk8Fu6QBUKtSLf7Z4Z8Mnzuc2X/O5vyGRTOhZFgLpKtVcf+xLrclVGqq6ujRqpUquT6z3bqpvvKK6oIFqr/9prp0qeonn+jvN4xQUB3DYNWaNVWHDFHdsMHnJWfOdJf6793zVEVUb7ihUN/Shj0btNKzlfT8Cedreka6z2P+/tt1sW3bVvXIkQLcbOJE92aXLCnARQJz9dWqTZqE5tqHjhzSM185U5u93kwPpx0OzU3CCAF2uS3yD/pg/ljSKIGmT3efPuA6/69Zk+OhTz2lKpKhmz/4XvWKK9wYhYgI9/zTT1VTUjKPTU9K1vpVk7Q7M10SKsRxGQcOH9DW/26tVZ6rohuTNvo8JilJtU0bN2Tkt98KeMMtW9zvb+zYAl7If0884XLxoUOhuf5Xa79ShqPPzHsmNDcII5Y0TMmwdq0bxQaqTZuqfv11rodnZKg2a6bauXOWjZs3qw4b5kodoBoZqRobqxoXp1qpko7gMQXVDatC9MnmQ1p6ml720WUa8VSEfvn7lz6P2b1b9ZxzXLhf+j4kcE2bqvboEaSL5e3jj92vfOnS0N3j6slXa7mR5XT9nvWhu0kYCDRp2NxT5uSyb59r+W3Rwg24GDPGdRft0SPX05Yvh99+c2tJZIqKcn1Ut2xxM9cOHuy2RUVBv34MmNabiAh456PyoX1PHqrKgzMeZNraabza41UuPevSE46ZOdONx/v1V5gyBS498ZD8ueQSt0bIgQNBumDuQtWDKquXe7xMZEQkt0+73WbBDaZAMkxx/7GSRhjbs0f1ySdVq1Z1X1EHDFDdutXv04cMUS1VSnXHjsBue8klqnXqqB49Gth5gUo9mqo3fnqjMhwd/M3gE/YvXKjas6d769HRqosXBzmA2bPdxadODfKFfTtyxJWUHn44tPcZv3i8Mhx98+c3Q3ujkxgBljTEnRMe4uPjdXEhLypjQmzXLjeb4Guvwf79rlvT448HNCAhPR3OOAPi4lwv3BP2Z6SzYscKlm5dyrrd69iwdwO7D+0mKTWJ3Us78de/Xyb6vv8jptOfRFWOonmN5rQ+vTUxNWMoU8p3D6xAbN2/leumXMcPm37g6fOf5tFOj2Z2r124EJ58Er79FqpVg4ceggcfhHInLtZXMEePQq1acPnlrvdZIejc2RVsli4N3T1UlYv/czELNi9gxZ0raHhqw9Dd7CQlIktU1e8Rq0XR5daYvG3bBmPHwr/+5cZJXHONmygqNjbgS82c6Rb9efVV9zo1LZWfEn/ix00/8uOmH1mweQH7j+wHoHREaRpWbZxxd0EAABVBSURBVEjNijU5vdLpRHVKZNvHu9g66yqONr6Bab9PIyUtBYDIiEha1mpJ69qtaRfVjvZR7Wleo7nPhZF8STmawosLX+S5H58jLSONj676KHMlvo0b4Z//dFVQNWvCCy/AoEFuudaQiIx0VVRfflloXW979nTTtmzbBrVrh+YeIsLbl79NzJsx3Pz5zczqP4tSEaF/b+HMShqmeDl4EJ5/3n1KHj7sGiEefdStx51P116rfPd9OsM+e5NZm75m7sa5pKSlIAgxNWPoWK8j551xHufUPYeGpzakdMTx36Wee859uK1YAdEtMtiQtIGlW5eydOtSlmxdwpK/l5CUmgTAKWVPoV1UO9rVbUfz6s05q9pZ1KhYg4qRFUnXdDYnb2bdnnV8ufZLvlz7JcmHk+nVrBfPX/g8jU9rjKr7on/ffW6d7yFDXOmiUNbxnjzZTcY4b16hLDCVkABnnw3vvQcDBoT2Xu8nvM+AzwcwvMtwnuz6ZGhvdpIJtKRhScMUH7Nnu0+PTZvcxFBPP+1GgeVDUkoS3234js9/nc/EAc9D/JvQ80GaVmvKRY0u4qJGF9GxXkdOLX9qntfas8e1jfftm2WUeBaqytrda1mYuJCFmxeyMHEhK3esRMn5/1a18tW4rOll3BJ3C53quw/o1FS4/Xb4z3/ceLv333fVaoVm3z63BGy/fjB+fMhvp+oGpHfq5PJVqN30v5v4cNmHzOo/ywb9ZWFJw5LGySc11X2Vf+kllyTefddNpBeAtIw0ft7yMzP+mMHMDTP5ecvPZGgG5X4ZSur0Z3nsoync2jOeBlUb5CvEu+5yYf31l6v6z0vK0RTWJ61n3e51JKUmceDIASIkgqgqUZxxyhnE1oo9rkSzbRv06uVWXx0xwhWuIoqib+OAAfDpp7B1a6EUb265BT77DHbuhNIhriw/cOQA8ePj2Xd4H0sGLuH0yqeH9oYnCUsaljROLr/+6mbaW72axP7D+Lb9E6zZUJbkZKha1TVhXHihq9fPLjUtlZnrZzJl9RSmrZ3G3tS9REgEbeu05eJGF9M1qif9urSjSRNh7tyChblunashu+uuY20jwZKQ4Nqfd+92s69fdVVwrx+QH35wLdSFUWeEmzzy2mvhxx+hY8eQ347l25fT/p32tKjRgjkD5uQ4a3BJEmjSKPJussH8sS63J5HUVNURIzSlVEX9T9W79LwWe9RVWKiWLevG20VGauaYu2uuUU1IcFNETF09VftO6auVnq2kDEdPHXWq3vTZTTp55WTdfWh35i1eesmdP2dOcEK+807V0qWDMAI7i88+U61QQTUqKrQD3fyWkaF61lmq551XKLdLSnKj2m++uVBup6qqn//2ucpw0asnX53jFC0lCTYi3BRrR45oxsSP9JeoK3UwY7RamWQF1caNVZ99VnXZMve55TlUly5Vvfu+VK1Y+bAi6Vq6zQTln9W12uhqetvnt+k3677RI2knTrx08KBqrVpuBpBg+f/2zj04qjrL45/TCR1CugkhgRiQBRKQDCKIIogzAypag49RXFnUcX2M4yDOIuOgzqyrWzXjaI2IouPorjvq6mo5iyVaiI9BhhEKFogiLzFICG8IkHQTEkiEPLrP/vG7SBtC+qah00nn96m61ffxu/eeU7e7v/f3Oqe8XNXvPzOpSsNh4y+ojhlj4ki1G2bPNoYVFbXJ7aZPNy8Ge/a0ye1UVXXOyjnKb9EZH8/Q8PEvXCeltaJhm6c6O6GQ6XgOBk2yiUOHzJj91FQzDNPrNRmLMjNPfGZkmPPq603ZhgazXl8PdXVoXT37DngoP5hCoDKVYHkjgR21BDYf5OsSD5/VjWAffUlNCXP9JA/TpsHll3+3DX/fkX0sKFnA/M3z+XTHpzTUZpC+6gnqVk4lwxfmmdkp3P2zFJqLFq4KU6eagLfLl7e6e6RFZs0yE9Lnzm0+iZMbKitNy88HH5zoXE9vm0nn7qipgfx8GDHCTBCJMzt3wqBBcP/9ZoJ/W6CqPLDoAZ4tepYZo2fw3MTnThl6PtmxfRpWNJqnshJKSqCkBN1cwv4NFVRsqaJq92G0sREv9eSxn76U4aWehhRo8EBDCjR6ICUM3hB0CUOXEET+vA7Rg88ZzWeM4TPG8DmjCXJyFsUUGhmYEWDUhWGuvC2PSf/oId1/lP01+9lauZXNwc2sO7COlXtWfptEpyCrgEmFk5hUOImxZ4+lZHMK99xj2sDHjTPTOIYO/e595swxw1QfeSTGHBMt0NAAl10GGzbA6tVQWNi68xcuNCOkysvNNJTp02lW+BLOs8+asCpLlpxI1BRHbr0VFiwwApKdHffbAUY4Hlz0IHOK5nDvqHt5/qrnTxpu3RmwotGZRKOuzsxaKyszn4GAWYJBCAYJBSrYV3WIoorurDmaz1fhEZQ2jmDP0eEcbWzhl9ktAJm7IHsLZJdCz1LICID3CIS8cCwLT+VgPOXDCZeNJnzwHHOehOnSq5T0futIP/sr0jIr6JJRhddfjbdnLSHfEer1GHWNddSF6qipr6Gm/ruxjnK65XBJv0v4fr/vc83gaxjaa+hJb4DhsOmnfeghUzm69lrTkez3w1tvmTf4yZPNMM54jEAqKzPzC7KzYfFiM2w0Gvv3mzkXb75pRO6NN+DCC8+8bWeMo0fNSLacHJNDPC0trrcrLjYz9idPNnmz2gpV5eG/P8ysFbOYMHACcyfPJadb50oSakUjGUVD1TQhrV8P69dzbMMadm9by66aMnZnYpY0P7vDBexvzCdYdw7VVcOpDwyH4BBQ5+2pSy303ghnbcDbp4TuOTV07650T88gQ7Lx1JxNuLoPxw5lU30gi8q92Rw6kIlq86/C3bIOkztkF70Gb6fnoFK6528htVstoXCIkIYIhUOENUxIzac3xUtaShppqWmkpaSR0SWD3hm9yfXlUpBVwJCcIeRm5LpuJggE4MUXTeruQMDsy86G++6DX/86vk0+y5aZYIE9epjQJKeaqF5ZaUZbPfOMab176CETBSXO/8Fnhg8/hB//2NQ4nnkm7rd74gl49FEj/D/5Sdxv9x1eW/ca9350L7m+XN684U3G9R/XtgYkECsaHVw0GgJV7Pm0lC2rdrC6+ADFB+rZU6tUeNKp8vg4go+6sA/q/XAsE+oyzWd99+9cx9erkryCAAMLD/O9YfWMPF84rzCDs7r3IqdbjuuYSXV1psng4EET+iktzbzR5+dDVvR5cW1CQ4OpaAWDprmoTWZPYzT86quhogLuuANuv93cv6bGpBN/7z0zpLSmxgyjnTXLtN13KO67D154wVTbpkyJ660aG2H8eJMGdulSU5trS1aXrWbKvCnsrNrJ7SNuZ/aVs+md0cxY7yTDikZHEQ1Vwjt2sfbdr5m3OMiK7RlsOTiAYPVQwuEm0ehS6vB4a/F2rSPdFyLDB927p5CTlUrvnl05K7sb/c5OIT8fCgrMH3pmZmLc6myUl5swIy+9ZAQ2Er/fiMXMmTGFzGofHD1qJsqsWmXa1OKcI333bjNDvKbGRGofMSKutzuJ2vpaHl/2OE+vehpvipe7zr+LmWNnJnWgQysa7VU0VKn88nP++vYi5i9XVu/6HnsD4wkdc95k0oN4zlpPr7ytDBxwmMLCHoy84B+45Nz+nHd2PmmpHaE9o/MSDMKaNWasgc8HAwfC2LFxiEabCGpqTDPV0qVmduNTT8W1Ordtm+l7P3zYDHRo66YqgJJgCU+ueJK3vnyLxnAj4/qP46Zzb+Lac66lX2a/tjcojnQI0RCRicAfgRTgFVV9ssnxNOAN4ELgIHCTqu6Mdt32IhoNoQZKD27hy/WLWL7ka1auzaR0zwhqy8bDYfOF8/jLyBv0BSMvCnD1pF5cMfpcCnrmu46QarG0KUePmuFozz1npufff79pk8uLTyiOXbuMWKxcCddfb/o6Rrmfs3zGKDtcxitrX2Fu8Vw2BzcDZkTfD/v/kPN6n8ew3sMY1nsYeb68Djtkt92LhoikAFuAK4G9wGrgFlXdFFHmF8BwVZ0mIjcDN6hq1FHxZ0o0VPXbjtzIz/pQPUfqjlD9TRXBqiqClVUc2F/Ott0V7Nj3DWXlDRw42JVgcCDhwFCoGAZ1PQDo2q2CwYNLGT9BuOXWwYwd2at9DrW0WFpi5Ur43e9MvHmAiy6Ciy82qfj69DFLbq6pYqWlmaVLl5jGFTc2mkrNU09BdbVpqpo40aRSGTTIaFd2dtvMcVFVigPFLN6+mCU7l1C0t4iK2opvj2d1zWJQz0Hk+fPo4+tDnj+PPF8eWelZ+L1+/Gl+fF4fPq8Pb4qXVE8qqZ5UUiTl2/XjS1uLT0cQjbHAb1X1R872wwCq+oeIMp84ZVaJSCpwAOilUYyNVTR6z+7NkfojNIYbCYVDJjrpMT/MOtT8CdpyPP6u3kP0z97F8EH1jL+8DxNu6suQQrEiYUkeNm0yPf2LFpksSrW1LZf3eo1wRC5gJtxE6fE+fNhMgHz/fVixwohJJJWViRmUEagNUBwopriimI0VG9lZtZP9NfvZd2QfwW+CMV/XIx484kGQbwVEnJlRp9rOzchl+y+3x3S/jiAak4GJqnq3s30bMEZVp0eU+cops9fZ3uaUOelJiMhUYKqzOQQoibMLbsgBYv/WtB+SwY9k8AGSw49k8AGSw49IH/qr6smzcU9Bh5/+qKp/BuIf/L8ViMgXrVHu9koy+JEMPkBy+JEMPkBy+HE6PiSi17UMiBx+cLazr9kyTvNUJqZD3GKxWCwJJBGisRoYLCIDRcQL3AwsaFJmAXCHsz4Z+DRaf4bFYrFY4k+bN0+paqOITAc+wQy5/W9VLRaRxzAhehcArwJvishWoBIjLB2JdtVcdhokgx/J4AMkhx/J4AMkhx8x+5BUk/ssFovFEl/sTDKLxWKxuMaKhsVisVhcY0XjNBCRiSJSIiJbReRfmzmeJiJvO8c/E5EBbW9ly7jwYZyIrBWRRmeOTbvEhR8zRWSTiHwpIn8Xkf6JsLMlXPgwTUQ2ish6Efk/ERna3HUSTTQ/IsrdKCIqIu1u+KqLZ3GniAScZ7FeRO5OhJ3RcPMsRGSK89soFpG/RL1oa3LD2uXEgunE3wbkA15gAzC0SZlfAC856zcDbyfa7hh8GAAMx8QCm5xom0/Dj8uAbs76vR30WXSPWL8OWJhou2PxwynnB5YBRcCoRNsdw7O4E3gh0baeAT8GA+uALGe7d7Tr2ppG7IwGtqrqdlWtB+YC1zcpcz3wP876PGCCtK+oZlF9UNWdqvolEE6EgS5x48cSVf3G2SzCzA9qT7jx4XDEZgbQHkexuPldAPwemAUca0vjXOLWh/aOGz9+DryoqocAVLWCKFjRiJ2+wJ6I7b3OvmbLqGojUA20UQZkV7jxoSPQWj9+Bvw1rha1Hlc+iMi/OGF1ngJmtJFtrSGqHyJyAdBPVT9qS8Nagdvv041Oc+c8EWmP8dLd+HEOcI6IrBCRIicCeYtY0bB0KkTkn4FRwOxE2xILqvqiqhYAvwEeTbQ9rUVEPMAc4IFE23KafAAMUNXhwN840aLQ0UjFNFFdCtwCvCwiPVo6wYpG7CRDOBQ3PnQEXPkhIlcAjwDXqWpd0+MJprXPYi4wKa4WxUY0P/zAMGCpiOwELgYWtLPO8KjPQlUPRnyHXsHk/mlvuPlO7QUWqGqDqu7ApK0Y3NJFrWjETjKEQ3HjQ0cgqh8iMhL4L4xgRG23TQBufIj8MV8DlLahfW5p0Q9VrVbVHFUdoKoDMP1L16lq4rOnncDNs4jMPnUd8HUb2ucWN7/v+ZhaBiKSg2muajnGeqJ7+DvyAlyNUeZtwCPOvscwPwKArsA7wFbgcyA/0TbH4MNFmLeRWkwtqTjRNsfox2KgHFjvLAsSbXMMPvwRKHbsXwKcm2ibY/GjSdmltLPRUy6fxR+cZ7HBeRaFibY5Rj8E01y4CdgI3BztmjaMiMVisVhcY5unLBaLxeIaKxoWi8VicY0VDYvFYrG4xoqGxWKxWFxjRcNisVgsrrGiYemwiEjIiTD6lYi8IyLdopT/ONps1xbOfeV0o8pGREZdJyKlIvKJiFwScfwxZwLiqc6f1F4j21o6D3bIraXDIiI1qupz1t8C1qjqnFZeQzC/g7gHZBSROzFzEqY725cB/wtcpqpRJ4eJyOvAh6o6L552WiwtYWsalmRhOTAIQETmi8gaJz/A1OMFRGSniOSIyAAnx8AbwFfAbSIyxynzSxHZ7qzni8gKZ32piIwSkRQRed2p3WwUkV85xwtEZKFz3+UiUhjNYFVdgsnVPNW5xuvi5CwRkSflRP6Pp50ayXXAbKd2VSAiPxeR1SKyQUTePV7Tcq7zvIisFJHtEpEHRUR+49i9QUSejNV2S+clNdEGWCynixPX6ypgobPrLlWtFJF0YLWIvKuqTWN+DQbuUNUiETkLmO7s/yFwUET6OuvLmpx3PtBXVYc59z7e3PVnYJqqlorIGOA/gMtdmL8WuKeJP9nADZhZxioiPVS1SkQWEFHTEJEqVX3ZWX8cE733T85l8oAfAIWY0BHzROQqTGjsMar6jYj0PE3bLZ0QKxqWjky6iKx31pcDrzrrM0TkBme9H0YgmorGLlUtAlDVAyLiExG/U/4vwDiMaLzX5LztQL6I/An4CFgkIj7gEuAdOZEuJc2lD83lV6nG5Jl4VUQ+BD48xbnDHLHoAfiATyKOzXea3DaJSK6z7wrgNXXyijjCejq2WzohVjQsHZmjqnp+5A4RuRTz5zjWeZteiokB1pTaJtsrgZ8CJRgBugsYS5MQ3qp6SERGAD8CpgFTgPuBqqa2uGQkTYLdqWqjiIwGJmACXU6n+Tf/14FJqrrB6S+5NOJYZBTflhJ/eYjddksnxPZpWJKNTOCQIxiFmNDbblgOPIhpjlqHSQ9bp6rVkYWcSKAeVX0Xk8/iAjUZ9XaIyD85ZcQRlhYRkfGY/oyXm+z3AZmq+jHwK+D4tY5gQosfxw/sF5EuwK0ufPwb8NOIvo+esdpu6bzYmoYl2VgITBORrzG1hiKX5y3HNE0tU9WQiOwBNjdTri/wmphkQgAPO5+3Av8pIo8CXTD5LjY0c/5NIvIDoBuwA7ixmZFTfuB9EemKqSXMdPbPxSTJmYGpgfw78BkQcD79tICqLhSR84EvRKQe+Bj4t1bYbrHYIbcWi8VicY9tnrJYLBaLa6xoWCwWi8U1VjQsFovF4horGhaLxWJxjRUNi8VisbjGiobFYrFYXGNFw2KxWCyu+X+uvCGt0ZLKWAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "colors = {'Tree':'blue', 'Hamming':'red', 'Jukes-Cantor 69 (JC69)':'green'}\n", "for k in sorted(pairwise_distances.keys()):\n", " distplot(pairwise_distances[k], kde=True, hist=False, color=colors[k], label=k);\n", "plt.title(\"Pairwise Distance Distributions\");\n", "plt.xlabel(\"Pairwise Distance\");\n", "plt.ylabel(\"Kernel Density Estimate\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also numerically compare these distributions by looking at their averages." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Average Pairwise Hamming Distance: 0.264545\n", "Average Pairwise Jukes-Cantor 69 (JC69) Distance: 0.333181\n", "Average Pairwise Tree Distance: 0.324325\n" ] } ], "source": [ "for k in sorted(pairwise_distances.keys()):\n", " print(\"Average Pairwise %s Distance: %f\" % (k, sum(pairwise_distances[k])/float(len(pairwise_distances[k]))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ideally, for two entities *u* and *v*, we would want the pairwise distance between *u* and *v* to be the same whether we are comparing their sequences or comparing them on the true phylogeny. As expected, because of loss of information due to sites that mutate multiple times, the Hamming distance between any two points is significantly lower than the true phylogenetic distance, and although the JC69 correction brings these pairwise sequence distances closer to their true phylogenetic values, it's not perfect.\n", "\n", "We can observe this more closely by looking at a plot of pairwise sequence distances vs. pairwise tree distances." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VMXawH+TDiEQmiAECCA1AQIEAZEiIIIK4idCIoigooA0r5RruRoQr9JEmgW5igoSiuUiotJBUS41QEDpLUgnCQnJpu18f8zZzSbZJJu6SZjf8+TJnpk5M++p75n3nXlHSCnRaDQajSY7XJwtgEaj0WhKNlpRaDQajSZHtKLQaDQaTY5oRaHRaDSaHNGKQqPRaDQ5ohWFRqPRaHLkjlMUQojOQohjhVRXvBCiQWHUpXEMIcQRIUQ3Z8tRFhBC/CSEeKYQ6nlNCLGkMGQqLQghugkhjjhbjuJClNZ5FEKIs0ANIA24DfwEjJFSxjtTroIihHgMmAo0AJKBQ8BzUsozThWsiBFC+ANnUNcS4/8eYJ6UcmM+63KXUqYWmpAlECHEMOA/QCJgBk4Db0gp1zlTroIghPgYGGJsegACSDK2f5VS9ini9qcDUwCTkfQ38Avwbynl5XzU5SelHFaoQhYzpb1H0VdKWQFoAwQDbxSkMiGEW6FIlf/27wG+BF4BKgH1gUUoZXin4Gtc01bARuA742WoyZ4/jHPmi1Iaq4QQlQtSoTOfBSnlSCllBeOY/g2stGzbUxJFJOtyKaUPUBV4AqgD7BVC1CiCtko+UspS+QecBXrabM8C1hm/hwN/AnGoL6wXbcp1A6Iy1TMF9eWeBIwAfrDJPwGsttm+AAQZvyVwj/H7YeCo0eZFYKLNPo8CEUAM8DvQMptjGgBE5HDMLsA/gVPADWAVUMUm/2ngnJH3uu05ApYC03M4D7WAb4BrqK/xcTZ5YUZbXxrHdwQItsmvA3xr7HsDWGiT96xxLaJRX2X1sjk2f+N8umVKnwhcAVwyX3fgXmAvcMso876Rft6oK9746wg0BLYY8l0HlqOUku19MNG4D2KBlYCXTf5jxjW8ZZz/3kZ6JdTL+ZJx3acDrnaOrxbqq9/2erU2ZHEH7gG2G21fR70cHXkOhgG/2Wx7G8ceDFQG1hnXJdr47WdTdhvwvE09O4G5xjmajrqX2hr5g416A4zt54Dvbe6PZcZvL2CZUUcMqldYIy/nKtPxWeu2SbvHkGW4ca23GOmdgF1GuxFAF5t9fIHPjbajgGkY95SdNqcDSzOluQGRwHvGdk/grE3+a6iexy3gL9Tz9SjKKpCCug/3GWWfJ/39dMpyDWzrBSYb1+1vYKhNfnnjGp037pUdgGdux1/QP6e/8PMteMYXRh3Uy+ttY/sR1ItBAF2BBKCNkdeNrIoiwqijHMrkE4N6KdcyHpYoo2wD1ANneWnZKopLQGfjd2Wb9loDV4H2gCvwjNGmp51jaoDq7s4FHgAqZMofb9wIfoAn8AmwwshrbtyMXYy894FUHFAUxrHuA95EdfUboBTsQzYPqwmlDF2Bd4FdRp4rcNCQ2Rv1orjfyHsMOAk0Qz1obwC/Z3M9/bGvKBoY6c3sXPc/gKeN3xWADtnVhXq5PGicm+qoB+yDTPfBbuOaV0E9yCONvHtRD+WDxrmqDTQ18r4zroM3cJdRx4vZHOMWYITN9izgY+P3CpRyd7E9hw48B8MwFIVxjsejXkCVSP8aLg/4AKsxXu5G+W1kVBSpwFijnnIYvVsjfzHqpTbK2P4SeNnm/rAoiheBH4w2XYG2QMW8nisbGa11Z7qWEvXiL2/IWgelnB4yzmFvlMKtauzzA/ChUb4G6n5/Lps2sygKI/3fwE7jt1VRAAGo90RNY7s+0CC7uoC+qPtaAN1RHxAtbepNBd5CfUD0Q5lhLefwE2AzcLdxfu83yuV4/AV+3xZGJc74Qz3Y8aiX+jnjJiiXTdnvgfHG725kVRTPZip/AWXOCjEekN1AU9QXzFqbcraK4rzxkFTMVNdHGArMJu0Y0DUbWTugvt6voV7OSzEUBurl1cOm7N2orxU31Es+3CbPG/U144iiaA+czyTHq8DnNg/rJpu85kCi8bujIaubnWP5CZuH0biBE7DTqyB7ReFlpHeyuV6WY9qB8udUc6SuTGX6Awcy3QdDbLZnkv4S/wSYa6eOGqheaDmbtFBgazZtPk/6168w7rMuxvaXxr3ml53M2dQ5DPViiUG9GHZh09POVDYIiLbZ3kZGRZH5HngO43437r3nLfcY6pmzfAyFka4onsVOrzmv58qmjLVumzSLoqhrk/a65X61SduM6gnVRr2MPW3yngY2ZtNmdopiDPCn8dtWUTRB9Wh72Ll/7daVqcw64CWbeuOx6WkBN1E9RFfjHAbYqSPb48/L/ZTdX2n3UfSXUvpKKetJKUdLKRMBhBB9hBC7hBA3hRAxqC/hajnUcyHT9nbUi7SL8XsbqmfS1di2xxNGO+eEENuFEB2N9HrAK0KIGMsfSvvXsleJlHKXlHKglLI60NmQ4XWbur6zqedPlP+ihlHfBZt6bqO+MByhHlArk4yvGfVasHXiJQBehm24DnBO2nca1wPm2dR5E/WCrO2gXNiUvWkn7zmgMfCXEGKPEOLR7CoRQtQQQoQLIS4KIW6hzCOZ74nMx1jB+F0H9TWdmXqor7lLNsf4Cepr2R7fAB2FEHejrqsZ+NXIm4w6N7uNkV3PZncsdthlPAfVpJQdpJSbjGMuL4T4RAhxzjjmHYCvEMI1m3rsPQedDXldUR8wnYzBApVQPfHMfIUyMYYLIf4WQswUQriT93PlCLby1gNCM93DHVDPRT1UT/KKTd4iMt7fjlAbO/ehlPIYyq84DbgqhFghhKiZXSVCiEeFEP+zeT/1IuO9eF1KaeuXtNyLNVA9/uzuxeyOv8A41XlbFAghPFEP5FDgv1LKFCHE96iHMDtkpu3tqO5hfVR3Mwb1ZdIRWGi3Ain3AI8ZD8UY1ENVB3UzvyOlfCevxyKl3COE+BYINJIuoHo/OzOXFUJcQpl4LNvlUaYHC7dR3W4LtjfyBeCMlLJRXmU09q0rhHCzoywsx748H/VaeBxlussypFlKeQL1cLgA/wesEUJUJev1BHUdJdBCSnlTCNGfbK6lHS6gTJn20pNQPZpcR1dJKaOFEBuAQahrFS6NTz+pRtOMABBC3A9sEkLskFKedFBGe7yC+tptL6W8LIQIAg6Q/bOQ4bxJKU8KIRJQ5qgdUspbQojLwAsoc5fZzjGmoHp5Uw2Fsh517daTh3PlCJZzZ3AB9UU9KnM5IUQd1Mu2ij2ZHcFQrn1RX//2ZFkGLBNCVAI+RZlnh5PpnAohygFrUNaKH4330zpyfj9ZuIKyEjREmdptyfb4C4PS3qOwhwfq6+EakCqE6IPS2HlhO8pHUE5KGYX66uuNevEeyFxYCOEhhBgshKhkPCi3UF+LoG6akUKI9kLhLYR4RAjhY6ee+4UQI4QQdxnbTVE2yl1GkY+Bd4QQ9Yz86sZwWlA336NGHR4Yzjqb6iOAh4UQVYyvnQk2ebuBOCHEFCFEOSGEqxAiUAjRzoFztRvln3nPODYvIUQnG3lfFUIEGPJWEkI86UCdlh7AGJSt9lV7D7gQYogQorqRF2Mkm1HX3oyyA1vwQXXpY4UQtYFJjshh8B9guBCihxDCRQhRWwjRVEp5CdgAzBFCVDTyGgohuuZQ19eoj5gBxm/LsTwphPAzNqNRL5h8vdRs8EGZXGKEEFVQ5zKvbEd9+Fh60tsybWdACPGAEKKF8WK9hTKNmvN5rvLCV8DjQogHjfvXy5CllpTygiHvbJu27xFCdMmtUiGEuxCiORCO8l19YKdMM6MtT9T5tgxVBvVy9xdCWBSBJ+oddQ1IE6oX3MORAzR6GUuBD4QQNY3j7GR8nGZ7/I7UnRtlTlFIKeOAcagv+mjgKWBtHus4jnqp/Gps30I5d3dm6hLa8jRwVqgu/khUDwQp5V7Ul+JCQ56TKHuwPWJQiuGwECIe+BnlAJxp5M8zjmWDECIOpUDaG+0cAV5CvXwuGW1F2dT9FcrpfBb1wK60Od401AiNINSIp+vAEpR5IUeMffui7MbnjTYHGXnfATNQZohbqFEjuY2BjxFC3AYOo0x5T0opP8umbG/giHGu5gEhUspEKWUC8A6w0+iGd0B95bZBOaV/RI3Scggp5W7U1+FcY//tqK4+qJe+B2rEWzRKYd+dQ3VrgUbAZSnlQZv0dsD/jGNZi/KpnQbrJMPBjsprwwcoR6/Fd/FzPurYjlI4O7LZzkxN1Dm4hTKNbkfde5D3c+UwUsqzqN7nv1Av4fOoHpXlHTcE5beztL2ajL3qzAw2nrFo4L+oF36wtD+PwhP1jF5HmS8rk24uXok65ptCiN1SyhjgZdRzfRP1wZCXOS8vo87rPmP/f6Pmw+V2/AWi1E640+SOUJMSn7fYrDUajSY/lLkehUaj0WgKF60oNBqNRpMj2vSk0Wg0mhzRPQqNRqPR5Eipm0dRrVo16e/v72wxNBqNplSxb9++68ZE3jxT6hSFv78/e/fudbYYGo1GU6oQQpzL777a9KTRaDSaHNGKQqPRaDQ5ohWFRqPRaHKk1Pko7JGSkkJUVBQmkyn3whpNLnh5eeHn54e7u7uzRdFoSgRlQlFERUXh4+ODv78/6bG3NJq8I6Xkxo0bREVFUb9+fWeLo9GUCIrM9CSE+EwIcVUIEZlNvhBCzBdCnBRCHBJCtMlvWyaTiapVq2oloSkwQgiqVq2qe6cajQ1F6aNYiorumR19UFE0G6Hi239UkMa0ktAUFvpe0mgyUmSmJynlDmPhkux4DPjSWHxklxDCVwhxtxG3XqPRaDSFwNatW7l0qWCvVWeOeqpNxqUMo8hmiUwhxAtCiL1CiL3Xrl0rFuHySoUKFTJsL126lDFjxhRb+3///TcDBgwotvY0Gk3J5sKFCwwcOJDu3bszalTBFr4rFcNjpZSLpZTBUsrg6tXzNQO9zFOrVi3WrFnjbDE0Go2TMZlMTJ8+nSZNmrB69WoAAgMDc9krZ5w56ukiak1pC35GWv5JTYWoqNzL5Qc/P3DL3+n64YcfmD59OsnJyVStWpXly5dTo0YNwsLCOHPmDKdPn+b8+fPMnTuXXbt28dNPP1G7dm1++OEH3N3d8ff3JzQ0lJ9++gk3NzcWL17Mq6++ysmTJ5k0aRIjR47k7NmzPProo0RGRrJ06VLWrl1LQkICp06d4vHHH2fmTLVI3n/+8x9mzJiBr68vrVq1wtPTk4ULHV06WqPRlGTWrVvH+PHjOX36NAA1a9Zk5syZDBkyBBeX/PcLnKko1gJjhBDhqOU8Ywvsn4iKgqIa0njmDOQQjDAxMZGgoCDr9s2bN+nXrx8A999/P7t27UIIwZIlS5g5cyZz5swB4NSpU2zdupWjR4/SsWNHvvnmG2bOnMnjjz/Ojz/+SP/+/QGoW7cuERERvPzyywwbNoydO3diMpkIDAxk5MiRWeSJiIjgwIEDeHp60qRJE8aOHYurqytvv/02+/fvx8fHh+7du9OqVatCPEkajcaZbNiwgdOnT+Pm5saECRP417/+RcWKFQtcb5EpCiHECqAbUE0IEYVa2N0dQEr5MbAetSbySSABtSZxqaVcuXJERERYt5cuXWoNXhgVFcWgQYO4dOkSycnJGcbn9+nTB3d3d1q0aEFaWhq9e6uBYi1atODs2bPWchal06JFC+Lj4/Hx8cHHxwdPT09iYmKyyNOjRw8qVVJLXjdv3pxz585x/fp1unbtSpUqVQB48sknOX78eOGeCI1GU2zEx8fj6elpnRw6bdo0rl69SlhYGE2bNi20dopy1FNoLvkSeKlQG/XzU1/+RYGfX753HTt2LP/4xz/o168f27ZtIywszJrn6ekJgIuLC+7u7tahmS4uLqSmptotZ/ltr1zm8gCurq52y2g0mtKJlJIVK1YwadIkJk2axIQJEwDw9fUlPDy80NsrEzOzrbi55WgechaxsbHUrq0GdH3xxRdOk6Ndu3ZMmDCB6OhofHx8+Oabb2jRooXT5NFoNHnn4MGDjB07ll9//RWAGTNmMGrUqAwfh4VNqRj1VNoJCwvjySefpG3btlSrVs1pctSuXZvXXnuNe++9l06dOuHv7281T2k0mpLNzZs3GTNmDG3atLEqif/7v//jjz/+KFIlAaguTGn6a9u2rczM0aNHs6Rp7BMXFyellDIlJUU++uij8ttvv3WyRCUTfU9pSgqpqanyk08+kVWrVpWABGSzZs3kxo0b81QPsFfm872rexR3GGFhYQQFBREYGEj9+vWto6o0Gk3JJDU1lTlz5nDjxg18fHyYM2cOBw8epGfPnsUmQ9nyUWhyZfbs2c4WQaPR5EJiYiLlypUD1MCUDz74gPDwcGbMmEHNmjWLXR7do9BoNJoSQkpKCnPmzKFOnTr89ddf1vQ+ffrwxRdfOEVJgFYUGo1GUyLYuHEjLVu2ZOLEidy4cYO33nrL2SJZ0YpCo9FonMjZs2d54okn6NWrF3/99RcuLi689NJLfPRRgVZeKFS0j0Kj0WicQGJiIrNmzeLdd9+1LpR1//33s2DBggzhgEoCukdRiGQONZ4Zf39/rl+/XqhtXr58mZCQEBo2bEjbtm15+OGH8x2WY+nSpfz9998FlmnVqlU0b96cgIAAnnrqKWv6lClTCAwMJDAwkJUrV2a7/4QJE9ixYwcA3bp1s4ZCiY+P58UXX7Qea7du3fjf//4HQExMDAMGDKBp06Y0a9aMP/74A1CTkzp27EiLFi3o27cvt27dAuDw4cMMGzaswMeq0eSX48ePM3XqVEwmE7Vq1WL58uXs2LGjxCkJQM+jKEy8vb1zzK9Xr568du1aobVnNptlhw4d5EcffWRNi4iIkDt27MhXfV27dpV79uzJ0z6pqakZto8fPy6DgoLkzZs3pZRSXrlyRUop5bp162TPnj1lSkqKjI+Pl8HBwTI2NjZLfdevX5ft27e3K9OgQYPkP//5T5mWliallPL06dNy3bp1Ukophw4dKj/99FMppZRJSUkyOjpaSillcHCw3LZtm5RSyv/85z/yjTfesNbdo0cPee7cObvHVVLuKU3ZIiUlJcP2+PHj5ZQpU+StW7eKvG0KMI/C6S/+vP7lpChS0lLkmegzRfKXkpaSpd3MeHt7y61bt8pHHnnEmvbSSy/Jzz//XEqZrigSEhJk79695eLFi6WUUn711VeyXbt2slWrVvKFF16QqampMjU1VT7zzDMyICBABgYGyvfffz9Le5s3b5adO3e2K0tcXJzs3r27bN26tQwMDJTff/+9lFLKM2fOyKZNm8rnn39eNm/eXD744IMyISFBrl69Wnp7e8vGjRvLVq1ayYSEBLlp0yYZFBQkAwMD5fDhw6XJZLIex+TJk2Xr1q3lihUrMrQ7adIk6wvblpkzZ8pp06ZZt5999lm5cuXKLOU++eQT+dZbb1m3LYri5MmT0t/fP4tiklLKmJgY6e/vL81mc5a8ihUrWtPPnz8vmzVrZs374IMP5IwZM+yeP60oNIVJbGysnDhxomzXrp3de7g4KIiiKFM+iqhbUdSfVzRhxs+MP4O/r3+B64mPjyckJIShQ4cydOhQ/vzzT1auXMnOnTtxd3dn9OjRLF++nICAAC5evEhkZCSA3QixkZGRtG3b1m47Xl5efPfdd1SsWJHr16/ToUMHawTaEydOsGLFCj799FMGDhzIN998w5AhQ1i4cCGzZ88mODgYk8nEsGHD2Lx5M40bN2bo0KF89NFH1uBjVatWZf/+/VnatZi9OnXqRFpaGmFhYfTu3ZtWrVoxdepUXnnlFRISEti6dSvNmzfPsv/OnTvtrtR35MgRgoKCcHV1zZJ35swZqlevzvDhwzl48CBt27Zl3rx5eHt7ExAQwH//+1/69+/P6tWruXAhfVHF4OBg3nvvPSZPnmz3HGo0BUVKyfLly5k0aRKXL18GVLy3Z5991smS5Q3toyhmHnvsMYYPH87QoUMB2Lx5M/v27aNdu3YEBQWxefNmTp8+TYMGDTh9+jRjx47l559/znNMeSklr732Gi1btqRnz55cvHiRK1euAFC/fn2rHbRt27YZwplbOHbsGPXr16dx48YAPPPMM1a/AcCgQYPstpuamsqJEyfYtm0bK1asYMSIEcTExNCrVy8efvhh7rvvPkJDQ+nYsaPdl/6lS5fI6yqGqamp7N+/n1GjRnHgwAG8vb157733APjss8/48MMPadu2LXFxcXh4eFj3u+uuuwrFJ6PR2OPAgQN07tyZp59+msuXL1OuXDnefvvtDH670kKZ6lH4VfTjzPiiCTPuV9GxMONubm6YzWbrtmU0g4VOnTrx888/89RTTyGEQErJM888w7vvvpulroMHD/LLL7/w8ccfs2rVKqZOnUrfvn0BGDlyJAEBAdkuf7p8+XKuXbvGvn37rKvkWWTJHII8MTHRoWOzxdvb2266n58f7du3x93d3apoTpw4Qbt27Xj99dd5/fXXAXjqqaesSsiWcuXKZTlnAAEBARw8eJC0tLQsCsbPz8/aLsCAAQOsiqJp06Zs2LABUL2dH3/80bqfyWSyzn7VaAqLGzdu8MYbb/DJJ58o+z5q7ZfZs2dTt25dJ0uXP8pUj8LNxQ1/X/8i+XNzcUyn1qtXj6NHj5KUlERMTAybN2/OkD9t2jQqV67MSy+ppTh69OjBmjVruHr1KqAiRFoWGTKbzTzxxBNMnz6d/fv3U6dOHSIiIoiIiGDkyJF0796dpKQkFi9ebK3/0KFD/Prrr8TGxnLXXXfh7u7O1q1bOXfuXK6y+/j4EBcXB0CTJk04e/YsJ0+eBOCrr76ia9euudbRv39/tm3bBsD169c5fvw4DRo0IC0tjRs3blhlPHToEL169cqyf7Nmzaxt2tKwYUOCg4N56623rA/f2bNn+fHHH6lZsyZ16tTh2LFjgOqlWcxalvNqNpuZPn16htUAjx8/XuC1hDWazCxZsoSPP/4YKSUBAQFs3ryZVatWlVolAWVMUTiT1NRUPD09qVOnDgMHDiQwMJCBAwfSunXrLGXnzZtHYmIikydPpnnz5kyfPp1evXrRsmVLHnzwQS5dusTFixfp1q0bQUFBDBkyxG6PQwjBd999x6ZNm2jYsCEBAQG8+uqr1KxZk8GDB7N3715atGjBl19+6dBqV8OGDWPkyJEEBQUhpeTzzz/nySefpEWLFri4uNhdcjUzDz30EFWrVqV58+Y88MADzJo1i6pVq5KSkkLnzp1p3rw5L7zwAsuWLcPNzhrkjzzyiFXR2J5XUA/glStXuOeeewgMDGTYsGHcddddACxYsIDBgwfTsmVLIiIieO211wBYsWIFjRs3pmnTptSqVYvhw9MXUty6dSuPPPJIrsek0eSG5eMFYPz48bRu3Zq5c+dy4MABunfv7kTJCon8esGd9VdSh8dGRETIdu3aOVuMMkGnTp1kdHS0NJlM0s/PT8bExBR6GyaTSbZv3z7LcEULJeGe0pR8Ll68KIcMGSJff/31DOmWIdwlCXSYcefy8ccfExoayvTp050tSplgzpw5HDp0iKCgIEaPHl0kiyudP3+e9957z26vRqPJjeTkZGbNmkWTJk1YtmwZs2bNyjAoxMWlbL1ahbTpMpUGgoODpWWmroU///yTZs2aOUkiTVlE31Oa7Pjll18YP3681SdWvXp13n33XYYPH16iFYQQYp+UMjg/+5bco9JoNJoSxJkzZ+jfvz+9e/fm2LFjuLq6Mm7cOI4fP85zzz1XopVEQdH9bo1Go3GA4cOHs337dkDFIJs/fz4tWrRwslTFQ9lVgRqNRlOIzJo1i7p16xIeHs6WLVvuGCUBWlFoNBpNFo4ePUqvXr2sUYgB2rVrx8mTJxk0aBBCCCdKV/w4ZHoSQtQDGkkpNwkhygFuUsq4ohWt9HDjxg169OgBqLDfrq6u1jAUu3fvzhA2QqPRlFxiY2OZOnUqCxYsIDU1lZs3b7J7926r/8Hd3d3JEjqHXBWFEGIE8AJQBWgI+AEfAz2KVrTSQ9WqVYmIiAAgLCyMChUqMHHixAxlrOORy7DDS6MprZjNZr766iumTJlijYl2zz33MHXqVP3M4liP4iXgXuB/AFLKE0KIu4pUqgJiL8idLdWrV7fGKkpLS8sQUdQe/v7++ZLj5MmT9OvXj9atW3PgwAF++uknWrVqZY0EGx4ezqZNm6wzjkeNGsX58+dxcXFh/vz5dOjQIV/tajQax9m3bx9jxoxh165dAJQvX5433niDf/zjHxniot3JOKIokqSUyRabnBDCDSjRky/q18851PiaNWt44oknABW+O7fyBZlr8tdff/Hll18SHBxMampqtuXGjRvH5MmT6dChA2fPnuXRRx+1hhjXaDRFw82bN+nSpQsJCQkAhISEMGvWLPz8HAsCeqfgiKLYLoR4DSgnhHgQGA38ULRilR0swexyY9OmTdYJPADR0dEkJibq6KYaTRFSpUoVXnnlFb7//nsWLFjgUODLOxFHFMU/geeAw8CLwHpgSVEKVVDOnMk51Ljtege+vr65li8ItuG4XVxcMvRObMNpSym141ujKWJ27NjB7NmzWblypfUj7PXXX+fNN98sk+FczNLMraRbxCfHF6geR85MOeAzKeWnAEIIVyMtoUAtFyF58Sm4urrm2weRV1xcXKhcuTInTpygYcOGfPfdd1al1bNnTxYtWsTLL78MQERERMlcZF2jKYVERUUxadIkwsPDAZg5cyZvvfUWQJn0Q6SkpRBjiiE2KRazNOPuUrDRWo648zejFIOFcsAmRyoXQvQWQhwTQpwUQvzTTn5dIcRWIcQBIcQhIcTDjoldepkxYwYPPfQQ9913XwY76KJFi9i5cyctW7akefPmfPrpp06UUqMpGyQlJfHee+/RtGlTq5Jo3749Dz9cNl81CSkJXLx1kTMxZ4g2RWOW5tx3coBcgwIKISKklEG5pdnZzxU4DjwIRAF7gFAp5VGbMouBA1LKj4QQzYH1Ukr/nOrVQQE1xYG+p0o/69evZ8KECZw4cQI1Eq0UAAAgAElEQVRQS9/OmDGDoUOHlqkhrxbzUowphuS0ZLtlTt08xcONH853UEBHTE+3hRBtpJT7AYQQbQFH1s68FzgppTxt7BcOPAYctSkjActi0JUAvYCxRqMpMOvXr7cuSmUJ3vfWW28VSch6Z5HZvJSZpNQkfj71M+GR4ey/tL9AbTmiKCYAq4UQfwMCqAkMcmC/2oDtBIUooH2mMmHABiHEWMAb6GmvIiHEC6hJf6V6OUGNRlM8PPTQQ7Rp04bKlSszf/5869K4ZYGElARiTDHZOqjPxpxlZeRKvv3zW2KSYgqlzVwVhZRyjxCiKdDESDompUwplNYhFFgqpZwjhOgIfCWECJQyo3qUUi4GFoMyPWUj5x0Xf0VTNJS2NVrudKSUrF69msjISKZNmwaoXsTGjRupXLlymXgvSCm5lXSLaFO0XfNSSloKW85sIfxIOL9f+N2a7ubiRs8GPRnSYghDwobku31Hx4O1A/yN8m2EEEgpv8xln4tAHZttPyPNlueA3gBSyj+EEF5ANeCqg3IB4OXlxY0bN6hatWqZuCk0zkNKyY0bN/Dy8nK2KBoHiIyMZNy4cWzduhUhBP369bPOW6pSpYqTpSs4KWkpxCbFEmuKJU2mZcm/FHeJVUdXsfrIaq4lXLOm1/KpxZPNn2RA8wHc5X1XgUc9ORLr6StUjKcIwCKpBHJTFHuARkKI+igFEQI8lanMeVTMqKVCiGaAF3CNPOLn50dUVBTXruV5V40mC15eXnpmbgknJiaGsLAwFi5cSFqaei317duXatWqOVmywiExJZFoU7Rd81KaOY3fzv/GisgVbD+33eqfEAi6+XcjJDCEznU74+riWmjyONKjCAaayzz2x6WUqUKIMcAvgCtqLsYRIcQ01CLfa4FXgE+FEC+jlM+wvLYDKqJjbmE4NBpN6cdsNrN06VJeffVVrl5VhodGjRoxb948+vTp42TpCobFvBRjiiEpLSlL/vWE63xz9BtWHlnJxbh040y18tUY0HwAA5sPpHbF2nbrdnct4h4FEIlyYF/Ka+VSyvWomdy2aW/a/D4KdMprvRqN5s5k8uTJzJkzB1BRD/71r38xYcKEUj1pLtWcqkYv2TEvSSnZ8/ceVkSuYOOpjaSY093DHfw6EBIYQo/6PfBwzRrRwUW44OPhQyWvSni5FcyU6oiiqAYcFULsBqxqTkrZr0AtazQaTR558cUXWbBgAQMGDGDmzJnUrm3/C7o0kJiSaB29JDPFWY01xfL9se8JjwzndPRpa3olz0r0b9qfkMAQGlRuYLfe8u7lqehZER8Pn0Lz2TqiKMIKpSWNRqPJA6mpqXz44Yc0adKEhx56CFBmppMnT1KnTp1c9i6ZSCmJS44jxhSDKdWUJe/w1cOsiFzB+hPrM+QH1QgitEUove/pbbd34ObiRkXPilTyrFRgM5M9HBkeu73QW9VoNJoc2LZtG2PHjiUyMpKGDRsSGRlpHYlWGpVETual28m3WXd8HeFHwjl6LX0+cnn38vRr0o/QwFCaVmuapU6BoIJHBSp6VsTbwztLfmHiyKinDsACoBnggXJM35ZSVsxxR41Go8kjFy5cYOLEiaxatQoAIQQ9e/YkOTm5VA5ZNqWaiE6MtmteOnb9GOFHwvnvX//ldspta3qTqk0IbRFK38Z9qeBRIUudnq6eVPSsSEXPioU6siknHDE9LUQNbV2NGgE1FGhclEJpNJo7C5PJxJw5c/j3v/9tXUTovvvuY8GCBbRp08bJ0uWNnMxLlrAaKw6v4MDlA9Z0D1cPHr7nYUICQwiqGZTFt1CYjun84NCEOynlSSGEq5QyDfhcCHEAeLVoRdNoNHcCUko6d+6MJdhnjRo1mDlzJkOGDClVwftSzanEmmKJTYol1ZxxNcuzMWdZecQIq2FKD6vh7+tPaGAo/Zv2x9fLN0udReGYzg+OKIoEIYQHECGEmIkaJlt6rp5GoynRCCF45plniIiIYPz48bz55ptUrFh6LNumVBMxphjikuIymJdS0lLYenYrKyJX2A2rERoYSvva7bMogKJ2TOcHR8KM1wOuoPwTL6OivC6SUp4qevGyYi/MuEajKT3Ex8cze/Zsxo0bZw2zkZqayunTp2ncuHRYtXMyL+UUVmNgwEAGNBtAde/qGfYpDse0EKJIw4z3l1LOA0zAVKPB8cC8/DSo0WjuTKSUhIeHM2nSJC5evMj169dZuHAhAG5ubqVCSaSZ06yhvW3NS2nmNH678BvhkeFsO7stQ1iNrvW6EhIYQpd6XbI4n53hmM4PjiiKZ8iqFIbZSdNoNBq7HDp0iLFjx7Jjxw5AxdOqUaOGk6VynOzMS/kJq+EiXKzKwRmO6fyQraIQQoSigvjVF0KstcmqCNwsasE0Gk3pJzo6mjfffJMPP/wQs1l9ZT/++OO8//77xbZWfX6RUhKfHE+0KTqDeSm3sBqDAgbRs0HPLGE1SopjOj/k1KP4HeW4rgbMsUmPAw4VpVAajab0c/78edq2bcv169cBaNq0KfPnz+fBBx90smQ5k515KaewGo83fZxBgYOyhNUoiY7p/JCtopBSngPOCSF6AolSSrMQojHQFDhcXAJqNJrSSZ06dQgKCmLXrl289dZbjBs3Dg+PrMHrSgpJqUlEm6IzmJfyE1ajOGdMFxeO+Ch2AJ2FEJWBDah1JgYBg4tSMI1GU7q4cuUKq1evZsyYMYAa9vrpp5/i6enJ3Xff7WTp7GMxL8WYYkhMTbSm306+zY8nfmRF5AqHw2qUFsd0fnBEUQgpZYIQ4jngQynlTCFERFELptFoSgcpKSksXLiQsLAwbt26RePGjenVqxdAifVDpJnTiE2KJcYUk8G8lNewGqXRMZ0fHFIUxnrWg1FLl4KK96TRaO5wtmzZwtixYzl6VH11V6lShdjYWCdLlT1JqUnEmGK4lXTLal6yhNUIjwxn/6X91rIerh70uacPoYGhWcJqlGbHdH5wRFFMQIXr+M5Yoa4BsLVoxdJoNCWZc+fOMXHiRNasWQOAi4sLL774Im+//TZVq1Z1snRZiUuKy2JeymtYDTcXNyp5VqKiZ8VS7ZjOD46GGd9us30aGFeUQmk0mpLL2rVrCQkJITFRvXTvv/9+FixYQFBQkJMly4g981JuYTVCAkPoULuDtZdgcUxX8qpEeffyTjmOkkBO8yg+kFJOEEL8AGSJ86FXuNNo7kzuvfde3NzcuPvuu5k1axZPPfVUiTK/2DMv5RRW48nmTzKg+QDu8r7Lmu7p6kklr0r4ePiUOcd0fsipR/GV8X92cQii0WhKJseOHePcuXNWB3XNmjX58ccfCQoKwsfHx8nSKeyNXjJLM7+e/9VuWI0u9boQGhiaIazGneKYzg85zaPYZ/zfLoSobvy+ll15jUZTtoiLi2P69OnMnTuXypUrc/z4cSpVqgRA586dnSydwp556XrCdb7981vCI8MdCqtxpzmm80OOPgohRBgwBhVWXAghUoEFUsppxSCbRqNxAlJKvv76ayZNmsSlS5cA8Pb25ty5c7Rs2dLJ0ikym5eklOy+uJvwI+F2w2qEBIbQo34Pa1iNsjJjurjIyUfxD6AT0E5KecZIawB8JIR4WUo5t5hk1Gg0xURERARjx47lt99+A6BcuXK8+uqrTJo0qUQsRWoxLyWkqFXw8hJWQyDw9vCmkmelMjNjurjIqUfxNPCglPK6JUFKeVoIMQQ1Q1srCo2mDDF16lSmTZtmDd43YMAAZs+eTb169Zwql1maiTUp81KKOSXPYTU8XD2sw1q1Yzp/5KQo3G2VhAUp5TUhhO6raTRljCZNmmA2m2nevDnz58+nR48eTpUnOS2Z6MRo4pLjMEtzrmE1QgJCaFa9GeD8NabLGjkpiuR85mk0mlLA77//Tp06dahTpw4AgwYNQkrJgAEDcHd33rfg7eTbRJuireal4zeOsyJyhd2wGiGBIfRr0s8aVqOcWznrsFbtmC48clIUrYQQt+ykC0CraI2mlHLp0iUmT57MsmXLePLJJ1m1ahWggviFhoY6RabM5qWcwmo8fM/DhLYIpVWNVgghtGO6GMhpeKw25mk0ZYjk5GTmzZvHtGnTiI+PB+DMmTPcvn0bb2/nOHeT05Kto5fM0uxwWA3tmC5eHIn1pNFoSjkbNmxg3LhxHDt2DIBq1arx7rvv8uyzz+Li4lLs8tial/ISVsPD1cPae9CO6eJDKwqNpgyTkpJCSEgI3377LaCC940ePZpp06ZRuXLlYpXFLM3cSrpFdGI0KeYUh8NqCAQ+nj5U8qxEOfdyxSqzRlGkikII0RuYhwpLvkRK+Z6dMgOBMFQ8qYNSyqeKUiaN5k7C3d0dT09PALp06cKCBQuKfdKcrXkpJS2F3y78ZjesRtd6XQltEUrnup1xdXHFy82LSp6V8PH0wUUUf6+nzGA2g8mUe7kccEhRCCHqAY2klJuEEOUANyllXC77uAKLgAeBKGCPEGKtlPKoTZlGqBDmnaSU0UKIu+zXptFoHEFKyZYtW+jevbt11M+sWbPo27cvISEhxToS6HbybWJMMdxOuc31hOt8c/QbVh5ZmWNYDVfhao235OnmWWyyljlSUyE+Hm7fhoQEcCtYnyDXvYUQI4AXgCpAQ8AP+BjIbZD1vcBJIyw5Qohw4DHgqE2ZEcAiKWU0gJTyal4PQKPRKP766y/GjRvHxo0bWbZsGYMHq9WKa9euXWyjmSzmpRhTDEmpSez5ew8rIlfYDasxKGAQPRv0xMPVg/Lu5ankWYkKHhX0sNb8kpSklEN8vPpty82bBaraETXzEuql/z8AKeUJB7/8awMXbLajgPaZyjQGEELsRJmnwqSUP2euSAjxAkpZUbduXQea1mjuHG7dusW0adOYN28eqakqMN6WLVusiqI4sDUvRSdGZxtWo3/T/oQEhtCgcgPcXdyVY9qrEm4u2l2aZ6SExMR05ZCavqQrZjP8+Sds367+Dh4sUFOOXJ0kKWWydSEPIdywsz5FAdpvBHRD9VR2CCFaSCljbAtJKRcDiwGCg4MLq22NplRjNptZtmwZU6ZM4fLlywA0bNiQDz74gEcffbRYZLCYl+KT47MNq9GqRitCA0Pp06gP5dzK6YWACkJamjInWf6McCuAUhY7d8K2bfDrr3Ct8IJ9O6IotgshXgPKCSEeBEYDPziw30Wgjs22n5FmSxTwPyllCnBGCHEcpTj2OFC/RnPHcvbsWQYPHszvv6vhpOXLl+e1117jlVdeKfLgfbbmpejEaIfCauiFgApASkp6r8FkUj0JUP9PnVI9hm3bYP/+jL0KgLp1oWtX6NEDhg3LtwiOKIp/As8Bh4EXgfXAEgf22wM0EkLURymIECDziKbvgVDgcyFENZQp6jQajSZHqlSpwunT6lEZOHAgs2fPtobiKCpS0lKIMcUQmxTLn9f+JPxIuN2wGqEtQunbuK/VKV3Js5J2TOcVkyldOSTbRExKTIT//S/dpHQx07e3uzu0a6eUQ5cuUL8+CKHSC4AjiqIc8JmU8lOwjmYqByTktJOUMlUIMQb4BeV/+ExKeUQIMQ3YK6Vca+T1EkIcBdKASVLKG/k/HI2mbJKamsrhw4dp3bo1ABUrVuTTTz/F29ubBx54oEjbTkhJIMYUw42EG/x86mdWHF7BgcsHrPmZw2pYZkxrx3QekFKNTrIoh7S09LwLF1SPYft2pSSSM4Xaq1kzXTF06AAVKqTnubpC+fIZ0/KBkDJnk78QYhfQU0oZb2xXADZIKe8rUMv5JDg4WO7du9cZTWs0TuHXX39l7NixnDp1imPHjlGrVq0ib9MszcQlxRFtiub4jeO5htWoVr6aNZS3jrfkIGlp6UNYb99ONyklJ8Pevem9hjNnMu7n6gqtW6crhyZNVK8B1H8vL/D2VgrCxgwphNgnpQzOj6iO9Ci8LEoCQEoZL4TQXiiNpoi5ePEikydP5uuvvwZU0L6ffvqJ5557rsjatJiXridcZ/OZzXbDavSo34PQFqF0rN1RzZjWjmnHSU5OVw6JienpV66kK4bff1e9C1uqVIHOnaFbN+jUCYwlaQFlVrIohvLloQhCsjiiKG4LIdpIKfcDCCHaAom57KPRaPJJUlISH3zwAW+//Ta3byv7/7333svChQtp165dkbSZkJJAdGI0J2+eZPXR1aw6sirbsBp1KtbRjum8YBnCevt2utkoNRUiImDHDqUc/vor634tWqgeQ7duEBiYrgBcXJRCsCiHYggJ74iimACsFkL8jQoxXhMYVKRSaTR3KDt27GDEiBEcP34cgOrVq/Pee+8xbNiwQg/eZxm9dCPhBlvObskxrEa3et3wLeerHdOOIKVSChblYPE33LyZrhh++w1uZVrFwccH7r9fmZQ6d4Zq1dLzMpuTitn3k6uikFLuEUI0BZoYSceM4awajaaQSUlJ4fjx47i6ujJmzBjCwsLw9fUt3DYM89Kp6FOsObqG8MjwbMNqNKraSDumHcHib4iPV2YjKdUchyNH0k1Khw+n+yEsNG6sFEO3bhAUlB5qw80toznJ1bk9N0enQ7YD/I3ybYQQSCm/LDKpNJo7hISEBG7evImfnx8APXr0YPr06Tz22GMEBgYWbluGeckS0ju7sBp97umjnNN6xnTOWPwNlvkNoHoJO3eqnsOOHXA902rS5cpBx45KOXTtCnffrdKFSFcK3t7g4VG8x5ILjsR6+goV4ykCNYQV1MxsrSg0mnwipWTNmjW88sorNGzYkC1btli/2F9//fVCbedW0i3OxZ5j1ZFVhEeGcyr6lDW/kmclHm/6OCGBIbSs0VI7pnMjs79BSjhxIr3XsH9/xqGtAP7+6SOU7r03XQl4eqYrhnLlit2clBcc+VwIBprL3MbRajQahzhy5Ajjxo1jy5YtAFy9epWjR48SEBBQaG2kpKUQnRjNb+d/Y3nk8ixhNYJqBBESGEL/pv2pUaGGdkxnh9mcPr/B4m9ISIBdu9Tchh074NKljPu4u0P79unKwd9fpbu6ZjQnFTCia3HiiKSRKAf2pdwKajSa7ImNjSUsLIwFCxaQZnx19uvXj7lz59KgQYNCaSMxJZGoW1GER4YTfiTcbliN0MBQOvh10I7p7MgcoltKOHtW9Rh27FCT3lIyuWnvvjvdnNShg1IEQqiegqXX4Fl6z7UjiqIacFQIsRuwxq6VUvYrMqk0mjLGV199xcSJE7l6VUXSb9SoEfPmzaNPnz4FrltKSVxyHLsu7OKLQ1/YDasREhhCSGAItX1qa8e0PSwhum/fVv6G5GTYvTtdOZw9m7G8qyu0aZOuHBo1UorBwyOjOckJy8wWBY4oirCiFkKjKevs3buXq1ev4u3tzb/+9S8mTJhgXXkuv6SaU7kcf5mVkSv5OvJr9l/ab83zcPWgzz19GNJyCF3rdcXXy1fPmLbFEjLDMow1NVWZkHbsUCalXbuyTnqrVk2Zkrp0UZPeKlZ0ypwGZ+DI8NjtmVa4K4+K3aTRaLLh2rVr+Pr64m68OKZOnYrJZOLNN9+kdu3aBarblGpi/9/7+SziM77585uMYTUq+RMSGMLgloOp71sfbw/vArVVpsgcojspSU16sziijbkrVoSAli2VYujaFQIClGKwzGnw9s4QIqMsk58V7mrj2Ap3Gs0dR2pqKh999BFvvvkmb731FhMmTADA19eXTz75JN/1Sim5mXiTNX+u4auDX7Hzwk5rnpuLGz0b9GRIiyH0vqc3vl6+2jFtwRKi2xIy49o1tVaDZdJbXKYVnStWzDjprWpV1Uuw7TWUEXNSXijKFe40mjuK7du3M3bsWA4fPgzAvHnzGDNmDG4FGN2Sak7lyNUjLNm/hJVHVmYJqzEwYCDDWw2ncbXGeLndGV+3uWIJ0W1RDpGRGSe9ZaZJEzXhrUsXNenNw0P5FyyKoYTNaXAGzl7hTqMp9URFRTFx4kRWrlwJqOB9zz//PO+8806+lcTt5Nt8/9f3fB7xOVvPbs0SVmNoq6H0b9ofXy9f7Zi2DdF9+zbcuKEmvVkc0ZnXiy5fXk16syiHmjXViCSLYijhcxqcQVGucKfRlGmSk5OZM2cO06dPJ8FwfHbo0IEFCxYQHJz3aM5SSk5Hn2bxvsUsP7w8S1iNJ5s/yYg2Iwi4KwAP1zv8K9fib7Aohz//THdEHziQcYlQUHMZunVTf23bZuwxeHs7PURGSacoV7jTaMo8S5cuJSEhgRo1ajBjxgyefvrpPAfvS0lLYf3J9SzZt4RfTv2SMaxG7Q48E/QMAwMGUtmr8p3de0hOTlcO16/DH3+km5SuXMlY1tNTzYK29Brq1SszcxqcgSOjnszAp8afRnNHExMTYw3S5+Hhwfz589mwYQNvvvkmlWzXCHCAv2/9zX8O/IelB5dyOjp9BeBKnpUY0HwAI9qMoG2ttnd2vKXExHTlcOxYumLYsyfrpLfatdPDcrdvr9ZssHVC38lKtoA4ssLdGez4JKSUhTOVNI/oFe40zuD27du88847LFiwgD179tC0adN81WM2m9l2dhsf7/uYH47/kCWsxvDWwxnSYghVylcpLNFLF7Yhum/eTO817NgB585lLOvmlj7prVs3NenNMmzV27tUhcgoDop6hTvbir2AJ1FDZTWaMo+UkpUrVzJx4kQuGgvZv/vuu3zxxRd5qifGFMPnBz5nyYElWcJq/F/T/+OFti/QqW4nXMSdN/SS1NR05WAbYG/XroyrwAFUr56+0tt996ntO2xOgzNwxPR0I1PSB0KIfcCbRSOSRlMyOHz4MGPHjmX79u0AeHp6MmXKFKZMmeJwHXsv7mXhnoWsObomQ1iNplWbMrz1cJ4NepZq3tVyqKGMYgnRHR2tlv609BpOnMhYTgg1ZNWiHFq2hAoV7ug5Dc7AkQl3bWw2XVA9DN2n05RZoqOjCQsLY9GiRdbgff379+f999+nfv36ue6fmJLIskPLWLxvMXsvpZtJPVw96Nu4LyODR9Ldv3uhr1hX4rGEzDhzBjZvVsph506lMGzx9U2f9Nali/I96DkNTsWRF/4cm9+pwFlgYJFIo9GUAM6dO8fChQsxm800btyY+fPn89BDD+W631/X/2LR7kUsP7ycaFO0Nb2+b32ebf0sI9qMoEaFGkUpesnCbFaK4dYt5WvYulUphyNHspZt3jw9VEb79mpZUIs5STuhnY4jpqcHikMQjcaZmEwmvAwbd1BQEBMnTqRatWqMHz8ejxy+YlPSUlhzdA0f7f2IX8//ak13c3Gj9z29ebHtizzS6JE7Z1irJUR3VBRs3KiUw6+/QkxMxnLe3srH0LUrdO8O9eun9xr0nIYShyOmp3/klC+lfL/wxNFoipcrV67w6quvcuDAAfbs2WOdST1jxowc9zsfc56FexbyxcEvuHr7qjW9lk8thrUaxqjgUfhV8itS2UsMSUkqZtLu3bBpk+o1RERknfTWsGH68NXOnZWJqXx5PaehFODoqKd2wFpjuy+wGziR7R4aTQknJSWFRYsWERYWRmxsLADh4eEMGTIk233M0sz6E+tZtGcRG05tyBBWo3v97oxuN5rHmjxW9gPySalGI12+DL/8km5Suno1YzlPT7WIT5cu8OCD0LRpqVj2U5MVRxSFH9BGShkHIIQIA36UUmb/RGk0JZgtW7Ywbtw4jhi28ipVqvDOO+8QGhpqt/zV21f5ZO8nLDmwhPOx563p1cpX45lWzzCm3Rj8K/sXh+jOw2xWJqWICPj5ZxUqY+9e+5PeLOakBx5QaziUsmU/NVlx5OrVAJJttpONNI2mVHH+/HkmTpzI6tWrARW878UXX2T69OlUrVo1Q1kpJdvPbmfhnoWsPbY2Q1iN++vez6jgUQxoPqBsx1xKSVGhMjZtgg0blHKIispYxt1dxU7q2hUeekgNX9VzGsocjiiKL4HdQojvjO3+QN5mG2k0JYBly5ZZlUSnTp1YsGABrVu3zlAmOjGapQeX8vHejzl+I30hG19PXwa3HMyYe8fQtFr+ZmWXCkwmOHoUfvxRDWHdtUv5IGy56y5lTurZU/3VrFmmlv3UZCXXEB5gnUvR2djcIaU8UKRS5YAO4aFxFCklaWlpVge1yWSiV69evPDCCwwePNg6EklKyZ6/97Bo9yJWHV2VIaxG27vb8lK7lwgJDKGcezmnHEeRIqUakbRlC/z0k+o1nDqVsYyLi5r01q2b6jXce6/qNZTRZT/LKkUdwgOgPHBLSvm5EKK6EKK+lPJMfhrUaIqD48ePM378eNq0acM777wDgJeXF9u3b7cqiPjkeL4+/DWL9izi0JVD1n293b0JDQzlpXtfIqhmkFPkL1LS0pQyWLdOmZR++03Nd7DF1ze919C7N9Spo81JdzCOBAV8CzXyqYmUsrEQohawWkrZqTgEzIzuUWhyIi4ujunTpzN37lxSUlLw8PDg1KlT+PmlD1WNvBrJh3s+ZNmhZcQlpy+FGVA9gFHtRjG05VB8PH2cIX7RYTKpkUnr1yufw9GjWcsEBCgHdJ8+0KmTmvSmzUllhqLuUTwOtAb2A0gp/xZCOPQUCSF6A/MAV2CJlPK9bMo9AawB2kkptRbQ5BkpJV9//TWTJ0/m77//BsDf35+5c+dSu3ZtTKkmvjn6DR/u/ZDfL/xu3c/T1ZMnmj/BmHZj6ODXoWxNjIuKUr2Gn39WSsLepLf774deveCRR6BBA21O0tjFEUWRLKWUQggJIITwdqRiIYQrsAh4EIgC9ggh1kopj2Yq5wOMx1iTW6PJKwcPHmTs2LH8+quaGe3l5cWrr77KpEmTiEqIYtLGSXwe8Tk3E9OXxLynyj2MCh7FsKBhVClXRoIhm80qwJ7FpHTwYNZJb40aqaGrvXur/xUrOkdWTanCEUWxSgjxCeArhBgBPItjixjdC5yUUp4GEEKEA48Bmfu8bwMzgEkOS63RGEgpGT16NL//rnoJT8+vauYAACAASURBVDzxBO/NfI9DSYfou7ovm89stpZ1c3GjX5N+jA4eTff63ctG7+HGDTVCaf16NfEt86Q3Ly8VKuPBB6FvX2jWTJuTNHnGkVhPs421sm8BjYE3pZQbHai7NnDBZjsKaG9bwBhNVUdK+aMQIltFIYR4AXgBoG7dug40rSnLSCmtL3khBB988AHDhw/ntXde41jFY3T+rjOX4y9by9epWIcRbUbwfJvnudvnbmeJXThIqXoKa9cqk9KePSq+ki116ijF0KePGqXkU8b8LZpix6FRT1LKjUKI/UAX4GZu5R1BCOECvA8Mc6D9xcBiUM7swmhfUzr5448/GDduHPPnz6djx46YpZkblW/Q4NUGPH3w6QxhNXrf05vR7UbT554+pTusxu3bKlSGxaRkLKBkxd1dRVx96CHo1w8CA3WvQVOoZKsohBDrgH9KKSOFEHejnNl7gYZCiMVSyg9yqfsiUMdm289Is+ADBALbjK/DmsBaIUQ/7dDWZOby5ctMmTKFL7/8EoDRY0Yz6INBLN6/mDMx6SO17/K+i+dbP8+ItiPw9/V3krSFwMmT8P33am7Dzp1ZJ73VrAk9eign9MMPq/WhNZoiIqceRX0pZaTxeziwUUo51HA+7wRyUxR7gEZCiPooBRECPGXJlFLGAtalvYQQ24CJWklobElJSWH+/PlMnTqVuDg1lLVyg8ocbnOYiC0R1nIP+D/AyOCR9G/av3SG1UhOVpPd1q5VvYeTJzPmu7ioUBm9eytfQ3CwDqynKTZyUhS20b56YDiwpZRxQgiz/V3SkVKmCiHGAL+ghsd+JqU8IoSYBuyVUq7NuQbNnc7GjRsZN24cf/31FwCu3q6kPZBGdJtocIHKXpV5ptUzjAweSZNqTZwsbT64eFEphh9/VEoi86S3KlVUr+Hhh5VyyBSPSqMpLnJSFBeEEGNRTug2wM8AQohygEODraWU64H1mdLsrrUtpezmSJ2aO4MrV67Qt29fkpKSQADBkPZAGpSH9rXbMyp4FAMDBpausBqpqSp20n//q0xK9lZ6a9UqvdfQsaP2NWhKBDkpiueAaUBPYJCU0jJbpwPweVELprkziU+OZ8XhFXy872OSOiSphXf7QIW6FRjSYggvBr9YusJqXLumegzr1qkZ0cbaF1YqVlSzoR99VCmHGjows6bkka2ikFJeBUbaSd8KbC1KoTR3FlJK5i2dx+y5s7k14BZx0gir0QVaDGrB6HajGdxicOkIq2E2w/79yqS0fr36nTlMTrNm6b2G++/Xs6E1JR69mojGaZhSTcz/cT7vvv4uMUeMDmsV8OzpycCAgYwKHlU6wmrExqphqxZH9LVrGfPLlVORVx99VP3puUCaUoZWFJpi5+TNk8zfPp9P536KaacJjKER5ZuXZ9yocUzsO5Gq5Uuw41ZKFVRv3Tr1t2tX1klv9eunO6G7ddPrQmtKNVpRaIqFlLQUfjj+Ax/t+YhN32+CjUC8yvOu4c0/3/4nrz33Gi4l1Xl7+7Zas2HdOmVSyrzSm4eHMiM9+qia29C4sXPk1GiKgFwVhRCiMfARUENKGSiEaAn0k1JOL3LpNKWeC7EXWLJ/CUsOLOHvuL9VpC9jrUR3T3f+MfkfhL0WhldJXOvg5EmlFNatgx07sk56q1UrvdfQvTtUqOAcOTWaIsaRHsWnqIB9nwBIKQ8JIb4GtKLQ2MUszWw4tYGP937MD8d/yBhWo29vjh85TutmrZkzZ07Jit2VlKQUgiXI3okTGfNdXdWQ1UceUX+BgXrSm+aOwBFFUV5KuTuTQzE1u8KaO5ert6/y2YHPWLzPCKthBvZB+fjyvPzmy4xoM4J6vvWI7x9PhZLy9X3hgprTsG6dMi1lnvRWvXr6CKX/b+++w6usssWPf1cICSAhlIjSe1HpEFCuZUQRCM2KIKIgKOUn19HrjOi9zMjYr3cQzTCiFLu0kVGcgOgoTQQkNAkCUgRBQUBIA0IK6/fHfpOThBAOgXNSWJ/n4eE9bzv77CfJOnvvd6/dvbtb+c2Yi4w/geKwiDQBstejuBPYH9BSmVJDVVm2ZxlT4qfw0ZaPyDjlTejfA5W/rEzqT6mckBPcWuVWGlRtAFC8QSIzE1audK2GuDhISMh7XMSlx+jTx2Vf7djRJr2Zi54/geL/4TK3thSRn4EfgXsDWipT4iWmJfLuxnd5Pf51th7emrM/Mj2SWqtqsfWrraR6o9VDhw6lXr16Z7pV4B086FJyx8W5x1fzT3qrWtVlXu3d2/1fs2bxlNOYEsqf9Sh2ATd7K9uFqGrK2a4xZZOqEv9LPH+P/zuzE2ZzIvNEzrHoy6JptK0RC6YtYGuqCxzR0dHExsbSpUuXM90yME6dgrVr3TjDggVuzYb8k95at3athpgYuPpqCLUHAI05E3+eenoEl7IjBZjqLTY0TlU/D3ThTMlwLP0YMxNmMnnNZDYc8GVsrRxWmcGtBzOy40je/MubTJkyBYCoqChefPFFhg0bFrzHXRMT3aS3BQvcmEP+ld4uuQRuvtkFhpgYqFs3OOUypgzw52vUA6r6qoj0AGoAQ4D3AAsUZdymXzcxJX4K7333HinpvoZk65qtGd1pNPe2uTcnrcZjjz3Gu+++y4gRI5gwYQJVAz3oq+rGFxYscF1K33wDWVl5z2ne3Ldew3XX2aQ3Y4rIn0CR/bhTDPCulyrcngkso9Iy0/jo+4+YvGYyK/etzNkfXi6cu666izGdxtCmehteeukllh9aTkxMDADNmjVjz549REVFnenW5y81Fb780tellH/SW3i4S7AXE+MGops2DVxZjLmI+BMo1orI50Aj4Elv4aKzrkdhSpcdR3bw+prXeXvj2xw54Vvttmm1poyOHs39be+nesXqfPTRRwz8r4H89NNPNGnShG7duuVMlgtIkPjhB19gWLrULfCTW/36vlZDt25QqdKFL4MxFzl/AsVwoB2wS1WPi0gN3Ip3ppTLyMpg/rb5TF4zmcW7fQmBQ0NC6deiH2M6jaFbo26ICN9//z13/+fdfPnllwCEhYVx9913o/kHic9XWpoLCNnBIf9Kb6GhLlVG9ljDlVfapDdjAqywNbNbqupWXJAAaGw9TmXD3qS9TImfwowNMziQeiBnf90qdXmww4M82OFBakXUAiApKYkJEyYQGxtLppf4rk+fPrzyyis0vVBdOz/95AsMX34Jx4/nPX7ZZb7A0L27rQ9tTJAV1qJ4DHgI+GsBxxToFpASmYA4padYuH0hk9dMZtHORXnSanRv0p2xncfSq2kvyoWU811z6hTXXHMNW7ZsAaBp06a8+uqrOeMSRZaR4Qafs4NDQZPeOnf2dSm1b2+T3owpRoUtXPSQ9/+NwSuOudAOHjvIG/FvMHXdVPYm783ZH1UpiuHthzOq0ygaVm1Y4LUhISGMGTOGJ554gvHjx/Poo48SXtQnhw4ccJPeFixwj7Hmn/RWrZpLlRET4ya9XXpp0d7HGHPBydn6mEXka2ApsBxYUdwT7jp16qTx8fHFWYQST1VZvHsxk7+dzKc/fOpLqwFcX/96RkeP5vYrbiesXFie6w4fPswzzzzDn//8Z6pXrw5AZmYmv/76K3Xq1Dm3QmRlQXy87/HVtWtPP6dtW1+roUsXm/RmTACJyFpV7VSUa/35zRwCXAfcAbwsIieB5ar6aFHe0ATO0RNHmbZuGlPXTWX7EV/m08jwSIa0GcLDnR+mRVSL067LzMxkypQpjB8/nsTERDIzM5k8eTIAoaGh/geJI0dcayEuzrUeDh/Oe7xyZTfG0Lu3az2ca/AxxhQLf1J4/CgiaUC69+9G4IpAF8z4R1VZuW8lk9dMZt6WeaRlpuUc61irI2OixzCo1SAqlq9Y4PXLli1j7NixfPfddwBERkZy1VVX+fvmsHGjb6xh5UqXPiO3li19rYZrr3UL/BhjShV/UnjsBA4DHwLTgbGqavMoilnKyRTe2vAWU9dNJeGgbzC4UvlKDLhyAGO7jKVDrQ5nvP7nn3/mD3/4AzNnzgRARBg+fDjPP/88lxY2PpCSAv/+ty84/PJL3uMVKrhJb717u0lvjRuf1+c0xhQ/f7qeXgOuBQYB7YGlIrJMVXcGtGSmQOt+Wcff1vyNOZvncCzDt3ZCy6iWjOw4kgfaP0CV8CqF3mPHjh20a9eOY97aC126dCE2Npbo6OjTT1b1TXqLi3ML+2Rk5D2nQQPfYj6/+51NejOmjPGn6+lV4FURqYybaPc0UBcoV9h15sI5ln6MD777gGnrp7HmlzU5+8PKhdGvRT/Gdh7LdfWvw995Lk2aNKFr165s3LiRl156ifvuuy9v8r4TJ9ykt+w1G378Me8NQkNd7qTsLqWWLW3SmzFlmD9dT3/FtSgqA98Af8I9AWUCbNOvm3g9/nVmJswkMS0xZ3/Dqg0Z3n44IzuO5NJLzv4Y6c6dO5k7dy7jxo0DXDfTjBkziIiIIDJ78tru3b5Ww+LFLljkdvnlLij07u2ysFYpvNVijCk7Cg0UXvK/DcD/quqvwSnSxe1ExgnmbJ7DtHXT+Hrv1zn7Q0NC6dmkJyM7jaR3s95+tR6OHTvGCy+8wMsvv0x6ejrt27enR48eANS97DL4+mtfcPAm1eUICXGT3rLXbGjXzloNxlykCg0Uqqoi8kdVfS9YBboYqSrbftvGG/Fv8MGmDzh0/FDOsdoRtRnSZgijO43OWUrUn/vNmTOHxx9/nH1ehtV69eq5x1dnzHCB4Ysv3MB0btWru8luffq4/2vUuGCf0RhTevkzmL1ORKJVdc3ZTzXnIi0zjU+2fsLUdVNZvHtxnrQaNzS8gZEdRnLHlXdQvlx5v++ZkJDA2LFjWbJkCQDhYWH8sXNnxiUnU+mee06/oG1bFxj69IHoaChnQ0/GmLz8CRRdgMEisgc4hlufQlW1TUBLVkad0lPsOrKL6eun8/6m99mX7FtToUbFGgxqPYhRHUdxVU0/5zLkMnPmTIYMGUKWt4BPv/LleSU9ncZf+7qwiIhwYwx9+7pJb7VqnfdnMsaUbf4Eih5FvbmI9ARexT0hNU1VX8x3/DFgBJAJHMKtprenqO9Xkh1PP87nuz5n2rppfL7z8zxpNa6uezVD2w5lcJvBVA6rfG43VoUNG+Bf/+LGf/6TSllZ1MJVes/sx1hbtnTjDH37wn/8B5T3v4VijDGFpRmvoqrJuLWyz5mIlAMmA92BfcAaEZmvqt/nOm090Mlb52I08L/A3UV5v5Io61QWe5P28tbGt/jguw/YedQ39SQyPJLbr7idkR1HEl0nmhA5h+yoycnwxRd8O2MGB7/5hj6J7omoy4EvgbYVKhB2ww0uMPTuDQ0bXtDPZYy5uBTWovgQ6AOsxaUVz/3IiwJnm3LbGdihqrsARGQW0B/ICRSqujjX+auAe/0ueQmlqqSmp7JszzLe2vAWcdvj8qTVaHdZO+5rex/3trnXr0dbvZvC1q0wfz7ExXHwm294MiuLGcClwA9A1QYNoFcvovv3hxtugIoFp+wwxphzVVia8T7e/42KeO86wN5cr/fhxjvOZDiwsKADIvIQbm0M6tevX8TiBFZ6Vjq/pPzCh5s+5MNNH7L50OacY5XKV6Jf836M6DCC6xpcd1rW1gIdPw5ffQWffgoLF8LevWQCf8dNZMlO0l27Vi1+nT6dqj172uOrxpiA8Cuvs4hUA5oBFbL3qeqyC1UIEbkX6ATcUNBxVX0TeBNcmvEL9b7n65SeIjU9lW/3fcvbG9/mk22fkJqemnO8eY3mDG49mPvb3k/dKnXzLApUoB9/hI8/dnMbli+HkydzDi0BxoaGkuCtMletalWefe45Ro4cSTl7UskYE0D+zMweATyCS9uxAbgaWMnZV7j7GaiX63Vdb1/++98M/Ddwg6qezH+8JErLTONg6kHmbpnLzE0zWbvft9ZCWLkwejbtybB2w7ip0U1EhEec+Ubp6bBkietSWrTo9PWhQ0IgOppHRZi0ahVkZiIiPPTQQzz77LNERUUF5gMaY0wu/rQoHgGigVWqeqOItASe9+O6NUAzEWmECxADgTwP8otIe+ANoKeqHjynkgdZ1qkskk8ms/HXjbz/3fvM2zKPo2lHc443jGzIwFYDua/tfTSu1pjw0DOsBLdvH3zyiWs1LF0Kx47lPV69unt8tXdvN7ehenWunj0bBg6ka9euxMbG0qHDmbPCGmPMheZPoEhT1TQRQUTCVXWriJy++k0+qpopIg8Di3CPx85Q1c0i8hcgXlXnAy/jckjN9VJS/KSq/Yr+cS4sVeV4xnEOHz/Mpz98yqyEWazYuyLneGhIKDc1uol729xLTNMYIitEnt69lJXlAsK//uUW9dm8mdO0bQu33AL9+6NduhD32We0adOG+t4qcwMGDCAiIoJevXr5nfjPGGMuFH8CxT4RqQp8DHwhIkcBv+Y6qOoCYEG+fX/KtX3zOZQ1aNKz0klKS2Lbb9uYlTCLOZvn5EmrUatyLQZcNYAhbYbQIqrF6XMffvnFBYaFC13XUmJi3uMRES4dd69e0L8/1K4NwPbt23mkXz8WLlzIXXfdxZw5cwCXxC8mJiaAn9gYY87MnzTjt3mbT4vIYiAS+CygpSoGp/QUKSdTOJp2lH/v+jczE2ayZPeSPGk1rm9wPfe0uoc+zftQo1INX/dSZiasWuWCw6JFbtW3/GuRN2/uWg19+7qFfXJNektNTeW5555j4sSJpKenA5CYmEhaWhoVKlTAGGOKU2ET7ioAo4CmwCZguqouDVbBguVExgmSTiaxJ3EPc7+fy+zNs/Ok1YiqFMWdV9zJPa3voVXNVkRWiHST4w4edOMMCxe6tNyHDuW9ccWKbhZ0z56u1dC06WnvrarMnj2bxx9/nJ9/duP8DRo04JVXXuHWW2+1biZjTIlQWIviHSADt/ZEL+BK3MB2qZc9MJ2YlsiKvSuYlTDrtLQaXep0YWCrgfRr0Y/LLrmMS0LCYe1a31jDunWuJZFbgwZuIDomxmVfveSSM5bh5MmT9OjRg6VLXeytUKECTzzxBH/84x+pZCvEGWNKkMICxZWq2hpARKYD3wanSIFzLP0YSSeT2J+yn39u/SezEmadllbj1pa3MqjVINpf3o6qqZmELVzqWg1ffXX6+tDly7s1G3r0gH79oE0bvye9hYeHU6dOHQBuu+02Jk6cSENLtWGMKYEKCxQ5X6+9J5iCUJwLLyMrg6STSSSlJbH+wHpmJsxkwfYFp6XVGNhqIP0a96T2T4lEzFxKyBfPwerVbq5Dbpdf7sYYevVy//ycy5CVlcUnn3zCbbfdltOl9PLLLzN06FC6d+9+wT6vMcZcaIUFirYikuxtC1DRe52dZrzEroWpqqSkp5B8MpnDxw8T90McMxNmnpZWo2+zvgyqH0PXbcepOm0F5ZdOcUuC5lauHHTo4OtS6twZwvxIwZHLypUrGTt2LGvXruX9999n8ODBANSuXZva3hNPxhhTUhWW66nU5YVIy0wj+WQyySeT2Xp4K7MTZvPxto/zptWo3ozBUTcxeFt56vx9A+GrRyH514euUQOuv951KfXqBXXrulnS5+jAgQOMGzeOd955J2ffhg0bcgKFMcaUBn7leirJsk5lkZKeQlJaEsknk1m0cxGzEmblTasREkafSm15aHskN3ywkwrbp5x+o9atfV1K117rBqKL2N2WkZFBbGwsEyZMIDnZNco6dOhAbGwsXbt2LdI9jTGmuJTaQHE84zhJaUmkpqeyO3E3szbPYt6WeSSm+Sa3NaYaw3+sxogFv1LzUL6VXKtUcQHh5pvdI6yNG0P4GdJunIOtW7dy++23s2XLFgCqV6/O888/z4gRIyx5nzGmVCp1gSLzVCa7ju7iRMYJFu9efHpaDRX67ruEMYtTuWnXUQRfPiZatHBrNdxyi5sZHRkJoRe2CurUqUNiYiIhISGMGjWKZ555hupeKg5jjCmNSl2gSMtMY+KK/2Puln9wMFdajXopwkNrlOHrlFqpbkxCK1WCa65xXUo9e7pAUalSkcYbzlietDTWrVuX06UUERHB22+/Tc2aNWnXrt0Fex9jjCkuovlTTZRwUluUkd62Qsx2GBUPvbZDOYVTDRsQcsPvoFs3uOkml401AKu9qSrz58/n0Ucf5dChQ2zbts2eYDLGlFgislZVOxXl2lLXogC4LBWGr4MH10H94+XJ6NQBvfNG6BlDSKtWULlynlxKF9q2bdt45JFHWLRoEQDly5fn66+/ZsCAAQF7T2OMKS6lrkXRpKLod9WiyOp6NSE3dqNSjz6E1LzMPaUU4MHilJQUnnnmGSZNmkRGhpuP2LNnTyZNmkSLFmfNvG6MMcXmompRRNRpBEv+TZXql7supSDNGI+Li+PBBx9k//79ADRq1IhJkybRt29fS95njCnTSl2gCK1anUvqNg76+1asWJH9+/dTsWJFnnrqKR5//HFLAW6MuSiUukARLL/99hspKSk5ifq6devGK6+8wu233079+vWLt3DGGBNEF+450TIiKyuLKVOm0Lx5c4YNG0buMZzf//73FiSMMRcdCxS5rFixgujoaEaPHs2RI0dYt24dO3fuPPuFxhhThlmgAPbv38+QIUO49tprWb9+PQBDhw7lhx9+oGkBK9MZY8zF5KIfo5g0aRLjx48n1ZvN3bFjR2JjY7nmmmuKuWTGGFMyXPQtit27d5OamkpUVBRTp07l22+/tSBhjDG5XHQtij179lC7dm3KezO3n376acLDwxk3bhzVqlUr5tIZY0zJc9G0KE6cOMHTTz9Ny5Yt+dvf/pazv2rVqrz00ksWJIwx5gzKfKBQVebNm8cVV1zBhAkTSEtLY9q0aWRlZRV30YwxplQo04Fiy5Yt3HLLLdxxxx3s2bOHsLAwnnzySVavXm2LCBljjJ/K5BjF8ePHGT9+PK+99hqZmZkAxMTEMGnSJJo1a1bMpTPGmNKlTAaK0NBQ4uLiyMzMpEmTJkyaNIk+ffoUd7GMMaZUKjOB4sCBA1x++eUAhIWFERsby5o1a3jssccseZ8xxpyHUj9GcfjwYUaOHEnDhg3ZsmVLzv7u3bvz1FNPWZAwxpjzVGoDRWZmJpMnT6Z58+a8+eabnDx5ktdee624i2WMMWVOQAOFiPQUkW0iskNExhVwPFxEZnvHV4tIQ3/uu3z5cjp27MjDDz/M0aNHqVq1Kq+99hqxsbEX+iMYY8xFL2BLoYpIOeAHoDuwD1gDDFLV73OdMwZoo6qjRGQgcJuq3l3YfWvUqKFHjhzJvp4HHniAF154gUsvvTQgn8MYY8qC81kKNZAtis7ADlXdparpwCygf75z+gPveNv/AG6Ss6wrmh0kOnfuzOrVq5k2bZoFCWOMCaBABoo6wN5cr/d5+wo8R1UzgSSgRv4bichDIhIvIvFVqlRhxowZrFy5kujo6AAV3RhjTLZSMZitqm+qaidV7dSsWTOGDRtGSEipKLoxxpR6gfxr+zNQL9frut6+As8RkVAgEvgtgGUyxhhzjgIZKNYAzUSkkYiEAQOB+fnOmQ/c723fCXylgRpdN8YYUyQBm5mtqpki8jCwCCgHzFDVzSLyFyBeVecD04H3RGQHcAQXTIwxxpQgAU3hoaoLgAX59v0p13YacFcgy2CMMeb82IiwMcaYQlmgMMYYUygLFMYYYwplgcIYY0yhApbrKVBEJAXYVtzlKCGigMPFXYgSwurCx+rCx+rCp4WqRhTlwtK4cNG2oia2KmtEJN7qwrG68LG68LG68BGR+KJea11PxhhjCmWBwhhjTKFKY6B4s7gLUIJYXfhYXfhYXfhYXfgUuS5K3WC2McaY4CqNLQpjjDFBZIHCGGNMoUpsoBCRniKyTUR2iMi4Ao6Hi8hs7/hqEWkY/FIGhx918ZiIfC8i34nIlyLSoDjKGQxnq4tc590hIioiZfbRSH/qQkQGeD8bm0Xkw2CXMVj8+B2pLyKLRWS993sSUxzlDDQRmSEiB0Uk4QzHRURe8+rpOxHp4NeNVbXE/cOlJd8JNAbCgI3AlfnOGQNM8bYHArOLu9zFWBc3ApW87dEXc11450UAy4BVQKfiLncx/lw0A9YD1bzXNYu73MVYF28Co73tK4HdxV3uANXF9UAHIOEMx2OAhYAAVwOr/blvSW1RdAZ2qOouVU0HZgH9853TH3jH2/4HcJOISBDLGCxnrQtVXayqx72Xq3CrCZZF/vxcADwDvASkBbNwQeZPXTwITFbVowCqejDIZQwWf+pCgSrediTwSxDLFzSqugy3ts+Z9AfeVWcVUFVEap3tviU1UNQB9uZ6vc/bV+A5qpoJJAE1glK64PKnLnIbjvvGUBadtS68pnQ9VY0LZsGKgT8/F82B5iKyQkRWiUjPoJUuuPypi6eBe0VkH26NnLHBKVqJc65/T4DSmcLDnIGI3At0Am4o7rIUBxEJASYCQ4u5KCVFKK776Xe4VuYyEWmtqonFWqriMQh4W1X/KiLX4FbWbKWqp4q7YKVBSW1R/AzUy/W6rrevwHNEJBTXnPwtKKULLn/qAhG5GfhvoJ+qngxS2YLtbHURAbQClojIblwf7PwyOqDtz8/FPmC+qmao6o/AD7jAUdb4UxfDgTkAqroSqIBLGHix8evvSX4lNVCsAZqJSCMRCcMNVs/Pd8584H5v+07gK/VGa8qYs9aFiLQH3sAFibLaDw1nqQtVTVLVKFVtqKoNceM1/VS1yMnQSjB/fkc+xrUmEJEoXFfUrmAWMkj8qYufgJsAROQKXKA4FNRSlgzzgfu8p5+uBpJUdf/ZLiqRXU+qmikiDwOLcE80zFDVzSLyFyBeVecD03HNxx24wZuBxVfiwPGzLl4GKgNzvfH8n1S1X7EVOkD8rIuLgp91sQi4RUS+B7KAP6hqmWt1+1kX/wVMFZFHcQPbQ8viF0sRmYn7chDljcf8GSgPoKpTCapQ/AAABLtJREFUcOMzMcAO4DgwzK/7lsG6MsYYcwGV1K4nY4wxJYQFCmOMMYWyQGGMMaZQFiiMMcYUygKFMcaYQlmgMEEjIlkiskFEEkRkrohUOsv5C0SkahHfa5qIXFm0koKIDPPKukFE0kVkk7f9YlHvWch77fPun+Blef2LiIR7x+qJyOxCrq0uIqMudJmMyc0ejzVBIyKpqlrZ2/4AWKuqE8/xHoL7uQ1a6gVvlncnVT1cwLFQL9fY+dx/H9BKVRNFpAowFUhV1eF+XNsU+IeqtjufMhhTGGtRmOKyHGgKICIfi8ha79v0Q9kniMhuEYkSkYbeWgPvAgnAEBGZ6J3ziIjs8rYbi8gKb3uJiHQSkXIi8rb3bX2TN+EKEWkiIp9577tcRFr6W3AReVZE3vXe620RCRWRiSLyrZfjf0Suc8fl2v+ns91bVZOBh4ABIhIpIk1FZIN3r9YissZr2XwnIo2BF4EW2a0dEakiIl+JyDrvnD7etU29Opju1fNCEangHWvuXbPRu65hUcpuyq4SOTPblG3icnP1Aj7zdj2gqkdEpCKwRkQ+KmAGcTPgflVdJSKXAw97+68DfhOROt72snzXtQPqqGor772zu7LeBEap6nYR6QL8Heh2Dh+jJXC9qqaJyBjgoKp29rqMVonI57i8U/WBLrj8/wtEpKuqflPYjVU1SUT24AJpUq5DY4D/U9XZ3vsIMA5omt2iEJHywK2qmiwiNYEVwL+861sAg1R1k4jMA27FpeSeCTytqp96wSNE3MI+51x2UzZZoDDBVDH72zGuRTHd2/5PEbnN266HCwr5A8UeL38+qnpARCqLSIR3/oe4BVuuA+blu24X0FhEYoE44HMRqQx0xZfyBCD8HD/LJ6qavd7FLcAVIpKdRibS+wy34ALiem9/ZVy+JX/+2Ba0tso3wP+IW8FwnqrukNOXYBHgRRG5FjgF1BOX5wncmg2bvO21QEMRqQZEqeqnANmfSUTOp+ymjLFAYYLpRP6+dBH5HXAzcI2qHheRJbiEbfkdy/f6G1yemm24oPMAcA0up08OVT0qIm2BHsAoYADweyDxPPv1c5dHgDGq+mXuE0SkH/Csqk7nHIhIJC4AbgdqZu9X1fdEZCXQG/hMRB7g9AV47sMFqg5eDqR9+Oozd1bhLAr//ZeilN2UTTZGYYpbJHDUCxItcanB/bEceBzX1bQetxzsSVXN3VWTnTU1RFU/Av4H9wc0GfhRRO7yzhEvmBTVImCM16WGiLTwutEWAcNF5BJvf91c3+4L5LWSXgfmeuXMfayxqu5Q1Vdx3UltgBRcevVskbhusEwR6c5ZFqXxVr87JCJ9vfeoIO5ptHMuuym7rEVhittnwCgR2YJrHazy87rluG/dy1Q1S0T2AlsLOK8O8Ja4RY0AnvT+Hwy8LiL/g8uuOQu31nJRvIHrz9/gdQUdBPqr6gIv+K3y9qcA9wCnPT0FLPfOCcF1nz1bwDn3iMggIAPXknjae1JqrYhswnWtTQQ+9V5/i2uVnM1g4A0ReQ5IB+44x7KbMs4ejzXGGFMo63oyxhhTKAsUxhhjCmWBwhhjTKEsUBhjjCmUBQpjjDGFskBhjDGmUBYojDHGFOr/Az7GdX6NqR1TAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for k in ['Hamming','Jukes-Cantor 69 (JC69)']:\n", " ax = regplot(x=pairwise_distances['Tree'], y=pairwise_distances[k], color=colors[k], label=k, scatter=False);\n", "regplot(x=pairwise_distances['Tree'], y=pairwise_distances['Tree'], color='black', label='True', scatter=False, line_kws={'linestyle':'--'})\n", "plt.title(\"Pairwise Sequence Distance vs. Pairwise Tree Distance\");\n", "plt.xlabel(\"Pairwise Tree Distance\");\n", "plt.ylabel(\"Pairwise Sequence Distance\");\n", "ax.legend(loc=\"best\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "**Author:** [Niema Moshiri](https://niema.net/), UC San Diego, October 2, 2018\n", "\n", "---" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }