project2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project 2: Sequence Alignment and Phylogeny"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "   <div class=\"alert alert-block alert-danger\">\n",
    "    <center>Due: <b>Wednesday, 21 September, 8:59pm</b>.</center> \n",
    "   </div>\n",
    "   \n",
    "   <div class=\"alert alert-block alert-warning\">\n",
    "   <center>\n",
    "       <b>Collaboration and Resource Policy</b>\n",
    "    </center>\n",
    "    For this assignment, you are encouraged to work with one other person. Your team must satisfy these constraints:\n",
    "    \n",
    "   1. You **did not work together on Project 1**.\n",
    "   2. You and your partner have a **total number of siblings that is divisible by two** (e.g., if you have one sibling, you need to find a partner with 1, 3, 5, or 7 siblings. If anyone has more than 7 siblings, they can partner with anyone!)\n",
    "    \n",
    "We expect most students will have the best learning experience on this assignment by working with a partner, but if you prefer to work alone it is permissible to do this assignment on your own.\n",
    "    \n",
    "You are encouraged to discuss these problems with anyone you want, including other students in the class. If you do discuss the specific questions in the assignment with anyone other than your assignment partner and the course staff, though, you should list them in the _External resources used_ section below.\n",
    "    \n",
    "You are welcome to use any resources you want for this assignment, other than ones that would defeat the purpose of the assignment. This means you should not look at answers or code from any other students in the class (other than your collaboration with your partner) or from previous offerings of this course, and if you find code that implements the problem you are being asked to do for the assignment, you should not use that code. \n",
    "\n",
    "You should document all external resource you use that are not part of the course materials in the _External resources used_ section below.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Team submitting this assignment:**  \n",
    "<div class=\"alert alert-success\">\n",
    "    <b><em>list each member of your team here, including both your name and UVA computing id</em></b>\n",
    "</div>\n",
    "\n",
    "**External resources used:** \n",
    "<div class=\"alert alert-success\">\n",
    "<em>It is not necessary to list the course materials, but if you used any other resources, including discussing problems with students not on your team, list them here.</em>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "<b>Submission</b>: Please submit the code you wrote to generate your answers for all parts using this form: https://forms.gle/gv144kv3KRo67uUX7. Your answers should be in the Jupyter Notebook, along with your code. Before submission, you should make a copy of your notebook file with the name uvaid1\\_uvaid2.ipynb (where uvaidn is each teammates UVA id) so the submitted file identifies you. You and your partner should submit a single file once together. Submission is due 8:59 pm on Wednesday, 21 September."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting Started"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install basic required packages, should be run only once. You may need to restart the jupyter python kernel (under the Kernel menu) after this. (You can execute this directly in the notebook but running the command below.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import blosum as bl\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "import utils\n",
    "from itertools import chain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Global Sequence Alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we provide the sequence alignment code from [Class 6](https://computingbiology.github.io/class6/). You are welcome to use and modify this code however you want in your solution, but should answer the questions below based on this provided code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simpleMatch(a, b):\n",
    "    return 1 if a == b else -1\n",
    "\n",
    "def distanceMatch(a, b):\n",
    "    return 0 if a == b else -1\n",
    "\n",
    "def linearGap(n):\n",
    "    return -1 * n\n",
    "\n",
    "def alignmentScore(s1, s2, gapPenalty, match):\n",
    "    if not s1 or not s2:\n",
    "        return gapPenalty(len(s1)) + gapPenalty(len(s2))\n",
    "    else:\n",
    "        return max(gapPenalty(1) + alignmentScore(s1, s2[1:], gapPenalty, match), \n",
    "                   gapPenalty(1) + alignmentScore(s1[1:], s2, gapPenalty, match),\n",
    "                   match(s1[0], s2[0]) + alignmentScore(s1[1:], s2[1:], gapPenalty, match)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def alignmentScoreDP(s1, s2, gapPenalty, match):\n",
    "    m = np.zeros((len(s1) + 1, len(s2) + 1))\n",
    "    m[0, 0] = 0\n",
    "    for i in range(1, len(s1) + 1):\n",
    "        m[i, 0] = gapPenalty(i)\n",
    "    for j in range(1, len(s2) + 1):\n",
    "        m[0, j] = gapPenalty(j)\n",
    "    for i in range(1, len(s1) + 1):\n",
    "        for j in range(1, len(s2) + 1):\n",
    "            m[i, j] = max(gapPenalty(1) + m[i, j - 1],  \n",
    "                          gapPenalty(1) + m[i - 1, j],    \n",
    "                          match(s1[i - 1], s2[j - 1]) + m[i - 1, j - 1]) \n",
    "    return m\n",
    "    \n",
    "def readAlignment(s1, s2, m, gapPenalty, match):\n",
    "    i = len(s1)\n",
    "    j = len(s2)\n",
    "    s1a = \"\"\n",
    "    s2a = \"\" \n",
    "    score = 0\n",
    "    while i > 0 or j > 0:\n",
    "        if i > 0 and j > 0 and m[i, j] == m[i - 1, j - 1] + match(s1[i - 1], s2[j - 1]):\n",
    "            i = i - 1\n",
    "            j = j - 1\n",
    "            score += match(s1[i], s2[j])\n",
    "            s1a = s1[i] + s1a\n",
    "            if s1[i] == s2[j]:\n",
    "                s2a = s2[j] + s2a\n",
    "            else:\n",
    "                s2a = s2[j].lower() + s2a\n",
    "        elif i > 0 and m[i, j] == m[i - 1, j] + gapPenalty(1):\n",
    "            i = i - 1\n",
    "            score += gapPenalty(1)\n",
    "            s1a = s1[i] + s1a\n",
    "            s2a = '-' + s2a\n",
    "        elif j > 0 and m[i, j] == m[i, j - 1] + gapPenalty(1):\n",
    "            j = j - 1\n",
    "            score += gapPenalty(1)\n",
    "            s1a = '-' + s1a\n",
    "            s2a = s2[j] + s2a\n",
    "        else:\n",
    "            assert False\n",
    "    return (s1a, s2a, score)\n",
    "\n",
    "def showAlignment(s1, s2, gapPenalty, match):\n",
    "    m = alignmentScoreDP(s1, s2, gapPenalty, match)\n",
    "    r = readAlignment(s1, s2, m, gapPenalty, match)\n",
    "    print (r[0] + \"\\n\" + r[1] + \"\\n\" + str(r[2]))\n",
    "    return (m, r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "G-ATT\n",
      "GCA-T\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "# Example\n",
    "r = showAlignment(\"GATT\", \"GCAT\", linearGap, simpleMatch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the version that supports affine gap penalties (from Class 6):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def alignmentScoreDPG(s1, s2, gapPenalty, match):\n",
    "    m = np.zeros((len(s1) + 1, len(s2) + 1))\n",
    "    m[0, 0] = 0\n",
    "    for i in range(1, len(s1) + 1):\n",
    "        m[i, 0] = gapPenalty(i)\n",
    "    for j in range(1, len(s2) + 1):\n",
    "        m[0, j] = gapPenalty(j)\n",
    "    for i in range(1, len(s1) + 1):\n",
    "        for j in range(1, len(s2) + 1):         \n",
    "            m[i, j] = max(chain((gapPenalty(g) + m[i, j - g] for g in range(1, j + 1)),\n",
    "                                (gapPenalty(g) + m[i - g, j] for g in range(1, i + 1)),   \n",
    "                                [(match(s1[i - 1], s2[j - 1]) + m[i - 1, j - 1])]))\n",
    "    return m\n",
    "    \n",
    "def readAlignmentG(s1, s2, m, gapPenalty, match):\n",
    "    i = len(s1)\n",
    "    j = len(s2)\n",
    "    s1a = \"\"\n",
    "    s2a = \"\"\n",
    "    score = 0\n",
    "    while i > 0 or j > 0:\n",
    "        if i > 0 and j > 0 and m[i, j] == m[i - 1, j - 1] + match(s1[i - 1], s2[j - 1]):\n",
    "            i = i - 1\n",
    "            j = j - 1\n",
    "            s1a = s1[i] + s1a\n",
    "            s2a = (s2[j] if s1[i] == s2[j] else s2[j].lower()) + s2a\n",
    "            score += match(s1[i], s2[j])\n",
    "        else:\n",
    "            foundit = False\n",
    "            for g in range(1, i + 1):\n",
    "                if m[i, j] == m[i - g, j] + gapPenalty(g):\n",
    "                    s1a = s1[i - g:i] + s1a\n",
    "                    s2a = ('-' * g) + s2a\n",
    "                    i = i - g\n",
    "                    score += gapPenalty(g)\n",
    "                    foundit = True\n",
    "                    break\n",
    "            if not foundit:\n",
    "                for g in range(1, j + 1):\n",
    "                    if m[i, j] == m[i, j - g] + gapPenalty(g):\n",
    "                        s1a = ('-' * g) + s1a\n",
    "                        s2a = s2[j - g:j] + s2a\n",
    "                        j = j - g\n",
    "                        score += gapPenalty(g)\n",
    "                        foundit = True\n",
    "                        break\n",
    "            assert foundit\n",
    "    return (s1a, s2a, score)\n",
    "\n",
    "def showAlignmentG(s1, s2, gapPenalty, match):\n",
    "    m = alignmentScoreDPG(s1, s2, gapPenalty, match)\n",
    "    r = readAlignmentG(s1, s2, m, gapPenalty, match)\n",
    "    print (r[0] + \"\\n\" + r[1] + \"\\n\" + str(r[2]))\n",
    "    return (m, r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def affineGap(n, gp = -1, gn = -0.2):\n",
    "    return gp + (n - 1) * gn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AAAGAATTCA\n",
      "AAA----TCA\n",
      "4.4\n"
     ]
    }
   ],
   "source": [
    "# Example\n",
    "s1 = \"AAAGAATTCA\"\n",
    "s2 = \"AAATCA\"\n",
    "r = showAlignmentG(s1, s2, affineGap, simpleMatch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 1 (a).** Run the given algorithm to find a global sequence alignment for the OCA2 genes (a key gene for the production of melanin) for humans and mice with the following gap penalties (still using simpleMatch as the match score function):\n",
    "\n",
    "  a. `linearGap` penalty\n",
    "    \n",
    "  b. `affineGap` penalty, with $gp=-0.2$\n",
    "\n",
    "  c. `affineGap` penalty, with $gp=-0.1$\n",
    "    \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "human_oca2, mouse_oca2 = utils.load_oca2_sequences()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 1 (b).** Use the given function to convert these sequences to their amino-acid sequences, and then re-run alignment for all sequences with the default parameters for `affineGap`.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'KCGV'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Convert sequence of nucleotides to amino acids using codon table lookup\n",
    "# Example\n",
    "utils.convert_to_amino(\"AAATGCGGCGTA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Alignment with Amino-Acids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The PAM<em>n</em> matrix (to be covered in [Class 6](https://computingbiology.github.io/class6/)) represents the likelihood of the occurrence of each tranformation during a time period where there are _n_ total mutation events per 100 amino acids."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 2 (a)** What would a negative value of an entry in a PAM 1 matrix $M$ indicate? Explain in terms of evolution and functionality of the proteins. Note that $M_{ij} = log(\\frac{q_{ij}}{p_ip_j})$ where $q_{ij}$ indicates the frequency of amino acids $i$ and $j$ observed to align in related sequences, and $p_i$ and $p_j$ represent the frequencies of occurrence of $i$ and $j$.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "**Problem 2 (b).** The BLOSUM<em>x</em> matices are created by clustering sequences with more than <em>x</em>% similarity into one single sequence and comparing sequences with more than <em>x</em>% divergence. Therefore, BLOSUM matrices are based on local alignments. Which of BLOSUM 50 and 60 contain more evoluationary divergence? \n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 2 (c).** Use the BLOSUM62 matrix as your scoring function to perform global alignment on the amino-acid sequences using `linearGap` (default parameters).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "blosum_matrix = bl.BLOSUM(62)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 2 (d).** How do your results for Problem 2c differ from the earlier ones of Problem 1a (with `linearGap`)? Which one would you say is more biologically plausible?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 2 (e).** We discussed in class that the PAM matrices follow the Markov property and a mismatch at any site depends only on the amino acid at that site and the transition probability. Is this a suitable representation of evolution? Think about if replacements are equaly likely to occur over entire sequences. It may help to consider the difference between PAM and BLOSUM matrices.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Local Sequence Alignment\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "<b>Problem 3 (a).</b> Implement local alignment (for both the normal and affine-gap penalties) using the Smith-Waterman algorithm. Feel free to re-use and modify the given Needleman–Wunsch algorithm. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def showAlignmentLocal(s1, s2, gapPenalty, match):\n",
    "    # Although it is often useful to return all high scoring local alignments for an input pair, \n",
    "    # it is sufficient if your algorithm just returns the single highest-scoring local alignment \n",
    "    # (as shown in the examples below).\n",
    "    \n",
    "    # Your code here (implement)\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've included some assert statements that can help you check the correctness of your algorithm. As with any algorithm, correctness on these test inputs does not guarantee algorithmic correcntess, but can be useful to debug."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GTTGAC\n",
      "GTT-AC\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "# Example expected output\n",
    "# Taken from https://en.wikipedia.org/wiki/Smith–Waterman_algorithm)\n",
    "r = showAlignmentLocal(\"GGTTGACTA\", \"TGTTACGG\", linearGap, simpleMatch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GTTGAC\n",
      "GTT-AC\n",
      "4\n",
      "GG\n",
      "GG\n",
      "2\n",
      "G\n",
      "G\n",
      "1\n",
      "TA-CGG\n",
      "TATCGG\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "# First assert\n",
    "r = showAlignmentLocal(\"GGTTGACTA\", \"TGTTACGG\", linearGap, simpleMatch)\n",
    "assert (r[1][2] == 4 and \"GTTGAC\" in r[1] and \"GTT-AC\" in r[1])\n",
    "\n",
    "# Second assert\n",
    "r = showAlignmentLocal(\"GGACTTAAATAGA\", \"TGTTGGTGATCCACGTGG\", linearGap, simpleMatch)\n",
    "assert (r[1][2] == 2 and \"GG\" == r[1][0] and \"GG\" == r[1][1])\n",
    "\n",
    "# Third assert\n",
    "r = showAlignmentLocal(\"TTGA\", \"GGCC\", linearGap, simpleMatch)\n",
    "assert (r[1][2] == 1 and \"G\" == r[1][0] and \"G\" == r[1][1])\n",
    "\n",
    "# Fourth assert\n",
    "r = showAlignmentLocal(\"TACGGGCCCGCTAC\", \"TAGCCCTATCGGTCA\", linearGap, simpleMatch)\n",
    "assert (r[1][2] == 4 and \"TA-CGG\" in r[1] and \"TATCGG\" in r[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "<b>Problem 3 (b).</b> Align the provided hemoglobin genes for:\n",
    "    <ol>\n",
    "        <li>`polar bears` & `black bears`,</li>\n",
    "        <li>`humans` & `chimps`,</li>\n",
    "        <li>`polar bears` & `humans`, and</li>\n",
    "        <li>`black bears` & `chimps`.</li>\n",
    "    </ol>\n",
    "\n",
    "Use `linearGap`.\n",
    "    \n",
    "Take note of the scores you get. What do you notice?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "polar_bear, black_bear, human, chimp = utils.get_hemoglobin_sequences()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "<b>Problem 3 (c).</b> Use <a href=\"https://blast.ncbi.nlm.nih.gov/\">BLAST</a> for the above pairs of sequences. Carefully inspect the returned results to see if they are similar to the alignments you obtained above.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "<b>Problem 3 (d).</b> Could you run an affine-gap-loss version of your local-alignment algorithm for the given sequences? How much time did BLAST take?\n",
    "Can you think of any optimizations you could make to make the affine-gap-loss version run faster- perhaps utilizing parallel processing or GPUs?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4: Phylogenetic Tree Reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this part, we'll briefly enter a fictional setup where you want to trace the evolution of Pokémon. The data is in the format of a two lists: one each for the sequences themselves, and names of the Pokémons."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "**Problem 4 (a).** Implement an algorithm for Phylogenetic Tree Reconstrution using the neighbor joining algorithm. Color intermediate nodes different from leaf nodes. Use given names as node labels in your visualization.\n",
    "    \n",
    "For computing the distances matrix, use affine-based gap-loss in your alignment score computations.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can either label intermediate nodes in the Phylogenetic tree such that they start with \"intermediate_\" and use the given functions below, or use your own nomenclature/way of handling those node, and modify the given helper functions accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've provided a helper function to plot a given Phylogenetic tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_alignment(dist, names):\n",
    "    # Your code here (implement)\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_graph_nice(G):\n",
    "    \"\"\"\n",
    "        Helper function to plot a given Phylogenetic tree.\n",
    "        Assumes intermediate node names start with 'intermediate_'\n",
    "    \"\"\"\n",
    "    nodes = list(G.nodes)\n",
    "    # Plot intermediate nodes smaller\n",
    "    sizes = [10 if \"intermediate_\" in x else 2000 for x in nodes]\n",
    "    labels = {}    \n",
    "    for node in nodes:\n",
    "        if not node.startswith(\"intermediate_\"):\n",
    "            labels[node] = node\n",
    "    fig, ax = plt.subplots(figsize=(15,15))\n",
    "    nx.draw_planar(G, node_size=sizes, with_labels=True, node_color = \"#ADD8E6\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the visualization for the given example on Wikipedia to get a sense of what the output should look like. We use `networkx` for creating and managing the graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wikipedia example: https://en.wikipedia.org/wiki/Neighbor_joining\n",
    "distances = np.array([\n",
    "    [0, 5, 9, 9, 8],\n",
    "    [5, 0, 10, 10, 9],\n",
    "    [9, 10, 0, 8, 7],\n",
    "    [9, 10, 8, 0, 3],\n",
    "    [8, 9, 7, 3, 0]\n",
    "], dtype=float)\n",
    "\n",
    "seq_names = [\"a\", \"b\", \"c\", \"d\", \"e\"]\n",
    "G = construct_alignment(distances, seq_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAAM9CAYAAAB5Rim2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACKXklEQVR4nOzdd3iUVf7+8Xtm0gskdJDOJJQkdKSJ9GJZ1l5WcVVW17au3VUJ1fWrYsFekLWtdXVdK4SiICJIJwVIofcWkpCezDy/P9T5ibFAmORMeb+uay+Kk2duQjbDPZ/znGOzLMsSAAAAAOCk2E0HAAAAAAB/RJkCAAAAgFqgTAEAAABALVCmAAAAAKAWKFMAAAAAUAuUKQAAAACoBcoUAAAAANQCZQoAAAAAaoEyBQAAAAC1QJkCAAAAgFqgTAEAAABALVCmAAAAAKAWKFMAAAAAUAuUKQAAAACoBcoUAAAAANQCZQoAAAAAaoEyBQAAAAC1QJkCAAAAgFqgTAEAAABALVCmAAAAAKAWKFMAAAAAUAuUKQAAAACoBcoUAAAAANQCZQoAAAAAaoEyBQAAAAC1QJkCAAAAgFqgTAEAAABALVCmAAAAAKAWKFMAAAAAUAuUKQAAAACoBcoUAAAAANQCZQoAAAAAaoEyBQAAAAC1QJkCAAAAgFqgTAEAAABALVCmAAAAAKAWKFMAAAAAUAuUKQAAAACoBcoUAAAAANQCZQoAAAAAaoEyBQAAAAC1QJkCAAAAgFqgTAEAAABALVCmAAAAAKAWKFMAAAAAUAuUKQAAAACoBcoUAAAAANQCZQoAAAAAaoEyBQAAAAC1QJkCAAAAgFqgTAEAAABALVCmAAAAAKAWKFMAAAAAUAshpgMAACBJVS63CiqqVFBepfyyKlW53XJbks0mOWw2xYaFqFFkqOLCQxUV6pDNZjMdGQAQ5GyWZVmmQwAAglNJVbW2HC3R7qJyVbjccthtcrstuX/l8SF2m75/2bKpcWSoEhrFqFlUGMUKAGAEZQoAUK8sy9KBkgrl5Bcrv7xKliXV9oUoxGaTw26TMz5a7eOiFO5g9ToAoP5QpgAA9aa4slrf7T2q4kqXXF58+bHbJLts6tG8gdo2iGRSBQCoF5QpAECdsyxLeUdLtPHwMbnq8FXHYbMpPiJU/VrFKTLEUXdPBACAKFMAgDpWUe3Wst1HdKyyuk6L1I9s+r5U9WsVp5YxEXX/hACAoEWZAgDUmbIqlxbvPKzyanet74uqLYdN6tm8odo1jKrnZwYABAvu1AUA1Inyape+MlSkJMllSesPFGpnYamBZwcABAPKFADA66rcbi3ZeUQVhorUj1yWtO5AofaXlBtMAQAIVJQpAIDXbThQpLJql9Ei9SOXJa3aW6CK6l87vQoAgNqhTAEAvOpgSYX2HCuT2xea1A9cbktr9xeYjgEACDCUKQCA11S53Vq1r6Bedu07GW5JB0srtecYy/0AAN5DmQIAeM3GQ8dU5fbN5XQu6/vplMuXRmYAAL9GmQIAeEW129L2Qt9a3vdzbkvafazMdAwAQICgTAEAvMIfSorLspSTX2I6BgAgQFCmAACnzLIs5RwplssPzoEvrarW0fIq0zEAAAGAMgUAOGWFFdUq85Otx92WtPUo0ykAwKmjTAEATtmRskrJJ06V+n2WpENllaZjAAACAGUKAHDKjpRV+tx26L+lrMrFrn4AgFNGmQIAnLJ8P7sHyWG3qbDCvzIDAHwPZQoAcEpcbktlVS6vXzf/wH49+re/6JqBybpxZH99/sYrXru2ZVkqoEwBAE4RZQoAcEoqXG7Zbd69ptvt1v/d+Ge179JNLy9Zqymvva/P3nhF65Yu9sr1XZbqpAACAIILZQoAcErcliWbvNum8jLWq+joEV1y8x0KDQtTizbtNPriK7Tsi/957TmquWUKAHCKQkwHAAD4t7roJIf27lb+wQOa0K+L5/fcLpe69u3vteew/OBMLACAb6NMAQBOicNm83qhatKilZq1bqvn0pZ5+cr/X4i31yYCAIIOy/wAAKckxG6T28tTHmf3XoqKjtFHs59VRXmZXC6XduZsVl7Geq9c3yYpzMFLIADg1PBKAgA4JWEOu0K9POVxOBy674XXtW1Tlm4aNUDXDEzW86l3qeRYkXeub7epYXioV64FAAheNotF4wCAU7R01xEdKq00HeOE2SWd5WyucKZTAIBTwKsIAOCUNY0M8/J+fnUrxGGjSAEAThmvJACAUxYfGSaHH23oEMcSPwCAF1CmAACnrGlUmOkIJyzEZlOHuCjTMQAAAYAyBQA4ZXabTR0bRvnFi4rNJrWMiTAdAwAQAPzhdQ8A4Ac6xkfL12+cstukjnHRstt8PCgAwC9QpgAAXhEV6lCTSN9f7teRJX4AAC+hTAEAvKZ7swby1X0oHLbvi1RkqMN0FABAgKBMAQC8pkF4qLo0ipHDB5fRhTscSmrSwHQMAEAAoUwBALwqsXGMonxs+mO3Sae3ivOr7dsBAL6PMgUA8Cq7zab+reJ8ZjrlsEnO+Gg18oP7uQAA/oUyBQDwugbhoRpwWrwchvuUwya1iI5QUpNYs0EAAAGJMgUAqBPNo8PVr5W5QuWw2dQsOlz9WsXJ5iNTMgBAYLFZlmWZDgEACFwHSyq0fM9RuS1L9fWC47DZdFpshPq0aEiRAgDUGcoUAKDOlVRWa+W+AhVVVMtVhy87Nn1fpHo1b6DWDSIpUgCAOkWZAgDUC8uytKWgVFmHjtVJobJZbjWNjlDflnGKCPGt3QQBAIGJe6YAAPXCZrPJGR+tUe2bqG2DCNltOuUd/76fREmVxwq08PUXNbh1I4oUAKDeMJkCABhR5XJrR2Gpco+WqNLlls1mU7X791+S7JLsdpvclqXWsZFyxkcrQi4lJibq/fff14ABA+o+PAAAokwBAAyzLEvFlS4drahSfmmlDpdXqqSyWq6fvDrZJIWH2BUfEaqmkeGKiwhVXESIQuz/f4HFyy+/rA8++EDz58+v/z8EACAoUaYAAD7rx5eoE9lIoqqqSp07d9Zrr72mM888s66jAQDAPVMAAN9ls9lOeEe+0NBQTZ48WZMmTRLvEwIA6gNlCgAQMK688kodOHBACxcuNB0FABAEKFMAgIAREhKiqVOnKjU1lekUAKDOUaYAAAHl0ksvVUlJib744gvTUQAAAY4yBQAIKHa7XdOmTWM6BQCoc5QpAEDAOf/88yVJH330keEkAIBAxtboAICA9Pnnn+vee+/Vhg0b5HA4TMcBAAQgJlMAgIB09tlnKyYmRu+//77pKACAAMVkCgAQsBYuXKibb75ZWVlZCgkJMR0HABBgmEwBAALWyJEj1aJFC/373/82HQUAEICYTAEAAtrXX3+tq6++WtnZ2QoNDTUdBwAQQJhMAQAC2plnnimn06lXX33VdBQAQIBhMgUACHjfffedLrroIuXm5ioiIsJ0HABAgGAyBQAIeP3791fPnj01e/Zs01EAAAGEyRQAICisW7dO55xzjvLy8hQVFWU6DgAgADCZAgAEhV69emngwIF6/vnnTUcBAAQIJlMAgKCRmZmpkSNHKi8vT7GxsabjAAD8HJMpAEDQSE5O1siRI/X000+bjgIACABMpgAAQSUnJ0eDBw9Wbm6u4uLiTMcBAPgxJlMAgKCSmJioc889V0888YTpKAAAP8dkCgAQdLZt26a+ffsqJydHjRs3Nh0HAOCnmEwBAIJOhw4ddPHFF2vmzJmmowAA/BiTKQBAUNq1a5d69OihTZs2qXnz5qbjAAD8EGUKABC0/v73v8tut+vJJ580HQUA4IcoUwCAoLV//35169ZN6enpat26tek4AAA/Q5kCAAS1u+++WyUlJXr++edNRwEA+BnKFAAgqB06dEhdunTRmjVr1L59e9NxAAB+hDIFAAh6kyZN0r59+zRnzhzTUQAAfoQyBQAIekePHlVCQoKWL1+uhIQE03EAAH6Cc6YAAEEvPj5ef//73zVt2jTTUQAAfoTJFAAAkoqKiuR0OrV48WJ169bNdBwAgB+gTAEA8INHH31Uq1ev1vvvv286CgDAD1CmAAD4QUlJiZxOp+bOnauePXuajgMA8HHcMwUAwA+io6N17733asqUKaajAAD8AJMpAAB+ory8XE6nUx999JH69etnOg4AwIcxmQIA4CciIiL0wAMPKDU11XQUAICPo0wBAPAzEydOVHZ2tr755hvTUQAAPowyBQDAz4SFhSk1NZXpFADgN1GmAAD4BVdddZV2796tL7/80nQUAICPokwBAPALQkJCNHXqVKWmpoq9mgAAv4QyBQDAr7jssstUUFCgefPmmY4CAPBBlCkAAH6Fw+HQtGnTmE4BAH4RZQoAgN9wwQUXqLq6Wh9//LHpKAAAH8OhvQAA/I5PPvlEkyZN0vr162W38z4kAOB7vCIAAPA7/vCHPygiIkIffPCB6SgAAB/CZAoAgBOQlpam2267TZmZmXI4HKbjAAB8AJMpAABOwJgxY9SkSRO9/fbbpqMAAHwEkykAAE7Q4sWLNXHiRG3evFmhoaGm4wAADGMyBQDACRo2bJjat2+v119/3XQUAIAPYDIFAMBJ+Pbbb3X55ZcrJydH4eHhpuMAAAxiMgUAwEkYNGiQkpKS9Morr5iOAgAwjMkUAAAnac2aNRo/frzy8vIUGRlpOg4AwBAmUwAAnKQ+ffro9NNP1wsvvGA6CgDAICZTAADUQnp6usaMGaO8vDzFxMSYjgMAMIDJFAAAtdC9e3cNGzZMzz77rOkoAABDmEwBAFBLmzdv1pAhQ5SXl6eGDRuajgMAqGdMpgAAqKUuXbrorLPO0qxZs0xHAQAYwGQKAIBTsGXLFvXv3185OTlq1KiR6TgAgHrEZAoAgFPQqVMnnX/++XrsscdMRwEA1DMmUwAAnKKdO3eqV69e2rRpk5o1a2Y6DgCgnlCmAADwgltuuUXh4eF6/PHHTUcBANQTyhQAAF6wd+9eJScnKzMzU61atTIdBwBQDyhTAAB4yZ133qnKyko988wzpqMAAOoBZQoAAC85ePCgunbtqnXr1qlt27am4wAA6hhlCgAAL7r//vt1+PBhvfzyy6ajAADqGGUKAAAvys/PV2Jior777jt16tTJdBwAQB3inCkAALyoUaNGuuWWWzR9+nTTUQAAdYzJFAAAXlZYWCin06mlS5eqS5cupuMAAOoIkykAALysYcOGuuOOOzR16lTTUQAAdYjJFAAAdaC4uFhOp1Pz589X9+7dTccBANQBJlMAANSBmJgY3XPPPZoyZYrpKACAOsJkCgCAOlJWVian06lPPvlEffr0MR0HAOBlTKYAAKgjkZGRuv/++zV58mTTUQAAdYDJFAAAdaiiokKJiYl69913NXDgQNNxAABexGQKAIA6FB4ertTUVKWmppqOAgDwMsoUAAB17M9//rO2b9+uxYsXm44CAPAiyhQAAHUsNDRUkydPVmpqqlhdDwCBgzIFAEA9uOKKK3T48GEtWLDAdBQAgJdQpgAAqAcOh0NTp07VpEmTmE4BQICgTAEAUE8uvvhilZeX67PPPjMdBQDgBZQpAADqid1u1/Tp05Wamiq32206DgDgFFGmAACoR3/84x8VEhKi//73v6ajAABOEYf2AgBQz+bOnau77rpL6enpcjgcpuMAAGqJyRQAAPVs3Lhxatiwod59913TUQAAp4AyBQAIOElJST59QK7NZtOMGTM0depUVVdX/+rjrr76ak2aNEmStHTpUnXu3Lm+IgIATgBlCgAQcLKysjRs2LDffVz79u21cOHCug/0C0aMGKHWrVvrjTfeOKHHDxkyRNnZ2Sf02Ndee01nnHHGqcSTJG3cuFF9+/ZVfHy84uPjNWrUKG3cuPGUrwsAgYIyBQBALViWdUo78tlsNk2fMUOL1mVp1Z58HSqt8GI672jVqpU++OAD5efn6/Dhwxo/frwuu+wy07EAwGdQpgAAAefHidPUqVN1ySWX6KqrrlJsbKySkpK0evVqSdKECRO0c+dO/eEPf1BMTIweffRRSdKKFSs0aNAgxcXFqUePHsctFxw2bJgeeOABDR48WFFRUdq6datsNpuef/55JSQkKDY2VqmpqdqyZYsGDhyoBg0a6JJLLlFlZaXnGp999pl69uypuLg43XrHneoxdJR2FVfo2935+nrFSvXu3VuxsbG69NJLVV5e7vm4xYsXq3Xr1p5fP/zww+rUqZNiY2PVrVs3ffTRR5KkTZs26YYbbtDy5csVExOjuLg4SVJFRYXuuusutW3bVs2bN9cNN9ygsrKy3/w8xsXFqX379rLZbLIsSw6HQ3l5eaf0dwMAgYQyBQAIaJ988okuu+wyFRQUaPz48brlllskSW+++abatm2rTz/9VMXFxbrnnnu0Z88enXPOOZo0aZLy8/P12GOP6cILL9ShQ4c813vzzTf18ssv69ixY2rXrp0kad68eVqzZo1WrFihRx99VNdff73eeust7dq1S5mZmXrnnXckSWvXrtW1116rl156SUeOHNGYi6/QY3+/XlWVFaqqrNTlF1+kCRMmKD8/XxdffLE+/PDDX/1zderUSUuXLlVhYaGmTJmiK6+8Uvv27VPXrl314osvauDAgSouLlZBQYEk6d5771VOTo7Wr1+vvLw87dmzR9OnTz+hz2FcXJwiIiL0t7/9Tffff39t/hoAICBRpgAAAe2MM87Q2WefLYfDoQkTJmjDhg2/+th///vfOvvss3X22WfLbrdr9OjR6tu3r7744gvPY66++molJSUpJCREoaGhkr4vKg0aNFBSUpKSk5M1ZswYdezYUQ0bNtRZZ52ldevWSZJmz56tv/71r+rfv7/mz5+v/Qf2KyQ0TNnr1mjjmlVyu6p12223KTQ0VBdddJH69ev3q1kvvvhitWrVSna7XZdeeqkSEhK0cuXKX3ysZVmaPXu2nnzySTVq1EixsbG6//77T3g3wYKCAhUWFurZZ59Vr169TuhjACAYhJgOAABAXWrRooXn51FRUSovL1d1dbVCQmq+BO7YsUP/+c9/9Omnn3p+r6qqSsOHD/f8uk2bNjU+rnnz5p6fR0ZG1vj1/v37Pdd//fXX9cgjj6i6ulpRUVFyud2q2L9T7732ilq1bCmbzeb52B8nX7/kjTfe0BNPPKHt27dLkoqLi3X48OFffOyhQ4dUWlqqPn36eH7Psiy5XK5fvf7PRUdH64YbblDTpk21adMmNWvW7IQ/FgACFWUKABC0flpcpO+L0oQJEzR79uwT/pgTZVmWqqqq5HA4dP311+vBBx9UTEyM578v++JjrV27VpZleZ5j586d6tSpU41r7dixQ9ddd50WLVqkgQMHyuFwqGfPnrIs6xczNmnSRJGRkcrKytJpp51Wq/yS5Ha7VVpaqj179lCmAEAs8wMABLHmzZtr69atnl9feeWV+vTTT5WWliaXy6Xy8nItXrxYu3fvPqXnKS0t1QUXXKC8vDxFRUXp8ssvV3R0tEpKSvT555/r2LFjmjVrlkpKSvToo4+qurpa//3vf3912V5JSYlsNpuaNm0qSXr11VeVmZl53J9r9+7dno0v7Ha7rrvuOt1+++06ePCgJGnPnj1KS0v7zdwLFizQunXr5HK5VFRUpDvuuEPx8fHq2rXrKX0+ACBQUKYAAEHrvvvu04MPPqi4uDg99thjatOmjT7++GM99NBDatq0qdq0aaOZM2fWegt0y7K0fv16ffTRR0pOTtbmzZv1+uuv65ZbblF8fLycTqdee+01SVLPnj01ZswYzZo1S/Hx8Xrvvfd0wQUX/OJ1u3XrpjvvvFMDBw5U8+bNlZGRocGDB3v++4gRI5SUlKQWLVqoSZMmkqRHHnlETqdTAwYMUIMGDTRq1KjfPbeqoKBAl19+uRo2bKhOnTopLy9P8+bNU0RERK0+HwAQaGzWj2sCAACA12zdulXXX3+9CgoKNGfOHPXo0eN3PyY3N1cDBw5Ubm6u4uPj6yElAOBUMJkCAMCLXC6XZs2apdNPP11jx47VihUrTqhISVJCQoL++Mc/6oknnqjjlAAAb2AyBQCAl2zcuFETJ05UWFiYZs+ercTExJO+xvbt29WnTx9lZ2d7lujVpYceekgPPfRQjd8fMmSI5s6dW+fPDwD+jDIFAMApqqys1COPPKKnn35aM2bM0PXXXy+7vfaLP2666SbFxMTo0Ucf9WJKAIC3UaYAADgFq1ev1sSJE3XaaafppZde+sVzqE7Wnj17lJKSoo0bNx53ThYAwLdQpgAAqIXS0lJNnTpVr7/+up544gn96U9/qvUZVL/ktttuk2VZeuqpp7x2TQCAd1GmAAA4SUuWLNFf/vIX9e3bV0899VSdHGC7f/9+JSUlaf369V6ZdgEAvI8yBQDACSoqKtK9996rTz/9VM8//7zGjx9fp8937733qrCwUC+++GKdPg8AoHbYGh0AgBPw+eefKzk5WS6XS5mZmXVepCTp7rvv1n/+8x9t27atzp8LAHDymEwBAPAbDh8+rNtuu03Lly/X7NmzNWLEiHp9/smTJ2vXrl169dVX6/V5AQC/j8kUAAC/wLIsvfvuu0pOTlbz5s2Vnp5e70VKku644w599tlnysnJqffnBgD8NiZTAAD8zJ49e3TTTTcpLy9P//rXv9S/f3+jef75z39q48aNeuutt4zmAAAcj8kUAAA/sCxLs2fPVs+ePdWrVy+tXbvWeJGSpFtvvVULFy5UZmam6SgAgJ9gMgUAgKQtW7bouuuuU3FxsebMmaOUlBTTkY4zc+ZMfffdd/rggw9MRwEA/IDJFAAgqLlcLj3xxBPq37+/zjnnHC1fvtznipQk3Xzzzfr222+1bt0601EAAD9gMgUACFqZmZmaOHGioqKiNHv2bDmdTtORftPTTz+tBQsW6NNPPzUdBQAgJlMAgCBUWVmpadOmafjw4Zo4caIWLVrk80VKkq6//nqtX79e3333nekoAABRpgAAQWblypXq06ePVq9erXXr1un666+X3e4fL4cRERGaNGmSUlNTTUcBAIgyBQAIEqWlpbrzzjs1fvx4PfDAA/rkk0/UunVr07FO2jXXXKPc3FwtXbrUdBQACHqUKQBAwPvqq6+UkpKi/fv3KzMzU5dddplsNpvpWLUSFhamKVOmaNKkSeK2ZwAwiw0oAAABq7CwUHfffbfmzp2rF154Qeeee67pSF5RXV2tpKQkPffccxo1apTpOAAQtJhMAQAC0qeffqrk5GTZ7XZlZmYGTJGSpJCQEE2dOlWpqalMpwDAICZTAICAcujQId16661atWqVXnnlFQ0bNsx0pDrhdrvVvXt3PfLIIzrnnHNMxwGAoMRkCgAQECzL0ltvvaXk5GS1bt1a6enpAVukJMlut2vatGmaPHky0ykAMITJFADA7+3atUs33nijdu7cqTlz5qhfv36mI9ULt9utvn37KjU1Veeff77pOAAQdJhMAQD8ltvt1osvvqjevXurf//+Wr16ddAUKen76dT06dM1efJkud1u03EAIOgwmQIA+KXc3Fxdd911Ki8v15w5c5SUlGQ6khGWZWngwIG67bbbdNlll5mOAwBBhckUAMCvVFdXa+bMmRo4cKDOO+88LVu2LGiLlCTZbDbNmDFDU6ZMUXV1tek4ABBUKFMAAL+xYcMGDRgwQGlpaVq5cqVuu+02ORwO07GMGzVqlJo3b6633nrLdBQACCos8wMA+LyKigo9+OCDevHFF/Xwww/r2muvlc1mMx3Lp3z99de6+uqrlZ2drdDQUNNxACAoMJkCAPi05cuXq1evXsrIyNCGDRs0ceJEitQvOPPMM9WpUye9+uqrpqMAQNBgMgUA8EklJSV64IEH9N577+mpp57SxRdfTIn6HStWrNAll1yinJwcRUREmI4DAAGPyRQAwOcsXLhQKSkpys/PV2Zmpi655BKK1AkYMGCAunfvrtmzZ5uOAgBBgckUAMBnFBQU6M4779TChQv14osv6qyzzjIdye+sXbtW5557rvLy8hQVFWU6DgAENCZTAACf8L///U9JSUmKiIhQRkYGRaqWevfurYEDB+qFF14wHQUAAh6TKQCAUQcOHNDf/vY3rV+/Xq+88orOPPNM05H8XmZmpkaOHKm8vDzFxsaajgMAAYvJFADACMuy9MYbb6h79+7q2LGjNmzYQJHykuTkZI0YMULPPPOM6SgAENCYTAEA6t3OnTv117/+Vfv27dOcOXPUp08f05ECTnZ2ts444wzl5uYqLi7OdBwACEhMpgAA9cbtduu5555T7969NWTIEK1atYoiVUc6d+6sc845R08++aTpKAAQsJhMAQDqRXZ2tv7yl7/I5XJpzpw56tq1q+lIAW/r1q06/fTTlZ2drcaNG5uOAwABh8kUAKBOVVVV6eGHH9bgwYN1ySWXaOnSpRSpetKxY0ddeOGFmjlzpukoABCQmEwBAOrMunXrNHHiRDVp0kQvv/yy2rdvbzpS0Nm1a5d69OihTZs2qXnz5qbjAEBAYTIFAPC68vJyPfDAAxo7dqxuvfVWpaWlUaQMadOmja688ko98sgjpqMAQMBhMgUA8Kply5Zp4sSJSkpK0rPPPquWLVuajhT09u3bp6SkJGVkZOi0004zHQcAAgZlCgDgFcXFxbr//vv1wQcf6JlnntGFF15oOhJ+4u6771Zpaamee+4501EAIGCwzA8AcMrmz5+v5ORkFRUVKTMzkyLlg+655x69++672rFjh+koABAwmEwBAGotPz9fd955p7766iu99NJLGjt2rOlI+A0PPPCADhw4oFdeecV0FAAICEymAAC18uGHHyo5OVkxMTHKyMigSPmBu+66S//73/+Ul5dnOgoABAQmUwCAk7Jv3z7dcsstysrK0iuvvKIzzjjDdCSchOnTpys3N1dvvvmm6SgA4PeYTAEATohlWXrttdfUo0cPdenSRevXr6dI+aHbbrtNaWlp2rhxo+koAOD3mEwBAH7X9u3bdf311+vw4cP617/+pZ49e5qOhFPwyCOPaM2aNXr//fdNRwEAv8ZkCgDwq9xut5555hn17dtXI0aM0HfffUeRCgC33HKLvv76a23YsMF0FADwa0ymAAC/aNOmTfrLX/4iu92uV155RZ07dzYdCV40a9YsLV68WP/73/9MRwEAv8VkCgBwnKqqKv3zn//UkCFD9Kc//UlLliyhSAWgG264QatXr9aqVatMRwEAv8VkCgDgsXbtWl177bVq0aKFXnrpJbVr1850JNSh559/Xp9++qnmzp1rOgoA+CUmUwAAlZWV6R//+IfOOuss3XHHHZo7dy5FKghMnDhRmzZt0rJly0xHAQC/RJkCgCC3dOlS9ezZU1u3blV6erquuuoq2Ww207FQD8LDw5WamqrU1FTTUQDAL7HMDwCCVFFRke677z7973//07PPPqvzzz/fdCQYUFVVpW7duumll17SiBEjTMcBAL/CZAoAgtDcuXOVkpKi8vJyZWZmUqSCWGhoqKZMmaLU1FTx/ioAnBwmUwAQRI4cOaLbb79dS5cu1ezZszVq1CjTkeADXC6XUlJS9MQTT2jcuHGm4wCA32AyBQBBwLIsvf/++0pOTlajRo2UkZFBkYKHw+HQtGnTmE4BwEliMgUAAW7v3r26+eablZ2drTlz5mjgwIGmI8EHud1u9erVS9OnT9cf//hH03EAwC8wmQKAAGVZlubMmaOePXsqJSVF69ato0jhV9ntds2YMUOTJ0+W2+02HQcA/AKTKQAIQFu3btX111+vgoICzZkzRz169DAdCX7Asiz1799fd911ly655BLTcQDA5zGZAoAA4nK5NGvWLJ1++ukaO3asVqxYQZHCCbPZbJo+fbqmTp0ql8tlOg4A+LwQ0wEAAN6RlZWliRMnKjw8XN9++60SExNNR4IfGjt2rBo1aqR33nlHV155pek4AODTmEwBgJ+rrKzUjBkzNGzYMF199dX66quvKFKoNZvNphkzZmjq1KmqqqoyHQcAfBplCgD82KpVq9S3b1+tWLFCa9eu1Q033CC7nW/tODXDhw9X27Zt9cYbb5iOAgA+jQ0oAMAPlZaWaurUqXr99df1xBNP6E9/+pNsNpvpWAggy5Yt05/+9Cfl5OQoPDzcdBwA8Em8fQkAfmbx4sXq0aOHdu3apYyMDF1xxRUUKXjd4MGD1a1bN82ZM8d0FADwWUymAMBPFBYW6t5779Vnn32m559/XuPHjzcdCQFu9erVOu+885Sbm6vIyEjTcQDA5zCZAgA/8Pnnnys5OVlut1uZmZkUKdSLvn37qm/fvnrxxRdNRwEAn8RkCgB82KFDh3TbbbdpxYoVmj17tkaMGGE6EoJMenq6xowZoy1btig6Otp0HADwKUymAMAHWZald999VykpKWrRooXS09MpUjCie/fuGjp0qJ599lnTUQDA5zCZAgAfs2fPHt14443aunWr5syZo/79+5uOhCC3adMmDR06VHl5eWrQoIHpOADgM5hMAYCPcLvdevnll9WzZ0/17t1ba9eupUjBJ3Tt2lXjxo3TrFmzTEcBAJ/CZAoAfEBeXp6uu+46lZSUaM6cOUpJSTEdCThOXl6eBgwYoJycHDVq1Mh0HADwCUymAMAgl8ulxx9/XAMGDNC5556r5cuXU6Tgk5xOp8477zw9/vjjpqMAgM9gMgUAhmRmZuraa69VdHS0Zs+eLafTaToS8Jt27Nih3r17a/PmzWratKnpOABgHJMpAKhnlZWVmjp1qoYPH66//OUvWrRoEUUKfqFdu3a67LLL9Mgjj5iOAgA+gckUANSjlStX6tprr1WHDh30wgsvqHXr1qYjASdl7969Sk5OVlZWllq2bGk6DgAYRZkCgHpQUlKiyZMn66233tKsWbN06aWXymazmY4F1Modd9yhqqoqPfPMM6ajAIBRlCkAqGNffvmlrrvuOg0YMEBPPfWUmjRpYjoScEoOHjyoLl26aP369Wrbtq3pOABgDGUKQMCwLEsuy5Lbkuw2mxw2GZ3+FBQU6O6779a8efP0wgsv6NxzzzWWBfC2++67T/n5+XrppZdMRwEAYyhTAPxSlcut/PIqHS2r1KGyShWWV6nSbckmyWaTLEuyJEWE2BUfEaomkWGKjwhVfESYHPa6L1iffPKJbrrpJp177rl65JFH1LBhwzp/TqA+HTlyRImJiVq1apU6duxoOg4AGEGZAuBXjpZXKTe/WHuLy2W32eRyWzqRb2J22/fTKllS+7hIdYqLVnRYiNfzHTx4ULfeeqtWr16tV155RcOGDfP6cwC+YurUqdq+fbtee+0101EAwAjKFAC/sPdYubIOH1NpVbXcP0ydauvH6VWjiFAlN22gRpFhp5zPsiy9/fbbuuOOO/TnP/9ZU6dOVVRU1ClfF/BlhYWFcjqdWrp0qbp06WI6DgDUO8oUAJ9WXu3S2v2FOlRaIVcdfLdy2KT2DaOU3LRBrZf/7dq1SzfccIN27dqlf/3rX+rbt6+XUwK+66GHHlJGRobeeecd01EAoN5xaC8An7X7WJnmbzukAyV1U6QkyWVJ2wtLNX/bQR0pqzypj3W73XrhhRfUq1cvDRgwQKtXr6ZIIejceuut+uqrr5SRkWE6CgDUOyZTAHyOZVlaf6BIO4vK5KrHb1EOm9S9WQN1iIv+3cfm5ubqL3/5iyoqKjRnzhwlJSXVQ0LANz3++ONatmyZ/vvf/5qOAgD1iskUAJ9iWZZW7SvQzqLSei1S0vdTqvSDRco5Uvyrj6murtajjz6qgQMH6vzzz9eyZcsoUgh6N954o1asWKE1a9aYjgIA9YrJFACfYVmW1h4o1O6isjpb1nciHDYppWkDdYw/fkK1YcMGTZw4UXFxcXr55ZfZDhr4iWeffVZz587V559/bjoKANQbJlMAfEbu0RLjRUr6fkKVcahIB0oqJEkVFRVKTU3VqFGjdOONN2rBggUUKeBnrrvuOmVkZGj58uWmowBAvWEyBcAnFFVU6asdh40XqZ8Ks9sUf2S7rrv2WiUmJur5559Xq1atTMcCfNbs2bP1/vvva8GCBaajAEC9oEwBMM6yLC3aflhFldWmoxzHVV2lNV8tUJ/mDXTRRRfJZqvd1ulAsKiqqlKXLl30r3/9S0OHDjUdBwDqHGUKgHHZR4q1+UhxvW84cSLskga0jleL6AjTUQC/8Prrr2vOnDlasmQJb0AACHjcMwXAqCq3W5uPHPPJIiVJbknrDxSJ952AE3PFFVfo4MGDWrhwoekoAFDnKFMAjNpZWGY6wu+qqHYrv7zKdAzAL4SEhGjq1KmaNGkSb0IACHiUKQDGWJalnPxin9p04pe4fsgJ4MRccsklKi0tZZt0AAGPMgXAmCNllar09Sb1gwMlFSqvdpmOAfgFu92u6dOnKzU1VW6323QcAKgzlCkAxuwsKvPZe6V+ziab9hVXmI4B+I3zzjtPdrtdH330kekoAFBnKFMAjDlSVlmn13/mH7fp7VmPeOVaLsvSkTLKFHCibDabZsyYoSlTpsjlYqoLIDBRpgAY4bYsFVf61z+wjpSxCQVwMs466yzFxsbqvffeMx0FAOoEZQqAEccqq+XwszNoSqtccvvJskTAF9hsNj344IOaOnWqqqt961BuAPAGyhQAIworqiUvd6mtGzN01wVjdEXvBD1++19VVVnu1es7bDYdq+QfhMDJGDFihFq1aqU333zTdBQA8DrKFAAjql1ur055qior9cgt12ro+Iv0+ncbNWjsH7Ri/hdeu74kySZV+8nug4Cv+PHeqenTp6uysm7vkwSA+kaZAmCEy7IkL/aSnA1r5aqq1rl/vk4hoaEaOO5cOZN7eO8JfuAvuw8CvmTIkCFKTEzUq6++ajoKAHgVZQqAETYv3y919OB+NWre4rjrNmnV2qvPIUl+dpsX4DNmzJihBx98UOXl3l1+CwAmUaYAGOGw2bxaTOKbNlf+gf2yfjI5Orxvj/ee4Af+tmkG4CtOP/109ezZUy+//LLpKADgNZQpAEZEhti9Op1K7NlHjhCHPn9zjlzV1Vox/wvlZaz32vWl77dzjwhxePWaQDCZPn26/u///k+lpaWmowCAV1CmABgRFxHq1Q0oQsPCdPfTc7T4o/d01eldtWzux+o/+iyvXV+S7DabIkP4tgnUVq9evTR48GA999xzpqMAgFfYLIu7qQGY8WnuflW5/edbUKPIUA1r28R0DMCvZWVlacSIEcrLy1NsbKzpOABwSniLFYAxDcNDTUc4KU0jw0xHAPxeUlKSRo0apaeeesp0FAA4ZUymABiTfaRYm44ckz8Mp0JsNvU/LV7No8NNRwH8Xk5OjgYNGqS8vDzFxcWZjgMAtcZkCoAx7RpGmo5wwux2qVkUkynAGxITEzV+/Hg98cQTpqMAwClhMgXAqOW787WvpMJ0jN9kt0ldGseoS2Pu7wC8Zdu2berbt6+ys7PVpAn3IgLwT0ymABiV0CjGL85u6tAwynQEIKB06NBBl1xyiWbOnGk6CgDUGpMpAEZZlqUF2w+puNJlOsovsklqFRuh/q3iTUcBAs7u3bvVvXt3bdy4US1atDAdBwBOGpMpAEbZbDad3jJedh8dTjnsNvVs1sB0DCAgtW7dWhMmTNDDDz9sOgoA1AqTKQA+IfNQkbYcLZHLh74jOWw29W3ZUKfF+s9GGYC/2b9/v5KSkrRhwwa1bt3adBwAOClMpgD4hG5NYhUR4jAdw8MmqWlUGEUKqGMtWrTQxIkT9c9//tN0FAA4aUymAPiMwooqLdlxRNWGvy3ZJIWH2DWyfVOFO3jPCahrhw8fVufOnbV69Wp16NDBdBwAOGH8KwGAz2gYHqrBbRoZ390v1GHTsLZNKFJAPWnSpIluuukmzZgxw3QUADgpTKYA+JzDpZVatjtfrnr+9mSTFO6wa2jbxooOC6nX5waC3dGjR5WQkKDly5crISHBdBwAOCGUKQA+qaC8Sst256vK7Za7Hr5LOWxSTFiIBrVupEgfuncLCCYPPvigNm/erH//+9+mowDACaFMAfBZ1W630g8WaVdRWZ3u8uewSV2bxCohPlo2PzhAGAhUx44dU6dOnfTVV18pKSnJdBwA+F2UKQA+73Bphb7bW6Bqt+XVpX8hNpuiwhzq3ypesSzrA3zCzJkztXLlSv3nP/8xHQUAfhdlCoBfcLkt7S0uV05+sYorq+W2pNp887JLkk1qHBmmxEYxahYVxjQK8CGlpaXq1KmT5s6dq549e5qOAwC/iTIFwO8Ullcp92iJDpZWqKLaLYfdJrfbkvsXHuuwSTabTS63pchQh06LiVCn+GhFhXJfFOCrnnrqKS1atEiffPKJ6SgA8JsoUwD8WpXbrcLyKh0tr1JRRbWq3ZbcliWH3aYQu01xEaGKCw9Vw/BQOexMoAB/UF5eroSEBH344Yc6/fTTTccBgF9FmQIAAD7nxRdf1EcffaS0tDTTUQDgV3EiJQAA8DnXXnutcnJytHTpUtNRAOBXUaYAAIDPCQsL0+TJk5WamioW0QDwVZQpAADgkyZMmKC9e/fqyy+/NB0FAH4RZQoAAPikkJAQTZ06lekUAJ9FmQIAAD7r0ksvVWFhoebOnWs6CgDUQJkCAAA+y+FwaNq0aZo8eTLTKQA+hzIFAAB82gUXXCCXy6WPP/7YdBQAOA7nTAEAAJ/36aef6v7779eGDRtkt/NeMADfwHcjAADg884991xFRUXpP//5j+koAODBZAoAAPiF+fPn69Zbb1VmZqZCQkJMxwEAJlMAAMA/jB49Wk2bNtXbb79tOgoASGIyBQAA/MiSJUt07bXXavPmzQoNDTUdB0CQYzIFAAD8xtChQ9WhQwe99tprpqMAAJMpAADgX5YvX65LL71Uubm5Cg8PNx0HQBBjMgUAAPzKwIEDlZKSotmzZ5uOAiDIMZkCAAB+Z82aNRo/frzy8vIUGRlpOg6AIMVkCgAA+J0+ffqof//+euGFF0xHARDEmEwBAAC/lJGRoVGjRmnLli2KiYkxHQdAEGIyBQAA/FJKSoqGDx+uZ555xnQUAEGKyRQAAPBbmzdv1pAhQ5SXl6eGDRuajgMgyDCZAgAAfqtLly46++yz9eSTT5qOAiAIMZkCAAB+bcuWLerfv7+ys7PVuHFj03EABBEmUwAAwK916tRJF1xwgR577DHTUQAEGSZTAADA7+3cuVO9evXSpk2b1KxZM9NxAAQJyhQAAAgIf/vb3xQWFqbHH3/cdBQAQYIyBQAAAsK+ffuUlJSkzMxMtWrVynQcAEGAMgUAAALGXXfdpfLycj377LOmowAIApQpAAAQMA4dOqQuXbpo7dq1ateunek4AAIcZQoAAASU+++/X4cOHdLs2bNNRwEQ4ChTAAAgoOTn5ysxMVErVqyQ0+k0HQdAAOOcKQAAEFAaNWqkv/3tb5o+fbrpKAACHJMpAAAQcAoLC5WQkKAlS5aoa9eupuMACFBMpgAAQMBp2LCh7rjjDk2dOtV0FAABjMkUAAAISMXFxXI6nZo/f766d+9uOg6AAESZAgAAAevJJ5/U119/rY8++sh0FAABiDIFAAACVllZmZxOpz7++GP17dvXdBwAAYZ7pgAAQMCKjIzU/fffr8mTJ5uOAiAAMZkCAAABraKiQomJiXrnnXc0aNAg03EABBAmUwAAIKCFh4crNTVVqamppqMACDCUKQAAEPD+/Oc/a8eOHfrqq69MRwEQQChTAAAg4IWGhmrKlClKTU0VdzgA8BbKFAAACAp/+tOfdOTIEc2fP990FAABgjIFAACCgsPh0LRp05hOAfAayhQAAAgaF110kSoqKvTpp5+ajgIgAFCmAABA0LDb7Zo+fbomT54st9ttOg4AP0eZAgAAQWX8+PEKDQ3Vhx9+aDoKAD/Hob0AACDozJ07V3feeacyMjLkcDhMxwHgp5hMAQCAoDNu3DjFx8fr3XffNR0FgB9jMgUAAILSl19+qb/+9a/atGmTQkJCTMcB4IeYTAEAgKA0YsQItWnTRm+88YbpKAD8FJMpAAAQtJYtW6YrrrhCOTk5CgsLMx0HgJ9hMgUAAILW4MGD1aVLF82ZM8d0FAB+iMkUAAAIaqtWrdL555+v3NxcRUZGmo4DwI8wmQIAAEGtX79+6tOnj1566SXTUQD4GSZTAAAg6G3YsEHjxo1TXl6eoqOjTccB4CeYTAEAgKDXo0cPDRkyRM8++6zpKAD8CJMpAAAASRs3btSwYcOUl5enBg0amI4DwA8wmQIAAJDUrVs3jR07Vk899ZTpKAD8BJMpAACAH+Tl5WnAgAHKzc1VfHy86TgAfByTKQAAgB84nU798Y9/1OOPP246CgA/wGQKAADgJ7Zv364+ffpo8+bNatq0qek4AHwYZQoAAOBnbrrpJkVHR2vmzJmmowDwYZQpAACAn9mzZ49SUlKUlZWlli1bmo4DwEdRpgAAAH7B7bffLpfLpaefftp0FAA+ijIFAADwCw4cOKCuXbtqw4YNatOmjek4AHwQZQoAAOBX/OMf/1BBQYFefPFF01EA+CDKFAAAwK84cuSIEhMTtWrVKnXs2NF0HAA+hnOmAAAAfkXjxo11yy23aMaMGaajAPBBTKYAAAB+Q0FBgRISEvTNN9+oc+fOpuMA8CFMpgAAAH5DXFycbrvtNk2bNs10FAA+hskUAADA7zh27JicTqcWLVqk5ORk03EA+AgmUwAAAL8jNjZWd999t6ZMmWI6CgAfwmQKAADgBJSWlsrpdOqzzz5T7969TccB4AOYTAEAAJyAqKgo3XfffZo8ebLpKAB8BJMpAACAE1ReXq7ExES9//77GjBggOk4AAyjTAEAAJyEl19+WR988IHmz59vOgoAw1jmBwAAcBKuueYa5eXl6euvvzYdBYBhlCkAAICTEBoaqsmTJ2vSpEligQ8Q3ChTAAAAJ+nKK6/UgQMHtHDhQtNRABhEmQIAADhJISEhmjp1qlJTU5lOAUGMMgUAAFALl156qUpKSvTFF1+YjgLAEMoUAABALdjtdk2bNo3pFBDEKFMAAAC1dP7550uSPvroI8NJAJhAmQIAAD4rKSlJixcvNh3jV9lsNs2YMUOTJ0+Wy+X61cddffXVmjRpkiRp6dKl6ty5c31FBFCHKFMAAMBnZWVladiwYb/7uPbt2xvbWe/ss89WTEyM3n///RN6/JAhQ5SdnX1Cj33ttdd0xhlnnEo8SdKKFSs0evRoNWrUSE2bNtXFF1+sffv2nfJ1gWBHmQIAAEHNsiy53e5af7zNZtODDz6oqVOnqrq62ovJvOfo0aO6/vrrtX37du3YsUOxsbG65pprTMcC/B5lCgAA+KwfJ05Tp07VJZdcoquuukqxsbFKSkrS6tWrJUkTJkzQzp079Yc//EExMTF69NFHJX0/jRk0aJDi4uLUo0eP45YLDhs2TA888IAGDx6sqKgobd26VTabTc8//7wSEhIUGxur1NRUbdmyRQMHDlSDBg10ySWXqLKy0nONzz77TD179lRcXJwmT56ssVddr49zD2jhtkNavmq1evfurdjYWF166aUqLy/3fNzixYvVunVrz68ffvhhderUSbGxserWrZvn/qtNmzbphhtu0PLlyxUTE6O4uDhJUkVFhe666y61bdtWzZs31w033KCysrLf/DyeddZZuvjii9WgQQNFRUXplltu0bJly07p7wYAZQoAAPiJTz75RJdddpkKCgo0fvx43XLLLZKkN998U23bttWnn36q4uJi3XPPPdqzZ4/OOeccTZo0Sfn5+Xrsscd04YUX6tChQ57rvfnmm3r55Zd17NgxtWvXTpI0b948rVmzRitWrNCjjz6q66+/Xm+99ZZ27dqlzMxMvfPOO5KktWvX6tprr9VLL72kI0eO6KKrrtG7z89SdXWVjhSX6oLzz9eECROUn5+viy++WB9++OGv/rk6deqkpUuXqrCwUFOmTNGVV16pffv2qWvXrnrxxRc1cOBAFRcXq6CgQJJ07733KicnR+vXr1deXp727Nmj6dOnn9Tn8uuvv1ZSUtJJfQyAmihTAADAL5xxxhk6++yz5XA4NGHCBG3YsOFXH/vvf/9bZ599ts4++2zZ7XaNHj1affv2Pe5MqKuvvlpJSUkKCQlRaGiopO+LSoMGDZSUlKTk5GSNGTNGHTt2VMOGDXXWWWdp3bp1kqTZs2dr4sSJKiws1D333KP1m7IVEhqqnPVrlbN+jaqqqnTbbbcpNDRUF110kfr16/erWS+++GK1atVKdrtdl156qRISErRy5cpffKxlWZo9e7aefPJJNWrUSLGxsbr//vv17rvvnvDnMT09XdOnT9fMmTNP+GMA/LIQ0wEAAABORIsWLTw/j4qKUnl5uaqrqxUSUvOfMzt27NB//vMfffrpp57fq6qq0vDhwz2/btOmTY2Pa968uefnkZGRx/06IiJC2dnZmjVrlj744AMdPnxYM2fOVGhoqEJCQlRV7dLRg/tUWVmp1qedJpvN5vnYHydfv+SNN97QE088oe3bt0uSiouLdfjw4V987KFDh1RaWqo+ffp4fs+yrN/cSfCn8vLydNZZZ+mpp57SkCFDTuhjAPw6yhQAAPB7Py0u0vdFacKECZo9e/YJf8wvKS0t1Ycffqi0tDS99957crlcio+PV69evXT66afrwQcf9DzWbVkqqXLp7NEjtW1LnizL8jzHzp071alTpxrX37Fjh6677jotWrRIAwcOlMPhUM+ePT2HAP88Y5MmTRQZGamsrCyddtppv5v/5881atQopaamasKECSf1sQB+Gcv8AACA32vevLm2bt3q+fWVV16pTz/9VGlpaXK5XCovL9fixYu1e/fu37yOy+XSypUrNWPGDKWnp+uOO+7Q7Nmz1bVrV11++eW69NJLNXv2bD300EN6/fXX9d1338myLJWUlGjuF19IFWV6fOZMFRcX67HHHlN1dbX++9///uqyvZKSEtlsNjVt2lSS9OqrryozM/O4P9fu3bs9G1/Y7XZdd911uv3223Xw4EFJ0p49e5SWlvabf649e/ZoxIgRuvnmm3XDDTf8/icUwAmhTAEAAL9333336cEHH1RcXJwee+wxtWnTRh9//LEeeughNW3aVG3atNHMmTN/cQv0vXv36rXXXpMkDRgwQNdcc42OHj2qtm3b6umnn9a8efN0++23q0mTJp5JUd++fTV79mzdcsstio+Pl9Pp9Fzj9NNP19ChQzVr1izFx8frvffe0wUXXPCLubt166Y777xTAwcOVPPmzZWRkaHBgwd7/vuIESOUlJSkFi1aqEmTJpKkRx55RE6nUwMGDFCDBg00atSo3z236pVXXtHWrVs1bdo0xcTEeP4H4NTYrB/nyAAAAEGgoqJC33zzjebNm6e0tDTt3r1bo0aN0tixYzVmzJhfvJfqZGVmZmrkyJHKy8tTbGysF1ID8EWUKQAAENAsy1Jubq6nPC1dulTdunXTuHHjNHbsWPXr1+8XN7E4VX/605+UlJSkBx54wOvXBuAbKFMAACDgFBUVadGiRUpLS1NaWpoqKys95WnUqFFq1KhRnWfIycnR4MGDlZub6zlwty499NBDeuihh2r8/pAhQzR37tw6f34gGFGmAACA33O73Vq7dq2nPK1bt04DBw7U2LFjNXbsWCUlJZ3Q7n3eds0116hNmzYnfaguAP9AmQIAAH5p//79mj9/vtLS0rRgwQI1btzYU56GDh2qqKgo0xG1bds29e3bVzk5OWrcuLHpOAC8jDIFAAD8QmVlpZYtW+aZPm3fvl0jRozwFKjfOhjXpBtuuEFxcXF6+OGHTUcB4GWUKQAA4LPy8vI85WnJkiXq3Lmzpzz1799foaGhpiP+rl27dqlHjx7atGmTmjdvbjoOAC+iTAEAAJ9x7NgxffXVV54CVVJS4ilPo0eP9py15G9uvfVWORwOPfnkk6ajAPAiyhQAADDG7XZrw4YNnvK0evVqnX766Ro7dqzGjRunlJQUIxtHeNu+ffuUlJSk9PR0tW7d2nQcAF5CmQIAAPXq4MGDWrBggadANWzY0FOehg0bpujoaNMR68Tdd9+tkpISPf/886ajAPASyhQAAKhTVVVVWr58udLS0jRv3jzl5eVp+PDhnnOfOnToYDpivTh06JC6dOmiNWvWqH379qbjAPACyhQAAPC6bdu2ecrT4sWL1alTJ095GjhwoF9sHFEXJk2apH379mnOnDmmowDwAsoUAAA4ZSUlJVq8eLHmzZuntLQ0FRYWHrdxRLNmzUxH9AlHjx5VQkKCVqxYIafTaToOgFNEmQIAACfNsixlZGR4pk8rV65U3759PQWqR48estvtpmP6pBkzZignJ0dvvvmm6SgAThFlCgAAnJDDhw97No6YP3++IiMjPUv3hg8frtjYWNMR/UJRUZGcTqcWL16sbt26mY4D4BRQpgAAwC+qrq7WihUrPLvuZWdna+jQoZ7pE8vUau/RRx/V6tWr9f7775uOAuAUUKYAAIDHjh07POXpyy+/VPv27T3lafDgwQoLCzMdMSCUlJTI6XRq7ty56tmzp+k4AGqJMgUAQBArLS3VkiVLPAXqyJEjGj16tMaOHasxY8aoRYsWpiMGrFmzZumrr77Sxx9/bDoKgFqiTAEAEEQsy1JWVpanPC1fvly9evXyHJrbq1cvNo6oJ+Xl5XI6nfroo4/Ur18/03EA1AJlCgCAAJefn6+FCxd6ClRoaKhn6d6IESPUsGFD0xGD1gsvvKCPP/5Y8+bNMx0FQC1QpgAACDAul0srV670lKesrCwNGTLEM31KSEiQzWYzHROSKisr1blzZ7355ps644wzTMcBcJIoUwAABIDdu3d7ytPChQvVunVrz7blZ5xxhsLDw01HxK/417/+pTfffFNfffWV6SgAThJlCgAAP1RWVqalS5d6Ds09cOCARo0apXHjxmnMmDFq1aqV6Yg4QdXV1eratateeukljRgxwnQcACeBMgUAgB+wLEubN2/2lKdly5ape/funulTnz595HA4TMdELb311lt6/vnn9c0337AEE/AjlCkAAHxUQUGBFi1a5ClQkjzlaeTIkYqLizMbEF7jcrnUvXt3PfbYYzrrrLNMxwFwgihTAAD4CJfLpTVr1mjevHlKS0tTenq6zjjjDM/Oe126dGFqEcA++OADPfzww1q1ahV/z4CfoEwBAGDQ3r17NX/+fM2bN08LFy5UixYtPOVpyJAhioyMNB0R9cTtdqt3796aOnWqzjvvPNNxAJwAyhQAAPWooqJC33zzjWf6tHv3bo0aNUpjx47VmDFj1KZNG9MRYdAnn3yiSZMmaf369RyeDPgByhQAAHXIsizl5uZ6ytPSpUuVlJTkmT7169dPISEhpmPCR1iWpf79++uuu+7SJZdcYjoOgN9BmQIAwMuKioo8G0ekpaWpqqrKU55GjRqlRo0amY4IH5aWlqbbbrtNmZmZ7NAI+DjKFAAAp8jtdmvt2rWe8rRu3ToNHDhQY8eO1bhx49StWzc2FMAJsyxLQ4YM0V//+ldNmDDBdBwAv4EyBQBALezfv1/z589XWlqaFixYoMaNG3umT0OHDlVUVJTpiPBjixcv1sSJE7V582aFhoaajgPgV1CmAAA4AZWVlVq2bJln+rR9+3aNGDFC48aN05gxY9SuXTvTERFgRo4cqcsvv1x/+ctfTEcB8CsoUwAA/Iq8vDxPeVqyZIk6d+7sOTS3f//+bByBOvXtt9/q8ssvV05OjsLDw03HAfALKFMAAPzg2LFj+uqrrzwFqqSkxFOeRo0apSZNmpiOiCBz9tln65xzztHNN99sOgqAX0CZAgAELbfbrQ0bNnjK0+rVq3X66ad7ClRKSgobR8Co1atX649//KPy8vI4wBnwQZQpAEBQOXjwoBYsWKB58+Zp/vz5iouL82wcMWzYMEVHR5uOCBzn/PPP15lnnqnbb7/ddBQAP0OZAgAEtKqqKi1fvlxpaWmaN2+etmzZouHDh3sKVIcOHUxHBH5Tenq6xowZo7y8PMXExJiOA+AnKFMAgICzbds2T3lavHixnE6npzwNHDiQrabhdy677DL17NlT//jHP0xHAfATlCkAgN8rKSnR4sWLNW/ePKWlpamoqEhjxozR2LFjNXr0aDVr1sx0ROCUbN68WUOGDFFeXp4aNmxoOg6AH1CmAAB+x7IsZWRkeMrTypUr1bdvX8/0qUePHrLb7aZjAl511VVXqVOnTpoyZYrpKAB+QJkCAPiFw4cPa8GCBUpLS9P8+fMVFRXlKU/Dhw9XbGys6YhAndqyZYv69++vnJwcNWrUyHQcAKJMAQB8VHV1tVasWOHZtjw7O1tDhw71FCin02k6IlDvrrvuOjVt2lQPPfSQ6SgARJkCAPiQHTt2eMrTl19+qfbt22vs2LEaN26cBg0apLCwMNMRAaN27typXr16adOmTdwLCPgAyhQAwJjS0lItWbLEU6COHDmi0aNHa9y4cRo9erRatGhhOiLgc2655RaFh4fr8ccfNx0FCHqUKQBAvbEsS1lZWZ7ytHz5cvXq1Uvjxo3T2LFj1atXLzaOAH7H3r17lZycrMzMTLVq1cp0HCCoUaYAAHUqPz9fCxcu9BSo0NBQT3kaMWKEGjRoYDoi4HfuvPNOVVZW6plnnjEdBQhqlCkAgFe5XC6tXLnSU56ysrJ05plnejaOSEhIkM1mMx0T8GsHDx5U165dtW7dOrVt29Z0HCBoUaYAAKds9+7dnvK0cOFCtWnTxlOezjjjDIWHh5uOCASc++67T0eOHNHLL79sOgoQtChTAICTVlZWpqVLlyotLU3z5s3TgQMHNHr0aI0dO1ZjxozhPg6gHuTn5ysxMVHfffedOnXqZDoOEJQoUwCA32VZljZv3uwpT8uWLVOPHj0806c+ffrI4XCYjgkEnalTp2rbtm16/fXXTUcBghJlCgDwiwoKCrRo0SLNmzdPaWlpstlsnvI0cuRIxcXFmY4IBL3CwkI5nU4tXbpUXbp0MR0HCDqUKQCApO83jlizZo2nPGVkZGjw4MGeAtWlSxc2jgB80P/93/9pw4YNevfdd01HAYIOZQoAgtjevXuP2ziiRYsWGjt2rMaNG6chQ4YoIiLCdEQAv6O4uFhOp1MLFixQSkqK6ThAUKFMAUAQqaio0DfffOOZPu3evVujRo3SuHHjNGbMGLVu3dp0RAC18MQTT+ibb77Rf//7X9NRgKBCmQKAAGZZlnJzcz3laenSpUpKSvJMn/r168fGEUAAKCsrk9Pp1CeffKI+ffqYjgMEDcoUAASYoqIiLVq0yLN8r6qqylOeRo4cqUaNGpmOCKAOPPfcc/riiy/0+eefm44CBA3KFAD4ObfbrbVr13rK07p16zRo0CDPxhHdunVj4wggCFRUVCgxMVHvvvuuBg4caDoOEBQoUwDgh/bv36/58+crLS1N8+fPV9OmTT3l6cwzz1RUVJTpiAAMmD17tt577z0tXLjQdBQgKFCmAMAPVFZWatmyZZ7p0/bt2zVy5EhPgWrbtq3piAB8QFVVlbp06aI5c+Zo2LBhpuMAAY8yBQA+Ki8vz1OelixZoi5dunjKU//+/RUSEmI6IgAf9MYbb2j27Nn6+uuvWeIL1DHKFAD4iGPHjumrr75SWlqa5s2bp7KyMk95GjVqlJo0aWI6IgA/4HK5lJycrKeeekpjxowxHQcIaJQpADDE7XZrw4YNnunT6tWr1b9/f0+BSklJ4V1lALXy3nvv6fHHH9d3333H9xGgDlGmAKAeHTx4UAsWLNC8efM0f/58xcXFebYtHzp0qKKjo01HBBAA3G63evbsqX/+85/6wx/+YDoOELAoUwBQh6qqqrR8+XLPoblbtmzR8OHDPdOnDh06mI4IIED973//07Rp07RmzRrZ7XbTcYCARJkCAC/btm2b576nxYsXy+l0eqZPAwYMUGhoqOmIAIKAZVnq16+f/vGPf+iiiy4yHQcISJQpADhFJSUlWrx4sWf6VFRUpDFjxmjcuHEaPXq0mjZtajoigCA1d+5c3XXXXUpPT5fD4TAdBwg4lCkAOEmWZSk9Pd2zccTKlSvVt29fjRs3TmPHjlX37t1ZUgPAJ1iWpcGDB+vmm2/WFVdcYToOEHAoUwBwAg4fPqwFCxYoLS1N8+fPV1RUlKc8DRs2TLGxsaYjAsAvWrRokW688UZt3LiR8+kAL6NMAcAvqK6u1ooVKzzTp+zsbA0bNsyzcUSnTp1MRwSAE2JZlkaMGKGrrrpK11xzjek4QEChTAHAD3bs2OEpT19++aU6dOjgKU+DBg1SWFiY6YgAUCvffPONJkyYoOzsbL6XAV5EmQIQtEpLS7VkyRJPgTpy5IjGjBmjsWPHasyYMWrevLnpiADgNWPHjtX555+vG264wXQUIGBQpgAEDcuylJWV5SlPy5cvV+/evT3Tp169erFxBICAtXLlSl1wwQXKy8tTRESE6ThAQKBMAQho+fn5WrhwoadAhYWFecrTiBEj1KBBA9MRAaDejB8/XiNHjtTf//5301GAgECZAhBQqqurtWrVKs+huRs3btSZZ57pOTTX6XTKZrOZjgkARqxfv15nnXWW8vLyFB0dbToO4PcoUwD83u7duz2Tp4ULF6pNmzae8jR48GCFh4ebjggAPuPiiy9Wv379dM8995iOAvg9yhQAv1NWVqalS5d6pk8HDhzQ6NGjPRtHtGrVynREAPBZWVlZGj58uLZs2cIZecApokwB8HmWZWnz5s2e8vTtt9+qe/funulT79695XA4TMcEAL9x5ZVXqkuXLpo0aZLpKIBfo0wB8EkFBQXHbRxhs9k0btw4z8YRcXFxpiMCgN/Kzc3VwIEDlZubq/j4eNNxAL9FmQLgE1wul9asWaN58+YpLS1NGRkZGjx4sKdAde7cmY0jAMCLJk6cqFatWmnGjBmmowB+izIFwJi9e/cet3FEixYtPOVpyJAhnIMCAHVo+/bt6tOnj7Kzs9WkSRPTcQC/RJkCUG8qKio8G0ekpaVpz549GjVqlGfjiNatW5uOCABB5cYbb1RsbKweffRR01EAv0SZAlBnLMtSbm6uZ+ne0qVLlZyc7Dk0t1+/fmwcAQAG7d69Wz169FBWVpZatGhhOg7gdyhTALyqqKhIixYt8kyfqqurPeVp1KhR3OgMAD7mtttuk2VZeuqpp0xHAfwOZQpBybIslVe7VVBRpfJqt9yWJUuS3SaF2e1qGBGqmFAHGx6cALfbrbVr13rK07p16zRo0CBPgerWrRufRwDwYfv371dSUpLWr1+vNm3amI4D+BXKFIJGaZVLOwpLdaC0QkUV1XJbluw2myzL0o//L7DZJJtssiRZshQbFqLGkWFq1zBK8RGhRvP7kv3792v+/PlKS0vT/Pnz1bRpU095OvPMMxUVFWU6IgDgJNx7770qLCzUiy++aDoK4FcoUwholmXpYGmlcvKLdaSsUrIkdy2u47BJUaEOJTaKUevYSDnswTVpqays1LJlyzyH5u7YsUMjR470FKi2bduajggAOAWHDx9W586dtXr1anXo0MF0HMBvUKYQsA6XVmjVvgJVuSxVe+nLPOSH5WrJTWPVIS4qoJev5eXleZbuLVmyRF26dNHYsWM1btw4nX766QoJCTEdEQDgRZMnT9bu3bv1r3/9y3QUwG9QphBwqt1upR8s0q6iMrnq6KvbYbMpLiJE/VrGKyo0MHajO3bsmL766ivP9KmsrMwzeRo9erQaN25sOiIAoA4VFBQoISFBy5YtU2Jiouk4gF+gTCGg5JdVavmeo6pyu+Wu469smyS7zaYezRqofZz/3SPkdru1YcMGT3las2aN+vfv7zk0Nzk5OaAnbwCAmv75z39q48aNeuutt0xHAfwCZQoBY39xub7be7TOplG/xmGzKaFRtLo2jvH58nHw4EEtWLBA8+bN04IFC9SwYUNPeRo6dKiio6NNRwQAGHTs2DE5nU4tWrRIycnJpuMAPo8yhYCwr7hcKw0UqR85bFLHuGilNGtgJsCvqKqq0vLlyz2H5m7ZskXDhw/3FKj27dubjggA8DEzZ87Ud999pw8++MB0FMDnUabg9w6VVujb3fnGitSPHDapc+NYdWkcYzTH1q1bPRtHLF68WE6n01OeBgwYoNBQtngHAPy60tJSOZ1Off755+rVq5fpOIBPo0zBr1VUu5S27ZCq6/oGqRPksEmDWzdSk6jwenvOkpISLV682DN9KioqOm7jiKZNm9ZbFgBAYHj66ae1YMECffrpp6ajAD6NMgW/9u3ufB0oqZAvfRGHO+wa27GpQuz2Orm+ZVlKT0/3TJ9Wrlypfv36eQpU9+7dZa+j5wYABIfy8nIlJCTogw8+UP/+/U3HAXwWZQp+a8+xMq3eV2B8ed/P2W1S2waR6t0izmvXPHz4sBYsWKC0tDTNnz9f0dHRnvI0fPhwxcSYXVoIAAg8L730kj788EPNnz/fdBTAZ1Gm4JcqXW7N23rQZ5b3/ZzDJp3RprEaR4bV6uOrq6u1YsUKz/QpOztbw4YN8xSoTp06eTkxAADHq6ysVOfOnfXGG29oyJAhpuMAPokyBb+Um1+sjYeP+dxU6qeaRYXpjDYnftDtjh07POXpyy+/VIcOHTzladCgQQoLq10xAwCgtl577TW9+uqrWrx4sc8f/wGYQJmC37EsS3O3HFS5y206ym+y26QxHZopKtTxi/+9tLRUS5Ys8RSoI0eOaMyYMRo7dqzGjBmj5s2b13NiAACOV11draSkJD333HMaNWqU6TiAz6FMwe8cKKnQd3uOqtrHv3TtkpyNopXc9PuzpyzLUlZWlqc8LV++XL1799bYsWM1btw49ezZk40jAAA+5+2339Yzzzyjb7/9lukU8DOUKfidb3Yd0cHSStMxTojDJlWsX+opUGFhYZ7yNHz4cDVo4FuH/AIA8HMul0s9evTQo48+qrPPPtt0HMCnUKbgVyzL0se5++Wj+07UUFZSrHkvPKb+PVI0duxYOZ1O3tUDAPidDz/8UA899JBWr17N6xjwE5Qp+JXiymot2n5YLj/5snXYpB7NG6p9wyjTUQAAqDW3262+ffsqNTVV559/vuk4gM/gBg34lYLyKvnT+2EuSzpS5h9LEgEA+DV2u13Tp0/X5MmT5Xb79gZQQH2iTMGvHCmv8vmNJ34unzIFAAgA55xzjqKjo/X++++bjgL4DMoU/EpRRZXXr3l43x49+reJumZgsv7cP0mzp9/v1euXVvEOHgDA/9lsNs2YMUNTp05VdXW16TiAT6BMwa+4vLzzhMvl0kM3/FlNW7XWC4tWavbXa3TGOX/06nNwWyIAIFCMGjVKzZo101tvvWU6CuATKFPwK97exS8vfZ2OHtyvq+5OVURUlMLCI9S1T3+vPgdVCgAQKGw2mx588EFNmzZNVVXeXy0C+BvKFPyKt8+0Pbx/r5q2ai1HSIh3L/wT7CALAAgkZ555pjp16qRXX33VdBTAOMoU/Eqozbtfsk1atNKhfXvkqsO133baFAAgwMyYMUMPPvigKioqTEcBjKJMwa/ERXh3guTs3kvxTZvr34//U+WlpaqsKNfmtSu9+hzRoQ6vXg8AANMGDBig7t27a/bs2aajAEZRpuBX4iPDFGL33qTH4XDovhde076d2/XX4f10/dA+WvbFJ167viQ1iQzz6vUAAPAF06dP10MPPaTS0lLTUQBjbBZbjcGPlFW5lLbtoNc3oqgrITaberVoqDYNIk1HAQDA6y688EINGjRId955p+kogBGUKfgVy7L0ad4BVftJm3LYbBrZvoliwupugwsAAEzJzMzUyJEjlZeXp9jYWNNxgHrHMj/4FZvNphbR4aZjnLBQu417pgAAASs5OVkjRozQM888YzoKYASTKfido+VV+nrnYbl8/CvXYZO6NYlVQqMY01EAAKgz2dnZOuOMM5Sbm6u4uDjTcYB6xWQKfic+IlRRfjDtsSS1axhlOgYAAHWqc+fOOuecc/Tkk0+ajgLUOyZT8Es7Cku1/kCRXD785ds6NkKnt4o3HQMAgDq3detW9evXTzk5OWrcuLHpOEC9YTIFv9Q6NlJhDt89DNdhk7o24UZcAEBw6Nixoy666CLNnDnTdBSgXjGZgt86Ulapb3Yd8bl7pxw2qUvjWHVuzL1SAIDgsWvXLvXo0UObNm1S8+bNTccB6gWTKfitxpFhatcwSr42oIoOC1Fio2jTMQAAqFdt2rTRlVdeqUceecR0FKDeMJmCX3O5Lc3fdlBl1W7TUSR9P5Ua0b6pYjlXCgAQhPbt26ekpCRlZGTotNNOMx0HqHNMpuDXHHabzmjTWKF28+Mph03q1zKeIgUACFotW7bUtddeq4ceesh0FKBeMJlCQCgor9LXu46o2m3my9lhk3o2b8hW6ACAoHfo0CF16dJFa9euVbt27UzHAeoUZQoBo6iiSkt2fl+o6vOL2mGT+rSIU+sGkfX4rAAA+K4HHnhABw4c0CuvvGI6ClCnKFMIKKVVLq3ad1QF5VV1vsufwyaFOuzq3ypejSPD6vbJAADwI0ePHlVCQoJWrFghp9NpOg5QZyhTCDiWZWl7YanSDx6T26qbKZXDJrVrGKWUpg3k8IH7tQAA8DXTp09Xbm6u3nzzTdNRgDpDmULAKq1yad3+Ah0qq5QsyRv7/TlsUlRoiHo3b6jGUUyjAAD4NUVFRXI6nVqyZIm6du1qOg5QJyhTCHilVS5tLSjR1oJSyZKqT/JL3mGTLEmtYiKU0ChG8RGhdRMUAIAA88gjj2jt2rV67733TEcB6gRlCkHDbVnaX1yhw2WVOlJWqaKKKkmS3WY7fimgJbksSzFhDjWKDFOTyDC1jIlQmIOTBAAAOBklJSXq1KmT0tLS1KNHD9NxAK+jTCFoWZal0iqXKlxuuSxLlvX9uVWhdptiwkJkt3EvFAAAp2rWrFlavHix/ve//5mOAngdZQoAAAB1pry8XE6nUx999JH69etnOg7gVaxbAgAAQJ2JiIjQ/fffr8mTJ5uOAngdkykAAADUqYqKCnXu3FlvvfWWBg8ebDoO4DVMpgAAAFCnwsPDlZqaqtTUVNNRAK+iTAEAAKDOXXXVVdq5c6e+/PJL01EAr6FMAQAAoM6FhoZq6tSpSk1NFXeZIFBQpgAAAFAvLr/8ch09elRpaWmmowBeQZkCAABAvXA4HJo2bRrTKQQMyhQAAADqzYUXXqjKykp98sknpqMAp4wyBQAAgHpjt9s1ffp0TZ48WW6323Qc4JRQpgAAAFCvxo8fr7CwMH344YemowCnhEN7AQAAUO/mzZunO+64QxkZGXI4HKbjALXCZAoAAAD1buzYsWrUqJHeeecd01GAWmMyBQAAACO++uorXXfdddq0aZNCQ0NNxwFOGpMpAAAAGDF8+HC1bdtWb7zxhukoQK0wmQIAAIAxy5Yt0xVXXKHs7GyFh4ebjgOcFCZTAAAAMGbw4MHq2rWr5syZYzoKcNKYTAEAAMCo1atX67zzzlNubq4iIyNNxwFOGJMpAAAAGNW3b1/17dtXL774oukowElhMgUAAADj0tPTNWbMGG3ZskXR0dGm4wAnhMkUAAAAjOvevbuGDh2qZ5991nQU4IQxmQIAAIBP2LRpk4YOHaq8vDw1aNDAdBzgdzGZAgAAgE/o2rWrxo0bp1mzZpmOApwQJlMAAADwGXl5eRowYIBycnLUqFEj03GA38RkCgAAAD7D6XTqvPPO0+OPP246CvC7mEwBAADAp+zYsUO9e/fW5s2b1bRpU9NxgF9FmQIAAIDPufnmmxUZGanHHnvMdBTgV1GmAAAA4HP27NmjlJQUZWVlqWXLlqbjAL+IMgUAAACfdMcdd6iqqkrPPPOM6SjAL6JMAQAAwCcdPHhQXbp00fr169W2bVvTcYAaKFMAAADwWffdd5/y8/P10ksvmY4C1ECZAgAAgM86cuSIEhMTtWrVKnXs2NF0HOA4nDMFAAAAn9W4cWPdcsstmj59uukoQA1MpgAAAODTCgoKlJCQoG+++UadO3c2HQfwYDIFAAAAnxYXF6fbb79dU6dONR0FOA6TKQAAAPi84uJiOZ1OLViwQCkpKabjAJKYTAEAAMAPxMTE6O6779aUKVNMRwE8mEwBAADAL5SWlsrpdOqzzz5T7969TccBKFMAAADwH88++6zmzZunzz77zHQUgDIFAAAA/1FRUaGEhAS99957GjhwoOk4CHKUKQAAAPiV2bNn6/3339eCBQtMR0GQYwMKAAAA+JWrr75aW7du1ZIlS0xHQZCjTAEAAMCvhIaGavLkyUpNTRWLrGASZQoAAAB+54orrtDBgwe1cOFC01EQxChTAAAA8DshISGaOnWqJk2axHQKxlCmAAAA4JcuueQSlZaW6vPPPzcdBUGKMgUAAAC/ZLfbNX36dKWmpsrtdpuOgyBEmQIAAIDfOu+882S32/XRRx+ZjoIgxDlTAAAA8GtffPGF7rnnHm3YsEEOh8N0HAQRJlMAAADwa2eddZZiY2P13nvvmY6CIMNkCgAAAH5v4cKFuummm7Rx40aFhISYjoMgwWQKAAAAfm/kyJFq1aqV3nzzTdNREESYTAEAACAgLF26VFdddZWys7MVFhZmOg6CAJMpAAAABIQhQ4YoMTFRr776qukoCBJMpgAAABAwVq5cqQsvvFC5ubmKiIgwHQcBjskUAAAAAsbpp5+unj176uWXXzYdBUGAyRQAAAACyrp163TOOecoLy9PUVFRpuMggDGZAgAAQEDp1auXBg0apOeee850FAQ4JlMAAAAIOFlZWRoxYoTy8vIUGxtrOg4CFJMpAAAABJykpCSNGjVKTz31lOkoCGBMpgAAABCQcnJyNGjQIOXl5SkuLs50HAQgJlMAAAAISImJiRo/fryeeOIJ01EQoJhMAQAAIGBt27ZNffv2VXZ2tpo0aWI6DgIMkykAAAAErA4dOuiSSy7RzJkzTUdBAGIyBQAAgIC2e/dude/eXRs3blSLFi1Mx0EAoUzhF1mWpSq3JZdlyW1ZsttscthsCrXbZLPZTMcDAAA4KX//+99ls9k0a9Ys01EQQChTkCRVVLt1uKxC+WVVOlxWqaKKalmWpe970/flyW1ZcthsahAeoiZRYWoUEaYmUWEKc7BaFAAA+Lb9+/crKSlJGzZsUOvWrU3HQYCgTAUxy7KUX16l3Pxi7S+pkN1mU7X7xL8cQmw2uWWpVUyEEhpFKz4irA7TAgAAnJp77rlHx44d0wsvvGA6CgIEZSoIWZal3cfKtPFwscqr3XJ54UvAYZMiQx1KbtpArWIivJASAADAuw4fPqzOnTtr9erV6tChg+k4CACUqSBTXu3S6n0FOlJW5ZUS9XMOm03NosLUu2Wcwln+BwAAfExqaqr27Nmjf/3rX6ajIABQpoKEZVnaVVSm9QeK5LIs1eVful2Sw25T7xZxOi2WKRUAAPAdR48eVUJCgpYvX66EhATTceDnKFNBwLIsbThQpB1FpXLV49+2wyZ1io9WUpNYdgAEAAA+48EHH9TmzZv173//23QU+DnKVICzLEur9xVob3F5vRapHzlsUtsGUerZvAGFCgAA+IRjx46pU6dO+uqrr5SUlGQ6DvwYN7UEMMuytO5AobEiJUkuS9pZVKrMQ8fMBAAAAPiZ2NhY3XXXXZo6darpKPBzTKYCWF5+sbIOHzNWpH7KYZN6NGug9nHRpqMAAACopKRETqdTc+fOVc+ePU3HgZ9iMhWgiiurfaZISd9PqDYcPKbSKpfpKAAAAIqOjtY//vEPTZ482XQU+DHKVACyLEsr9x71mSL1I7dladW+o2IYCgAAfMFf//pXrVu3TitXrjQdBX6KMhWA8o6W6FhltekYNViSCsqrtL2w1HQUAAAARURE6IEHHlBqaqrpKPBTlKkA43Jb2ni42OemUj9yWVLmoWNyM50CAAA+4Nprr1VOTo6++eYb01HghyhTAWb3sTL5+gbkliXtKy43HQMAAEBhYWGaPHmyJk2axK0IOGmUqQCTk1+sah//RlBtWcrJLzEdAwAAQJI0YcIE7d27V19++aXpKPAzlKkAcrS80m92yyusqNKxCt+7rwsAAASfkJAQTZ06VampqUyncFIoUwFkd5G5w3lPlmVJe4rLTMcAAACQJF166aUqLCzU3LlzTUeBH6FMBZDDZZVevd4NI07Xhm+/9uo1f2RJOlzq3bwAAAC15XA4NG3aNE2ePJnpFE4YZSpAWJalIj9bNldQUWU6AgAAgMcFF1wgl8uljz/+2HQU+AnKVIAoq3bp+3mP/6hyWapwuU3HAAAAkCTZ7XZNnz5dqampcrv5Nwp+H2UqQBRWVMtm8/6m6FsyNujv5wzVVad31bP33abKCu9tae6w21TEdAoAAPiQc889V1FRUfrPf/5jOgr8AGUqQFS762Yq9fVn/1XqK2/rufnfau/2rfrghae8ev26yg0AAFAbNptNM2bM0JQpU1Rd7V+3UKD+UaYChNuyVBf3Sp51xTVq0vI0xcbF68Ib/q5vPv+f9y5ufZ8bAADAl4wePVpNmzbV22+/bToKfBxlKkB4f4Hf95q0aOX5edNWrXX04AHvXdwm2eosOQAAQO3YbDY9+OCDmjZtmqqquCUBv44yFSAc9rqpJYf37/3/P9+3R/HNmnv1+g47ZQoAAPieoUOHqkOHDnrttddMR4EPo0wFiMgQR52Mp+a99ZqO7N+rYwVH9d+Xntbgs8Z77dqWJUWG8CUIAAB804wZMzRjxgxVVFSYjgIfxb9kA0TD8FC56mAzhyHnnq/pEy/XTaMHqnnrdrroxr977doVFRX604Xn67777tPbb7+tjIwMVVZykC8AAPANAwcOVEpKimbPnm06CnyUzeKI54Axb+tBlVa5TMc4YRFyq3LTKmVkZCg9PV0ZGRnasWOHEhISlJKSou7du3t+PO200+pk63cAAIDfsmbNGo0fP155eXmKjIw0HQc+hjIVQFbuPardx7x3DlRd69AwSr1aNDzu98rKyrRp0yZPufrxx4qKiuMKVkpKipKTk9WgQQND6QEAQLC44IILdMYZZ+iOO+4wHQU+hjIVQHYUlmr9gSK5/OCvNMRuU58WDXVa7Im9w3Pw4EFlZGR4/peenq6NGzeqWbNmNUpWYmKiQkJC6vhPAAAAgkVGRoZGjRqlLVu2KCYmxnQc+BDKVACpdlv6PO+AX5SpULtN5ziby34KS/dcLpe2bt3qmV79WLL27Nmjzp071yhZLVu2ZKkgAAColcsuu0w9evTQfffdZzoKfAhlKsCs21+o7YWl8uW/VLskZ6NoJTetmyV6JSUl2rhxY42SZVmWp1j9WLKSkpJ4hwkAAPyuzZs3a8iQIcrLy1PDhg1//wMQFChTAeZYZbUWbT+kOtjYz2vsNmlMh2aKCnXU23NalqUDBw4cdx9Wenq6Nm/erFatWtUoWU6nUw5H/eUDAAC+789//rM6dOigqVOnmo4CH0GZCkBf7zyiI2WVPjmdsklqHh2uQa0bmY4iSaqurlZeXl6NKdaBAwfUtWvXGksFmzf37qHFAADAf2zZskX9+/dXdna2GjdubDoOfABlKgCVVlVrwbbDPnnvVIjNpjEdmyoixLenPseOHVNWVlaNkhUaGlpjitWtWzdFRUWZjgwAAOrB9ddfr8aNG+v//u//TEeBD6BMBaitR0uUceiYTxUqh82mXs0bqG1D/ywelmVp7969xy0VzMjIUHZ2ttq2bVujZHXs2FF2O+diAwAQSHbu3KlevXpp06ZNatasmek4MIwyFaAsy9KSnUd0tLzKJ5b72SQ1jQrT4NaNAm5HvaqqKuXk5NQoWYcPH1ZSUlKNktWkSRPTkQEAwCn429/+prCwMD3++OOmo8AwylQAK61y6csdh1TpMv9XHBFi14h2TXx+eZ83FRYWKjMz87iClZGRocjIyOPuw+revbu6du2qiIgI05EBAMAJ2Ldvn5KSkpSZmalWrVqZjgODKFMB7lhltRbvOKwqg9v7hTlsGtGuiaJCOUjXsizt2rWrxhQrLy9PHTp0qDHFateuHUsFAQDwQXfddZfKy8v17LPPmo4CgyhTQeBYRbWW7DqsKpdVr0v+bJLCHHYNa9tY0WEUqd9SWVmpzZs31yhZhYWFSk5OrlGy4uPjTUcGACCoHTp0SF26dNHatWvVrl0703FgCGUqSJRWVWvZ7nyVVrnrZVMKh02KCQvRoNaNFBlES/u87ejRo8ftJpiRkaHMzEw1aNCgxlLBLl26KCwszHRkAACCxv33369Dhw5p9uzZpqPAEMpUEHFblnKOFCs7v1h1eRuVwyZ1axIrZ3x0wG024Qvcbrd27NhRo2Rt27ZNnTp1qlGy2rRpw98DAAB1ID8/X4mJiVqxYoWcTqfpODCAMhWEiiqq9N3eApVWubw6pXLYbIoJc6h/q3jFsKyv3pWXl2vTpk01lgqWlpYqOTn5uJKVkpKihg0bmo4MAIDfmzZtmrZs2aI33njDdBQYQJkKUpZl6VBppXLyi3W4rFKyJHctrmP/YeDRLCpciY2i1TgyjCmIjzl8+HCNKVZWVpYaN2583H1YKSkp6ty5s0JDQ01HBgDAbxQWFiohIUFLlixR165dTcdBPaNMQaVVLm0rKNHe4gqVVFZ7ytAvTa1CbDZZkixZigkNUavYCHWMiwqqLc8Dgdvt1tatW2uUrF27dikhIaHGUsFWrVpRkgEA+BUPP/yw1q1bp/fee890FNQzyhSO47YsHausVkF5lY5VVqvaZcllWXLYbQqx29QgLERxEaGKDQvhH9cBqLS0VBs3bqxRsqqqqmpMsZKTkxUbG2s6MgAAxhUXF8vpdGr+/Pnq3r276TioR5QpAL/rwIEDxx08nJ6erk2bNql58+Y1SlZCQoJCQrhnDgAQXJ588kl9/fXX+uijj0xHQT2iTAGoFZfLpS1btnimVz/+uHfvXnXp0qVGyWrRogXTTABAwCorK5PT6dTHH3+svn37mo6DekKZAuBVxcXFysrKqrFU0Gaz1Th8OCkpSdHR0aYjAwDgFc8995w+//xzffHFF6ajoJ5QpgDUOcuytH///uO2bE9PT1d2drZOO+20GlOsTp06yeFgUxMAgH+pqKhQYmKi3nnnHQ0aNMh0HNQDyhQAY6qrq5Wbm1ujZB08eFDdunWrUbKaNWtmOjIAAL/plVde0TvvvKNFixaZjoJ6QJkC4HOKioqUlZVVo2SFh4fXWCrYrVs3RUZGmo4MAIAkqaqqSl27dtXs2bM1fPhw03FQxyhTAPyCZVnas2dPjYKVm5urdu3a1ShZHTp0kN1uNx0bABCE3nzzTb300ktaunQpmy8FOMoUAL9WWVmpnJycGiXr6NGjSkpKqlGyGjdubDoyACDAuVwuJScna9asWRo7dqzpOKhDlCkAAamgoECZmZnHlayMjAxFR0d7itWPP3bt2lXh4eGmIwMAAsj777+vxx57TN999x3TqQBGmQIQNCzL0s6dO4/bsj0jI0NbtmxRx44da0yx2rVrxwsgAKBW3G63evXqpRkzZmj8+PGm46COUKYABL2Kigpt3rz5uJKVnp6uY8eOeQrWT0tWXFyc6cgAAD/w8ccfa8qUKVq7di338QYoyhQA/Ir8/PwaU6zMzEzFxcUdt2V79+7d1blzZ4WFhZmODADwIZZl6fTTT9c999yjiy++2HQc1AHKFACcBLfbre3bt9coWdu3b5fT6axRslq3bs1SQQAIYnPnztWdd96pjIwMDqQPQJQpAPCCsrIybdq0qUbJKi8vr7FUMDk5WQ0aNDAdGQBQDyzL0hlnnKGbbrpJV1xxhek48DLKFADUoUOHDh23ZXtGRoY2btyoJk2a1JhiJSYmKiQkxHRkAICXffnll/rrX/+qTZs28X0+wFCmAKCeuVwubd26tUbJ2r17txITE2uUrJYtW7JUEAD83IgRI3TllVfq2muvNR0FXkSZAgAfUVJSoo0bN9bYVdDtdh+3m2BKSoqSk5MVExNjOjIA4AQtW7ZMV1xxhXJyctiwKIBQpgDAh1mWpQMHDtSYYm3atEktW7Y8rmR1795dTqeTG5wBwEeNGzdOf/zjH3XjjTeajgIvoUwBgB+qrq5WXl5ejZK1f/9+denSpcZSwebNm5uODABBb9WqVTr//POVm5uryMhI03HgBZQpAAggx44dU1ZWVo2S5XA4aiwVTEpKUlRUlOnIABBU/vjHP2r48OG67bbbTEeBF1CmACDAWZalffv2Hbdle3p6unJyctS6desaJatjx44sFQSAOrJhwwaNGzdOeXl5io6ONh0Hp4gyBQBBqqqqSjk5OTWmWIcPH1a3bt1qlKymTZuajgwAAeGSSy5Rnz59dO+995qOglNEmQIAHKewsFCZmZk1SlZkZORx92GlpKSoW7duioiIMB0ZAPzKxo0bNWzYMOXl5XGIu5+jTAEAfpdlWdq9e3eNpYJ5eXlq3759jSlW+/btZbfbTccGAJ81YcIEJSYmKjU11XQUnALKFACg1iorK5WdnV2jZBUWFiopKalGyWrUqJHpyADgE/Ly8jRgwADl5uYqPj7edBzUEmUKAOB1R48eVWZm5nGHD2dmZqpBgwY1lgp26dJF4eHhpiMDQL2bOHGiWrZsqQcffNB0FNQSZQoAUC8sy9KOHTtqTLG2bdumTp061ShZbdu2lc1mMx0bAOrM9u3b1adPH23evJlNfvwUZQoAYFR5ebk2b95co2SVlJR4CtZPS1bDhg1NRwYAr7npppsUHR2tmTNnmo6CWqBMAQB80pEjR47bTTAjI0OZmZlq3LhxjYLVuXNnhYaGmo4MACdtz549SklJUVZWllq2bGk6Dk4SZQoA4Dfcbre2bdtWo2Tt2LFDiYmJNUrWaaedxlJBAD7v9ttvl8vl0tNPP206Ck4SZQoA4PfKysq0cePG4wpWRkaGKioqjttNsHv37kpOTlZsbKzpyADgceDAAXXt2lUbNmxQmzZtTMfBSaBMAQAC1sGDB2scPrxx40Y1a9asRslKSEhQSEiI6cgAgtQ//vEPFRQU6MUXXzQdBSeBMgUACCoul0tbtmypsVRwz5496ty5s6dk/fhjixYtWCoIoM4dOXJEiYmJWrVqlTp27Gg6Dk4QZQoAAEklJSXKyso6rmSlp6dLUo0pVlJSkqKjow0nBhBopkyZop07d+rVV181HQUniDIFAMCvsCxL+/fvr7FUcPPmzWrVqlWNktWpUyc5HA7TsQH4qYKCAiUkJOibb75R586dTcfBCaBMAQBwkqqrq5Wbm1ujZP14E/nPS1azZs1MRwbgJ/75z38qKytLb7/9tukoOAGUKQAAvOTYsWPKzMyscT9WaGjocfdhde/eXd26dVNkZKTpyAB8zLFjx+R0OrVo0SIlJyebjoPfQZkCAKAOWZalPXv21Jhi5eTkqG3btseVrJSUFHXs2FF2u910bAAGPfbYY1q+fLk+/PBD01HwOyhTAAAYUFVVpezs7Bol68iRI0pKSjquYKWkpKhJkyamIwOoJ6WlpXI6nfrss8/Uu3dv03HwGyhTAAD4kIKCAs9SwZ+WrOjo6BpTrK5duyoiIsJ0ZAB14JlnnlFaWpo+++wz01HwGyhTAAD4OMuytGvXruO2bM/IyNCWLVvUoUOHGiWrXbt2LBUE/Fx5ebkSExP1/vvva8CAAabj4FdQpgAA8FMVFRXavHlzjSlWUVGRkpOTa5Ss+Ph405EBnISXX35ZH3zwgebPn286Cn4FZQoAgACTn5/vKVg/lqzMzEzFxcUdt2V7SkqKunTporCwMNORAfyCqqoqde7cWa+99prOPPNM03HwCyhTAAAEAbfbrR07dhy3ZXt6erq2b98up9NZY4rVpk0b2Ww207GBoPf6669rzpw5WrJkCf+f9EGUKQAAglh5ebk2bdpUo2SVlZXVmGIlJyerYcOGpiMDQaW6ulpJSUl69tlnNXr0aNNx8DOUKQAAUMOhQ4eOWyqYkZGhrKwsNWnSpEbJSkxMVGhoqOnIQMB655139NRTT2n58uVMp3wMZQoAAJwQt9utrVu3HrfZRXp6unbt2qXOnTvXKFmtWrXiH36AF7jdbvXo0UMPP/ywzjnnHNNx8BOUKQAAcEpKS0u1cePGGksFq6urj7sPq3v37kpOTlZMTIzpyIDf+e9//6sHH3xQa9as4U0KH0KZAgAAdeLAgQPHTbEyMjK0ceNGtWjRokbJcjqdCgkJMR0Z8FmWZalPnz6aNGmSLrjgAtNx8APKFAAAqDcul0t5eXk1StbevXvVtWvXGksFmzdvzrvwwA8+//xz3XvvvdqwYYMcDofpOBBlCgAA+IDi4mJlZWXVuB/LbrfXmGIlJSUpKirKdGSg3lmWpYEDB+rvf/+7Lr/8ctNxIMoUAADwUZZlad++fTWmWJs3b1br1q1rlKyOHTvybj0C3sKFC3XzzTcrKyuLpbE+gDIFAAD8SlVVlXJzc2uUrIMHD6pbt241SlbTpk1NRwa8xrIsDRs2TNdcc42uvvpq03GCHmUKAAAEhKKiImVmZtYoWeHh4TUKVteuXRUZGWk6MlArX3/9ta6++mplZ2dzxpthlCkAABCwLMvS7t27j9uyPSMjQ7m5uWrXrl2NktW+fXvZ7XbTsYHfNWbMGF100UW6/vrrTUcJapQpAAAQdCorK5WdnV1jinX06FElJSUdV7JSUlLUuHFj05GB43z33Xe66KKLlJubq4iICNNxghZlCgAA4AdHjx6tsVQwMzNTMTExx23ZnpKSoq5duyo8PNx0ZASxP/zhDxozZoz+9re/mY4StChTAAAAv8GyLO3YsaPGUsGtW7eqY8eOx5Ws7t27q23btpyNhXqxbt06nXPOOcrLy+O4AEMoUwAAALVQUVGhTZs21ShZxcXFSk5OrrFUMC4uznRkBKALL7xQAwcO1F133WU6SlCiTAEAAHjRkSNHPAXrx5KVlZWl+Pj4GksFO3furLCwMNOR4ccyMzM1cuRI5eXlKTY21nScoEOZAgAAqGNut1vbtm2rMcXasWOHEhISapSs1q1bs1QQJ+xPf/qTkpKS9MADD5iOEnQoUwAAAIaUlZVp06ZNnnL1448VFRXHbdmekpKi5ORkNWjQwHRk+KCcnBwNHjxYubm5LCetZ5QpAAAAH3Pw4MEaSwU3btyoZs2a1ShZiYmJCgkJMR0Zhl1zzTVq06aNpk+fbjpKUKFMAQAA+AGXy6WtW7cedy5Wenq69uzZo86dO9coWS1btmSpYBDZtm2b+vbtq5ycHM5Fq0eUKQAAAD9WUlKijRs31ihZlmXVKFhJSUmKiYkxHRl15IYbblBcXJwefvhh01GCBmUKAAAgwFiWpQMHDhx3H1Z6ero2b96sVq1a1ShZTqdTDofDdGycol27dqlHjx7atGmTmjdvbjpOUKBMnQDLslRa5dLR8irll1fpSFmlqt1uudySzSbZbTZFhTjUJCpU8RFhiosIVZjDbjo2AADAcaqrq5WXl1djinXgwAF17dq1RsniH+T+59Zbb5XD4dCTTz5pOkpQoEz9hsLyKuUeLdGeY+WSJJuk6t/4dNkkOew2udyWwkPs6hQXrfZxUQqnWAEAAB927NgxZWVl1ShZISEhx23Z3r17d3Xr1k1RUVGmI+NX7Nu3T0lJSUpPT1fr1q1Nxwl4lKmfcVuW9hwrV3Z+sUoqq+W2pNp+guw/3PPZMiZCiY1iFB8R6rWcAAAAdcmyLO3du/e4pYIZGRnKzs5W27Zta0yxOnbsKLudN5B9wd13362SkhI9//zzpqMEPMrUTxRVVOm7vQUqrXLJ5eVPi8MmtWkQqe7NGiiEbzQAAMBPVVVVKScnp0bJOnz4sJKSkmqUrCZNmpiOHHQOHTqkLl26aM2aNWrfvr3pOAGNMqXvp1E5R4q1Ob9Y7jr8bNhtUqjdrv6t4tQkKrzunggAAKCeFRYWKjMzs8ZSwaioqBpLBbt27aqIiAjTkQPapEmTtG/fPs2ZM8d0lIAW9GWqotqtb3YfUXGl96dRv8ZhkzrFRyupSSznPwAAgIBlWZZ27dpVY4qVl5en9u3b1yhZ7dq1Y6mglxw9elQJCQlasWKFnE6n6TgBK6jLVFm1S0t2HFFZtavW90XVlsNmU+vYCPVu0ZBCBQAAgkplZaU2b95co2QVFhYqKSnJU7J+/DE+Pt50ZL80Y8YM5eTk6M033zQdJWAFbZmqqHbpyx2HVV7trvci9SOHTWrdIFK9m1OoAAAAjh49etwSwYyMDGVmZqpBgwbHTbFSUlLUpUsXhYdz28RvKSoqktPp1OLFi9WtWzfTcQJSUJYpl9vSoh2HVFJZ/xOpn3PYJGd8tJKaNjCcBAAAwPe43W7t2LGjRsnatm2bOnXqVGOpYJs2bXiT+iceffRRrV69Wu+//77pKAEpKMtU+sFCbSsolctH/uQOmzSkTWM1igwzHQUAAMAvlJeXa9OmTTWWCpaWlio5ObnGJKthw4amIxtRUlIip9OpuXPnqmfPnqbjBJygK1P5ZZVauuuIzxSpH0WFODS6Q1M57LyTAgAAUFuHDx+uMcXKyspS48aNj7sPKyUlRf+vvXsPjrJc7Dj+e3c3CYEkbLglmAtgEgMJAxTEw23aSMMI4SDeEMvF4bRkytj8gYqKgiOlDKAyjMeOiJLB1GAtcxSKTB04g5ICo9w8VTRwQpCbhEsCIUBINnt7+wclBReO4XXDbna/nxlHNrt53mf3r/3med/nzc3NVUxM5N8H9K233tL27du1adOmUE8l4kRVTPn8pv54rFbNXn+opxLAbkj9nJ01qFd0/tUEAACgvfj9fh09ejQgsn766Sfl5OQERFZaWlpEnSrocrmUnZ2tjRs3avjw4aGeTkSJqpiqunBFf77QGHarUtfZDKmwb08lxDpCPRUAAICI19TUpIMHDwZElsfjCbj58MCBA5WYmBjqKVv27rvvatOmTdqyZUuopxJRoiamTNPUf/14Tu5wLSlJhqR7nZ01OIXVKQAAgFA5d+5ca2Bdj6xDhw4pJSUlYBUrJydHDkf4/yHc7XYrNzdX5eXlGjNmTKinEzGiJqbONLq073SDvGH+du2God9mp3DtFAAAQBjx+Xz68ccfW1evrv//9OnT6t+/f0Bkpaamht2pgmvXrlV5ebm2b98e6qlEjKiJqf8+eV4Xmj2hnsYvshuGhqQkqU/XzqGeCgAAAH5BY2OjKisrA04VNAwj4FTB/Px8denSJWRz9Xq9GjBggN577z2NHTs2ZPOIJFERUy6vT1uO1srfQd5p1ziH/rZvz1BPAwAAABaYpqmzZ8/etGX7gQMHVFVVpbS0tIDIysrKkt1uvytz++ijj7Rq1Srt2rUr7FbOOqKoiKnTjS7tP9MgbwepKZshTc4Jv6VhAAAAWOf1elVdXR0QWbW1tcrLyws4VbBXr15Bn4PP59OgQYO0YsUKTZgwIejjR5uoiKnKusuqqr8a6mm0md0w9GCf7kqKi/z7HgAAAESS/Px8vfPOOyooKGjz71y+fFmVlZUBkRUXFxewipWXl6f4+PhfNcdPPvlEy5cv1759+277x/tZs2YpPT1dS5Ys0c6dOzV79mxVVVX9quNGoqiIqY5yvdR1jv+7biqT66YAAAAiUt++fVVaWqrCwsJbPm+apmpqagICq7q6Wn369AmIrH79+slms7Xp2H6/X0OHDtWiRYv0yCOP3PI1N8bUnSgrK1Npaal27dp1R7/3c263W9OmTdP+/ft14sQJbd++/Y4C9W4J/30cg+ByizfoY254/1+17Q//rksXzqtH73s0be58/WZccJZKvaapC80eZbJDOgAAQFQyDEPp6elKT09XUVFR68/dbreqqqr03XffqbKyUqWlpTpw4IDq6+s1cODAgMjq3r17wNg2m02LFy/WwoUL9dtJk+S4S9dr3akxY8Zo7ty5mjJlSqincltty9cOztMO10qlZvbVknUbVb6/Sk/+03P6/Yslulh7Lmjju7y+oI0FAACAu6Nv377atm2bFi1apCeffFJPP/20EhMTlZ+fr/3790uSZs6cqZMnT2rSpElKSEjQG2+8IUnavXu3Ro0aJafTqcGDB6uioqJ13IKCAi1YsEAPPvigHnjgAY0YMULLly9XUVGR4uLiZBiG+vfvr/T0dC1dulSFhYVKTU1V7969NX78eL3wwgt67rnnlJubK6fTqWXLlmvC7BJtqj6nPx6r1e59+zV06FAlJiZq6tSpcrlcrceuqKhQenp66+Ply5crKytLiYmJysvL08aNGyVJhw4d0pw5c/T1118rISFBTqdTktTS0qJ58+YpMzNTKSkpmjNnjpqbm//i5xgbG6u5c+dqzJgxd21zDisiPqb87XQW46jxk9QtJVU2m02jiyard59+qv7+f4I2frjfDwsAAAB/2WeffaannnpKDQ0Nevjhh1VSUiJJKi8vV2ZmpjZv3qzGxka9+OKLqqmp0cSJE7Vw4ULV19drxYoVevzxx1VXV9c6Xnl5ud5//31duXJFffr0kSRt2bJF33zzjXbv3q2PP/5YO3bs0FdffaULFy4oJydH8+bNU0lJiZqamrRq1Sp5PB65XC45nN30weuL5fW41dDYrEcffVQzZ85UfX29pkyZok8//fS27ysrK0s7d+7UpUuX9Nprr2nGjBk6c+aMBgwYoNWrV2vkyJFqbGxUQ0ODJOmll17S4cOH9e233+rIkSOqqanR4sWL2++Dv4ui4jS/9lDxn3/Q5rL3VFtzSpLkarqqyxfrgzb+tm3b9Nd//1TQxgMAAMDdMW7cuNZ/T5w48abnbtzw4cbX3e71P9/Rb+DAgTc93rx5s7p2/f9rQ7788ktlZWW1Pp43b95Nrz927JgkKTk1TY5DB3X42z9JMuX1ejV37lwZhqEnnnhCK1euvO37u/G0u6lTp2rZsmXau3evJk+eHPBa0zS1Zs0aHThwQN26dZMkvfLKK5o2bZqWLVt222N0FBEfU7Z22F68tuaU3n31BS0qW6/7htwvu92u5x8plIK4mjR+3Dj9M6tTAAAAHcr1jSV27dqlI0eOaN26dZKk48ePq1+/fvJ4PHI4HAEbUDzzzDNau3atOnXq1DqWx+PRq6++qvnz56ugoEDTp09XcXFx6/OGYai6ulrZ2dmSrl1jNHv2bM2aNUuStHDhQp09e1alpaUqKipSRUWFYmNjJUk7Nq6Xq8WthrqzstvsykxPuyn0rq983cqHH36olStX6vjx45Ku3bj4/Pnzt3xtXV2dmpqaNGzYsNafmaYpny8yLmmJ+JiSJIfNCOo9plqam2QYhpKSr13Q9+Wn/6GT1cHdKjLOEfFnYAIAAEStn29JnpGRoZkzZ2rNmjVt/p07kZGRoQULFmjBggWtPzNNUy0+v77etVPr3jgt0zRbj3Hy5MmbVriuO3HihIqLi/XFF19o5MiRstvtGjJkiK5vEP7zOfbo0UPx8fGqrKxUWlqa5fmHq6j4xp4UG9xmzMi+T5N+94965e8e1j+MHqQTh/+s/n81PGjj2w1D3eJjgzYeAAAAwktKSoqOHj3a+njGjBnavHmztm7dKp/PJ5fLpYqKCp06dSooxysuLtbq1au1Z88emaapq1ev6vPPP5enuUmjR42Sw+HQ22+/La/Xqw0bNmjv3r23HOfq1asyDEM9e/aUJH3wwQf64Ycfbnpfp06dktvtlnRt58Di4mI9++yzqq2tlSTV1NRo69atvzjnlpaW1o0w3G63XC6Xwu2uTlERU93bIUymPztf/7bnoMp2V+p3Ly/Sv6zboMIp04MytiEpmRv2AgAARKyXX35ZS5YskdPp1IoVK5SRkaFNmzZp6dKl6tmzpzIyMvTmm2/K7/cH5Xj333+/1qxZo5KSEiUnJys7O1tlZWWSru2ct2HDBpWVlSk5OVnr16/XY489dstx8vLy9Pzzz2vkyJFKSUnR999/r9GjR7c+P3bsWOXn5ys1NVU9evSQJL3++uvKzs7WiBEjlJSUpMLCwjbdADg3N1fx8fGqqanRQw89pPj4eJ04ceLXfxhBFBU37a250qxvzl4K6ql+7cmQNPm+1Ha53gsAAABAcETFylRyp9h22yK9PSTEOggpAAAAIMxFRUx1jrEH/bqp9mI3pHudnUM9DQAAAKBdLV26VAkJCQH/TZgwIdRTa7OoOM1Pkk5dadafzlwK+5vh2g2pKDtFMbao6FwAAACgw4qab+z3JHRSRzhzLj0pnpACAAAAOoCo+dZuMwzd6+wsWxgHld2QcpK7hHoaAAAAANogamJKku7rnhC2qz4249qqVBJbogMAAAAdQniWRTuJsdk0vLdT9jBcnYqx2TS4V1KopwEAAACgjaIqpiSpV5c4pSXGh9XpfnZDeuAepxxhumoGAAAAIFBUfnsfnJKkOLtN4dBTdkPKTOqsnp3jQj0VAAAAAHcgKmMqxmbT32T2UEyIz/ezGddWyoakcHofAAAA0NFEZUxJ127kW5DZQ7EhWqGyG4Z6dY7Tb+5JltER9mwHAAAAcJOouWnv7TR5fNpx8oJafD757tInYTcMpSV20rDUroQUAAAA0EFFfUxJks9v6oe6yzp+qaldg8omyW4zNDS1q9IS49vvQAAAAADaHTF1gwvNbu05fVEenz/oUWU3DPXsHKthqU7FOaL27EoAAAAgYhBTP+PzmzrWcFWHL16Vx2fK9ys+HkOSYUjdOsUot3uiUrqwYx8AAAAQKYip2zBNU3VNbh2ub1Rdk1t2myG/acr/C5+WwzBk6lpI9XXGK8vZRV1iHXdjygAAAADuImKqDXx+U5daPLro8uh8s1sNLo+8flN+05QhyWYYinPY1D0+Vt3jY+XsFKOEGDubSwAAAAARjJgCAAAAAAvYCQEAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMACYgoAAAAALCCmAAAAAMCC/wWzvTtYz9jl0wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1080x1080 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_graph_nice(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get sequences\n",
    "sequences, seq_names = utils.get_sequences_for_ancestry()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "**Problem 4 (b).** Given $n$ sequences each of roughly the same length $m$, what would the time complexity be for constructing such a phylogenetic tree? Can you think of any algorithms or heuristics that might make the process faster? \n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 4 (c).** Assume a direct correlation between the distance between any two nodes and the number of years (in millions) between their evolution. Assuming `Grumpig` was the first Pokémon to evolve, when did life first come to be in the fictional scenario?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "For this part, feel free to use any of `networkx`'s in-built functions (or any graph-specific library you may have chosen for Problem 4).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Life evolved ??? million years ago in the Pokémon world\n"
     ]
    }
   ],
   "source": [
    "how_long_ago = \"???\" # Replace with your answer\n",
    "print(f\"Life evolved {how_long_ago} million years ago in the Pokémon world\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to test the robustness of such a tree reconstruction algorithm is to consider collection of nodes independently and see if the recontructed sub-trees match the bigger tree.\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "**Problem 4 (d).** Find an  edge between intermediate nodes with the largest weight in the phylogenetic tree and remove that edge- this will produce two disjoint cluster of nodes. Re-run your tree reconstruction algorithm on these two sets of Pokémons. Do your reconstructed tree match the larger phylogenetic tree?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "For this part, feel free to use any of `networkx`'s in-built functions (or any graph-specific library you may have chosen for Problem 4).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Problem 4 (d).** Can you inspect the matrix of distances between the Pokémons and predict whether the reconstructed trees would always be unique? Why/why not?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 5: Tracing Evolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    This problem is a \"Challenge Problem\". This means it is a problem of unknown difficulty that might be quite challenging (unlike the earlier problems, we don't have a reference solution for this one, or a clear idea how hard it might be). We do hope all students will at least attempt this and that more ambitious students will work hard to solve it and learn interesting things by the attempt (whether or not it is successful), but not get frustrated if you can't get to the desired answer.  As a \"Challenge Problem\" it means that you shouldn't be worried if you are not able to solve this, though, and you can get full expected credit on this assignment without answering it.\n",
    "</div>\n",
    "\n",
    "\n",
    "Now that we can construct Phylogenetic trees using sequence alignment, we can attempt to construct these trees for different organisms and trace their evolution through time. You're given reads processed from a FASTA file for Hemoglobin Beta Proteins, which can be used to then trace evolution based on how similar their sequences are across organisms from different kingdoms. Each record has the following relevant information in Tuple format:\n",
    "\n",
    "`((uniprot identifier, full name, shortened name, group), (sequence))`\n",
    "\n",
    "As you may notice, running our nearest-neighbor reconstruction algorithm on this data will give a Phylogenetic tree that does not fully correspond to what we know about the evolution of these species."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "**Problem 8 (Challenge).** Construct a Phylogenetic Tree using the given sequences. Implement and use any tree-reconstruction method of your choice, and see if it works better than the nearest-neighbor method with a linear-gap penalty.\n",
    "    \n",
    "For visualization, use the short name to display in the evolution tree.\n",
    "    \n",
    "</div>\n",
    "    \n",
    "This is an open-ended question, and is inspired by https://www.mimuw.edu.pl/~lukaskoz/teaching/sad2/lab6/readme.html. You are free to use any approach to deal with the issue. Make sure you provide your code, along with any assumptions you may have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(('P01941.1', 'Tupaia glis', 'Tgli', 'Mammalia'), 'VLSPGDKSNIKAAWGKIGGQAPQYGAEALERMFLSFPTTKTYFPHFDMSHGSAQIQAHGKKVADALSTAVGHLDDLPTALSALSDLHAHKLRVDPANFKLLSHCILVTLACHHPGDFTPEIHASLDKFLANVSTVLTSKYR')\n"
     ]
    }
   ],
   "source": [
    "sequences = utils.get_sequences_for_tree()\n",
    "print(sequences[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Write a description of your algorithm, and things you learned from working on this here._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    \n",
    "Is this (using Hemoglobin proteins) the best way to trace and visualize evolution? Why do you think it is useful, and what could the possible downsides of this be?\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Type your answer here_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "   <div class=\"alert alert-block alert-danger\">\n",
    "    <center>\n",
    "      \n",
    "**End of Project 2!**\n",
    "        \n",
    "Remember to follow the submission directions above to submit your assignment.\n",
    "    \n",
    " </center>\n",
    " </div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}