Scistree2_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scistree2 Tutorial\n",
    "\n",
    "This tutorial consists of two parts:\n",
    "\n",
    "- **Part I** introduces the usage of Scistree2 with a toy example.\n",
    "- **Part II** guides you through the entire analysis of a high-grade serous ovarian carcinoma (HGSOC) dataset, as shown in the paper *\"ScisTree2: An Improved Method for Large-scale Inference of Cell Lineage Trees and Genotype Calling from Noisy Single Cell Data.\"*, accepted to be presented in the RECOMB 2025.\n",
    "\n",
    "## Resources:\n",
    "\n",
    "- Scistree2 binary file: [GitHub - ScisTree2](https://github.com/yufengwudcs/ScisTree2)\n",
    "- HGSOC dataset: [Zenodo - HGSOC Dataset](https://zenodo.org/records/3445364)\n",
    "- Other utilities: Get `popgen` from [GitHub - popgen](https://github.com/haotianzh/popgen)\n",
    "- Paper: https://www.biorxiv.org/content/10.1101/2024.11.08.622704v1\n",
    "\n",
    "## Setup Instructions:\n",
    "\n",
    "Once you download the dataset, place the uncompressed files, Scistree2 binary, `popgen` source folder, and this tutorial in the same directory.\n",
    "\n",
    "This dataset has also been analyzed using HUNTRESS. If you're interested, you can compare our results with theirs.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part I: Toy example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here, we provide a Python interface for conveniently calling Scistree2.\n",
    "import multiprocessing as mp\n",
    "import uuid\n",
    "import numpy as np\n",
    "import subprocess as sp\n",
    "import popgen\n",
    "\n",
    "\"\"\"\n",
    "Scistree2 Caller.\n",
    "\n",
    "Args:\n",
    "    bin_path: Path to binary executable file.\n",
    "    threads: Number of threads in use. Default uses all.\n",
    "    nj: Call NJ only, serves when M is very big.\n",
    "    spr: Enable SPR local search. Default is on.\n",
    "    nni: Enable NNI local search. Default is off.\n",
    "    iterative: Enable iterative optimization. Default is off.                                                                           \n",
    "    verbose: Show outputs.\n",
    "\"\"\"\n",
    "class ScisTree2():\n",
    "    def __init__(self, bin_path='/home/haz19024/softwares/Scistree2/ScisTree2/ScisTree2-source-code/scistree ',\n",
    "                 threads=-1, nj=False, spr=True, nni=False, iterative=False, verbose=True):\n",
    "        self.nj = nj\n",
    "        self.spr = spr\n",
    "        self.nni = nni\n",
    "        self.cmd = self.build_cmd(bin_path, threads, nj, nni, verbose)\n",
    "        \n",
    "    def build_cmd(self, bin_path, threads, nj, nni, verbose):\n",
    "        if threads == -1:\n",
    "            threads = mp.cpu_count()\n",
    "        cmd = [bin_path, '-T', str(threads)]\n",
    "        if verbose:\n",
    "            cmd.append('-v')\n",
    "        if nj:\n",
    "            cmd.append('-n')\n",
    "        if nni:\n",
    "            cmd.append('-q')\n",
    "        return cmd \n",
    "    \n",
    "    \"\"\"\n",
    "    IO: Write the genotype matrix to a file.\n",
    "    \"\"\"\n",
    "    def write_to_scistree(self, genotype_matrix):\n",
    "        nsite, ncell = genotype_matrix.shape\n",
    "        prefix = uuid.uuid4()\n",
    "        output = f'{prefix}.scistree.out'\n",
    "        with open(output, 'w') as out:\n",
    "            out.write(f'HAPLOID {nsite} {ncell}')\n",
    "            for i in range(ncell):\n",
    "                out.write(f' {i}')\n",
    "            out.write('\\n')\n",
    "            for i in range(nsite):\n",
    "                out.write(f's{i}')\n",
    "                for j in range(ncell):\n",
    "                    prob = genotype_matrix[i, j]\n",
    "                    out.write(f' {prob:.5f}')\n",
    "                out.write('\\n')\n",
    "        return output\n",
    "    \n",
    "    \"\"\"\n",
    "    Get the genotype matrix from the outputs of Scistree2.\n",
    "    \"\"\"\n",
    "    def get_scistree_genotype(self, prefix):\n",
    "        geno_file = f'{prefix}.genos.imp'\n",
    "        genotypes = []\n",
    "        with open(geno_file, 'r') as f:\n",
    "            for line in f.readlines():\n",
    "                if line.startswith('Site'):\n",
    "                    line = line.strip()\n",
    "                    genos = line.split('\\t')[1].split()\n",
    "                    genos = list(map(int, genos))\n",
    "                    genotypes.append(genos)\n",
    "        return np.array(genotypes)\n",
    "    \n",
    "    \"\"\"\n",
    "    Run Scistree2 local search.\n",
    "\n",
    "    Args:\n",
    "        geno [numpy.ndarray]: Genotype probability matrix.\n",
    "\n",
    "    Returns:\n",
    "        imputed_genotype [numpy.ndarry]: Imputed genotype. 0 as wild type, 1 as mutation.\n",
    "        nwk [str]: Newick string of the optimal tree.\n",
    "        ml [float]: Log Likelihood of the optimal tree.\n",
    "        time [int]: Running time in seconds.\n",
    "    \"\"\"\n",
    "    def infer(self, geno):\n",
    "        output = self.write_to_scistree(geno)\n",
    "        cmd = self.cmd + [f'{output}']\n",
    "        cmd = ' '.join(cmd)\n",
    "        res = sp.run(cmd, shell=True, stdout=sp.PIPE, encoding='utf-8').stdout.strip().split('\\n')\n",
    "        if not self.nj:\n",
    "            time = float(res[-1].split('=')[1].split('seconds')[0])\n",
    "            nwk = res[-2].split(':')[1].strip() + ';'\n",
    "            ml = float(res[-4].split(',')[0].split(':')[1])\n",
    "            imp_geno = self.get_scistree_genotype(output)\n",
    "            return imp_geno, nwk, ml, time\n",
    "        else:\n",
    "            time = float(res[-1].split('=')[1].split('seconds')[0])\n",
    "            nwk = res[-2].split(':')[1].strip() + ';'\n",
    "            return nwk, time\n",
    "    \n",
    "    \"\"\"\n",
    "    Evaluate a tree given genotype probabilities and return the imputated genotype.\n",
    "\n",
    "    Args:\n",
    "        geno [numpy.ndarray]: Genotype probability matrix.\n",
    "        nwk [str]: Newick string of the optimal tree.\n",
    "\n",
    "    Returns:\n",
    "        imputed_genotype [numpy.ndarry]: Imputed genotype. 0 as wild type, 1 as mutation.\n",
    "        ml [float]: Log likelihood of the optimal tree.\n",
    "    \"\"\"\n",
    "    @staticmethod\n",
    "    def evaluate(geno, nwk):\n",
    "        assert isinstance(nwk, str), 'tree should be a newick string.'\n",
    "        tree = popgen.utils.relabel(popgen.utils.from_newick(nwk), offset=-1)\n",
    "        traveror = popgen.utils.TraversalGenerator(order='post')\n",
    "        max_mls = np.zeros(geno.shape[0])\n",
    "        max_ml_nodes = [None] * geno.shape[0]\n",
    "        g = np.log(1-geno) - np.log(geno) # in log space to avoid numerical overflow\n",
    "        for node in traveror(tree):\n",
    "            if node.is_leaf():\n",
    "                likelihood = g[:, int(node.name)]\n",
    "            else:\n",
    "                likelihood = node.get_children()[0].likelihood + node.get_children()[1].likelihood\n",
    "            for i, l in enumerate(likelihood):\n",
    "                if l > max_mls[i]:\n",
    "                    max_mls[i] = l\n",
    "                    max_ml_nodes[i] = node\n",
    "            node.likelihood = likelihood\n",
    "        max_mls += np.log(geno).sum(axis=1)\n",
    "        imputed_genotype = np.zeros_like(geno, dtype=int)\n",
    "        for i, ml_node in enumerate(max_ml_nodes):\n",
    "            inds = [int(leaf.name) for leaf in ml_node.get_leaves()]\n",
    "            imputed_genotype[i, inds] = 1\n",
    "        return imputed_genotype, sum(max_mls)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "imputed_genotype: \n",
      " [[1 0 1 0 0]\n",
      " [0 1 0 1 0]\n",
      " [1 0 1 0 0]\n",
      " [0 0 0 0 1]\n",
      " [1 0 1 0 0]\n",
      " [1 1 1 1 0]]\n",
      "newick of inferred tree:  (((1,3),(2,4)),5);\n",
      "maximum likelihood:  -6.27126\n",
      "time:  0.0\n"
     ]
    }
   ],
   "source": [
    "\"\"\" \n",
    "Next, we invoke Scistree2 with SPR local search on a toy example. The input should be formatted with rows representing SNPs and columns representing cells.\n",
    "# Toy example:\n",
    "#       cell1   cell2   cell3   cell4   cell5\n",
    "# snp1  0.01    0.6     0.08    0.8     0.7\n",
    "# snp2  0.8     0.02    0.7     0.01    0.3\n",
    "# snp3  0.02    0.8     0.02    0.8     0.9\n",
    "# snp4  0.9     0.9     0.8     0.8     0.02\n",
    "# snp5  0.01    0.8     0.01    0.8     0.9\n",
    "# snp6  0.05    0.02    0.7     0.05    0.9\n",
    "\"\"\"\n",
    "genotype_probability = np.array([[0.01, 0.6, 0.08, 0.8, 0.7],\n",
    "                                 [0.8, 0.02, 0.7, 0.01, 0.3],\n",
    "                                 [0.02, 0.8, 0.02, 0.8, 0.9],\n",
    "                                 [0.9, 0.9, 0.8, 0.8, 0.02],\n",
    "                                 [0.01, 0.8, 0.01, 0.8, 0.9],\n",
    "                                 [0.05, 0.02, 0.7, 0.05, 0.9]]) \n",
    "caller = ScisTree2(threads=8) # use 8 threads\n",
    "imputed_genotype, tree, ml, time = caller.infer(genotype_probability) # run Scistree2 inference\n",
    "print('imputed_genotype: \\n', imputed_genotype)\n",
    "print('newick of inferred tree: ', tree)\n",
    "print('maximum likelihood: ', ml)\n",
    "print('time: ', time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "imputed_genotype: \n",
      " [[1 0 1 0 0]\n",
      " [0 1 0 1 0]\n",
      " [1 0 1 0 0]\n",
      " [0 0 0 0 1]\n",
      " [1 0 1 0 0]\n",
      " [1 1 1 1 0]]\n",
      "newick of inferred tree:  (((1,3),(2,4)),5);\n",
      "maximum likelihood:  -6.27126\n",
      "time:  0.0\n"
     ]
    }
   ],
   "source": [
    "\"\"\" \n",
    "One can also replace SPR local search with NNI local search.  \n",
    "Note: There is no difference between SPR search and NNI search.\n",
    "\"\"\"\n",
    "caller_nni = ScisTree2(threads=8, nni=True)\n",
    "imputed_genotype, tree, ml, time = caller_nni.infer(genotype_probability)\n",
    "print('imputed_genotype: \\n', imputed_genotype)\n",
    "print('newick of inferred tree: ', tree)\n",
    "print('maximum likelihood: ', ml)\n",
    "print('time: ', time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "imputed_genotype: \n",
      " [[1 0 1 0 0]\n",
      " [0 1 0 1 0]\n",
      " [1 0 1 0 0]\n",
      " [0 0 0 0 1]\n",
      " [1 0 1 0 0]\n",
      " [1 1 1 1 0]]\n",
      "newick of inferred tree:  (((1,3),(2,4)),5);\n",
      "maximum likelihood:  -6.271255186813891\n",
      "time:  0.0\n"
     ]
    }
   ],
   "source": [
    "\"\"\" \n",
    "You can also invoke Scistree2 with Neighbor Joining to obtain only the initial tree.\n",
    "In NJ mode, no further steps occur after invoking NJ, meaning there are no outputs for either the imputed genotype or the likelihood.\n",
    "Note: There is no difference between the initial tree and the optimal tree.\n",
    "\"\"\"\n",
    "caller_nj = ScisTree2(threads=8, nj=True)\n",
    "tree, time = caller_nj.infer(genotype_probability)\n",
    "imputed_genotype, ml = ScisTree2.evaluate(genotype_probability, tree)\n",
    "print('imputed_genotype: \\n', imputed_genotype)\n",
    "print('newick of inferred tree: ', tree)\n",
    "print('maximum likelihood: ', ml)\n",
    "print('time: ', time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   ca53fe9761                                      \n",
      " ┌─────┴───────────────────────┐                   \n",
      " 5                         33c374f53d              \n",
      "                   ┌───────────┴───────────┐       \n",
      "               a3348a947d              84d0f52026  \n",
      "             ┌─────┴─────┐           ┌─────┴─────┐ \n",
      "             4           2           3           1 \n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Visualize the optimal tree.\n",
    "\"\"\"\n",
    "t = popgen.utils.from_newick(tree)\n",
    "t.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tree: ((((1,2),3),4),5); \n",
      "ml: -10.835603378281727\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "One can also evaluate a tree using genotype probabilities. Here, we use the same example as before to evaluate another tree,\n",
    "we can observe that the maximum likelihood (ML) of this tree is lower than the optimal one (-6.27), as expected.\n",
    "\"\"\"\n",
    "nwk = '((((1,2),3),4),5);'\n",
    "_, ml = ScisTree2.evaluate(genotype_probability, nwk)\n",
    "print('tree:', nwk, '\\nml:', ml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part II: HGSOC\n",
    "\n",
    "## Inputs preparetion and Excution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/haz19024/.venv/work/lib/python3.10/site-packages/matplotlib/projections/__init__.py:63: UserWarning: Unable to import Axes3D. This may be due to multiple versions of Matplotlib being installed (e.g. as a system package and as a pip package). As a result, the 3D projection is not available.\n",
      "  warnings.warn(\"Unable to import Axes3D. This may be due to multiple versions of \"\n"
     ]
    }
   ],
   "source": [
    "# load some libraries for data manipulation and visualization\n",
    "import pandas as pd \n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "from datetime import date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/haz19024/.local/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3460: DtypeWarning: Columns (0) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Load SNV counts from a CSV file\n",
    "\"\"\"\n",
    "df = pd.read_csv('./ov2295_snv_counts.csv') # load SNV counts \n",
    "df = df.astype({'chrom': 'str', \n",
    "                'coord': 'int32', \n",
    "                'ref': 'str',\n",
    "                'alt': 'str',\n",
    "                'ref_counts': 'int32',\n",
    "                'alt_counts': 'int32'})\n",
    "df['snp_id'] = df['chrom'] + ':' + df['coord'].astype(str) + ':' + df['ref'] + ':' + df['alt'] # create a unique identifier for each SNV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Calculate the posterior probability of each SNV being a true positive.\n",
    "\n",
    "We use GATK likelihood with ADO included as described in CellCoal manual (https://dapogon.github.io/cellcoal/cellcoal.manual.v1.1.html#537_genotype_likelihoods)\n",
    "    P(D|G={g_1, g_2}) = (1-ado)\\prod_{i=1}^{r}P(b_i|G={g_1, g_2}) + 0.5*ado[\\prod_{i=1}^{r}P(b_i|G={g_1}) + \\prod_{i=1}^{r}P(b_i|G={g_2})]\n",
    "\"\"\"\n",
    "\n",
    "def posterior_prob_with_ado_seqerr_gt(ref_counts, alt_counts, prior_ref=0.5, ado=0.2, seqerr=0.01, margin=1e-5):\n",
    "    p00, p01, p10, p11 = np.log(1-seqerr), np.log(seqerr), np.log(seqerr), np.log(1-seqerr)\n",
    "    l00 = np.exp(ref_counts*p00 + alt_counts*p01)\n",
    "    l01 = (1-ado)*np.exp(ref_counts*np.log(0.5*np.exp(p00)+0.5*np.exp(p10))+alt_counts*np.log(0.5*np.exp(p01)+0.5*np.exp(p11)))+(0.5*ado)*(np.exp(ref_counts*p00+alt_counts*p10)+np.exp(ref_counts*p10+alt_counts*p11))\n",
    "    l11 = np.exp(ref_counts*p10 + alt_counts*p11)\n",
    "    g00 = prior_ref**2 * l00\n",
    "    g01 = 2*(1-prior_ref)*prior_ref * l01\n",
    "    g11 = (1-prior_ref)**2 * l11\n",
    "    prob =  max(min(np.round(g00/(g00+g01+g11), 5), 1-margin), margin)\n",
    "    return prob\n",
    "\n",
    "\"\"\"\n",
    "Calculate the phred-scaled likelihood of each SNV, this is required by CellPhy GT10 model.\n",
    "\"\"\"\n",
    "def phred_likelihood_with_ado_seqerr_gt(ref_counts, alt_counts, ado=0.2, seqerr=0.01):\n",
    "    # Q-phred score = -10 * log10(p)\n",
    "    p00, p01, p10, p11 = np.log(1-seqerr), np.log(seqerr), np.log(seqerr), np.log(1-seqerr)\n",
    "    l00 = np.exp(ref_counts*p00 + alt_counts*p01)\n",
    "    l01 = (1-ado)*np.exp(ref_counts*np.log(0.5*np.exp(p00)+0.5*np.exp(p10))+alt_counts*np.log(0.5*np.exp(p01)+0.5*np.exp(p11)))+(0.5*ado)*(np.exp(ref_counts*p00+alt_counts*p10)+np.exp(ref_counts*p10+alt_counts*p11))\n",
    "    l11 = np.exp(ref_counts*p10 + alt_counts*p11)\n",
    "    q00 = int(-10 * np.log10(l00))\n",
    "    q11 = int(-10 * np.log10(l11))\n",
    "    q01 = int(-10 * np.log10(l01))\n",
    "    return q00, q01, q11\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SA921 616\n",
      "SA922 588\n",
      "SA1090 739\n"
     ]
    }
   ],
   "source": [
    "# There are 3 samples [SA921, SA922, SA1090]. Check the number of cells for each sample\n",
    "for sample in df['sample_id'].unique():\n",
    "    print(sample, df[df['sample_id'] == sample]['cell_id'].unique().shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional: according to preclustered clones, we filter out the cells that are not in the clone clusters.\n",
    "# Note: You can skip this step if you don't have clone information, since Scistree2 is a poweful tool designed to be run without preclustering. \n",
    "df_clones = pd.read_csv('./ov2295_clone_clusters.csv')\n",
    "filtered_cells = df_clones['cell_id'].unique()\n",
    "df = df[df['cell_id'].isin(filtered_cells)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter out non-informative sites, criterion: more than 650 cells are missed\n",
    "informative_df = df.groupby('snp_id').size().sort_values(ascending=False)\n",
    "informative_sites = informative_df[informative_df > len(filtered_cells) - 650].index.tolist()\n",
    "df = df[df['snp_id'].isin(informative_sites)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "891 789\n"
     ]
    }
   ],
   "source": [
    "# Preview on processed data\n",
    "cells = df['cell_id'].drop_duplicates().tolist()\n",
    "sites = df['snp_id'].drop_duplicates().tolist()\n",
    "cells = sorted(cells)\n",
    "sites = sorted(sites)\n",
    "print(len(cells), len(sites))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "891 789\n",
      "ov2295_789_891\n"
     ]
    }
   ],
   "source": [
    "# Check the number of cells and the number of sites again, there are 891 cells with 744 mutations in HUNTRESS, a slight difference.\n",
    "cells = df['cell_id'].drop_duplicates().tolist()\n",
    "sites = df['snp_id'].drop_duplicates().tolist()\n",
    "print(len(cells), len(sites))\n",
    "prefix = f'ov2295_{len(sites)}_{len(cells)}'\n",
    "print(prefix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dataframe of meta information for each SNP.\n",
    "df_site_meta = df[['snp_id', 'chrom', 'coord', 'ref', 'alt']].drop_duplicates()\n",
    "df_site_meta.sort_values(by='snp_id', inplace=True)\n",
    "df_site_meta['id'] = [f's{ix}' for ix in np.arange(df_site_meta.shape[0])]\n",
    "df_site_meta.to_csv(f'./{prefix}_meta_site.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dataframe of meta information for each cell.\n",
    "df_cell_meta = df[['cell_id', 'sample_id']].drop_duplicates().sort_values(by='cell_id')\n",
    "df_cell_meta['id'] = [ix for ix in np.arange(df_cell_meta.shape[0])]\n",
    "df_cell_meta = pd.merge(df_cell_meta, df_clones, on='cell_id')\n",
    "df_cell_meta.to_csv(f'./{prefix}_meta_cell.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Get the maximum-likelihood genotype for each site.\n",
    "\"\"\"\n",
    "def get_ml_gt(ref_counts, alt_counts):\n",
    "    a, b, c = phred_likelihood_with_ado_seqerr_gt(ref_counts, alt_counts)\n",
    "    min_phred = min(a, b, c) \n",
    "    if a == min_phred:\n",
    "        return 0\n",
    "    elif b == min_phred:\n",
    "        return 1\n",
    "    else:\n",
    "        return 2\n",
    "\n",
    "df['ml_gt'] = df.apply(lambda x: get_ml_gt(x['ref_counts'], x['alt_counts']), axis=1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Get the phred likelihood for each site.\n",
    "\"\"\"\n",
    "def get_phred_likelihood(ref_counts, alt_counts, ml_gt):\n",
    "    gts = ['0/0', '0/1', '1/1']\n",
    "    a, b, c = phred_likelihood_with_ado_seqerr_gt(ref_counts, alt_counts)\n",
    "    return f'{gts[ml_gt]}:{a},{b},{c}'\n",
    "\n",
    "df['phred'] = df.apply(lambda x: get_phred_likelihood(x['ref_counts'], x['alt_counts'], x['ml_gt']), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    " Calculate the allele frequency for each site.\n",
    "\"\"\"\n",
    "def af(gts):\n",
    "    counts = [0, 0, 0]\n",
    "    for gt in gts:\n",
    "        counts[gt] += 1\n",
    "    counts = np.array(counts)\n",
    "    freqs = counts / sum(counts)\n",
    "    return np.round(freqs[0] + 0.5 * freqs[1], 2)\n",
    "\n",
    "allele_freqs = df.groupby('snp_id')['ml_gt'].agg(func=af).to_frame().rename(columns={'ml_gt': 'af'})\n",
    "df = df.join(allele_freqs, on='snp_id')\n",
    "df_site_meta = df_site_meta.join(allele_freqs, on='snp_id')\n",
    "df_site_meta.to_csv('./meta_site.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the posterior probability for each site\n",
    "df['genotype_prob'] = df.apply(lambda x: posterior_prob_with_ado_seqerr_gt(x['ref_counts'], x['alt_counts'], prior_ref=x['af']), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "IO: Write the genotype matrix to a file.\n",
    "\"\"\"\n",
    "def write_to_scistree(genotype_matrix):\n",
    "    nsite, ncell = genotype_matrix.shape\n",
    "    output = f'ov2295_{nsite}_{ncell}.scistree'\n",
    "    with open(output, 'w') as out:\n",
    "        out.write(f'HAPLOID {nsite} {ncell}')\n",
    "        for i in range(ncell):\n",
    "            out.write(f' {i}')\n",
    "        out.write('\\n')\n",
    "        # for i in range(self.nsite):\n",
    "        # idx = 1\n",
    "        for i in range(nsite):\n",
    "            out.write(f's{i}')\n",
    "            # idx += 1\n",
    "            for j in range(ncell):\n",
    "                prob = genotype_matrix[i, j]\n",
    "                out.write(f' {prob:.5f}')\n",
    "            out.write('\\n')\n",
    "\n",
    "\"\"\"\n",
    "IO: Write the genotype matrix to a file.\n",
    "\"\"\"\n",
    "def write_to_huntress(matrix):\n",
    "    flags = [0, 1, 1, 3]\n",
    "    nsite, ncell = matrix.shape\n",
    "    output = f'ov2295_{nsite}_{ncell}.huntress'\n",
    "    with open(output, 'w') as out:\n",
    "        out.write(f'cellID/mutID')\n",
    "        for i in range(nsite):\n",
    "            out.write(f'\\ts{i}')\n",
    "        out.write('\\n')\n",
    "        for i in range(ncell):\n",
    "            out.write(f'c{i}')\n",
    "            for j in range(nsite):\n",
    "                val = flags[int(matrix[j, i])]\n",
    "                out.write(f'\\t{val}')\n",
    "            out.write('\\n')\n",
    "\n",
    "\"\"\"\n",
    "IO: Write the genotype matrix to a file.\n",
    "\"\"\"\n",
    "def write_to_cellphy(matrix):\n",
    "    nsite, ncell = matrix.shape\n",
    "    output = f'ov2295_{nsite}_{ncell}.vcf'\n",
    "    header = \\\n",
    "f'''##fileformat=VCFv4.3\n",
    "##fileDate={date.today().isoformat()}\n",
    "##source=ov2295\n",
    "##ncell={ncell}\n",
    "##nsite={nsite}\n",
    "##reference=NONE\n",
    "##contig=<ID=1>\n",
    "##phasing=NO\n",
    "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n",
    "##FORMAT=<ID=PL,Number=G,Type=Integer,Description=\"Phread-scaled genotype likelihoods\">\n",
    "'''\n",
    "    with open(output, 'w') as out:\n",
    "        out.write(header)\n",
    "        out.write('#CHROM\\tPOS\\tID\\tREF\\tALT\\tQUAL\\tFILTER\\tINFO\\tFORMAT')\n",
    "        for i in range(ncell):\n",
    "            out.write(f'\\t{i}')\n",
    "        out.write('\\n')\n",
    "        for i in range(nsite):\n",
    "            chrom, coord, ref, alt, af, snp_id = df_site_meta.iloc[i][['chrom', 'coord', 'ref', 'alt', 'af', 'id']]\n",
    "            out.write(f\"{chrom}\\t{coord}\\t{snp_id}\\t{ref}\\t{alt}\\t.\\tPASS\\tAF={af}\\tGT:PL\")\n",
    "            for j in range(ncell):\n",
    "                out.write(f'\\t{matrix[i, j]}')\n",
    "            out.write('\\n')                               "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare Scistree2's inputs\n",
    "scistree_matrix_df = df.pivot(index='snp_id', columns='cell_id', values='genotype_prob')\n",
    "scistree_matrix = scistree_matrix_df.to_numpy()\n",
    "genotype_matrix_scistree = np.nan_to_num(scistree_matrix, nan=0.5)\n",
    "write_to_scistree(genotype_matrix_scistree) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare HUNTRESS's inputs\n",
    "huntress_matrix_df = df.pivot(index='snp_id', columns='cell_id', values='ml_gt')\n",
    "huntress_matrix = huntress_matrix_df.to_numpy()\n",
    "huntress_matrix = np.nan_to_num(huntress_matrix, nan=3)\n",
    "write_to_huntress(huntress_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare Cellphy's inputs\n",
    "cellphy_matrix_df = df.pivot(index='snp_id', columns='cell_id', values='phred') \n",
    "cellphy_matrix_df.fillna('./.:0,0,0', inplace=True)\n",
    "cellphy_matrix = cellphy_matrix_df.to_numpy()\n",
    "write_to_cellphy(cellphy_matrix)\n",
    "# cellphy_matrix_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run Scistree2 with 30 threads\n",
    "input_scis = f'{prefix}.scistree'\n",
    "! /usr/bin/time -o {prefix}.scistree.time /path/to/scistree -v -T 30 -e -o {prefix}.prob.muttree.gml $input_scis > {prefix}.scistree.log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run HUNTRESS with 30 threads\n",
    "input_huntress = f'{prefix}.huntress'\n",
    "! /usr/bin/time -o {prefix}.huntress.time python /path/to/HUNTRESS/HUNTRESS.py --i $input_huntress --o {prefix}.SC --t 30 > {prefix}.HUNTRESS.log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run Cellphy with 30 threads\n",
    "input_cellphy = f'{prefix}.vcf'\n",
    "! /usr/bin/time -o {prefix}.cellphy.time /path/to/cellphy/cellphy.sh FAST -a -t 30 -r $input_cellphy > {prefix}.cellphy.log 2>&1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results analysis and visulization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt \n",
    "import pickle\n",
    "import popgen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you want to annotate the SNVs, you can use ANNOVAR. Installation can befoudn at https://annovar.openbioinformatics.org/en/latest/user-guide/download/\n",
    "\n",
    "# Annovar annotation\n",
    "humandb = '/path/to/annovar/humandb/humandb/'\n",
    "xref = '/path/to/annovar/example/gene_fullxref.txt'\n",
    "!perl /path/to/annovar/table_annovar.pl {prefix}.vcf $humandb -buildver hg19 -out {prefix} -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f -nastring . -vcfinput -polish -xref $xref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Chr</th>\n",
       "      <th>Start</th>\n",
       "      <th>End</th>\n",
       "      <th>Ref</th>\n",
       "      <th>Alt</th>\n",
       "      <th>Func.refGene</th>\n",
       "      <th>Gene.refGene</th>\n",
       "      <th>GeneDetail.refGene</th>\n",
       "      <th>ExonicFunc.refGene</th>\n",
       "      <th>AAChange.refGene</th>\n",
       "      <th>...</th>\n",
       "      <th>Otherinfo894</th>\n",
       "      <th>Otherinfo895</th>\n",
       "      <th>Otherinfo896</th>\n",
       "      <th>Otherinfo897</th>\n",
       "      <th>Otherinfo898</th>\n",
       "      <th>Otherinfo899</th>\n",
       "      <th>Otherinfo900</th>\n",
       "      <th>Otherinfo901</th>\n",
       "      <th>Otherinfo902</th>\n",
       "      <th>Otherinfo903</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10</td>\n",
       "      <td>36258731</td>\n",
       "      <td>36258731</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>PCAT5;ANKRD30A</td>\n",
       "      <td>dist=168883;dist=1156054</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>10</td>\n",
       "      <td>44790486</td>\n",
       "      <td>44790486</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>downstream</td>\n",
       "      <td>LINC02881</td>\n",
       "      <td>dist=389</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>1/1:40,10,8</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>1/1:20,3,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>10</td>\n",
       "      <td>47079935</td>\n",
       "      <td>47079935</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>GPRIN2;NPY4R2</td>\n",
       "      <td>dist=72054;dist=1885</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>1/1:20,3,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10</td>\n",
       "      <td>50851414</td>\n",
       "      <td>50851414</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>intronic</td>\n",
       "      <td>CHAT</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>10</td>\n",
       "      <td>79932700</td>\n",
       "      <td>79932700</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>intergenic</td>\n",
       "      <td>RPS24;LINC00856</td>\n",
       "      <td>dist=116136;dist=75682</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>784</th>\n",
       "      <td>8</td>\n",
       "      <td>97876720</td>\n",
       "      <td>97876720</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>intronic</td>\n",
       "      <td>CPQ</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>1/1:40,10,8</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>1/1:20,3,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>785</th>\n",
       "      <td>8</td>\n",
       "      <td>98325993</td>\n",
       "      <td>98325993</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>ncRNA_intronic</td>\n",
       "      <td>LOC101927066</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>0/0:0,6,20</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>786</th>\n",
       "      <td>8</td>\n",
       "      <td>98413399</td>\n",
       "      <td>98413399</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>ncRNA_intronic</td>\n",
       "      <td>LOC101927066</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/1:20,6,7</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>1/1:40,7,1</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>1/1:20,3,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>1/1:20,3,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>787</th>\n",
       "      <td>8</td>\n",
       "      <td>99166875</td>\n",
       "      <td>99166875</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>intronic</td>\n",
       "      <td>POP1</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>0/1:20,6,7</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>788</th>\n",
       "      <td>8</td>\n",
       "      <td>99770344</td>\n",
       "      <td>99770344</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>intronic</td>\n",
       "      <td>STK3</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>.</td>\n",
       "      <td>...</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "      <td>0/0:0,4,13</td>\n",
       "      <td>./.:0,0,0</td>\n",
       "      <td>0/0:0,2,6</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>789 rows × 972 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Chr     Start       End Ref Alt    Func.refGene     Gene.refGene  \\\n",
       "0     10  36258731  36258731   C   T      intergenic   PCAT5;ANKRD30A   \n",
       "1     10  44790486  44790486   C   A      downstream        LINC02881   \n",
       "2     10  47079935  47079935   G   T      intergenic    GPRIN2;NPY4R2   \n",
       "3     10  50851414  50851414   C   A        intronic             CHAT   \n",
       "4     10  79932700  79932700   A   T      intergenic  RPS24;LINC00856   \n",
       "..   ...       ...       ...  ..  ..             ...              ...   \n",
       "784    8  97876720  97876720   T   C        intronic              CPQ   \n",
       "785    8  98325993  98325993   A   T  ncRNA_intronic     LOC101927066   \n",
       "786    8  98413399  98413399   C   T  ncRNA_intronic     LOC101927066   \n",
       "787    8  99166875  99166875   G   A        intronic             POP1   \n",
       "788    8  99770344  99770344   A   C        intronic             STK3   \n",
       "\n",
       "           GeneDetail.refGene ExonicFunc.refGene AAChange.refGene  ...  \\\n",
       "0    dist=168883;dist=1156054                  .                .  ...   \n",
       "1                    dist=389                  .                .  ...   \n",
       "2        dist=72054;dist=1885                  .                .  ...   \n",
       "3                           .                  .                .  ...   \n",
       "4      dist=116136;dist=75682                  .                .  ...   \n",
       "..                        ...                ...              ...  ...   \n",
       "784                         .                  .                .  ...   \n",
       "785                         .                  .                .  ...   \n",
       "786                         .                  .                .  ...   \n",
       "787                         .                  .                .  ...   \n",
       "788                         .                  .                .  ...   \n",
       "\n",
       "    Otherinfo894 Otherinfo895 Otherinfo896 Otherinfo897 Otherinfo898  \\\n",
       "0      ./.:0,0,0   0/0:0,4,13    0/0:0,2,6   0/0:0,4,13    ./.:0,0,0   \n",
       "1      ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0   \n",
       "2      ./.:0,0,0   0/0:0,4,13   1/1:20,3,0    ./.:0,0,0    ./.:0,0,0   \n",
       "3      ./.:0,0,0   0/0:0,4,13    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0   \n",
       "4      0/0:0,2,6    ./.:0,0,0    ./.:0,0,0   0/0:0,4,13    ./.:0,0,0   \n",
       "..           ...          ...          ...          ...          ...   \n",
       "784  1/1:40,10,8    ./.:0,0,0    ./.:0,0,0   0/0:0,4,13    ./.:0,0,0   \n",
       "785   0/0:0,6,20    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0   0/0:0,4,13   \n",
       "786    ./.:0,0,0   0/1:20,6,7    0/0:0,2,6   1/1:40,7,1    ./.:0,0,0   \n",
       "787    ./.:0,0,0    0/0:0,2,6   0/0:0,4,13   0/1:20,6,7    ./.:0,0,0   \n",
       "788    ./.:0,0,0    ./.:0,0,0   0/0:0,4,13   0/0:0,4,13    ./.:0,0,0   \n",
       "\n",
       "    Otherinfo899 Otherinfo900 Otherinfo901 Otherinfo902 Otherinfo903  \n",
       "0      0/0:0,2,6    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0  \n",
       "1    1/1:40,10,8    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0   1/1:20,3,0  \n",
       "2      ./.:0,0,0    ./.:0,0,0    0/0:0,2,6    ./.:0,0,0    ./.:0,0,0  \n",
       "3      ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0  \n",
       "4      ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0  \n",
       "..           ...          ...          ...          ...          ...  \n",
       "784    ./.:0,0,0    ./.:0,0,0   1/1:20,3,0    ./.:0,0,0    ./.:0,0,0  \n",
       "785    0/0:0,2,6    ./.:0,0,0    0/0:0,2,6    0/0:0,2,6    ./.:0,0,0  \n",
       "786    ./.:0,0,0   1/1:20,3,0    ./.:0,0,0    ./.:0,0,0   1/1:20,3,0  \n",
       "787    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0    ./.:0,0,0  \n",
       "788    ./.:0,0,0    0/0:0,2,6   0/0:0,4,13    ./.:0,0,0    0/0:0,2,6  \n",
       "\n",
       "[789 rows x 972 columns]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# SNV annotation\n",
    "df_gene_annotation = pd.read_table(f'../{prefix}.hg19_multianno.txt')\n",
    "df_gene_annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>variant_id</th>\n",
       "      <th>node</th>\n",
       "      <th>loss</th>\n",
       "      <th>origin</th>\n",
       "      <th>presence</th>\n",
       "      <th>ml_origin</th>\n",
       "      <th>ml_presence</th>\n",
       "      <th>ml_loss</th>\n",
       "      <th>snp_id</th>\n",
       "      <th>chrom</th>\n",
       "      <th>coord</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>id</th>\n",
       "      <th>af</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>4.567698e-05</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>10</td>\n",
       "      <td>36258731</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>s0</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>1</td>\n",
       "      <td>4.536726e-05</td>\n",
       "      <td>6.020528e-98</td>\n",
       "      <td>3.097139e-07</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>10</td>\n",
       "      <td>36258731</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>s0</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>2</td>\n",
       "      <td>3.076417e-07</td>\n",
       "      <td>5.940235e-98</td>\n",
       "      <td>2.072204e-09</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>10</td>\n",
       "      <td>36258731</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>s0</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>3</td>\n",
       "      <td>2.072204e-09</td>\n",
       "      <td>3.062727e-106</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>10</td>\n",
       "      <td>36258731</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>s0</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>4</td>\n",
       "      <td>2.072204e-09</td>\n",
       "      <td>7.672154e-107</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>10</td>\n",
       "      <td>36258731</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>s0</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13408</th>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>12</td>\n",
       "      <td>3.417474e-07</td>\n",
       "      <td>2.554902e-65</td>\n",
       "      <td>2.345168e-09</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>8</td>\n",
       "      <td>99770344</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>s788</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13409</th>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>13</td>\n",
       "      <td>2.323599e-09</td>\n",
       "      <td>3.455902e-65</td>\n",
       "      <td>2.156852e-11</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>8</td>\n",
       "      <td>99770344</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>s788</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13410</th>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>14</td>\n",
       "      <td>2.128003e-11</td>\n",
       "      <td>5.044326e-65</td>\n",
       "      <td>2.877698e-13</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>8</td>\n",
       "      <td>99770344</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>s788</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13411</th>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>15</td>\n",
       "      <td>1.586328e-11</td>\n",
       "      <td>1.338978e-63</td>\n",
       "      <td>5.705658e-12</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>8</td>\n",
       "      <td>99770344</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>s788</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13412</th>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>16</td>\n",
       "      <td>2.323299e-09</td>\n",
       "      <td>3.504392e-65</td>\n",
       "      <td>2.186873e-11</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>8</td>\n",
       "      <td>99770344</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>s788</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>13413 rows × 15 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            variant_id  node          loss         origin      presence  \\\n",
       "0      10:36258731:C:T     0  0.000000e+00   4.567698e-05  0.000000e+00   \n",
       "1      10:36258731:C:T     1  4.536726e-05   6.020528e-98  3.097139e-07   \n",
       "2      10:36258731:C:T     2  3.076417e-07   5.940235e-98  2.072204e-09   \n",
       "3      10:36258731:C:T     3  2.072204e-09  3.062727e-106  0.000000e+00   \n",
       "4      10:36258731:C:T     4  2.072204e-09  7.672154e-107  0.000000e+00   \n",
       "...                ...   ...           ...            ...           ...   \n",
       "13408   8:99770344:A:C    12  3.417474e-07   2.554902e-65  2.345168e-09   \n",
       "13409   8:99770344:A:C    13  2.323599e-09   3.455902e-65  2.156852e-11   \n",
       "13410   8:99770344:A:C    14  2.128003e-11   5.044326e-65  2.877698e-13   \n",
       "13411   8:99770344:A:C    15  1.586328e-11   1.338978e-63  5.705658e-12   \n",
       "13412   8:99770344:A:C    16  2.323299e-09   3.504392e-65  2.186873e-11   \n",
       "\n",
       "       ml_origin  ml_presence  ml_loss           snp_id chrom     coord ref  \\\n",
       "0              0            0        0  10:36258731:C:T    10  36258731   C   \n",
       "1              0            0        0  10:36258731:C:T    10  36258731   C   \n",
       "2              0            0        0  10:36258731:C:T    10  36258731   C   \n",
       "3              0            0        0  10:36258731:C:T    10  36258731   C   \n",
       "4              0            0        0  10:36258731:C:T    10  36258731   C   \n",
       "...          ...          ...      ...              ...   ...       ...  ..   \n",
       "13408          0            0        0   8:99770344:A:C     8  99770344   A   \n",
       "13409          0            0        0   8:99770344:A:C     8  99770344   A   \n",
       "13410          0            0        0   8:99770344:A:C     8  99770344   A   \n",
       "13411          0            0        0   8:99770344:A:C     8  99770344   A   \n",
       "13412          0            0        0   8:99770344:A:C     8  99770344   A   \n",
       "\n",
       "      alt    id    af  \n",
       "0       T    s0  0.91  \n",
       "1       T    s0  0.91  \n",
       "2       T    s0  0.91  \n",
       "3       T    s0  0.91  \n",
       "4       T    s0  0.91  \n",
       "...    ..   ...   ...  \n",
       "13408   C  s788  0.93  \n",
       "13409   C  s788  0.93  \n",
       "13410   C  s788  0.93  \n",
       "13411   C  s788  0.93  \n",
       "13412   C  s788  0.93  \n",
       "\n",
       "[13413 rows x 15 columns]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the precomputed clonal tree with node information in the HGSOC dataset\n",
    "with open('./ov2295_tree.pickle', 'rb') as f:\n",
    "    tree = pickle.load(f)\n",
    "df_nodes = pd.read_csv('./ov2295_nodes.csv')\n",
    "df_mutations = pd.merge(df_nodes, df_site_meta, left_on='variant_id', right_on='snp_id') # merge the node information with the SNV information\n",
    "df_mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],\n",
       " 1: [2, 3, 4, 5, 6, 7],\n",
       " 2: [3, 4],\n",
       " 3: [],\n",
       " 4: [],\n",
       " 5: [6, 7],\n",
       " 6: [],\n",
       " 7: [],\n",
       " 8: [9, 10, 11, 12, 13, 14, 15, 16],\n",
       " 9: [10, 11],\n",
       " 10: [],\n",
       " 11: [],\n",
       " 12: [13, 14, 15, 16],\n",
       " 13: [14, 15],\n",
       " 14: [],\n",
       " 15: [],\n",
       " 16: []}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Obtain the Ancestor - Descendent relationships in the given HGSOC clonal tree\n",
    "orders_nodes = {}\n",
    "for node in tree.nodes:\n",
    "    orders_nodes[node.label] = []\n",
    "    for d in node.descendents:\n",
    "        orders_nodes[node.label].append(d.label)\n",
    "orders_nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>node</th>\n",
       "      <th>loss</th>\n",
       "      <th>origin</th>\n",
       "      <th>presence</th>\n",
       "      <th>ml_origin</th>\n",
       "      <th>ml_presence</th>\n",
       "      <th>ml_loss</th>\n",
       "      <th>snp_id</th>\n",
       "      <th>chrom</th>\n",
       "      <th>coord</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>id</th>\n",
       "      <th>af</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>variant_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>10:36258731:C:T</th>\n",
       "      <td>12</td>\n",
       "      <td>6.020910e-98</td>\n",
       "      <td>0.993219</td>\n",
       "      <td>0.006781</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>10:36258731:C:T</td>\n",
       "      <td>10</td>\n",
       "      <td>36258731</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>s0</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10:44790486:C:A</th>\n",
       "      <td>0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>10:44790486:C:A</td>\n",
       "      <td>10</td>\n",
       "      <td>44790486</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>s1</td>\n",
       "      <td>0.41</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10:47079935:G:T</th>\n",
       "      <td>0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>10:47079935:G:T</td>\n",
       "      <td>10</td>\n",
       "      <td>47079935</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>s2</td>\n",
       "      <td>0.77</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10:50851414:C:A</th>\n",
       "      <td>1</td>\n",
       "      <td>1.655850e-176</td>\n",
       "      <td>0.993264</td>\n",
       "      <td>0.006736</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>10:50851414:C:A</td>\n",
       "      <td>10</td>\n",
       "      <td>50851414</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>s3</td>\n",
       "      <td>0.83</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10:79932700:A:T</th>\n",
       "      <td>0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>10:79932700:A:T</td>\n",
       "      <td>10</td>\n",
       "      <td>79932700</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>s4</td>\n",
       "      <td>0.55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8:97876720:T:C</th>\n",
       "      <td>9</td>\n",
       "      <td>1.347097e-137</td>\n",
       "      <td>0.993095</td>\n",
       "      <td>0.006905</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>8:97876720:T:C</td>\n",
       "      <td>8</td>\n",
       "      <td>97876720</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>s784</td>\n",
       "      <td>0.86</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8:98325993:A:T</th>\n",
       "      <td>12</td>\n",
       "      <td>4.809655e-113</td>\n",
       "      <td>0.993210</td>\n",
       "      <td>0.006790</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>8:98325993:A:T</td>\n",
       "      <td>8</td>\n",
       "      <td>98325993</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>s785</td>\n",
       "      <td>0.87</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8:98413399:C:T</th>\n",
       "      <td>9</td>\n",
       "      <td>7.527466e-144</td>\n",
       "      <td>0.993211</td>\n",
       "      <td>0.006789</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>8:98413399:C:T</td>\n",
       "      <td>8</td>\n",
       "      <td>98413399</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>s786</td>\n",
       "      <td>0.86</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8:99166875:G:A</th>\n",
       "      <td>9</td>\n",
       "      <td>4.995262e-76</td>\n",
       "      <td>0.993197</td>\n",
       "      <td>0.006803</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>8:99166875:G:A</td>\n",
       "      <td>8</td>\n",
       "      <td>99166875</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>s787</td>\n",
       "      <td>0.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8:99770344:A:C</th>\n",
       "      <td>2</td>\n",
       "      <td>2.541829e-65</td>\n",
       "      <td>0.993219</td>\n",
       "      <td>0.006781</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>8:99770344:A:C</td>\n",
       "      <td>8</td>\n",
       "      <td>99770344</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>s788</td>\n",
       "      <td>0.93</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>789 rows × 14 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 node           loss    origin  presence  ml_origin  \\\n",
       "variant_id                                                            \n",
       "10:36258731:C:T    12   6.020910e-98  0.993219  0.006781          1   \n",
       "10:44790486:C:A     0   0.000000e+00  1.000000  0.000000          1   \n",
       "10:47079935:G:T     0   0.000000e+00  1.000000  0.000000          1   \n",
       "10:50851414:C:A     1  1.655850e-176  0.993264  0.006736          1   \n",
       "10:79932700:A:T     0   0.000000e+00  1.000000  0.000000          1   \n",
       "...               ...            ...       ...       ...        ...   \n",
       "8:97876720:T:C      9  1.347097e-137  0.993095  0.006905          1   \n",
       "8:98325993:A:T     12  4.809655e-113  0.993210  0.006790          1   \n",
       "8:98413399:C:T      9  7.527466e-144  0.993211  0.006789          1   \n",
       "8:99166875:G:A      9   4.995262e-76  0.993197  0.006803          1   \n",
       "8:99770344:A:C      2   2.541829e-65  0.993219  0.006781          1   \n",
       "\n",
       "                 ml_presence  ml_loss           snp_id chrom     coord ref  \\\n",
       "variant_id                                                                   \n",
       "10:36258731:C:T            1        0  10:36258731:C:T    10  36258731   C   \n",
       "10:44790486:C:A            1        0  10:44790486:C:A    10  44790486   C   \n",
       "10:47079935:G:T            1        0  10:47079935:G:T    10  47079935   G   \n",
       "10:50851414:C:A            1        0  10:50851414:C:A    10  50851414   C   \n",
       "10:79932700:A:T            1        0  10:79932700:A:T    10  79932700   A   \n",
       "...                      ...      ...              ...   ...       ...  ..   \n",
       "8:97876720:T:C             1        0   8:97876720:T:C     8  97876720   T   \n",
       "8:98325993:A:T             1        0   8:98325993:A:T     8  98325993   A   \n",
       "8:98413399:C:T             1        0   8:98413399:C:T     8  98413399   C   \n",
       "8:99166875:G:A             1        0   8:99166875:G:A     8  99166875   G   \n",
       "8:99770344:A:C             1        0   8:99770344:A:C     8  99770344   A   \n",
       "\n",
       "                alt    id    af  \n",
       "variant_id                       \n",
       "10:36258731:C:T   T    s0  0.91  \n",
       "10:44790486:C:A   A    s1  0.41  \n",
       "10:47079935:G:T   T    s2  0.77  \n",
       "10:50851414:C:A   A    s3  0.83  \n",
       "10:79932700:A:T   T    s4  0.55  \n",
       "...              ..   ...   ...  \n",
       "8:97876720:T:C    C  s784  0.86  \n",
       "8:98325993:A:T    T  s785  0.87  \n",
       "8:98413399:C:T    T  s786  0.86  \n",
       "8:99166875:G:A    A  s787  0.91  \n",
       "8:99770344:A:C    C  s788  0.93  \n",
       "\n",
       "[789 rows x 14 columns]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Double check that each mutation occurs exactly once\n",
    "(df_mutations.groupby(['variant_id']).sum()['ml_origin'] == 1).all()\n",
    "# Find the node that each mutation originates from\n",
    "df_origin = df_mutations.groupby(['variant_id']).apply(lambda x: x[x['ml_origin']==1])\n",
    "df_origin = df_origin.set_index('variant_id')\n",
    "df_origin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "\"\"\"\n",
    "Get the genotype matrix from the output of Scistree2\n",
    "\"\"\"\n",
    "def get_scistree_genotype(prefix):\n",
    "    geno_file = f'{prefix}.scistree.genos.imp'\n",
    "    genotypes = []\n",
    "    with open(geno_file, 'r') as f:\n",
    "        for line in f.readlines():\n",
    "            if line.startswith('Site'):\n",
    "                line = line.strip()\n",
    "                genos = line.split('\\t')[1].split()\n",
    "                genos = list(map(int, genos))\n",
    "                genotypes.append(genos)\n",
    "    return np.array(genotypes)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Get the genotype matrix from the output of HUNTRESS\n",
    "\"\"\"\n",
    "def get_huntress_genotype(prefix):\n",
    "    geno_file = f'{prefix}.SC.CFMatrix'\n",
    "    genotypes = pd.read_table(geno_file, index_col=0).values.T\n",
    "    return genotypes\n",
    "\n",
    "\"\"\"\n",
    "Get the mutation list from the output of Cellphy\n",
    "\"\"\"\n",
    "def read_mutation_list(file):\n",
    "    edges = {}\n",
    "    with open(file, 'r') as f:\n",
    "        for line in f.readlines():\n",
    "            items = line.strip().split()\n",
    "            if int(items[1]) > 0:\n",
    "                eid, num_muts, muts = items\n",
    "                edges[eid] = [int(_[1:]) for _ in muts.split(',')]\n",
    "            else:\n",
    "                eid = items[0]\n",
    "                edges[eid] = []\n",
    "    return edges\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Get the mutation tree from the output of Cellphy\n",
    "\"\"\"\n",
    "def read_mutation_tree(file):\n",
    "    with open(file, 'r') as f:\n",
    "        line = f.readline()\n",
    "        line = line.strip()\n",
    "    return popgen.utils.from_newick(line)\n",
    "\n",
    "\"\"\"\n",
    "Get the genotype matrix from the output of Cellphy\n",
    "\"\"\"\n",
    "def get_cellphy_genotype(prefix):\n",
    "    raxml_tree = f'{prefix}.vcf.Mapped.raxml.mutationMapTree'\n",
    "    raxml_list = f'{prefix}.vcf.Mapped.raxml.mutationMapList'\n",
    "    num_snp, num_cell = len(sites), len(cells)\n",
    "    genotype = np.zeros([num_snp, num_cell], dtype=int)\n",
    "    edges = read_mutation_list(raxml_list)\n",
    "    tree = read_mutation_tree(raxml_tree)\n",
    "    generator = popgen.utils.TraversalGenerator(order='post')\n",
    "    for node in generator(tree):\n",
    "        if not node.is_root():\n",
    "            eid = re.search(r'\\[(\\d+)\\]', node.branch).group(1)\n",
    "            muts = edges[eid]\n",
    "            # muts = [mut-1 for mut in muts] # to 0-indexed\n",
    "            leaves = [int(leaf.name) for leaf in node.get_leaves()] # to 0-indexed\n",
    "            genotype[np.ix_(muts, leaves)] = 1\n",
    "    return genotype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _assign(ids, x):\n",
    "    # print(x)\n",
    "    chr = x['chrom']\n",
    "    coord = x['coord']\n",
    "    k = f'{chr}:{coord}'\n",
    "    ids[k] = x['id']\n",
    "\n",
    "# conversion for mutationlist of cellphy\n",
    "def convert_to_id(mut_file, out_path=f'{prefix}.vcf.Mapped.raxml.mutationMapList'):\n",
    "    cellphy_ids = {}\n",
    "    df_origin.apply(lambda x: _assign(cellphy_ids, x), axis=1)\n",
    "    with open(out_path, 'w') as out:\n",
    "        with open(mut_file, 'r') as f:\n",
    "            for line in f.readlines():\n",
    "                items = line.strip().split()\n",
    "                ix, num = items[0], int(items[1])\n",
    "                if num != 0:\n",
    "                    muts = [cellphy_ids[mut] for mut in items[2].split(',')]\n",
    "                    muts = ','.join(muts)\n",
    "                else:\n",
    "                    muts = ''\n",
    "                new_line = f'{ix}\\t{num}\\t{muts}\\n'\n",
    "                out.write(new_line)\n",
    "    \n",
    "# convert_to_id(f'{prefix}.vcf.Mapped.raxml.mutationMapList.swp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(789, 891) (789, 891) (789, 891)\n"
     ]
    }
   ],
   "source": [
    "# grab all inferred genotypes\n",
    "cellphy_genotype = get_cellphy_genotype(prefix)\n",
    "huntress_genotype = get_huntress_genotype(prefix)\n",
    "scistree_genotype = get_scistree_genotype(prefix)\n",
    "print(cellphy_genotype.shape, huntress_genotype.shape, scistree_genotype.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "sites = sorted(sites)\n",
    "# Get all Ancestor-Descendent pairs in the clonal tree in the HGSOC dataset\n",
    "def precede(s1, s2):\n",
    "    n1, n2 = df_origin.loc[s1, 'node'], df_origin.loc[s2, 'node']\n",
    "    if n2 in orders_nodes[n1]:\n",
    "        return 0\n",
    "    if n1 in orders_nodes[n2]:\n",
    "        return 1\n",
    "    return 2\n",
    "\n",
    "true_mutations = {f's{_}': [] for _ in range(len(sites))}\n",
    "for i in range(len(sites)):\n",
    "    for j in range(i+1, len(sites)):\n",
    "        res = precede(sites[i], sites[j])\n",
    "        if res == 0:\n",
    "            true_mutations[f's{i}'].append(f's{j}')\n",
    "        elif res == 1:\n",
    "            true_mutations[f's{j}'].append(f's{i}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Check if clade1 covers clade2.\n",
    "\"\"\"\n",
    "def is_covered(clade1, clade2):\n",
    "    diff = clade1 - clade2\n",
    "    if 0 in diff and 1 in diff and -1 in diff:\n",
    "        return 0 # no relationship\n",
    "    if 1 in diff:\n",
    "        return 1 # clade1(snp1) happens before clade2(snp2)\n",
    "    else:\n",
    "        return -1 # clade2(snp2) happens before clade1(snp1)\n",
    "\n",
    "\"\"\"\n",
    "Get AD pairs from given genotype matrix.\n",
    "\"\"\"\n",
    "def get_ancestor_descendant_pairs(geno):\n",
    "    num_snp, _ = geno.shape\n",
    "    mutations = {f's{_}': [] for _ in range(num_snp)}\n",
    "    for i in range(num_snp):\n",
    "        for j in range(i+1, num_snp):\n",
    "            code = is_covered(geno[i], geno[j])\n",
    "            if code == 1:\n",
    "                mutations[f's{i}'].append(f's{j}')\n",
    "            elif code == -1:\n",
    "                mutations[f's{j}'].append(f's{i}')\n",
    "    return mutations\n",
    "\n",
    "\"\"\"\n",
    "Calculate the error rate for Ancestor-Descendant pairs.\n",
    "\"\"\"\n",
    "def ancestor_descendant_error(mutation1, mutation2):\n",
    "    count = 0\n",
    "    total = 0\n",
    "    for snp1 in mutation1:\n",
    "        for snp2 in mutation1[snp1]:\n",
    "            count += int(snp2 not in mutation2[snp1])\n",
    "            total += 1\n",
    "    return count / total\n",
    "\n",
    "\"\"\"\n",
    "Calculate the error rate for Different-Lineage pairs.\n",
    "\"\"\"\n",
    "def different_lineage_error(mutation1, mutation2):\n",
    "    count = 0\n",
    "    total = 0\n",
    "    muts = list(mutation1.keys())\n",
    "    num_mut = len(muts)\n",
    "    for i in range(num_mut):\n",
    "        for j in range(i+1, num_mut):\n",
    "            if muts[i] not in mutation1[muts[j]] and muts[j] not in mutation1[muts[i]]:\n",
    "                count += int(muts[j] in mutation2[muts[i]]) + int(muts[i] in mutation2[muts[j]])\n",
    "                total += 1\n",
    "    return count / total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "scistree_mutations = get_ancestor_descendant_pairs(scistree_genotype)\n",
    "huntress_mutations = get_ancestor_descendant_pairs(huntress_genotype)\n",
    "cellphy_mutations = get_ancestor_descendant_pairs(cellphy_genotype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ancestor-Descendent accuracy and precision\n",
    "scistree_ad_acc = 1 - ancestor_descendant_error(true_mutations, scistree_mutations)\n",
    "scistree_ad_pre = 1 - ancestor_descendant_error(scistree_mutations, true_mutations)\n",
    "huntress_ad_acc = 1 - ancestor_descendant_error(true_mutations, huntress_mutations)\n",
    "huntress_ad_pre = 1 - ancestor_descendant_error(huntress_mutations, true_mutations)\n",
    "cellphy_ad_acc = 1 - ancestor_descendant_error(true_mutations, cellphy_mutations)\n",
    "cellphy_ad_pre = 1 - ancestor_descendant_error(cellphy_mutations, true_mutations)\n",
    "\n",
    "\n",
    "# Different-Lineage accuracy and precision\n",
    "scistree_dl_acc = 1 - different_lineage_error(true_mutations, scistree_mutations)\n",
    "scistree_dl_pre = 1 - different_lineage_error(scistree_mutations, true_mutations)\n",
    "huntress_dl_acc = 1 - different_lineage_error(true_mutations, huntress_mutations)\n",
    "huntress_dl_pre = 1 - different_lineage_error(huntress_mutations, true_mutations)\n",
    "cellphy_dl_acc = 1 - different_lineage_error(true_mutations, cellphy_mutations)\n",
    "cellphy_dl_pre = 1 - different_lineage_error(cellphy_mutations, true_mutations)\n",
    "\n",
    "\n",
    "def f1_score(pre, acc):\n",
    "    return 2 / (1/pre + 1/acc)\n",
    "\n",
    "# F1 score\n",
    "scsitree_ad_f1 = f1_score(scistree_ad_pre, scistree_ad_acc)\n",
    "huntress_ad_f1 = f1_score(huntress_ad_pre, huntress_ad_acc)\n",
    "cellphy_ad_f1 = f1_score(cellphy_ad_pre, cellphy_ad_acc)\n",
    "scistree_dl_f1 = f1_score(scistree_dl_pre, scistree_dl_acc)\n",
    "huntress_dl_f1 = f1_score(huntress_dl_pre, huntress_dl_acc)\n",
    "cellphy_dl_f1 = f1_score(cellphy_dl_pre, cellphy_dl_acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Method</th>\n",
       "      <th>Measure</th>\n",
       "      <th>Value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>Ancestor-Descendent Accuracy</td>\n",
       "      <td>0.937556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>Ancestor-Descendent Precision</td>\n",
       "      <td>0.674161</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>Different-Lineage Accuracy</td>\n",
       "      <td>0.733573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>Different-Lineage Precision</td>\n",
       "      <td>0.956050</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>Ancestor-Descendent F1-Score</td>\n",
       "      <td>0.784335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>Different-Lineage F1-Score</td>\n",
       "      <td>0.830164</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>Ancestor-Descendent Accuracy</td>\n",
       "      <td>0.990460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>Ancestor-Descendent Precision</td>\n",
       "      <td>0.619685</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>Different-Lineage Accuracy</td>\n",
       "      <td>0.641850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>Different-Lineage Precision</td>\n",
       "      <td>0.997233</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>Ancestor-Descendent F1-Score</td>\n",
       "      <td>0.762383</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>Different-Lineage F1-Score</td>\n",
       "      <td>0.781014</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>Ancestor-Descendent Accuracy</td>\n",
       "      <td>0.171243</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>Ancestor-Descendent Precision</td>\n",
       "      <td>0.143840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>Different-Lineage Accuracy</td>\n",
       "      <td>0.419005</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>Different-Lineage Precision</td>\n",
       "      <td>0.472635</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>Ancestor-Descendent F1-Score</td>\n",
       "      <td>0.156349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>Different-Lineage F1-Score</td>\n",
       "      <td>0.444207</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Method                        Measure     Value\n",
       "0   ScisTree   Ancestor-Descendent Accuracy  0.937556\n",
       "1   ScisTree  Ancestor-Descendent Precision  0.674161\n",
       "2   ScisTree     Different-Lineage Accuracy  0.733573\n",
       "3   ScisTree    Different-Lineage Precision  0.956050\n",
       "4   ScisTree   Ancestor-Descendent F1-Score  0.784335\n",
       "5   ScisTree     Different-Lineage F1-Score  0.830164\n",
       "6   HUNTRESS   Ancestor-Descendent Accuracy  0.990460\n",
       "7   HUNTRESS  Ancestor-Descendent Precision  0.619685\n",
       "8   HUNTRESS     Different-Lineage Accuracy  0.641850\n",
       "9   HUNTRESS    Different-Lineage Precision  0.997233\n",
       "10  HUNTRESS   Ancestor-Descendent F1-Score  0.762383\n",
       "11  HUNTRESS     Different-Lineage F1-Score  0.781014\n",
       "12   CellPhy   Ancestor-Descendent Accuracy  0.171243\n",
       "13   CellPhy  Ancestor-Descendent Precision  0.143840\n",
       "14   CellPhy     Different-Lineage Accuracy  0.419005\n",
       "15   CellPhy    Different-Lineage Precision  0.472635\n",
       "16   CellPhy   Ancestor-Descendent F1-Score  0.156349\n",
       "17   CellPhy     Different-Lineage F1-Score  0.444207"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "methods = ['ScisTree'] * 6 + ['HUNTRESS'] * 6   + ['CellPhy'] * 6\n",
    "measures = ['Ancestor-Descendent Accuracy',\n",
    "            'Ancestor-Descendent Precision',\n",
    "            'Different-Lineage Accuracy',\n",
    "            'Different-Lineage Precision',\n",
    "            'Ancestor-Descendent F1-Score',\n",
    "            'Different-Lineage F1-Score'] * 3\n",
    "vals = [scistree_ad_acc, scistree_ad_pre, scistree_dl_acc, scistree_dl_pre, scsitree_ad_f1, scistree_dl_f1,\n",
    "        huntress_ad_acc, huntress_ad_pre, huntress_dl_acc, huntress_dl_pre, huntress_ad_f1, huntress_dl_f1,\n",
    "        cellphy_ad_acc, cellphy_ad_pre, cellphy_dl_acc, cellphy_dl_pre, cellphy_ad_f1, cellphy_dl_f1\n",
    "        ]\n",
    "df_result = pd.DataFrame({'Method': methods, 'Measure': measures, 'Value': vals})\n",
    "df_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAEYCAYAAABBfQDEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnuUlEQVR4nO3debwXdd3//8eLg4AKiiJqAYoK5nqJiphLqbkvuZdahpWh9lO7LOObXuZGuWXpRUaZe2q5pGmY5NJVeBm5Qa6ACxcUoKiEu2iAvX5/zBwa1vMBPucc4Dzut9u58ZmZ98y853M+vM/z8573zERmIkmSJKnQrrUrIEmSJC1PDMiSJElShQFZkiRJqjAgS5IkSRUGZEmSJKnCgCxJkiRVGJCl5VREfDki/tza9ZC0YomIKyPi7Mr01yPitYh4LyK6RcQuEfFSOX1oK1Z1iUXE2IjYvbXroZWfAVlExMiIeDMiOrZyHb5Wx+31jogs/wC8V/5x+F1E7F2vfaxIyveiT2vXQ9Kyi4i/RcQHEfFuRLwVEX+JiJMioh1AZp6Umd8ry64CXAbsk5mdM3MGMAT4STl9dwvX/YaI+H4TZc6LiJsXtiwzt8zMkc1SOanCgNzGRURv4FNAAge3bm2WThQW9VnumpmdgW2AB4G7IuLLLVY5SWoen83MLsCGwMXAd4BrF1JuPaATMLYyb8P5pmsWEe2XZj1pRWNA1kDgUeAG4LjGmeW3/GERcW/ZS/FYRGxSWb5lRDwYEW+UvbP/Vc5vFxFnRMT/RcSMiLg9ItYul3WKiJvL+W9FxBMRsV5EXEAR0n9S9vb+pCy/c1nm7fLfnSv7HxkRF0TEKGAmsPHiDjIzX83MocB5wCWNgToiPh4Rd0bE9IiYFBHfqOxjQESMjoh3ymO8rLJs17LX5q2ImNIYuiOiY0T8MCIml+tcGRGrlst2j4ipEXF6RLweEdMi4iuVbXaLiOHl/h4H5r7f5fLNKu/5CxHx+Vp+XxHxv2Wxp8v396jFvVeSVhyZ+XZmDgeOAo6LiK0ae2kjYlPghbLoWxHxx4j4P4r28p6yPegYEWtGxLVlm/RyuW4DzB3qNSoiLo+IGcB5S9vORcQJwBeB/1fu+54lPd4oes/3Kl+fV/6NubFs98ZGRP9K2aba90fKNnxaRPwkIjpUlu9TtrNvR8RPI+KhqJzljIivRsT4KM6+3h8RGy7psWj5ZkDWQOCX5c++EbFeZdnRwPnAWsAE4AKAiOgC/AG4D/g40Af4n3KdU4FDgd3KZW8Cw8plxwFrAr2AbsBJwAeZeRbwMHBKecrvlChC9b3Aj8uylwH3RkS3Sv2+BJwAdAH+XuPx/gZYF/hEGZLvAZ4GegB7AqdFxL5l2aHA0MxcgyKs3l4e/4bA74ErgO5AP+Cpcp2LgU3LeX3K7Z5T2f/65XvQAzgeGBYRa5XLhgEfAh8Dvlr+UO5zdYoe8F+V9T8a+GlEbFHZ9kJ/X5n56XL5NuX7e1uN75WkFURmPg5MpehsaJz3IrBlOdk1Mz+TmZsAkyl6oDtn5j8pOkjmULRZ2wL7ANUhbzsCEyl6oy9gKdu5zLyK4m/ND8p9f7YOh34wcCvQFRgONHawNNW+fwR8E1gH2Klc/v+V664D3AGcSfH35wWg2kFzCPBfwOEUfwMeBm6pw7FoOWJAbsMiYleKU223Z+YY4P+AL1SK3JWZj2fmHIpGrV85/yDg1cz8UWZ+mJnvZuZj5bKTgLMyc2rZ8J4HHBnFabnZFI1Nn8z8KDPHZOY7i6jegcBLmXlTZs7JzFuA54Fqg3pDZo4tl8+u8bBfKf9dG9gB6J6ZQzJzVmZOBK6mCJqU9e0TEetk5nuZ+Wg5/wvAHzLzlsycnZkzMvOpiAiKwP7NzHwjM98FLqxsr3GbQ8r1RgDvUYT1BuAI4JzMfD8znwN+UVnvIOBvmXl9ebxPAncCn6uUWdTvS1Lb8ApF21azslPkAOC0su15HbicedutVzLzirJt+ZClbOeW+qgW78+ZOSIzPwJuohhOB0207+Xfn0fL9vRvwM8pOnageD/GZuZvymP+MfBqZZ8nARdl5vhy+YVAP3uRVy6OJWrbjgMeyMx/lNO/KuddXk5XG4SZQOfydS+KML0wG1KM8/1XZd5HFD0PN5Xr3hoRXYGbKcL0wsLtx1mwV/jvFD0BjaY0voiI9yrzq72q82tc/w1ga+DjEfFWZXkDRW8AFD0fQ4DnI2IScH5m/o5FH393YDVgTJGVi6qV22w0o2xQGzW+r90p/j9OqSyrHv+GwI7z1bU9xXvaaFG/L0ltQw+Ktm1JbAisAkyrtFvtmLctqr5elnZuARHxRYpwCvBwZu6/hPWfv93rVHbIbMhi2vdy+MllQP/yeNoDY8pyH6dyzJmZETG1sp0NgaER8aPqoVC8/7WezdRyzoDcRpXjxT4PNEREYwPTEegaEdssek2gaDiOXsyyr2bmqEUsPx84P4qLA0dQnLq6luIiwapXKBqhqg0ohnU0mrtOeSHeXOX2F+Yw4PVyv12BSZnZd2EFM/Ml4JjyVN3hwB3lEI8pwICFrPIP4ANgy8x8eRH7X5TpFKc4e1H0lENxvI2mAA9lZpu8C4ekxYuIHSgC2p8phkTUagrwT2Cd+UJtVbV9XpZ2bv5tkZmNQ/zqbQqLad+BnwFPAsdk5rsRcRpwZLlsGtCzsWB5drBnZd0pwAVl3bWScohF23UoRc/uFhSn4vsBm1N8ux7YxLq/Az4WEaeVF2t0iYjGBvlK4ILGU00R0b0cr0VE7BERW5fDCd6hOA3X2NP8GvNeaDcC2DQivhAR7aO4sGyLct9LLIqLAU8BzgXOzMx/AY8D70bEdyJi1YhoiOIClx3KdY6NiO5l2bfKTf2LojHfKyI+X9atW0T0K8tdDVweEeuW2+hRGfO2SOXpwd9QXACzWjm2+LhKkd+V78eXImKV8meHiNi8xrdg/vdX0kogItaIiIMoxuHenJnPLsn6mTkNeAD4UbmtdhGxSUTstojyS93OlWpti9pFcWF348+S3oZ0se07xbUr7wDvRcRmwNcr694LbB0Rh5a90SdTjKtudCVwZkRsCRDFRY7V4W5aCRiQ267jgOszc3J5h4dXM/NVigscvshizi6UY872phgP/CrwErBHuXgoxYUSD0TEuxR3yGgMz+tTXPjwDjAeeIh/DxEYSjFW+c2I+HF5r86DgNOBGcD/Aw6qDAep1VsR8T7wLMW4ss9l5nXlcXxU7qMfMImiZ+QaiotLAPYDxpbDN4YCR2fmB5k5udzW6RSnM5/i3+PevkNxgdyjEfEOxcWMtY69O4XiNOSrFBfNXN+4oHzP96HouX+lLHMJRa9/Lc4DfhHFFdufb6qwpOXePWUbOwU4i2K4wFcWv8oiDQQ6AOMoLqy+g+Ji4UVZlnbuWmCLsi26ezHljqHoqW78WdSwvoWqoX3/NsX1JO9SBP7bKuv+g+L6jh9Q/P3ZAhhN0dNOZt5F0f7eWh7/c8CSDg3Rci4y5z+zLUmSJJh7R4ypwBcz80+tXR+1DHuQJUmSKiJi34joWg7t+C+Ki/AebWI1rUQMyJIkSfPaiWJYxz8ohhMempkftG6V1JIcYiFJkiRV2IMsSZIkVbTafZDXWWed7N27d2vtXpJa1ZgxY/6Rmd1bux62xZLaskW1xa0WkHv37s3o0aNba/eS1KoiYrl44pZtsaS2bFFtsUMsJEmSpAoDsiRJklRhQJYkSZIqDMiSJElSRZMBOSKui4jXI+K5RSyPiPhxREyIiGciYrv6V1OSJElqGbX0IN8A7LeY5fsDfcufE4CfLXu1JEmSpNbRZEDOzP8F3lhMkUOAG7PwKNA1Ij5WrwpKkjybJ0ktqR5jkHsAUyrTU8t5C4iIEyJidESMnj59eh12LUltxg14Nk+SWkSLXqSXmVdlZv/M7N+9e6s/QEqSVhiezZOkllOPgPwy0Ksy3bOcJ0lqOZ7Nk6Q6qcejpocDp0TErcCOwNuZOa0O211pTB6ydWtXoWYbnPNsa1dBUjPLzKuAqwD69++frVwdqU0wC6xYmgzIEXELsDuwTkRMBc4FVgHIzCuBEcABwARgJvCV5qqsJGmRPJsnSXXSZEDOzGOaWJ7AyXWrkSRpaXg2T5LqpB5DLCRJzcyzeZLUcgzIkrQC8GyeJLWcFr3NmyRJkrS8MyBLkiRJFQZkSZIkqcKALEmSJFV4kZ4kSVohbT/4xtauQs3u6tLaNdCSsAdZkiRJqlhhe5D91ihJkqTmYA+yJEmSVGFAliRJkioMyJIkSVKFAVmSJEmqMCBLkiRJFQZkSZIkqcKALEmSJFUYkCVJkqQKA7IkSZJUYUCWJEmSKgzIkiRJUoUBWZIkSaowIEuSJEkVBmRJkiSpwoAsSZIkVRiQJUmSpAoDsiRJklRhQJYkSZIqDMiSJElShQFZkiRJqjAgS5IkSRUGZEmSJKnCgCxJkiRVGJAlSZKkCgOyJEmSVGFAliRJkioMyJIkSVJFTQE5IvaLiBciYkJEnLGQ5RtExJ8i4smIeCYiDqh/VSVJkqTm12RAjogGYBiwP7AFcExEbDFfse8Ct2fmtsDRwE/rXVFJauvsrJCkllFLD/IAYEJmTszMWcCtwCHzlUlgjfL1msAr9auiJMnOCklqObUE5B7AlMr01HJe1XnAsRExFRgBnLqwDUXECRExOiJGT58+fSmqK0ltlp0VktRC6nWR3jHADZnZEzgAuCkiFth2Zl6Vmf0zs3/37t3rtGtJahPsrJCkFlJLQH4Z6FWZ7lnOqzoeuB0gMx8BOgHr1KOCkqSa2VkhSXVQS0B+AugbERtFRAeKcW3D5yszGdgTICI2pwjIdktIUv3YWSFJLaTJgJyZc4BTgPuB8RQXgIyNiCERcXBZ7HRgUEQ8DdwCfDkzs7kqLUltkJ0VktRC2tdSKDNHUIxnq847p/J6HLBLfasmSWqUmXMiorGzogG4rrGzAhidmcMpOiuujohvUlywZ2eFJC2FmgKyJKn12VkhSS3DR01LkiRJFQZkSZIkqcKALEmSJFUYkCVJkqQKA7IkSZJUYUCWJEmSKgzIkiRJUoUBWZIkSaowIEuSJEkVPklPkrTUth98Y2tXoWZjLh3Y2lWQtIKwB1mSJEmqMCBLkiRJFQZkSZIkqcKALEmSJFUYkCVJkqQKA7IkSZJUYUCWJEmSKgzIkiRJUoUBWZIkSaowIEuSJEkVBmRJkiSpwoAsSZIkVbRv7QpIK6rtB9/Y2lWoyZhLB7Z2FSRJWqHYgyxJkiRVGJAlSZKkCgOyJEmSVGFAliRJkioMyJIkSVKFAVmSJEmqMCBLkiRJFd4HWZLUJkwesnVrV6EmG5zzbGtXQWrz7EGWJEmSKgzIkiRJUoUBWZIkSapwDLK0kltRxl2CYy8lScuHmnqQI2K/iHghIiZExBmLKPP5iBgXEWMj4lf1raYkSZLUMpoMyBHRAAwD9ge2AI6JiC3mK9MXOBPYJTO3BE6rf1Ulqe2yo0KSWk4tQywGABMycyJARNwKHAKMq5QZBAzLzDcBMvP1eldUktqqSkfF3sBU4ImIGJ6Z4yplqh0Vb0bEuq1TW0la8dUyxKIHMKUyPbWcV7UpsGlEjIqIRyNiv4VtKCJOiIjRETF6+vTpS1djSWp75nZUZOYsoLGjosqOCkmqk3rdxaI90BfYHTgGuDoius5fKDOvysz+mdm/e/fuddq1JK306tZRAXZWSFJTagnILwO9KtM9y3lVU4HhmTk7MycBL1IEZklSy6ipowLsrJCkptQSkJ8A+kbERhHRATgaGD5fmbspGmUiYh2KnoyJ9aumJLVpdlRIUgtqMiBn5hzgFOB+YDxwe2aOjYghEXFwWex+YEZEjAP+BAzOzBnNVWlJamPsqJCkFlTTg0IycwQwYr5551ReJ/Ct8keSVEeZOSciGjsqGoDrGjsqgNGZObxctk/ZUfERdlRI0lLzSXqStAKwo0KSWk697mIhSZIkrRQMyJIkSVKFAVmSJEmqMCBLkiRJFQZkSZIkqcKALEmSJFUYkCVJkqQKA7IkSZJUYUCWJEmSKgzIkiRJUoUBWZIkSaowIEuSJEkVBmRJkiSpwoAsSZIkVRiQJUmSpAoDsiRJklRhQJYkSZIqDMiSJElShQFZkiRJqjAgS5IkSRUGZEmSJKnCgCxJkiRVtG/tCkiSJGn5scsVu7R2FWo26tRRzbJde5AlSZKkCgOyJEmSVGFAliRJkioMyJIkSVKFAVmSJEmqMCBLkiRJFQZkSZIkqcKALEmSJFUYkCVJkqQKA7IkSZJUYUCWJEmSKtq3dgUkSdK/7XLFLq1dhZqNOnVUa1dBahY19SBHxH4R8UJETIiIMxZT7oiIyIjoX78qSpLAtliSWkqTATkiGoBhwP7AFsAxEbHFQsp1Af4TeKzelZSkts62WJJaTi09yAOACZk5MTNnAbcChyyk3PeAS4AP61g/SVLBtliSWkgtAbkHMKUyPbWcN1dEbAf0ysx7F7ehiDghIkZHxOjp06cvcWUlqQ2rW1ssSVq8Zb6LRUS0Ay4DTm+qbGZelZn9M7N/9+7dl3XXkqTSkrTFdlZI0uLVEpBfBnpVpnuW8xp1AbYCRkbE34BPAsO9OESS6qpubbGdFZK0eLUE5CeAvhGxUUR0AI4GhjcuzMy3M3OdzOydmb2BR4GDM3N0s9RYktom22JJaiFNBuTMnAOcAtwPjAduz8yxETEkIg5u7gpKkmyLJakl1fSgkMwcAYyYb945iyi7+7JXS5I0P9tiSWoZPmpakiRJqjAgS5IkSRUGZEmSJKnCgCxJkiRVGJAlSZKkCgOyJEmSVGFAliRJkioMyJIkSVKFAVmSJEmqMCBLkiRJFQZkSZIkqcKALEmSJFUYkCVJkqQKA7IkSZJUYUCWJEmSKgzIkiRJUoUBWZIkSaowIEuSJEkVBmRJkiSpwoAsSZIkVRiQJUmSpAoDsiRJklRhQJYkSZIqDMiSJElShQFZkiRJqjAgS5IkSRUGZEmSJKnCgCxJkiRVGJAlSZKkCgOyJEmSVGFAliRJkirat3YFJEmS2rJ/dejCzG0H8dEavYBo7erw/XYrTv/p+PHjayrXqVMnevbsySqrrFJTeQOyJElSK5q57SDW3mgbuq7egYjWD8irNDS0dhVqttl6mzVZJjOZMWMGU6dOZaONNqppuyvOVwRJkqSV0Edr9FpuwvHKKCLo1q0bH374Yc3rGJAlSZJaVRiOm9mSvr8GZEmSJKmipjHIEbEfMBRoAK7JzIvnW/4t4GvAHGA68NXM/Hud6ypJbZbtsNR27Hz5mLpu7y/f3L7JMhcP/Tm33T2ChoZ2zGnXwHmXnsc2222zQLnnnnqO3/76t5x1wVkLLHtx/It855TvADDt5Wl07tKZLmt0oevaXbn+19cv+4G0oCYDckQ0AMOAvYGpwBMRMTwzx1WKPQn0z8yZEfF14AfAUc1RYUlqa2yHJTWnR0c/xe//8L88et+v6dixA3996x1mz5690LJb9duKrfpttdBlm26+KXf9z10AnPmNM9l9793Z97P7zlNmzpw5tG+//N8jopYaDgAmZOZEgIi4FTgEmNswZ+afKuUfBY6tZyUlqY2zHZbUbF59fTrd1u5Kx44dAFir21oAPPvks1x49oV8MPMDOnTowPV3XM/Yp8dy3c+u48qbr+TxvzzORWdfBBRjfG+6+yZW77z6AtsfeNhANttqM/762F858LADGbDzAC4+92JmzpzJWmuvxYVDL2Td9dZl8t8m870zv8cbM95g1VVXZcgPh7Bx341b7o2oqCUg9wCmVKanAjsupvzxwO8XtiAiTgBOANhggw1qrKIktXl1a4fBtljSvPbabRcuvPxKttr1QD7zqU+yy6EH0K9/P7514re47OeXsfW2W/Peu+/RsVPHeda7/mfXc/ZFZ7PdgO14//336dix4yL2ALNnzeaOB+5g9uzZDDxsIMNuGMba66zNiLtHMPSioVzw3xdw7rfP5dwfnEvvjXvz9F+fZsgZQ7jhzhua+egXrq593BFxLNAf2G1hyzPzKuAqgP79+2c99y1JarodBttiSfPqvPpqPHLf7fz5sTE89JfHOf3E0znxP0+k+3rd2XrbrYsyXTovsN52A7bjknMv4aAjDmLvA/Zm9Y8v2HvcaP9D9gdg0oRJvPT8Sxx/1PEAfPTRR3Rfrzvvv/8+T45+km8O+ubcdWbNmlXPw1witQTkl4Feleme5bx5RMRewFnAbpn5z/pUT5KE7bCkZtbQ0MBuOw9gt50HsO4Wm/Gr63/V5DqDTh3EbnvtxkN/eIgvfPYLXHPrNYscErHaaqvNfd3nE3249d5b51n+3rvv0WWNLnPHMLe2Wm7z9gTQNyI2iogOwNHA8GqBiNgW+DlwcGa+Xv9qSlKbZjssqdm8OGESEyb++6Y348eOZ+O+GzP9tek8++SzALz/3vvMmTNnnvUm/20ym26+KYNOHcTW/bZm4oSJTe6r9ya9eXPGmzw5+kkAZs+ezUvPv0TnLp3puUFP7ht+H1A8/e75sc/X6xCXWJM9yJk5JyJOAe6nuL3QdZk5NiKGAKMzczhwKdAZ+HV5I+bJmXlwM9ZbktoM22Gpbanltmz19N7MmXzruxfy1jvv0r59A+tvtCFDLh3C4UcfzvfP+j7//PCfdOzUket+fd0869141Y08Nuox2rVrR59P9OHTn/l0k/vq0KED/33Nf3PBdy/gvXfeY86cOQw8YSB9N+vLpcMu5fwzzufK/76SObPnsP+h+7PZlk0/Sro51DQGOTNHACPmm3dO5fVeda6XJKnCdlhSc9nuP7Zk5PBfzp2e1NAAFHezuG3EbfOUHbDLAAbsMgCA71743UVu86IfXzT39Y133TjPss232pyb7755gXV6btiTq2+5eskPoBks/zeiU4va5YpdWrsKNRl16qjWroIkSVpJ+ahpSZIkqcKALEmSJFUYkCVJkqQKA7IkSZJUYUCWJEmSKryLhSRJ0nLktWuOruv21vvarU2W6dZ3B2a89MTc6btuvYvnnn6Osy86mzO/cSa77707+35237nLt994e8ZMHMPLk19mrwF7cdb3z+LYrx0LwPfO/B5bbbMVzz71LE8+8SSzZ81m6pSpbLTJRgCcdNpJjHxwJE888gRd1uhCZvKd87/DTp/aCYCBhw1k+uvT6dSpEwAb9N6AodcOZdKESZw7+FzefeddZs2axfY7bs/tN93OzJkzGTRoEM888wyZSdeuXbnvvvvo3HnBx2PXyoAsabmxotxmELzVoCQ16rZON2665iY+P/DzdOjQYe78cy4ubtX+8uSXOelLJ83zGOmRD45k8DmD2fez+/LYnx/jnMHncP8j989dfumwS9mq31bz7OeC717AcScex5777QnAi+NfBGDo0KGst956PPts8dS/F154gVVWWWWZjskhFpIkSVpqa3Vbi09+6pPcffvdS7V+v/79eP3V15ssN/216az/sfXnTm+6+aYATJs2jR49esyd/4lPfIKOHTsuVV0a2YMsSZLUxn3w4T8ZsPcRAMyK4O233maPffaoef2vnfI1TvjCCRxxzBFLvO+H//Tw3F7hRoNPHjx3iMXOn96ZwecO5rgTjuPLR36Zbftvy86778zhRx8O68FXv/pV9tlnH+644w723HNPjjvuOPr27bvE9agyIEuSJLVxq3bqyOMP3gkUj5puHIMMEBELrjDfrF4b9uI/tv0Pfveb39W8z0uHXMrlF13Oa9Ne45bf3TLvsoUMsTj8mMPZdY9defiPD/PH+//I7TfdzvjnxtOvXz8mTpzIAw88wB/+8Ad22GEHHnnkETbffPOa6zI/h1hIkiRpkbqu1ZW333577vRbb77FWmuvtUC5E//zRK75yTVkZk3bHXzOYO77y32cftbpnPXNs2paZ9311+WILxzBsF8Mo6GhgeeeK0J8586dOfzww/npT3/Ksccey4gRI2ra3qIYkCVJkrRIO+y8A7//7e+ZNWsWAHffdjc77rLjAuU27rsxfTbtw8gHRy7R9r94/BfJfyV//tOfF1vu4T8+zOzZswGY/vp03nrzLXr06MGoUaN48803AZg1axbjxo1jww03XKI6zM8hFpIkScuRWm7L1pL22GcPxj0zjiP3OZKGhgZ69e7FeZect9CyJ552IofvdfgSbT8iOOm0k7h22LXsuseuwLxjkLuu3ZXrf309ox4axYVnXzj3ArzB5wxm/fXX54EHHuDrX/86mcm//vUvDjzwQI44YsnHQlcZkCVJktq46j2QAQ47+jAOO/qwudMnf/tkTv72yQus12ODHtzz0D1zpzfbcjPGTRu32DIAF/34onmm9zloH/Y5aB8AbrzrxoXW8Yzzz+CM889YYP7AgQMZOHDgQtdZWg6xkCRJkioMyJIkSVKFAVmSJEmqMCBLkiRJFQZkSZIkqcKALEmSJFV4mzdJkqTlyGfuPKGu2/vjEVc1WebV1//B4HMvZvTTz7HqmmvQbZ1unPm9M9lok40WWn77jbdnzMQxvDz5ZU760knc89A9PD7qcU7+8sn03KAns2bN4oBDDuDkb58897HVZ190dl2PqzkZkCVJktqwzOSo4/+TYz93MDf97IdMamjg+bHPM2P6jEUG5EXZfsftufLmK5n5/kwO2+swdt9n9+apdDNziIUkSVIbNnLU47RfpT2DBh41d95mW25G/0/259ph1/K5fT/HIXscwhU/uKLmba62+mps+R9bMnnSZABef+11Bh0ziH132pdLh1wKwJ2/upMLz75w7jq333w7F51z0UK319IMyJIkSW3YuBdeYrutt1hg/qiRo/j7pL9z+323c9f/3MXYZ8byxCNPLGQLC3rzjTd5eszT9PlEHwCef+55Lvv5ZQz/03B+/9vfM+3laex3yH6MfGAks2fPBuCuW+/iiGOW7RHR9eIQC0mSJC1g1MhRjBo5isP3OhyAme/P5O+T/s4OO+2wyHXGPDaGw/c6nGgXDDp1EH0368tzTz3HJz/1Sbqs0QWAPpv24ZWpr/CxHh9jx113ZOSDI9mk7ybMmT2HTTfftEWOrSkGZEmSpDZs80378Jt7H1xgfmZywjdO4KjK0IumNI5Bnl+HDh3mvm7X0I6P5nwEwJFfPJKrhl7FRn024rCjD1uK2jcPh1hIkiS1YXvsuiOzZs3impt/PXfeC+NeoHOXzvzmlt/w/vvvA/DatNeYMX1GXfe9zXbbMO2Vadx7170ceNiBdd32srAHWZIkaTlSy23Z6ikiuO2aoQw+9xJ+9NNradexIz169eDM751JlzW7cMyBxwDFhXc/GPYDunXvVtf973/w/ox/bjxrdl2zrttdFgZkSZKkNu7j66/LL3/+IwAmNTTMnd97494MHDRwgfJjJo4BoMcGPbjnoXsAGLDLAAbsMmCBsocdfdg8wyfmH4Ix5rExHHficct+EHXkEAtJkiS1uHfefof9dt6PTp06sdOndmrt6szDHmRJkiS1uDXWXIP7/nJfa1djoexBliRJalVJZrZ2JVZqS/r+GpAlSZJaUcM7U3jr/VmG5GaSmcyYMYNOnTrVvI5DLCRJklrRak9ezRsMYvoavYBo7erwj3YrTv9pvlHbl4pOnTrRs2fPmrdrQJYkSWpF7Wa9S+fHLmvtasx10lprtHYVajbq1FHNst2aviJExH4R8UJETIiIMxayvGNE3FYufywiete9ppLUxtkWS1LLaDIgR0QDMAzYH9gCOCYitpiv2PHAm5nZB7gcuKTeFZWktsy2WJJaTi09yAOACZk5MTNnAbcCh8xX5hDgF+XrO4A9I6L1B9FI0srDtliSWkgtY5B7AFMq01OBHRdVJjPnRMTbQDfgH9VCEXECcEI5+V5EvLA0lV7RbNg8m12H+d7ftiS+4d/8WjXT5w/8DC7rJpb0V2NbvIxsi+vPtrh2fv6aR3O1xS16kV5mXgW07APGV1IRMToz+7d2PdR2+RlccdkW14//D9Sa/Pw1n1qGWLwM9KpM9yznLbRMRLQH1gRm1KOCkiTAtliSWkwtAfkJoG9EbBQRHYCjgeHzlRkOHFe+PhL4Y3q3a0mqJ9tiSWohTQ6xKMexnQLcDzQA12Xm2IgYAozOzOHAtcBNETEBeIOi4Vbz8vSoWpufwRZkW7zc8v+BWpOfv2YSdi5IkiRJ/7biPEtQkiRJagEGZEmSJKnCgNxCIuLQiMiI2Ky166K2JSI+ioinIuK5iPh1RKxWh20OiYi9FrP8pIgYuKz7kerJdlitybZ4xeIY5BYSEbcBH6e4qvzcZtpHQ2Z+1Bzb1oorIt7LzM7l618CYzLzssry9pk5p9UqKLUQ22G1JtviFYs9yC0gIjoDuwLHU15VHhENEfHD8pvkMxFxajl/h4j4S0Q8HRGPR0SXiPhyRPyksr3fRcTu5ev3IuJHEfE0sFNEnBMRT5TbvarxMbMR0Sci/lBu968RsUlE3BgRh1a2+8uImP/RtVq5PAz0iYjdI+LhiBgOjCs/j5eWn51nIuLExhUi4jsR8Wz52bm4nHdDRBxZvr44IsaV6/2wnHdeRHy7fN0vIh4tl98VEWuV80dGxCXl5/zFiPhUS78Zajtsh7WcsS1ezrXok/TasEOA+zLzxYiYERHbAwOA3kC/8vZNa0dxb9PbgKMy84mIWAP4oIltrw48lpmnA0TEuMwcUr6+CTgIuAf4JXBxZt4VEZ0ovhxdC3wTuDsi1gR25t/3UNVKJooHR+wP3FfO2g7YKjMnRfHo4bczc4eI6AiMiogHgM0oPr87ZubMiFh7vm12Aw4DNsvMjIiuC9n1jcCpmflQFLckOxc4rVzWPjMHRMQB5fxFniqUlpHtsJYLtsUrBnuQW8YxwK3l61vL6b2AnzeeTsnMN4BPANMy84ly3js1nG75CLizMr1HRDwWEc8CnwG2jIguQI/MvKvc7oeZOTMzH6J48ED3sk53enpnpbRqRDwFjAYmU/xBBng8MyeVr/cBBpblHgO6AX0pPqfXZ+ZMmPs5rXob+BC4NiIOB2ZWF5Z/8LuWnzWAXwCfrhT5TfnvGIqgIjUX22G1NtviFYg9yM2s/Jb3GWDriEiKG/wnxVOxajWHeb/MdKq8/rBxvFvZI/FToH9mTomI8+YruzA3AsdSnHL8yhLUSSuODzKzX3VGecb3/eosip6F++crt+/iNlz2ug0A9qR4ctspFJ/3Wv2z/PcjbI/UTGyHtZywLV6B2IPc/I4EbsrMDTOzd2b2AiYBTwMnlqdaGhvwF4CPRcQO5bwu5fK/Af0iol1E9KI4LbgwjY3wP6IYb3ckQGa+C0xtHOcWER3j31fP3kB5iiUzx9XtqLWiuR/4ekSsAhARm0bE6sCDwFcaPy8LOa3XGVgzM0dQnCbepro8M98G3qyMafsS8BBSy7Id1orCtng54beE5ncMcMl88+4ENqc4xfJMRMwGrs7Mn0TEUcAVEbEqxbi3vYBRFI35OGA88NeF7Sgz34qIq4HngFeZt3fkS8DPy3FHs4HPARMz87WIGA/cXY+D1QrrGorTan+NoktjOnBoZt4XEf2A0RExCxgB/FdlvS7Ab8teswC+tZBtHwdcWTbsE7GHTC3PdlgrCtvi5YS3eWvjyv8ozwLbld8wJUktyHZYWv44xKINi+Lm4uOBK2yUJanl2Q5Lyyd7kCVJkqQKe5AlSZKkCgOyJEmSVGFAliRJkioMyJIkSVKFAVmSJEmq+P8Bk4w9RMCs0PIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plots\n",
    "import seaborn as sns\n",
    "\n",
    "fig, axes = plt.subplots(figsize=(10,4), nrows=1, ncols=2)\n",
    "\n",
    "g1 = sns.barplot(ax=axes[0], data=df_result[df_result['Measure'].isin(['Ancestor-Descendent Accuracy', 'Ancestor-Descendent Precision'])], x='Measure', y='Value', hue='Method')\n",
    "axes[0].set_title('Ancestor-Descendent')\n",
    "axes[0].set_xlabel('')\n",
    "axes[0].set_xticks([0, 1])\n",
    "axes[0].set_xticklabels(['Accuracy', 'Precision'])\n",
    "axes[0].set_ylabel('')\n",
    "g1.legend_.remove()\n",
    "g2 = sns.barplot(ax=axes[1], data=df_result[df_result['Measure'].isin(['Different-Lineage Accuracy', 'Different-Lineage Precision'])], x='Measure', y='Value', hue='Method')\n",
    "axes[1].set_title('Different-Lineage')\n",
    "axes[1].set_xlabel('')\n",
    "axes[1].set_xticks([0, 1])\n",
    "axes[1].set_xticklabels(['Accuracy', 'Precision'])\n",
    "axes[1].set_ylabel('')\n",
    "g2.legend(loc='lower right')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmk0lEQVR4nO3dfZwWdb3/8debBYRAwWQxAxRUBOWgEAuKaJJoYqKYt3CsE1pC9gPvypOmRz3oqYyO93QUC83UxEyTOBwxUyrxJhZF5UYU8QbQckNR8Y67z++PmaWL5drlWtxhZ+H9fDzmsTPf+c53PjPX7H52Zr7XjCICMzOzvGnW2AGYmZkV4wRlZma55ARlZma55ARlZma55ARlZma55ARlZma55ARlTZqkXSX9WdL7kv67seOx4iQdKmlRY8dhTUvzxg7Atj+SXgV2BdYBHwD/B4yNiFVb0Nxo4B/ATuEv9eVWRPwF6NHYcVjT4jMoayzHRkRb4AtABXBJfRZWohmwB7BgS5KTJP+DthV4P9uWcoKyRhURy0nOoP4FQNJBkh6XtFLSs5IGV9eVNFPSf0maBXwI3A58A/h3SaskHSFpB0nXSnojHa6VtEO6/GBJyyR9X9LfgFslXS7pN5LuSC8TPi9pH0kXSXpL0lJJXy6I4XRJC9O6SySNKZhX3f5302XflHR6wfzWkv5b0muS3pX0mKTWm9vumiR1kXSfpCpJKyTdmJY3k3RJ2v5bkm6X1C6d11VSpPEvlfSOpG9L6i/puXS9NxasY5SkWZJuTGN9QdKQeu6Hwv08WNKygjrfl7Q8XX5Rddslfn5F969tgyLCg4etOgCvAkek412A+cAVQCdgBfAVkn+ejkyny9O6M4HXgV4kl6dbALcBVxa0PR54EugIlAOPA1ek8wYDa4GrgB2A1sDlwMfAUWmbtwOvABen7Z8JvFLQ/jHAXoCAw0gS5RdqtD8+XfYr6fyd0/kT023oBJQBB6dx1LndNfZdGfAscA3QBmgFHJLOOwNYDOwJtAXuA36VzusKBHBTusyX0+3+XbqvOgFvAYel9Uel23Jeui2nAu8Cn63Hfijcz4OBZen8HsBS4PMFse1Vj8+v6P71sO0NjR6Ah+1vIElQq4CVwGvAz9I/Yt+v/oNaUHcG8I10fCYwvsb829g4Qb0MfKVg+ijg1XR8MLAaaFUw/3LgDwXTx6axlaXTO6Z/2NvXsi2/A84paP8joHnB/LeAg0gSz0fAAUXaqHO7a5QPBKoK11Ew74/AdwqmewBrSBJv13Q7OhXMXwGcWjD9W+DcdHwU8Aaggvl/Bb5e4n6ouZ8H888EtXe6X44AWtRoZ3OfX9H929jHtIdsBl/is8ZyfES0j4g9IuI7EfERyf2kk9PLTSslrQQOAXYrWG7pZtr9PEnSq/ZaWlatKiI+rrHM3wvGPwL+ERHrCqYhOSNB0tGSnpT0dhrfV4AOBcuviIi1BdMfpst2IDlzeblIzKVsd7UuwGs11lGt2LY3J+mQUtu21pxuWzC9PCIK7+1t2Jcl7Idi+xmAiFgMnEvyz8Fbku6WVP0Zbe7zq23/2jbICcryZCnJmUT7gqFNRPy4oM7mOkO8QfIHv9ruaVmpy9cqvRfyW+CnwK4R0R6YTnKZa3P+QXJJba8i80rZ7sK6u6t4x4Ni276WjZNQfXSSVLhtuwNvlLgf6tzPEXFXRBySxhsklwNr24Y3sO2SE5TlyR3AsZKOklQmqVV6Y7xzPdr4NXCJpHJJHYBL03YbQkuSeypVwFpJR5Pcy9msiFgPTAaulvT5dPsGpn/s67PdfwXeBH4sqU1ad1A679fAeZK6SWoL/BCYUsvZVik6AmdLaiHpZGBfkkS0xfsBQFIPSYen2/4xyZnb+oJtyOrzsybGCcpyIyKWAsOBH5D88VsKXED9jtMrgUrgOeB54Om0rCHiex84G7gHeAf4V2BqPZr4XhrTbOBtkrOGZvXZ7vTS47Ek93FeB5aRdGCAJAH+CvgzSUePj4Fx9dnGGp4CupOc/f0XcFJErGiA/bAD8OO03b+RJMKL0nmZfX7W9GjjS8xmZkk3c+Bb6WU4s0bhMygzM8slJygzM8slX+IzM7Nc8hmUmZnlUpN7iGOHDh2ia9eujR2GmZk1kDlz5vwjIsprlje5BNW1a1cqKysbOwwzM2sgkl4rVu5LfGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktOUGZmlktN7lFHZrZt6HfB7ZmvY86Ef8t8HZYdn0GZmVkuOUGZmVkuOUGZmVkuZZqgJA2VtEjSYkkXFpm/u6RHJT0j6TlJX8kyHjMzazoyS1CSyoCJwNHAfsBISfvVqHYJcE9E9AVGAD/LKh4zM2tasjyDGgAsjoglEbEauBsYXqNOADul4+2ANzKMx8zMmpAsu5l3ApYWTC8DDqxR53LgIUnjgDbAERnGY2ZmTUhjd5IYCdwWEZ2BrwC/krRJTJJGS6qUVFlVVbXVgzQzs60vywS1HOhSMN05LSv0TeAegIh4AmgFdKjZUERMioiKiKgoLy/PKFwzM8uTLBPUbKC7pG6SWpJ0gphao87rwBAASfuSJCifIpmZWXYJKiLWAmOBGcBCkt568yWNl3RcWu27wJmSngV+DYyKiMgqJjMzazoyfRZfREwHptcou7RgfAEwKMsYmho/n8zMLOGHxZrZNuv18b0zX8fulz6f+Tq2V43di8/MzKwoJygzM8slJygzM8slJygzM8sld5LYDvnGsZk1BT6DMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXHKCMjOzXMo0QUkaKmmRpMWSLiwy/xpJc9PhRUkrs4zHzMyajsyeZi6pDJgIHAksA2ZLmhoRC6rrRMR5BfXHAX2zisfMzJqWLM+gBgCLI2JJRKwG7gaG11F/JPDrDOMxM7MmJMsE1QlYWjC9LC3bhKQ9gG7AI7XMHy2pUlJlVVVVgwdqZmb5k5dOEiOAeyNiXbGZETEpIioioqK8vHwrh2ZmZo0hywS1HOhSMN05LStmBL68Z2ZmBbJMULOB7pK6SWpJkoSm1qwkqSewM/BEhrGYmVkTk1mCioi1wFhgBrAQuCci5ksaL+m4gqojgLsjIrKKxczMmp7MupkDRMR0YHqNsktrTF+eZQxmZtY05aWThJmZ2UacoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJecoMzMLJcyTVCShkpaJGmxpAtrqXOKpAWS5ku6K8t4zMys6WieVcOSyoCJwJHAMmC2pKkRsaCgTnfgImBQRLwjqWNW8ZiZWdOS5RnUAGBxRCyJiNXA3cDwGnXOBCZGxDsAEfFWhvGYmVkTkmWC6gQsLZhelpYV2gfYR9IsSU9KGlqsIUmjJVVKqqyqqsooXDMzy5PG7iTRHOgODAZGArdIal+zUkRMioiKiKgoLy/fuhGamVmjyDJBLQe6FEx3TssKLQOmRsSaiHgFeJEkYZmZ2XYuywQ1G+guqZuklsAIYGqNOr8jOXtCUgeSS35LMozJzMyaiMwSVESsBcYCM4CFwD0RMV/SeEnHpdVmACskLQAeBS6IiBVZxWRmZk1HZt3MASJiOjC9RtmlBeMBnJ8OZmZmGzR2JwkzM7OinKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXnKDMzCyXSkpQkg6RdHo6Xi6pW4nLDZW0SNJiSRcWmT9KUpWkuenwrfqFb2Zm26rmm6sg6TKgAugB3Aq0AO4ABm1muTJgInAksAyYLWlqRCyoUXVKRIzdgtjNzGwbVsoZ1FeB44APACLiDWDHEpYbACyOiCURsRq4Gxi+pYGamdn2pZQEtToiAggASW1KbLsTsLRgellaVtOJkp6TdK+kLsUakjRaUqWkyqqqqhJXb2ZmTVkpCeoeSTcD7SWdCTwM3NJA6/890DUi9gf+APyyWKWImBQRFRFRUV5e3kCrNjOzPKvzHpQkAVOAnsB7JPehLo2IP5TQ9nKg8Iyoc1q2QUSsKJj8OfCTEto1M7PtQJ0JKiJC0vSI6E1yhlMfs4HuaY+/5cAI4F8LK0jaLSLeTCePAxbWcx1mZraN2mwvPuBpSf0jYnZ9Go6ItZLGAjOAMmByRMyXNB6ojIipwNmSjgPWAm8Do+oXvpmZbatKSVAHAqdJeo2kJ59ITq7239yCETEdmF6j7NKC8YuAi+oVsZmZbRdKSVBHZR6FmZlZDZvtxRcRrwHtgWPToX1aZmZmlpnNJihJ5wB3Ah3T4Q5J47IOzMzMtm+lXOL7JnBgRHwAIOkq4AnghiwDMzOz7VspX9QVsK5gel1aZmZmlplSzqBuBZ6SdH86fTzwi8wiMjMzo4QEFRFXS5oJHJIWnR4Rz2QalZmZbfdKed3GQcD8iHg6nd5J0oER8VTm0ZmZ2XarlHtQ/wOsKphelZaZmZllppR7UEpftwFARKyXVMpyZmbbvEE31Pnu1gYxa9yszNeRR6WcQS2RdLakFulwDrAk68DMzGz7VkqC+jZwMMkTyZeTPJtvdJZBmZmZldKL7y2SV2WYmZltNbWeQUk6U1L3dFySJkt6N309+xe2XohmZrY9qusS3znAq+n4SOAAYE/gfOC6bMMyM7PtXV0Jam1ErEnHhwG3R8SKiHgYaJN9aGZmtj2rK0Gtl7SbpFbAEODhgnmtsw3LzMy2d3UlqEuBSpLLfFMjYj6ApMMosZu5pKGSFklaLOnCOuqdKCkkVZQeupmZbctq7cUXEdMk7QHsGBHvFMyqBE7dXMOSyoCJwJHAMmC2pKkRsaBGvR1J7nf50UlmZrZBnd+Dioi1NZITEfFBRKyqbZkCA4DFEbEkIlYDdwPDi9S7ArgK+LjEmM3MbDtQyhd1t1QnYGnB9LK0bIO0u3qXiPjfDOMwM7MmKMsEVSdJzYCrge+WUHe0pEpJlVVVVdkHZ2ZmjW6LEpSkniVUWw50KZjunJZV2xH4F2CmpFeBg4CpxTpKRMSkiKiIiIry8vItCdnMzJqYLT2DeqiEOrOB7pK6SWpJ8rikqdUzI+LdiOgQEV0joivwJHBcRFRuYUxmZrYNqbUXn6Tra5sFtN9cwxGxVtJYYAZQBkyOiPmSxgOVETG17hbMzGx7VtfDYk8nuT/0SZF5I0tpPCKmA9NrlF1aS93BpbRpZmbbh7oS1GxgXkQ8XnOGpMszi8jMzIy6E9RJ1PLdpIjolk04ZmZmibo6SbSNiA+3WiRmZmYF6kpQv6sekfTb7EMxMzP7p7oSlArG98w6EDMzs0J1JaioZdzMzCxzdXWSOEDSeyRnUq3TcdLpiIidMo/OzMy2W3W9bqNsawZiZmZWqNEeFmtmZlYXJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8slJygzM8ulTBOUpKGSFklaLOnCIvO/Lel5SXMlPSZpvyzjMTOzpiOzBCWpDJgIHA3sB4wskoDuiojeEdEH+AlwdVbxmJlZ05LlGdQAYHFELImI1cDdwPDCChHxXsFkG/xaDzMzS9X1uo1PqxOwtGB6GXBgzUqS/h9wPtASOLxYQ5JGA6MBdt999wYP1MzM8qfRO0lExMSI2Av4PnBJLXUmRURFRFSUl5dv3QDNzKxRZJmglgNdCqY7p2W1uRs4PsN4zMysCckyQc0GukvqJqklMAKYWlhBUveCyWOAlzKMx8zMmpDM7kFFxFpJY4EZQBkwOSLmSxoPVEbEVGCspCOANcA7wDeyisfMzJqWLDtJEBHTgek1yi4tGD8ny/WbmVnT1eidJMzMzIpxgjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1xygjIzs1zKNEFJGippkaTFki4sMv98SQskPSfpj5L2yDIeMzNrOjJLUJLKgInA0cB+wEhJ+9Wo9gxQERH7A/cCP8kqHjMza1qyPIMaACyOiCURsRq4GxheWCEiHo2ID9PJJ4HOGcZjZmZNSJYJqhOwtGB6WVpWm28C/1dshqTRkiolVVZVVTVgiGZmllfNGzsAAElfAyqAw4rNj4hJwCSAioqK2Iqhmdl2Yn3LHfmw75ms26kLoJKXu7JZ9n3NFi5cmPk6toZWrVrRuXNnWrRoUVL9LBPUcqBLwXTntGwjko4ALgYOi4hPMozHzKxWH/Y9k892O4D2bVoilZ6gWpSVZRhVoueuPTNfR9YighUrVrBs2TK6detW0jJZpv7ZQHdJ3SS1BEYAUwsrSOoL3AwcFxFvZRiLmVmd1u3Upd7JyUoniV122YWPP/645GUyS1ARsRYYC8wAFgL3RMR8SeMlHZdWmwC0BX4jaa6kqbU0Z2aWMTk5Zay++zfTe1ARMR2YXqPs0oLxI7Jcv5mZNV256CRhZpY3B18zp0Hbu/N7AzZb56Zrb2LafdMoKyujWbNmXD7hcg74wgGb1KusrOT222/n+uuv32Te888/z9e//nUAXn/9ddq1a0e7du3o0KEDDz/88KffkK3ICcrMLAeeqXyGmX+YyX1/uI+WO7TknRXvsGbNmqJ1KyoqqKioKDqvd+/ezJ07F4BRo0YxbNgwTjrppI3qrF27lubN8//n38/iMzPLgaq/V7HzZ3em5Q4tAdh5l53p+LmOPP/M8xx88MEccMABDBgwgPfff5+ZM2cybNgwAP70pz/Rp08f+vTpQ9++fXn//feLtj948GDOPfdcKioquO6665gzZw6HHXYY/fr146ijjuLNN98E4OWXX2bo0KH069ePQw89lBdeeGHr7IAi8p9Czcy2A4MGD+J/rv4fhh48lIGHDuTo4UfTp6IP5485n/vuvY/+/fvz3nvv0bp1642W++lPf8rEiRMZNGgQq1atolWrVrWuY/Xq1VRWVrJmzRoOO+wwHnjgAcrLy5kyZQoXX3wxkydPZvTo0dx00010796dp556iu985zs88sgjWW9+UU5QZmY50KZNG+596F7mPDmHp2Y9xXfHfJcx54yhfNdy+vfvD8BOO+20yXKDBg3i/PPP57TTTuOEE06gc+fanxh36qmnArBo0SLmzZvHkUceCcC6devYbbfdWLVqFY8//jgnn3zyhmU++aTxvp7qBGVmlhNlZWUMGDSAAYMGsM+++3DXrXdtdpkLL7yQY445hunTpzNo0CBmzJhBz57Fv9jbpk0bIPnSbK9evXjiiSc2mv/ee+/Rvn37DfewGpsTlGVi0A2DMl/HrHGzMl+H2dbyyuJXUDPRdc+uACycv5A9u+/JY48+xuzZs+nfvz/vv//+Jpf4Xn75ZXr37k3v3r2ZPXs2L7zwQq0JqlqPHj2oqqriiSeeYODAgaxZs4YXX3yRXr160a1bN37zm99w8sknExE899xzHHDApj0JtwYnKDOzIh4/r19J9V5poEcdffjBh1x58ZW8/+77lDUvY/duuzN+wnhOGHEC48aN46OPPqJ169abdBW/9tprefTRR2nWrBm9evXi6KOP3uy6WrZsyb333svZZ5/Nu+++y9q1azn33HPp1asXd955J2eddRZXXnkla9asYcSIEY2WoBTRtJ69WlFREZWVlY0dRmb6XXB75uu4f8cJma9j5M6bXitvaD6Datrydqy/e+S17LPH5+q9joZKUHXZFp7FV23hwoXsu+++G5VJmhMRm/SbdzdzMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJScoMzPLJX8PysysiL//fERJ9T5TYnsfjvnNZuv027Mfc5b88zUf9999P/Oencedk+8s+mTytm3bsmrVKl599VW6devG9ddfz7hx4wAYO3YsFRUVzJ49m1mzZrF69WpeeeUVevToAcAll1zCtGnT+NOf/kS7du2ICK6++mqGDBkCJA+XffPNNzd8MXjvvffm3nvvZdGiRYwZM4aVK1fyySefcOihhzJp0iQ+/PBDzjzzTJ577jkigvbt2/Pggw/Stm3bEvfQppygzMy2AR07duS6665jzJgxtGzZckP5xIkTAXj11VcZNmzYRo8xmjZtGhMmTOCkk07i0UcfZfTo0bz00ksb5t95552bvNbj7LPP5rzzzmP48OFA8v4pgOuuu45dd911w/SiRYto0aLFp9omX+IzM9sGlJeXM2TIEH75y19u0fIDBw5k+fLlm6335ptvbvRA2t69e28o79Sp04byHj16sMMOO2xRLNUyTVCShkpaJGmxpAuLzP+ipKclrZV0UrE2zMy2Fx9//DFfHfLVDcMNE26o1/Lf//73+elPf8q6devqve4HH3yQ448/fqOy0047bcO7pi644AIAzjvvPA4//HCOPvporrnmGlauXAnAGWecwVVXXcXAgQO55JJLNjoT21KZXeKTVAZMBI4ElgGzJU2NiAUF1V4HRgHfyyoOM7OmolWrVtz/x/s3TFffgwKQtEn9mmV77rknBx54IHfdtfmnoFe74IIL+MEPfsCyZcs2ebp5sUt8p59+OkcddRQPPvggDzzwADfffDPPPvssffr0YcmSJTz00EM8/PDD9O/fnyeeeGKTxxrVR5ZnUAOAxRGxJCJWA3cDwwsrRMSrEfEcsD7DOMzMmrxddtmFd955Z8P022+/TYcOHTap94Mf/ICrrrqKUp+zOmHCBF588UWuuuoqzjjjjJKW+fznP88ZZ5zBAw88QPPmzZk3L0mibdu25YQTTuBnP/sZX/va15g+fXpJ7dUmywTVCVhaML0sLas3SaMlVUqqrKqqapDgzMyaksGDBzNlyhRWr14NwG233caXvvSlTer17NmT/fbbj9///vf1an/s2LGsX7+eGTNm1FnvwQcfZM2aNQD87W9/Y8WKFXTq1IlZs2ZtSKCrV69mwYIF7LHHHvWKoaYm0YsvIiYBkyB5mnkjh2Nm24Fdv3V3SfW2xtPMAYYNG8acOXPo168fZWVl7LXXXtx0001F61588cX07du3Xu1L4pJLLuEnP/kJRx11FJDcg6ruZt6hQwcefvhhHnroIc4555wNr5afMGECn/vc53jooYc466yziAjWr1/PMcccw4knnvgptjjD121IGghcHhFHpdMXAUTEj4rUvQ2YFhH3bq5dv27j0/PrNiwP8nas+3UbW0deXrcxG+guqZuklsAIYGqG6zMzs21IZgkqItYCY4EZwELgnoiYL2m8pOMAJPWXtAw4GbhZ0vys4jEzs6Yl03tQETEdmF6j7NKC8dlA55rLmZmZ+UkSZmaWS05QZmaWS05QZmaWS03ie1BmZlvb4b8d3aDt/eKUX2y2TtVbVfzoP37EvLnz2LHdjuzSYRcuuuKiWruZF75uY9iwYcybN4+ZM2cyfPhwunXrxieffMKIESO47LLLuO2226isrOTGG29s0O3KkhOUmVkORATjTh/H8accz9U3Xw3AC/NfYEXVinq3deihhzJt2jQ++OAD+vTpw7HHHtvQ4W4VvsRnZpYDT816iuYtmjPiG/98UWLPXj2pOKiCCRMm0L9/f/bff38uu+yyktts06YN/fr1Y/HixQC88cYbDB06lO7du/Pv//7vAEyePJlzzz13wzK33HIL5513XsNs1KfkBGVmlgMvvfASvfbvtUn5rJmzeOmll/jrX//K3LlzmTNnDn/+859LanPFihU8+eST9OqVtDt37lymTJnC888/z5QpU1i6dCmnnHIKv//97zc8X+/WW28t+aGxWfMlPjOzHJs1cxaPPPTIhmfrrVq1ipdeeokvfvGLtS7zl7/8hb59+9KsWTMuvPBCevXqxezZsxkyZAjt2rUDYL/99uO1116jS5cuHH744UybNo19992XNWvWbHgJYWNzgjIzy4G9e+zNjGmbPkk8IrjooosYM2ZMyW1V34OqqfANt2VlZaxduxaAb33rW/zwhz+kZ8+enH766VsQfTZ8ic/MLAcOOuQg1nyyhnt+dc+GskULFtF2x7ZMnjyZVatWAbB8+XLeeuutBl33gQceyNKlS7nrrrsYOXJkg7b9afgMysysiEdOnFRSvYZ6mrkkbrj1Bn70Hz/i5zf+nJY7tKRTl05cdMVF7N15bwYOHAgkXcvvuOMOOnbs2CDrrXbKKacwd+5cdt555wZt99NwgjIzy4mOn+vINbdcs0n50IFDOeecczYprz6r6tq164a32g4ePJjBgwdvUnfUqFGMGjVqw3TNS4CPPfZYbnrvVfMlPjOz7djKlSvZZ599aN26NUOGDGnscDbiMygzs+1Y+/btefHFFxs7jKJ8BmVmBkCQ1RvGLVHf/esEZWYGlL23lJUfrHaSykhEsGLFClq1alXyMr7EZ2YGfOaZW3ibM6naqQugkpf7R7Ps/8+Pt7eNpNmqVSs6dy79HbVOUGZmQLPV79P2qavrvdy3d94pg2g2NmvcrMzXkUeZpn5JQyUtkrRY0oVF5u8gaUo6/ylJXbOMx8zMmo7MEpSkMmAicDSwHzBS0n41qn0TeCci9gauAa7KKh4zM2tasjyDGgAsjoglEbEauBsYXqPOcOCX6fi9wBBJpV/8NTOzbVaW96A6AUsLppcBB9ZWJyLWSnoX2AX4R2ElSaOB6tdbrpK0KJOItxN71H+RDtT4TPJAZ/t/Gaubj/Umo+hH1SQ6SUTEJKC0B2NZg5NUGREVjR2HWdZ8rOdLlpf4lgNdCqY7p2VF60hqDrQD6v9+YzMz2+ZkmaBmA90ldZPUEhgBTK1RZyrwjXT8JOCR8LfkzMyMDC/xpfeUxgIzgDJgckTMlzQeqIyIqcAvgF9JWgy8TZLELH98edW2Fz7Wc0Q+YTEzszzys/jMzCyXnKDMzCyXnKAyIul4SSGp51Za3w8aoI1XJT2fDgskXSmp9EcPZ0zSYEnTNl+z1uU/9T6yLSNpnaS5kuZLelbSdyU1S+dVSLo+Hd9B0sNp3VMlHZouM1dS6wzjGyXp87XM6yppXpHy8ZKOyComc4LK0kjgsfTn1lCvP75KFPv8vxQRvUmeBLIncHNDBJcTTlCN56OI6BMRvYAjSR6BdhlARFRGxNlpvb5pWZ+ImAKcBvwonf5ocyup47jenFFA0QRVm4i4NCIe3oJ1WakiwkMDD0Bbku947QMsSssGAzNJHun0AnAn/+yk0h94HHgW+CuwI0nPxwkk3fWfA8akdXcD/gzMBeYBhwI/BtalZXem9c5P588Dzk3LugKLgNuB+cAeNeJ+FehQML0T8C7w2XT6goJ4/jMtawP8bxr7PODULdimuvbN0LTsaeB6YFrBeienbT8DDE/LRwH3AQ8CLwE/Scs32UceturvxKoa03uSfOdR6ec/DegILE6PubnAGJLeva8UHNfFjsFNjus66i0EbknrPQS0JvmKy6q0jblA6xqxdgXmFdmm24CTCn53/jM9Tp8Hem7mOO0K/CWt/zRwcFreDPhZesz/AZhesI5+wJ+AOSS9o3dr7M818+OmsQPYFgeS//p+kY4/nh5Yg9NfvM7pQfgEcAjQElgC9E/r70TS/X80cElatgNQCXQDvgtcnJaXATum46sK1t8v/SVpQ5Is55P8Z9oVWA8cVEvcr1KQoNKyuSSPqPoySRdcpfFPA74InAjcUlC/3RZsU237phXJo7C6p+u9h38mqB8CX0vH2wMvpts7Kl13u3T514AuNfeRh63+O7HJvgdWArumn3/157phPJ2+reAPdG3H4EbH9WbqrQX6pPXuKTiGZgIVtcTeldIS1Lh0/DvAzzdznH4GaJWWdyf56g0kyXJ6GvfngHfSshYkf0vK03qnknx1p9E/2yyHJvGooyZoJHBdOn53Oj0N+GtELAOQNJfkwH8XeDMiZgNExHvp/C8D+0s6KW2nHcmBPBuYLKkF8LuImFtk/YcA90fEB2lb95GcaU0FXouIJ+uxLdUPAftyOjyTTrdN4/kL8N+SriL5w/IXSb3ruU2ra9k3q4BXIuKltPwO/vlMxi8Dx0n6XjrdCtg9Hf9jRLybLrOA5D/qwudCWtNU2zH4Ohsf13XVe6Xgd2YOyXHWUO4raPeEgliKHadvADdK6kNyZr9POv8Q4DcRsR74m6RH0/IewL8Af0ifp10GvNmAseeSE1QDk/RZ4HCgt6QgOZCC5DLYJwVV11H3/hfJf2Qziqzji8AxwG2Sro6I2+sRYnXSKiP5RQKYGhGXFlnPjiS/wC+m8fwoIja5JyXpC8BXgCsl/RG4vz7bJGkw9ds31W2dGBEbPThY0oFb0JZtZZL2JPls3gL2LXUxihyD6XvkPiixXs1jY5OOF+kxVL3spSSXCUtR3XbhMVfbcXo58HfgAJKzpY8307aA+RExsMRYtgnuJNHwTgJ+FRF7RETXiOhCcg390FrqLwJ2k9QfkqSQPpdwBnBWeqaEpH0ktZG0B/D3iLgF+DnwhbSdNdV1Sc5qjpf0GUltgK+mZRtExLpIbjz3qSU5tSW5Fv67iHgnjeeMtBxJnSR1THs+fRgRd5DcX/pCfbepjn35AtBV0l7pdGGHkxnAuOrXs0jqW0c71Qr3kTUSSeXATcCNkV6vKlHRY/BT1Cv0Psl9UiLiqYLfjZqPZ6uv2o7TdiRXGdYDXyf5RxZgFnCipGaSqi9/QvI7VS5pYNpOC0m9PmVsuef/LBveSDZ98eJvgbOAl2tWjojVkk4Fbki70X4EHEGSfLoCT6cHdxVwPMkBe4GkNSSXwP4tbWoS8JykpyPiNEm3kdyYheR6+DMq7Y3Fj6bra0ZyJnRFGudDkvYFnkh/11YBXwP2BiZIWg+sAc7agm0qKiI+VvKqlf+V9CFJkt0xnX0FcG26zc1I/gkYtplt22gflbAvrOG0Ti/dtiC5D/QroF7vV6/jGFy3JfVquA24SdJHwMDYtMdgD0nLCqbPKzHs2o7TnwG/lfRvJB16qs8AfwsMARaQXJZ+Gng3/Z06CbheUjuSv93Xktxf3mb5UUdmZjkiqW1ErJK0C8k/mYMi4m+NHVdj8BmUmVm+TJPUnqQ37BXba3ICn0GZmVlOuZOEmZnlkhOUmZnlkhOUmZnlkhOUmZnlkhOUmZnl0v8HBh22XXlIYj0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# f-score \n",
    "fig, ax = plt.subplots(figsize=(6,4), nrows=1, ncols=1)\n",
    "\n",
    "g1 = sns.barplot(ax=ax, data=df_result[df_result['Measure'].isin(['Ancestor-Descendent F1-Score', 'Different-Lineage F1-Score'])], x='Measure', y='Value', hue='Method', width=0.4)\n",
    "ax.set_title('Performance comparison')\n",
    "ax.set_xlabel('')\n",
    "ax.set_xticks([0, 1])\n",
    "ax.set_xticklabels(['Ancestor-Descendent', 'Different-Lineage'])\n",
    "ax.set_ylabel('F1 Score')\n",
    "ax.legend_.set_title('')\n",
    "ax.legend(loc='lower right')\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Method</th>\n",
       "      <th>Measure</th>\n",
       "      <th>Value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>User Time</td>\n",
       "      <td>549.92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ScisTree</td>\n",
       "      <td>Elasped Time</td>\n",
       "      <td>100.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>User Time</td>\n",
       "      <td>5960.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>HUNTRESS</td>\n",
       "      <td>Elasped Time</td>\n",
       "      <td>277.64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>User Time</td>\n",
       "      <td>208125.69</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>CellPhy</td>\n",
       "      <td>Elasped Time</td>\n",
       "      <td>6962.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Method       Measure      Value\n",
       "0  ScisTree     User Time     549.92\n",
       "1  ScisTree  Elasped Time     100.50\n",
       "2  HUNTRESS     User Time    5960.54\n",
       "3  HUNTRESS  Elasped Time     277.64\n",
       "4   CellPhy     User Time  208125.69\n",
       "5   CellPhy  Elasped Time    6962.00"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# running time\n",
    "def get_running_time(file):\n",
    "    with open(file) as f:\n",
    "        line = f.readline().strip()\n",
    "        user_time_str, _, real_time_str, *args = line.split()\n",
    "        user_time = float(user_time_str.split('user')[0])\n",
    "        real_time = .0\n",
    "        real_time_str = real_time_str.split('elapsed')[0].split(':')[::-1]\n",
    "        for i in range(len(real_time_str)):\n",
    "            real_time += 60**i * float(real_time_str[i])\n",
    "    return user_time, real_time\n",
    "\n",
    "scsitree_user_time, scistree_real_time = get_running_time(f'{prefix}.scistree.time')  \n",
    "huntress_user_time, huntress_real_time = get_running_time(f'{prefix}.huntress.time') \n",
    "cellphy_user_time, cellphy_real_time = get_running_time(f'{prefix}.cellphy.time')\n",
    "\n",
    "methods = ['ScisTree'] * 2 + ['HUNTRESS'] * 2 + ['CellPhy'] * 2\n",
    "measures = ['User Time', 'Elasped Time'] * 3\n",
    "vals = [scsitree_user_time, scistree_real_time, huntress_user_time, huntress_real_time, cellphy_user_time, cellphy_real_time ]\n",
    "df_time = pd.DataFrame({'Method': methods, 'Measure': measures, 'Value': vals})\n",
    "df_time\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_735248/3541172569.py:17: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.\n",
      "  axes[0].set_yticklabels(axes[0].get_yticklabels(), fontsize=12)\n",
      "/tmp/ipykernel_735248/3541172569.py:25: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.\n",
      "  axes[1].set_xticklabels(['User Time', 'Elapsed Time'], fontsize=16)\n",
      "/tmp/ipykernel_735248/3541172569.py:28: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.\n",
      "  axes[1].set_yticklabels(axes[1].get_yticklabels(), fontsize=12)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAE/CAYAAAC+ZjJjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAABQNElEQVR4nO3dd7wcVf3/8dc7CQTITSMEiEAIVSkS1KAUA+iXokBCEETpsRARUBB+UqQY+SJfQUFpAqGFGqWGBAhNBRIQpAhCUGLQ0CEJIZUSQj6/P85s2Gx2793bdu/ufT8fj33cu2fOzJydmZ39zJlzzigiMDMzM6t1XapdADMzM7O24KDGzMzM6oKDGjMzM6sLDmrMzMysLjioMTMzs7rgoMbMzMzqgoOaOiZpkKSQNLraZSmXpC6SRkv6j6QlkjrsmAOSjpT0L0kfZtt5ULXLZG1D0oOSZrTDckse35L2k/SspPez42lnSSNz/7dgXS2e16xW1U1QI2k1ScdKmixpjqSPJL0t6e7sy92t2mW0shwG/Bz4C/A94JDqFqc4SV8BLgb+BRxBKuesqhbKakHR41vSpsA4YB5wdJb+zyqV0fLkBYf7lZieu3i8qMT0oZJukPTfLGB9T9J0STdKGiFJBfl7SzpV0jOS5kpamM07XtL3S6xjK0nXSHo5u8iaI+kvkr4jqWsjn20dSedI+oekBdm8MyRdL+l/mrOdOoq6+KGXtDFwF7Ap8ADwf8BsYE1gF+BqYHPghGqVsUpeBlYFllS7IM2wK+nE/v3o2CND7pr9/W5EzKlqSayWlDq+dyadj4+NiKdziZKuA/4ALG7Bulozr7WSpC7ARcAPgdeAm4BpwFJgA+BrwO3Az0i/WUjqBTwBbAjcAlxF2n8bAl8GjgGuKFjPD4ELgXeBsaQLrb7Avtn8B0jaJyIWFcy3JymQ7g7cDIwB3gcGASOAByTtGRF3t8kGqZCaD2okrQrcSdrp+0bEbQVZzpa0DbBNxQtXJZJ6RsSC7KT5QbXL00xrA3M7YkCTXfF0j4j3SOWkrQOa7Hj+KCJqKRC18pU6vtfO/i53PEXEx8DHLVlRa+a1NnE6KaC5AfheRHxYMP1nWY3vp/LSDgc2IQW35xcuUNLaBe93JdUYPwf8T0TMzpv8G0lnAqcAlwCH5s23BSmQmQPsGhHL1QpKOh04iNr7/YCIqOkX8CMggF81c74RwCPAImBh9v/eRfLNAB4EBpNqgRYCM4FzSUHhKsBvgNdJB8DDwGYFyxiZlXEXYDSpBuVD4B/At4usczfgj8B/SJHzXOA+YKcieR/MypiL7Oek3RqQIu4ARhfM0w04EXghK/M7pCuGzxbkWzY/sBfpCuID4E3g10C3Zmzv7wNPZ59nXvZ5vpw3fedsXYWvsU0sN0hXJ7sAjwHvAW8B5wMNRfL3Bs4Gpmf7YBbpamXDRvbZacBLwEd56YWvB/Pm3Srbnu9k2+sFUi1h14J1jM3m7U+6onqbdBU3KNvmQaph/F22zd8D/gR8Opv/G3nbdAYwqsjn/RYwAXgl+7yzgfHAVo0c658h1XwuyPbVLcDaRfL3An5Juk2SO46mUHBMAwNIJ9VXSFedb5CuCtcs89j5YratpmXbYAHp+7pPkby5bdo7W+fMrGyPAF8qkr8vcHm2XRZln/8L2d8ZzTi+v5V99gVZGR8H9ivn+C6RPqPgONy5YH0rZ8fUM9n65gFPAkcXOYYL5+1Oqh2Ymm2bucBE4HMF+XJlHgl8J8v/Ien8dUKJ7fA50o/l21neV0nfr42yMs8CHikx70+z9e1YxvZu7nesrOOhxLpy23G/EtMHZdMvyktbk/S9/A/pQqjc4+jSbFlblpn/KdI5Y4sS00U6LwZ553fg1izta+WWrVZeNV9TA+Tuc44pdwZJR/JJe4gzsuSRwHhJP4iIwmWtC9xPCjRuIQUdx5Fu62xBusXzK2AN4P9ly9ksIpYWLOdsoAfw++z9d4BxklaJiLF5+UYCqwPXkqot1yEFBX+S9JWImFyw3AbgIdIX9RTSF6oxNwD7Z5/pEtJV4lHAXyUNjYi/F+TfAziS9IW7Ctg7+5zvAmc1sS4knU064fyNdDLtCYwC/iJp70jVm/8ktSM4hbQdf5LN/lJTywc+TzoOLidts68APwa2lLRrbj9I6g08CgzMPsdU0g/ukcDjkoZExMsFy/4NsFK27PnAv7NyjgKG8kmbn7ezdQwh7YuPSMfYW8Aw0r4fTLr6KXR/lu9/ScfHwrxp12TvzyIFP8cD90o6DTiHtP+uIrXPuEzSCxExJW/+o0kn/jHZOjbKyv6IpM9HxL8LyrIO6Qf9dtKPzGDgB6QAZrdcJkl9SD/iW5C+E5cAXUk/anuRbnsgaSDwV9IP2pWk/bkx6Qr2K9k2n1dkm+TbhxRo3UT6Qe1Haptym6SDIuLGIvPcS/oBPSPLfxxwl6QNImJBVraVsnzbkG7VPAZsTbp4eaeJMi2TdzV8DykAXpqV+WZJR0fExTR+fD9AClD3ydJns/wxULi+lbNy70y6OLie9EP92Ww5Rdt25H3me4Dts898EekH/3DSMbFjRDxZMNsRwFqk/TcXOJhUA/5a/raXtBfpx3IR6RbJdNK5ZXfSj/RLkq4Bjpf06Yh4sWA93wWmRcTDpcqfracl37Emj4cy9JS0RpH0vkXS9iRd8F4XK9bQNCZ3vvuOpBOjkRpbSRuQzn2PRMTUYnkiIiRdCXyJdGw8J2mVrHyvRsQ9zShbbah2VNXaF+nkM68Z+fuSThjTgV556b1IB9QCoE9e+gxSRPvNguXkIuQ7AOWl/zjLv3te2sgs7WWgd1567yxtDrBqXnqPIuVei3Syu7sg/cFs2WcWmWcQBTU1pHv6QQrQ8ss9mBSkTS4y/yJgUF66gOeBN8vY3p/OttMUYOW89E+RTpAzyLu6ovlXyLkr2xEF6edn6d8uSHsfGFyQd31SwDK2yD57EVityHrHktWIFaQ/km3HrfLSRPpBDlIV8XLLAK4vspzR2bSJJY6v+cB6een9ST9s4wqWU+xY2ox0Ff37gvQZ2bL3L0i/OEv/dF7a77O0YrVDXfL+v4N0dbxuQZ4h2XYaXTh/keUV+wyrZfvmhWL7pchn+2aW/oO8tFFZ2i8K8h5LXm1JE2X7fJb3rCLTxmf7qWdTx3fe/h5UkJ47DnfOSzuhkXV2aWLen1BwfsrSe5Fq0h7MS9s5y/sGy5+3ViMFCH8tkjYTWKdUuUjtHgM4p2D6Dll60RqgNviONXk8NLK+3HZs6pVfU3NullasNrE3KbDNvfrkTeub7YcgXSjdQqpV/3L+vs3yDsvyXVDmMXpL9v6z2fsJTX32WnzVQ++nXqRApFy7kq6GL4iI+bnE7P8LSLUeuxTM83pE3FyQNoX0RbowsiMlk6tF2aTIui+JvKvS7P9LSQfyznnpyxp0SWqQ1I90b/xxUsRdzG9KpBfaJ/v7y/xyR8SzpB/QL0vqXzDP+IiYkZc3SL031pbU0MT69iZtp3MiYlmDxYh4g9SAe33S1X1rvBgR4wvSfpX93Qcg62FwEOn24OuS1si9SEHbY+TVROS5JFIbmiZJWpN0BTwhIv6RS8+21y/zy1OgsX13QYnja0JEvJq3jlmkH/nljrvcsaSkV/Z5c3mLHUtvRMRNBWl/zv5uki2rC/Bt4J+xYq0msXzN2F6k218fFGzzGaQLi2LbvHB5+d+H1bLvw2pZuTbLGlcW+m1jnyEzgvS9Orcg7yWkYKQcB5F+IK7J/3zZZ5xAqpXcrsxllesgUi3pGYUTYsXa4UIHk2qonyoo68qkGsMvZ+268l1dcN56j/R9yd+Wu5N+oM+NiNdLlSsippFqWQ7V8j1Sv0cKVK5prPCt+I6Vczw05QzS70fh6+AieXPHZLHj6E+k72DutaxmNSLeJd3+PJt0S3Ff0rlsMvCSpPzvS24dTdV05srQu4yy1bx6uP00n3TiKNcG2d9i1XW5tA0L0v9bJO+7Jabl0vsVmadYF80XCtcpaSPSF3R3oE9B/mBFsyJibpH0YjYg1ZwUK8tU0ol+A5bvnvyfInlz1fP9aKSqnPK3d2GVd3Os8Fki4k1Jc/lku/YnlXU3Sne9LvaDMK0Z5Wjss/4zW37hsdXUOgq3fanjLjdt/fwESZ8j3dbamRTM5yu2jKb2NaQfr76k2xiN+TRp2IjvZa9iiq1vOdkP2ZmkALnYrdU+rHiCXm65EfFO1nM2/3u5Iam2cX5B3g8l/YfitxUKbUYK2v/VSJ61ylhOc2wCPBMRLWnEuRnpdnljww+sQWoLk1PqmMjflrng4O9llGEM6Rb4XqRb9T1Jt8PvjIi3m5i3pd+xco6HpjwXEQ8UJqr4+FS5Y6pYwH1kXvr1hROzC5STgJOyAH470vY5GLhd0uCImM6KwUophcFPbr7m/G7WjHoIap4HdpS0YUQ0eYJsocZ6EJSaphLpjcpqPh4m/QD9jtSqfQHpy3oy8NUis5VVk9AKjX3+Fn3OKsiV8wHSVVC52nvb5q58Sym17Zs87rL2LA+TTmL/S6qdWUQKjH9HqpUsd7nLLbtMufzXU/oK/P1GF5B+ee4j/RifTwp+55HK+R3gQIqMtxWp509jZWorIm3Pr1N62xVt71AlIp1TjmskT2HA09Y9qG4l1Yp/j3SL7luk890VjczTKhU8HnKez/5uTWqfll+Wvy1budRoYBoR75B6994p6VVSm8Rvk4L83Do+30RZctOfy/7+m3T7eesm5qtJ9RDU3ArsSGpI+7My8ucCny1I1YD5Ni/I09Y2I7UxaGyd/0Nqb/LdiLg6P2PWILG1/kP6EdiM1PuqWFmKXcG3Zn2Qtndho9+22t6bFSZIGkC6gs8texapDU+vYldbbSS33bYoMu0zpO3eXsdWMfuQApfhEfGX/AnZFWBzGjDmm02qFRrcRL7ppB/8lVuxzbfK1nNGRPw8f4JKDETWDP8BdpPUK7+2RlJ30tX+uyXn/MS/SeONvBIF3WLb0TTgM5K6R/MaoUIqb3/gz2XcqmpumSD9UN7XWMasJuxa4MeSPkUKbl6n6Zo/6HjfsVLuIrVxO0TSWS3YT8U8lv1dByAi/ivpGWD7rGPKCsdfdlGQqyW9PZvvA0l3A/tI2i0iGt1ftaYe2tRcQboC/X+S9i6WQdIXsh5PkO4bLwJ+lFV75vL0JHUPX5jlaQ8/zNoZ5NbZm9SzYC7pPjN8clVUOMrkbpRuT9Mc47O/J2cHfG75WwLDgSlZ9WdbmUD6Yftp1vMit74BpCvtlymvyroxn5Y0oiDtxOzveFh2T/8G4IsqPTJoU73GGhURM0m9q4Zl2zO3XJFq2aDgqq2dlTqWDueTcVGaLduW44DNJa1wWyl3XGVXmXcD35C0bbF8RdpvFSr1GbakeNuJ5riD1GPr+IL0H1L8tkEx12V/z1KRkVsltfWtJ0jHcV/g1CLra6rm4VrSvi9aU9OK8t5HCnaPz77bTZXrctK2PxvYltRIv8kaoQ74HSsqK+c5pNtlV2WBcjGFx/V2Sj0LixmR/X0hLy33mW/ILlQKjSZt3+vy2yCRxtB5H7hC0qeLFkw6UFKxOwMdWs3X1ETEe0pdCe8i3Z+9jxSUvEO6IvkKqW3KOVn+uZJOIPXoeFzS2GxRI0ldTX8QTXcxbanZ2TpzNTDfIXUv/n7eLYgppC6K52b3al8jXf0cQqo+/GxrChAR90u6iVSF2VfSnXzSpfsDUu+aNhMRL0r6NanHxsOS/sgnXbobgIPKOZk14TngekmXk65Ev0Lq4v0QqZdXzimkXhY3ZdvgMdK4KeuTuq0/RToOWuOYbL2TJeW6m+5FOgZvjIjC2sH2NIl0++w6pSHc3yV9/j1ItWat+f6fSroVekUWcOcazn8uW26uq/sPs2kPZ1fnfyddTG1IaiNzLenEW8o/SbdvTpCU6/G0Kamb+XOkRpUtdTXpODxdqXvsX7Pyf5Myt09EPKH0bLXRwDOSbib1FhqQlW0PUiPctnQ+qefLqUoDi95H+u5uQWrHVNjRoXDeXYFfZz9YfybdnhxIqiX+gPT9aZbsPPw9Um+d5yXlunT3Jx3755FXSx0R/5Q0hdROJEjDEpSrI33HGvMLUhuwI0hNJG4iHb+QhgkZTtrud+bNcxCpO/ddpCEwcm2X9iDtlxfI21YRcY+kH5P26z+z35YXSUHvN0iNqu8nfQ/Jm+95Sd8kXZw8m5XtcVKgsz7puzmYdFu1tlSz61Vbvki9IX5COoG+SxrD4G1SsHMIKw7KtA8p4l+UvR6loFtwlm8Ged0c89JHU7wL5iBW7EY9MkvbhXSg5wZCew44sMiytyJVxb5Lak/zIGlMlLEUdCOmkS7QxcqSpecG3/tnVo45pBqNkoPvlfv5G9k/h5N+0D4gnUTvB4YWyVfy85RYbvDJ4Hu5L+XbpGHDexbJvxppLJHnsrwLsu1wOXmDcVFi4LK86Svsi7xpg7PtOSfbvv+kkYHBSiyj7OOrsW1HujWbGxRuLun7sGWJvDMofqzvnK1zZEF6H9LFwnRScPgOqZdGYZfwNUiDNU7jk8HeniOdiDcvYx+vTxrQbRYpSPsb6fu7wjZqYpsGBYM5ksaDujIre27wvSEtOA73JI2Fktvnr5KCyiPKOb4b2d9Fj0PSGCinsPwAek8AR5YxbzfSxcsTfHL++zepBmi3pvZ7Y9uZNFDieNIF3Iekc90NFAxumeU9JFv+n8rdzm38HVvheCiRL7cdyx58r2D6Ttk2mJHtq/dJQfM4st6heXm3JLWXeYQUHC8mfXf/nh0jvUqs43OkWsNXs3nezY617xRuk4L51iF9N58j3aX4kHSL77rC46ZWXso+mLUjSSNJV4VfiYgHq1ua+qL0lONrImJktctiZuWTtD+pJvXAiBhX7fJYfaiHNjVmZlZ7jiLV6BQ+r8+sxWq+TY2ZmdWGrDH+/5Bup+8InBxt0zPIDHBQY2ZmlbM5cCOpDdClrDias1mruE2NmZmZ1QW3qTEzM7O64KDGzMzM6oKDGjMzM6sLDmrMzMysLjioMTMzs7rgoMbMzMzqgoMaMzMzqwsOaszMzKwuOKgxMzOzuuCgxszMzOqCgxozMzOrCw5qzMzMrC44qDEzM7O64KDGzMzM6oKDGjMzM6sLDmrqiKSRkiLvtUDSs5KOltStjdbRQ9J1kmZm6/hdWyzXzDqeIueUxZJeknSWpFWqXT5YroyDKrjOBwu2S6nXIEkzJI2tVNk6uzb5obMO55vAa0Cv7P8LgTWB09tg2UcBBwDfBaYBb7bBMs2sY8udU3oC+wAnZ///qJqFytwFbEdlz0VHks6vOacB2wDDC/K9Sdpe8ytUrk5PEVHtMlgbkTQSuBrYJCKm56X/Bfh8RPRuxbK7R8SHkq4G/iciBra6wHnLbYtlmVnbauSccj+wPdAzIpZWqXgdRlYTs0tErFvtsnR2vv3UOTwB9JK0pqTBkiZIelfS+5IekTQ0P7OksZJek7SdpEclvQ+cIymAkcB6edWrO2fzfFrS7ZLmZst9TNLXCpY7OptnS0n3SloI3JRNC0lnSjpe0suS3pN0V1bmNSXdJGmepFclnViw3P6SLpM0LZvvVUk3SlqnxPo3yZa9MFvX6ZK6FFnm77NlfZj9vU5S97w8TW5Lszr1NLAasAaApJ3zzwc5xW4NZbdjrpf0bUn/lLRI0pOSvlwwb+489DlJk7Pv9r8lHdFW68jyHpvl/0DS3yRt35a3jAqXlVfe7bPz2gJJb0s6OZv+NUl/z8r8hKQvFFnmN7Jz7HvZOfdmSW1yoVnrHNR0DhsAHwMbA48CqwOHA/sC7wAPFPni9Ab+AIwDvg7cSKrivRd4K/t/O+BpSZ8CpgCDgaOB/YG5wF2Svl6kPHcAD5Gqan+bl34I8FVS1e7RwFDgWuB24B9Zee8GfiVpj7z5Vgc+IFWJfw34KbAJ8IiK3/e/HfgzMAIYD/wCOCw3UVLfbDt9CzgP2AM4AVgJWDnL83nK35Zm9WYQMI90zLfEUOB40m2bbwFdgTsl9SnI14t07rke2Jt0gXaJpK+0xTokfZ90DnogW/7YbH2F5WgP1wDPkW5PjQfOknQ28Gvg7KzMPYDxklbOK/MRwK3AC8B+wA+ALYGHJPWsQLk7tojwq05epFqUAD5Nai/Vl3TAf0z60vwJ+Cewct48XbO08XlpY7Pl7F1kHdcDMwrSfgMsATYuWO6LwNN5aaOz5R5TZLlBaqPTLS/tvCz91Ly0bsBM4OpGtkNXYL1s3n2KrP87BfmfA+7Le39Gts0+18g6ytqWfvlVy68S55TvZt/3o/Py7Zzl27nE/IPy0mYA7wJ989KGZPkOzEvLnYe+kpfWnRRIjWntOkgX9a8CdxeU+RtZvrHN2E5jgddKTJuRv6y88p6el5Y7r30EbJCXPjzLu1P2voEUTF5VsI4NgMXAsdU+Zqr9ck1NffoX6csxB/g9cAOpge9OwM3AUkndlHpEiXSVsmPBMj4C7ixzfTsCj0XePfeI+JhUy7O1pF4F+W8vsZz7I2JJweeAVDuUW+4SYDopaFlG0g+VenotJJ1wX8kmfbrIeu4qeP88kF91uxvwRET8vVghJa1K87alWa3LP6dcCVwWERe1Ynl/jYh3894/l/0tvIXyXkT8JfcmUvu7aUXytWQd62avmwvmu4N0Dmlvk3L/5J3XpkXEf/Py5M6BufPddqTaqxty553s3PNqlrfTn3vc+6k+7UPqqbAAeDkiPsjal3QlVcWeVmwmSV3ik0Z/s7LApByrA8UCgLdIP/R9Wb71f6leCu8WvF/cSPqy20qSfgRcQKrZ+WmWvwvwWH6+PHMK3n9YkK8f8GyJMkL6vM3Zlma1LndO6Q8cBxwp6fGIuLaFy1vuOxipEwKs+H0t/O7Dit/Xlq5jQPZ3ZkG+jyXNLmP5rVXsvFbqHJgr85rZ3wfKXGan46CmPj2fX2uSmQssBS4mtVNZQcGPcHO6xc0B1i6Svna2nMIvWlt3ufs28KeIOD6XIGmDVixvNrBOI9Pn0rxtaVbrlp1TJP2Z1Mbt15JujYhFpDZtkLU5y9OvgmVsrtzF1Zr5iZK6kjWA7oBybZhGAlOLTF9QuaJ0TA5qOomIWCRpMqkx79Nt/KP7EHCspEERMQOWnRi+Bfw9Itp7jIbVWHEciO+0Ynn3AadKGhwRK9TYtPO2NOvQshqPn5Ju0xxJatj6cjZ5S9L3J2fPChevOV7LXt8kdVvPGUHH/W18lBS4bBwR11S7MB1RR91x1j6OAx4G7pV0JelKZQ3g80DXiDiphcv9LenK4X5JPycFGEcCm1KZk9o9wImSfgb8jdSDar9WLO+3wIGknkxnku7Fr0HqHXFERCyg/balWYcXERMkPQEcL+miiHhT0kPAydmtm5nAwcCGVS1oIyJiqaRfAJdLuoLUtmZD4CRSY9wOd7ESEfOzgPJiSf1J7XLmkWqWdwIejIgbq1nGanND4U4kIp4mjXr5DqkNyn3A+cBnST/QLV3uG8CXSdWhlwC3kNqd7BkR97Sy2OU4A7gM+AmpEfJWwO4tXVhEzAV2yJZ1EiloOpfUeHBxlqddtqVZDTkVWAvIjRtzMKkd2wWk3kCvAGdWpWRliogrSOeNXUk1T98jfY4gBQsdTkRcRuoV9WngOtIwF6NJlRTPVK1gHYRHFDYzM8tIGkIaD+fQiLiu2uWx5nFQY2ZmnVLWoeAoYDLptvlmwM9INbJbRsR7VSyetYDb1JiZWWf1Pqlx86GkoSfeJXWXPskBTW1yTY2ZmZnVBTcUNjMzs7rgoMbMzMzqQt21qVljjTVi0KBB1S6GWaf01FNPzY6I/tUuR3uSNAwY1rNnz8M33XTTahfHrFMqda6puzY1Q4YMiSeffLLaxTDrlCQ9FRFDql2OSvC5xqx6Sp1rfPvJzMzM6oKDGjMzM6sLDmrMzMysLjioMTMzs7rgoMbMzMzqgoMaMzMzqwsOaszMzKwu1N3ge2a16qLjJ7Z43qPPHdaGJbHOYIcLd2jxvI/86JE2LIlZ23FNjZlZM0gaJmnMvHnzql0UMyvgoMbMrBkiYmJEjOrdu3e1i2JmBSoe1EhaXdLtkhZJelnSgSXydZd0qaS3Jc2RNFHSOpUur5mZmdWGatTUXAwsBtYCDgIukbRFkXzHANsBWwGfAt4FLqxUIc3MzKy2VDSokdQD2Bc4LSIWRsQUYAJwSJHsGwD3RsTbEfEB8EegWPBjZmZmVvGamk2BJRExLS/tWYoHK1cCO0j6lKTVSLU6kypQRjMzM6tBle7S3QDML0ibB/QskvffwKvA68DHwHPA0cUWKmkUMApg4MCBbVVWMzMzqyGVrqlZCPQqSOsFLCiS92KgO9AP6AHcRomamogYExFDImJI//7927C4ZmZmVisqHdRMA7pJ2iQvbTAwtUjerYGxETEnIj4kNRL+oqQ12r+YZmZmVmsqevspIhZJug04Q9L3SYHL3sD2RbI/ARwq6UHgPeBI4I2ImF2h4taklo5K6xFpzcys1lXjMQlHAlcBM4F3gB9GxFRJQ4FJEdGQ5ft/wAWktjUrA88D+1ShvNYJvXLGZ1s038DTn2vjklglSBpEupDK1Rp/MyJmVa9EZtYSFQ9qImIOMKJI+mRSQ+Lc+3dIPZ7MzCrhoYjYr9qFMLOW82MSzMySHSRNlnSWJFW7MGbWfA5qzKxuSDpa0pOSPpQ0tmBaY49oeRPYGNgRWBP4RuVKbWZtpRptajqFlrbJALfLMGuFN4Azgd2BVQum5T+iZWvgLknPRsTUrIflhwBZZ4ZtgVsrVWgzaxuuqTGzuhERt0XEeFInhGWaekSLpPwBQIcC0ytTYjNrSw5qzKwzaOoRLV+W9JSkycA6wI3FFiJpVHZ768lZs9w5yqyj8e0nM+sMGn1ES0RMooxny0XEGGAMwJAhQ6KNy2hmreSaGjPrDJrziJZGSRomacy8efPapGBm1nYc1JhZZ9CcR7Q0KiImRsSo3r17t1nhzKxtOKgxs7ohqZukVYCuQFdJq0jqFhGLSA/FPUNSD0k7kB7Rcl01y2tmbctBjZnVk1OB94GTgIOz/0/Nph1J6uY9ExhH9oiW5q7At5/MOi4HNWZWNyJidESo4DU6mzYnIkZERI+IGBgRRXs4lbEO334y66Ac1JiZmVldcFBjZtYMvv1k1nFVPKhp4vkr+fkmSVqY91osyc8PMLOq8u0ns46rGoPvlXz+Sn6miPh6/ntJDwJ/rlAZzczMrMZUtKamqeevNDLfINLzWK5t90KamZlZTar07aemnr9SyqHA5IiY0V4FMzMrh9vUmHVclQ5qGn3+SiMOBcaWmuiHzJlZpbhNjVnHVemgptnPX5H0ZWBt4JZSeSJiTEQMiYgh/fv3b5OCmpmZWW2pdFDTkuevHAbcFhEL27VkZmZmVtMqGtQ09/krklYF9qeRW09mZpXkNjVmHVc1Bt8r+vwVSUMlFdbGjADmAn+paAnNzEpwmxqzjqvi49RExBxSsFKYPpnUkDg/bRwp8DEzMzNrlB+TYGZmZnXBQY2ZmZnVBQc1ZmbN4IbCZh2Xgxozs2ZwQ2GzjstBjZmZmdUFBzVmZmZWFxzUmJmZWV1wUGNmZmZ1wUGNmZmZ1QUHNWZmzeAu3WYdl4MaM7NmcJdus47LQY2ZmZnVBQc1ZmZmVhcqHtRIWl3S7ZIWSXpZ0oGN5P28pIclLZT0tqRjKllWMzMzqx3dqrDOi4HFwFrA1sBdkp6NiKn5mSStAdwD/AS4BVgZWLeyRTUzM7NaUdGaGkk9gH2B0yJiYURMASYAhxTJfhxwb0TcEBEfRsSCiPhnJctrZmZmtaPSt582BZZExLS8tGeBLYrk3RaYI+lRSTMlTZQ0sNhCJY2S9KSkJ2fNmtUOxTYzM7OOrtJBTQMwvyBtHtCzSN51gcOAY4CBwH+BccUWGhFjImJIRAzp379/GxbXzMzMakWl29QsBHoVpPUCFhTJ+z5we0Q8ASDpF8BsSb0jwqNemZmZ2XIqXVMzDegmaZO8tMHA1CJ5/wFE3vsoksfMrKI8orBZx1XRoCYiFgG3AWdI6iFpB2Bv4Loi2a8G9pG0taSVgNOAKa6lMbNq8ojCZh1X2UGNkuGSfiPpaknrZ+k7SfpUM9Z5JLAqMJPURuaHETFV0lBJC3OZIuLPwM+Au7K8GwMlx7QxMzOzzq2sNjWS+gJ3A18itX9pAC4EXgYOB+YAPy5nWRExBxhRJH1yttz8tEuAS8pZrpnVLkldgC4RsSQvbXdgS+DPEfH3qhXOzGpGuTU1vwbWA3YA+gHKm/YA8D9tXC4z61zGAVfl3kg6AphEOvc8JmmXahXMzGpHuUHN3sApEfFXVmyw+wop4DEza6ltSbXBOT8FrgB6k9rhnVKNQplZbSk3qGkAXi8xbRWWr7kxM2uuNcnOMZI2BjYALoqIBaROA5+tYtnMrEaUG9S8COxWYtpOwHNtUxwz66Tmk25tA+wMzI6If2TvPyZdPJmZNarcwfd+D1wkaR5wY5bWR9J3gKOBUe1RODPrNB4FTpK0BDiW5W9FbQy8Vo1CmVltKaumJiLGAOcBvwCmZ8n3A2OA30XEDe1TPDPrJE4g1dRMINXKjM6b9i3gr1Uok5nVmLIfkxARJ0m6BNiVdP/7HeD+iPhPexXOzDqHiPg3sImkfhHxTsHkY4C3KlEOSQcAF0SEHyJnVoOaDGokrUw6oYyMiAmkHglmZm2uSEBDRFSkzZ6krsA3gVcrsT4za3tNBjURsTi7z/1BBcpjZp2EpNObkT0i4n/brTDJAcDNwPHtvB4zayfl3n4aD+wH3Nd+RTGzTmZ0wfug+PAQubGxmgxqJB0NjCR1AR8XESPzpq0OXEnqyTkbODkibsymdQX2J4127qDGrEaVG9RMAi6QdAspwHmTgkH4smc1mZmVJSKWdVSQtDmpkfAY4A/A28BapNqTw4G9ylzsG8CZwO6kZ8zluxhYnC13a+AuSc9GxFTgYOCmiFgqedgts1pVblBza/b3G9krJ3dlFUDXNiyXmXUuFwFXRMQ5eWmvAGdnz4W6mDIexxIRtwFIGgKsm0uX1APYF9gyIhYCUyRNAA4BTgI2Bz4n6WBSg+ULIqKs59mZWcdRblDzlbZaYWNVwAX5RpOGRv8wL3kr97Yyq0tfAs4qMe0J4NRWLn9TYElETMtLe5Y0eCgRcWIuUdKTpQIaSaPIxuUaOHBgK4tkZm2trKAmIh5qw3U2VgVc6I8RcXAbrtvMOqZ5pOEiHigybbdsems0kEYtLlxnz8KMETGk1EKyMbvGAAwZMqTwOXhmVmVlj1MDy2pZtgNWB+YAf42IOc2Yv6kqYDPrnK4CTpbUQOqBlGtTsz+pZqRULU65FgK9CtJ6AQuauyBJw4BhG2+8cSuLZGZtrdxnPyHpTNID5yYA1wATgdclNaebZakq4C1K5B8maY6kqZJ+2Iz1mFltOZ0UuBwK/Al4Hvgz6YLnLFbsKdVc04BukjbJSxsMFKshblRETIyIUb17925lkcysrZVVUyPpWOBnpLYw15MG41ub1GPgZ5JmRcQFZSyq7Cpg4CZSNe/bpPvtt0qaGxHjipTP97nNalhELAVOk3QuqTv2AFIvy39ERNm3niR1I53XugJdJa1CupBaJOk24AxJ3yfd+t4b2L5tP4mZVVO5t5+OAM6PiJ/kpb0IPCRpIXAkUE5QU3YVcES8kPf2UUnnk8bKWSGo8X1us/oQEXOBya1YxKnAz/PeH0x6Zt1o0nnqKmAm6TEvPyzRlq9Rvv1k1nGVG9QMAu4qMe0uoNxbQ8uqgLNnvUD5VcClBuYyszqQdd3+IjCQ9FDL5UTEtU0tIyJGU+JWVdb+b0RrypgtZyIwcciQIYe3dllm1rbKDWreAbakeM+ELbLpTWpOFbCkvYGHgbnANsCPSbfAzKzOZIPvjQc2ovSowk0GNWbWuZUb1NwO/K+kd0hDjy/J7l1/EziD1HC4XEWrgCUNBSZFREOW79tZvu7Aa8DZEdGc9ZhZ7fg96Xy0P/Acy49P1aH49pNZx1VuUHMy6TbRNcBVkuaQunV3BabQjBqUUlXAETGZ1JA49/6AcpdpZjXv88DI3IjAHZlvP1lHsMOFO7R43kd+9EgblqRjKXfwvQWSdgT2BIbyyTg1D5FqV9w418xaYzZpUE4zsxYre/C9LHC5M3uZmbWl3wJHSZoUER9XuzBmVpvKHadmL2BQRFxUZNpRwH8j4u62LpyZdRr9gU8DL0i6n1QTnC8i4ucrzlZ5HalNzStnfLblM/ctHF3DrPaVW1NzGlDqXveq2XQHNWbWUvkPrNykyPRg+fFnqsZtasw6rnIfk/AZ4OkS054BNmuT0phZpxQRXZp4da12Gc2s4ys3qOlCXs+kAj2BldqmOGZmZmYtU+7tp2eBg0jj1RQ6CPhHm5XIzDqtrP3eTnzSw/LBiCg1mnlVdKQ2NWa2vHKDmnNJD5S8GbicNBjeOqSHSO5DGoTPrNNrzdgRB3BSG5aktkjqSepZORRYQhqYsx9wnKTJwF4RsbCKRVzGbWrMOq6ybj9FxO3AMcDuwCTSiJ/3Zu9/XAsDZplZh3YWaQC+Q4BVI2IAqRPCoVn6WVUsm5nViOaMU3OhpLGk5zT1Iw2W9WhHuXoys5q2L3BqRNyQS8jGq7lB0hrACaTnv5mZlVR2UANpZGFSDY2ZWVvqB7xQYtoL2XQzs0aVvP0kaQ1JWxVJ30zSTZKel3S/pK+1bxHNrBP4L7BXiWl7ZNPNzBrVWE3NmcB2pAdZAiBpAPAI0IvUI2owMFHSLhHxUHsW1Mzq2mXAuZIagBuAN4G1gW8D3weOq2LZluPeT2YdV2MNhbcHxhWk/QToDewTEV8ANgD+Dvy03BVKWl3S7ZIWSXpZ0oFN5F9Z0j8lvVbuOsystkTEb4FfkYaIuB94HvgTcBjwq4g4v4rFW05ETIyIUb179652UcysQGM1NeuSTiz59gD+lXVpJCIWSboQ+E0z1nkx6Wm8awFbA3dJejYippbI/1NgFmmQPzOrUxHxM0m/Brblk3FqHouId6tbMjOrFY0FNSsDi3JvJPUhPQ7hsoJ8M4A+5axMUg9SL4cts15TUyRNIHXjXGGQDkkbAAeTqp4vL2cdZla7sgBmUrXLYWa1qbHbTy+T154G2Dn7+3BBvj5AuVdSmwJLImJaXtqzwBYl8l8I/Ax4v8zlm1kNknRiVutbbNoFksq+xW1mnVdjQc0twEmS9pK0DekJuYtY8SpqB8rvmdAAzC9Im0eRW0uS9gG6ZgP/NUrSKElPSnpy1qxZZRbFzDqQ71D6cSvPZNPNzBrVWFDzG+AVYALwGKk25diImJfLIKkbqWHfn8tc30JSz6l8vYAF+QnZbapzKHOwrYgYExFDImJI//79yyyKmXUgA4F/l5j2H2D9CpalUZKGSRozb968pjObWUWVbFMTEQskbcsnD5d7OiIKa2R6AceSgp5yTAO6SdokInInsMFAYSPhTYBBwGRJkNr39Jb0FrBtRMwoc31mVhveIz1Prph1gQ8rWJZG+dlPZh1XoyMKR8RS4C+NTJ8D3FruyrLeUrcBZ0j6Pqn3096k7uP5ngfWy3u/PXAR6Rkwvr9kVn8mAz+VdEtELAtgJHUHjs+mm5k1qlmPSWgjRwJXATNJT+L9YURMlTQUmBQRDRGxBHgrN4OkOcDSiHir6BLNrNaNBh4Fpkm6HnidVHNzMOkRCSOrVjIzqxkVD2qy2p0RRdInkxoSF5vnQVIVtJnVoYh4VtJXSG35TiS191sKTAH2jYhnq1k+M6sN1aipMTNbQUT8DdhR0qpAX+DdiPBwDmZWtsZ6P5mZVUNXYCVgSbULYma1xUGNmXUI2ZhYT5PGrnoJ+GyWfkVTz4gzMwMHNWbWAUgaAdwBzOaTNjU5/yU92NLMrFGtDmok7Svp47YojJl1Wj8Hro6I3YDfFUx7Htiy4iUys5rjmhoz6wg2A/6Y/R8F094ldevuEDyisFnHVbL3k6RDy1zGNm1UFjPrvOYDa5SYNogONOimRxQ267ga69I9lnTFpDKWU3hlZWbWHPcDJ0uaxCfPgotsROGjWfFBumZmK2gsqJkDTATObGIZXwfOb7MSmVlndArwN+BF4G7ShdJJwFZAb4oM2GlmVqixoOYpYMOIeKmxBUh6s22LZGadTUTMkPR54BfA7sDHwI7APcDpEfFGNctnZrWhqaDm6DKWMQt4uG2KY2adVUS8Bnyv2uUws9pVsvdTRPwsIno1tYCIeDgivtK2xTKzzkxSb0lDJPmZb2ZWNnfpNrOqkLS7pF8VST8FmAk8Drws6UZJfk6dmTWpZFAj6auSij41uzUkrS7pdkmLJL1cavhzST+R9B9J8yW9Iem3PrGZ1ZUjgE3zEyTtCvwv8C/gWOAy4FvAMZUunJnVnsZqau4HNs+9kdRF0sOSNmnlOi8GFgNrAQcBl0jaoki+CcDns1tgWwKDgR+3ct1m1nF8DrirIO07wAfA7hFxYUQcSQps/OwnM2tSY0FN4fg0Ar4M9GzpyiT1APYFTouIhRExhRS8HFKYNyJeioi5eeteCmzc0nWbWYezJunBlfl2BaZExFt5aXdRUKPT1iStJelRSQ9J+rOkAe25PjNrH5VuU7MpsCQipuWlPQsUq6lB0oGS5pMecjeYdMVWLN8oSU9KenLWrA4z8KiZNW4B0CP3JqsF7gc8VpBvPtC1ncsyG/hyROwEXIt7YZnVpEoHNQ2kE1S+eZSo/YmIG7PbT5sClwJvl8g3JiKGRMSQ/v37t2V5zaz9/AvYO+/93qRB9+4ryLcBJb77bSUiPo6IpdnbnsDU9lyfmbWPpoKadSRtKGlDYMPCtPxXmetbCBR2E+/FJ8OiFxUR/yadZH5f5nrMrOP7LfB9SbdIupg08N5zwCMF+fYg1eg2SdLRWa3th5LGFkxrtJOCpK0lPU4an+vpFn4mM6uipnoT3VIkbXyJvOVUD08DuknaJAtUIN1WKueqqBuwURn5zKwGRMR4SccCxwOrk247HRERy54lJ2ltYBfgZ2Uu9g3So112B1YtmJbfSWFr4C5Jz0bE1Kw8zwBfkrQ/cDKpd5aZ1ZDGgprvtPXKImKRpNuAMyR9n3Ri2RvYvjBvNn1CRMyUtDnpJHNvW5fJzKonIi4ALmhk+luUfnp3sfy3AUgaAiwbuC+vk8KWEbEQmCIp10nhJEkrR8TiLPs84L3mfhYzq76SQU1EXNNO6zwSuIo0uNY7wA8jYqqkocCkiMiNjbMD8MtsrJxZwM3Aae1UJjOrb6U6KeyU/b+1pN+Qnjn1AfDdYguRNAoYBTBw4MD2K62ZtUjFB7OLiDkUeeJuREwmNSTOvW/zmiIz67Qa7aQQEX8jPUCzURExBhgDMGTIkGgiu5lVmB+TYGadQYs6KZhZbXFQY2adwbJOCnlp5XZSWI6kYZLGzJs3r80KZ2Ztw0GNmdUNSd0krULqjdlV0iqSukXEIiDXSaGHpB1InRSua+46ImJiRIzq3bt32xbezFrNQY2Z1ZNTgfeBk4CDs/9PzaYdSermPRMYR9ZJobkrcE2NWcfloMbM6kZEjI4IFbxGZ9PmRMSIiOgREQMj4sYWrsM1NWYdlIMaMzMzqwsOaszMzKwuOKgxM2sGt6kx67gc1JiZNYPb1Jh1XA5qzMzMrC44qDEzawbffjLruBzUmJk1g28/mXVcFX+gpaTVgSuB3YDZwMnFxouQ9FPgMGD9LN/vI+LXlSyrmZmt6KLjJ7Z43qPPHdaGJTFbXsWDGuBiYDGwFrA1cJekZ4uM7CngUOAfwEbAfZJejYg/VLKwZmZmVhsqevtJUg9gX+C0iFgYEVOACcAhhXkj4pyIeDoilkTEi8AdwA6VLK+ZWSG3qTHruCrdpmZTYElETMtLexbYorGZJAkYSgueqGtm1pbcpsas46p0UNMAzC9Imwf0bGK+0aSyXl1soqRRkp6U9OSsWbNaXUgzMzOrPZUOahYCvQrSegELSs0g6WhS25o9I+LDYnkiYkxEDImIIf3792+zwpqZmVntqHRQMw3oJmmTvLTBlLitJOm7wEnA/0TEaxUon5mZmdWoigY1EbEIuA04Q1IPSTsAewPXFeaVdBBwFrBrRPynkuU0MzOz2lONwfeOBFYFZgLjgB9GxFRJQyUtzMt3JtAPeELSwux1aRXKa2a2jHs/mXVcFR+nJiLmACOKpE8mNSTOvd+ggsUyMytLREwEJg4ZMuTwapfFzJbnxySYmZlZXXBQY2ZmZnXBQY2ZmZnVBQc1ZmZmVheq8UBLa8IOF7b8EVcHcFIblsTMCkkaBgzbeOONq10UMyvgmhozs2bws5/MOi7X1JiZmVXBK2d8tuUz9y184pCBa2rMzMysTjioMTMzs7rgoMbMzMzqgoMaMzMzqwsOaszMzKwuOKgxMzOzulDxoEbS6pJul7RI0suSDiyR7yuS/iJpnqQZFS6mmZmZ1Zhq1NRcDCwG1gIOAi6RtEWRfIuAq4CfVrBsZmaNkjRM0ph58+ZVuyhmVqCiQY2kHsC+wGkRsTAipgATgEMK80bE3yLiOuA/lSyjmVljPKKwWcdV6ZqaTYElETEtL+1ZoFhNjZmZmVnZKh3UNADzC9LmAT1bs1BJoyQ9KenJWbNmtWZRZmZmVqMqHdQsBAofWNELWNCahUbEmIgYEhFD+vfv35pFmZmZWY2qdFAzDegmaZO8tMHA1AqXw8zMzOpMRYOaiFgE3AacIamHpB2AvYHrCvNK6iJpFWCl9FarSFq5kuU1MzOz2lGNLt1HAqsCM4FxwA8jYqqkoZIW5uXbEXgfuBsYmP1/X6ULa2ZmZrWhW6VXGBFzgBFF0ieTGhLn3j8IqGIFMzMzs5rmxySYmZlZXXBQY2adnqQvSvqrpIcljZO0UrXLZGbN56DGzAxeBb4aETsCM0gdGMysxlS8TY2ZWUcTEW/mvV0MLK1WWcys5VxTY2Z1Q9LR2ejiH0oaWzBtdUm3S1ok6WVJBxaZf31gN2BihYpsZm3INTVmVk/eAM4EdicNHZHvYlItzFrA1sBdkp6NiKkAknqRxswaGREfVazEZtZmXFNjZnUjIm6LiPHAO/npknoA+wKnRcTCiJgCTAAOyaZ3A/4A/CIiXqxsqc2srbimxsw6g02BJRExLS/tWWCn7P8DgC8Bp0k6DbgkIv5YuBBJo4BRAAMHDmzfEpu1k4uOb/nd1aPPHdaGJWl7DmrMrDNoAOYXpM0DegJExHUUeVxLoYgYA4wBGDJkSLRxGc2slXz7ycw6g4VAr4K0XsCCKpTFzNqJgxoz6wymAd0kbZKXNhiY2twFSRomacy8efParHBm1jYc1JhZ3ZDUTdIqQFegq6RVJHWLiEXAbcAZknpI2oE0wF6Tt5wKRcTEiBjVu3fvti28mbWagxozqyenAu8DJwEHZ/+fmk07ktTNeyYwDvhhrju3mdWHigc15QyAleWTpLMlvZO9zpbkp3abWUkRMToiVPAanU2bExEjIqJHRAyMiBtbsg7ffjLruKpRU5M/ANZBwCWStiiSbxQwgnTfeytgGPCDCpXRzKwo334y67gqGtQ0NQBWgcOAcyPitYh4HTgXGFmxwpqZmVlNUUTlhlqQ9DngkYhYLS/t/wE7RcSwgrzzgN0i4vHs/RDgLxHRs8hylw2IBXwaqOcRQdcAZle7ENZi9b7/1o+I/tUuRHuSNIxUc/wt4N9VLk57qvdjtd7V+/4req6p9OB7jQ6AVSTvvIJ8DZIUBZFY/oBY9U7SkxExpNrlsJbx/qt9ETGR9MDLUU3lrWU+VmtbZ91/lW5T05wBsArz9gIWFgY0ZmZmZlD5oKY5A2BNzaY1lc/MzMysskFNMwfAuhY4TtI6kj4FHA+MrVhhO65OcZutjnn/Wa3wsVrbOuX+q2hDYUjj1ABXAbsC7wAnRcSNkoYCkyKiIcsn4Gzg+9msVwAn+vaTmZmZFVPxoMbMzMysPfgxCWZmZlYXHNR0QJIulxSSfltk2shsWu61SNKM7NET+/tREtXT3H2T5TmzWuW16pI0VtJrJabtnB0fu1S6XHllKDyeS71GZ6+arfZv4rPOLZJvUPVK2zrZOWlsI9MfLHO/D2pqWdVQ6XFqrAmSVgX2z94eKOmnEbGkSNZvAq8B3YGBwJ6kh/SNkjQsIt6vSIGtGO8bqwd3Advlvf886TE3PwaeyEvPBWb3VKhc7Sn33c1X7Pxbz45k+eFUTgO2AYYX5HsT2IcVx56rKgc1Hc8I0gF1N7AH8DXgziL5nomI6Xnvr5N0M3AzcA7wo3Yup5XmfWM1RVL3iPgwPy0iZgGz8vKskv37z4h4rMhiitY61ZjC726nExEv5L+XNAtYXGKf/70ypSqfbz91PIcB75Kec/V+9r4sEXErcAdwuKTVmspvleN9Y60l6UBJf5e0UNJ8Sc9J+kFBnp0k/UnSguz2572StizI86CkKdnTxv8u6UPS1XlryrbC7afc7VVJx0t6WdJ7ku6StGb2uknSPEmvSjqxyDI3kHSDpFmSPpT0jKR9WlPO9iDp25L+nJVzYbZNVzhvZ9vjl5JOkfSapPclPSxp64J8u0t6NNs2CyW9KOn0gjyDJU2Q9G62nEeUehAXrvOY7BbRB5KeLJanlZ99udtPebfnts/27wJJb0s6OZv+tWz7LJL0hKQvFFnmNyQ9lh0vcyXdLGlguWVyUNOBKI3Hswvwx+wqaTwwTFLfZizmbtJtj043PHYN8L6xFpH0ZeB64CFSbe5+wOVAn7w8ewJ/Io3GfjBwIOkRNJMlrVewyE2BC4ALgd2z+drDIcBXSUHT0cBQ0hhktwP/ID3g+G7gV5L2yPss6wGPkwZd/Qnp1sfTwK2SCm+DtKWukroVvJr6ndwQuAU4iLRvJgJXSDqiSN5DSTXwR5MuXNcC/qQ01AmSNiQ95Pm/pGeLDQfOA3rkFiDp88CjwOrA4aRt+A7wQH6QIOl7wO+Av2TlGku6Dd6c35OWugZ4jnR7ajxwlqSzgV+Thmr5FukzjZe0cl6ZjwBuBV4gHeM/ALYEHpJU7HFKK4oIvzrICzgBCGC77P3u2fsj8vKMzNI2LrGM3Dzfqvbn6Wyv5u6b7P8zq11uv6rzIv3IvFZi2s7Z8bFL9v7/AXOaWN504E8Fab1IDzX8XV7ag8BSYOtmlne5MhVMG51+TpZLC7JR5PPSzsvST81L6wbMBK7OS7uSdOurX8Ey7yfdImrrfZH77hZ73Vkk36ASy+mSfZ7LgWeLbI/ZQI+8tEHAR8D/Zu/3y/L1aqSsfwL+Caycl9Y1SxufV45XgXsK5v1WtvyxbXSczshfVt72Ob3I/v0I2CAvfXiWd6fsfe55j1cVrGMDYDFwbDnldU1Nx3IY8O+I+Gv2/gHgDZpxCwrI9bCp2Z4Idcz7xlrqCaCvpOsl7SWpT/5EpUfPbATckF/LALwH/BXYsWB5MyLimQqU+/5YvqPDv7K/9+YSsunTgfzapK+RanDmFXyee4HBkgqfIdhW9iE1is1/HdvYDJI2kTRO0uukH+6PSIPGfrpI9rsjjawPQETMAB7jkwbZz2Tz/0HSfpLWLFjXqsBOpPZ5S/O2i0i/F7n9vG72uqlg/bdSmYbPk3L/5O3faRHx37w8uWMht9+3IwXhhcfwq1newmO4KAc1HYSkIcDmwG2S+mQnrZ6kx0psK2nTMheVO0DebPtSWit531i+JaQr7GK65uUhIh4i9cxZj3TrZpakByRtleXL/fhdySc/rLnXXkC/guVX6hh8t+D94kbSV8l7vybpVk3hZ/l1Nr3w87SV5yPiyYJXyYbDkhpItUeDgZNIt9e2IY2a373ILG+XSFsHIFvX7qTf5uuAt7L2JTtleVcnHRunseK2OZoU+HYBBhRbXxZgvNP4JmgTxfZvqWMht99zx/ADrPjZPkuZ+9y9nzqOXG3Midmr0KHAqWUsZ0/gA+CpNiqXtR3vG8s3E1hD0soRsbhg2qeyv8t+lCLiFuCW7Id0Z1LbhHskrcsnP1Qnk34UChUuv6PXFr4DTCZ9xmLeqGBZGrMdsD4wNCKm5BKzGoZi1iqR9nruTUT8BfiLpO7ADsAZwF1KY+PMJd06vJjUNmkFEbFUUi5oXW59WbnaKyBsrdwxPJLiD69eUM5CHNR0AFlDqQNIDeNOKpLlt8Ahkk5rYjn7ku5Tnh8R77V5Qa3FvG+siL+QgpDhpIam+fYl1aa8WDhTRCwE7swalZ5P+pF6kdS+YYuI+FU7lrlS7iEFDFOjY4/rlOvJ+FEuIevYsXeJ/HtI6pG7BZUFKtsCK+yzSF3s/5wFsXeQ2qM8IWkyqWbo6YhYWmI9r5Fu2+xPqjXK2ZeO+7v/KClw2TgirmnpQjrqh+ts9iSdmI6PiAcLJ0q6DLiEdHWWs7WkNYCVSQO87UWqnr6fdKK06mnOvvmMpP2KLONPEVFYXWv15QHSMTFW0mdIFzU9gW+TfhS/k/vRknQG6ar7L6RainVJg+A9E6mnJJKOAu7ILpJuIjVKXQvYHnglIs6r4GdrrdOBvwEPS7qIFLD1JfWE2TAivttO6819dws9GcUHQX2UNPjcxZJ+TurRcypp2/cukv994D5JvybdnvpFNv9vYVnvnx1J7YleBdYgnTPeAJ7PlnEc8DBwr6QrScHvGqTBEbtGxElZbc0vSL2wrgb+AGxMumjuUIPl5UTEfEk/JW3L/qR2OfNIt+Z2Ah6MiBubWo6Dmo7hMFKEenOJ6eNIvQYOI/VcIC/vB6Rq7KdJJ8NbImsyblXTnH2zb/YqtA3wZLuV0KouIkLS3sAppNvLp5FuEz0DjIiIO/KyP04KYn5LalcxE7gvmye3vLsl7Zgt7wpgVeAtUkPUP7b352lLEfFK1s5wNHAW0J90e+J5Unfh9lLqHNyfFKgsJyJmKY2dcy6ptu0NUu3Z6sDPiyznWmARcBEpEHkC+HZEzMmmPwt8Hfg/UhuTOcAU4KBcjVVEPC1pm2z5F5CCp1mk88yleWW7MqvlOY50J+D57O/15WyIaoiIyyS9CvyUNCRBN9Ktucmk70WT/JRuMzOzdqY0OOEvI6KctpHWQu79ZGZmZnXBQY2ZmZnVBd9+MjMzs7rgmhozMzOrCw5qzMzMrC44qDEzM7O64KDGzMzM6oKDGjMzM6sLDmrMzMysLjioMTMzs7rgoMbMzMzqgoMaMzMzqwt+Sncn99RTT3Xv2rXrAx9//PH2OMg1s/qxtGvXro9+/PHHu3zhC1/4sNqFscpwUNPJdenS5X979Ojx+Y022qhLly6OacysPixdurTL9OnTh8yePftc4Ohql8cqw79inVxE/GDgwIGrOaAxs3rSpUsX1l9//VW6des2avjw4V+tdnmsMvxL1slFRM+VV1652sUwM2tzK6+8Ml26dFkJOHT48OEbVbs81v4c1JgkVbsMZmZtThLZ+S2AtapcHKsABzVm1qHMmDEDSSxZsgSAnXfemSuuuKLJ+R588EHWXXfd9i5ep7TFFlvw4IMPVrsYrdW12gWw9uegxjqsQYMG8cADDyyXNnbsWL785S+3etkjR47k1FNPbfVyrHE33ngjQ4YMoaGhgQEDBvD1r3+dKVOmtHh5Y8eOpWvXrjQ0NNCrVy+23npr7rzzzjYscf2bMmUK22+/Pb1792b11Vdnhx124Iknnmh0nqlTp7LzzjuXnD558mQaGhpoaGigR48eSFr2vqGhgVdeeaVNyj5z5kwOOOAAPvWpT9G7d2922GEHHn/88TZZttUH936yFbxyxmfbdfkDT3+uXZffFpYsWUK3brX39djhwh3adfmP/OiRsvOed955/OpXv+LSSy9l9913Z+WVV+aee+7hjjvuaFVgut122zFlyhSWLl3KxRdfzP7778/rr7/e4uVVwkXHT2zX5R997rCy8s2fP5+99tqLSy65hP3335/FixczefJkunfv3qr1Dx06lIULFwKppm2DDTZg7ty5Rb9DrfluLVy4kG222YbzzjuPNddckyuvvJI999yTGTNm0NDQ0KrPYPXBNTVWsyQxffr0Ze/za19ytyLOPfdc1lxzTQYMGMDVV18NwJgxY7jhhhs455xzaGhoYNiw9IMwaNAgzj77bLbaait69OjBkiVLeOyxx9h+++3p06cPgwcPXq4KfuzYsWy44Yb07NmTDTbYgBtuuAGA6dOns9NOO9G7d2/WWGMNvvWtb1Voi3Qc8+bN4/TTT+fiiy/mG9/4Bj169GCllVZi2LBh/PrXv2bp0qX86le/YqONNqJfv37sv//+zJkzp1nr6NKlC9/97nd5//33eemll5alF9vnTzzxBGuttRYff/zxsny33XYbgwcPbpsPXCOmTZsGwAEHHEDXrl1ZddVV2W233dhqq60AuPzyy9lss83o2bMnm2++OU8//TSwfK3p3/72N4YMGUKvXr1Ya621OO644xpd5+jRo9lvv/04+OCD6dWrF2PHjmXevHl873vfY8CAAayzzjqceuqpy+2bq666is0224y+ffuy++678/LLLwOw4YYbctxxxzFgwAC6du3KqFGjWLx4MS+++GKbbyurTQ5qrG699dZbzJs3j9dff50rr7ySo446infffZdRo0Zx0EEHccIJJ7Bw4UImTvzkKnrcuHHcddddzJ07l7fffps999yTU089lTlz5vCb3/yGfffdl1mzZrFo0SJ+/OMfM2nSJBYsWMCjjz7K1ltvDcBpp53Gbrvtxrvvvstrr73Gj370oyptger561//ygcffMA+++xTdPqFF17I+PHjeeihh3jjjTfo27cvRx11VLPWsWTJEq644goaGhrYZJNNgNL7fJtttqFfv37cd999y+a/7rrrOPTQQ1v+IWvQpptuSteuXTnssMOYNGkS77777rJpN998M6NHj+baa69l/vz5TJgwgX79+q2wjGOOOYZjjjmG+fPn89JLL7H//vs3ud477riD/fbbj7lz53LQQQcxcuRIunXrxvTp0/n73//Offfdt6zd1B133MFZZ53FbbfdxqxZsxg6dCgHHHBA0eU+88wzLF68mI033riFW8TqjYMa69BGjBhBnz59lr2OPPLIsuddaaWVOP3001lppZXYY489aGhoaPKK7sc//jHrrbceq666Ktdffz177LEHe+yxB126dGHXXXdlyJAh3H333UCqKXj++ed5//33GTBgAFtsscWy9b788su88cYbrLLKKm3SBqjWvPPOO6yxxholbzNceuml/PKXv2Tdddele/fujB49mltuuWVZ4+DGPPbYY/Tp04e1116bcePGcfvtt9O7d2+g8X1+2GGHcf311wMwZ84c7r33Xg488MA2+sS1oVevXkyZMgVJHH744fTv35/hw4fz9ttvc8UVV3DCCSewzTbbIImNN96Y9ddff4VlrLTSSkyfPp3Zs2fT0NDAtttu2+R6t9tuO0aMGEGXLl2YP38+d999N7/73e/o0aMHa665Jj/5yU/4wx/+AKRj4+STT2azzTajW7du/OxnP+OZZ55ZVluTM3/+fA455BB+/vOfL9v/Zg5qrEMbP348c+fOXfb6/e9/X/a8/fr1W+5HdbXVVlt237+U9dZbb9n/L7/8MjfffPNyQdWUKVN488036dGjB3/84x+59NJLGTBgAHvuuSf/+te/ADjnnHOICL74xS+yxRZbcNVVVzXzU9e+fv36MXv27JJByssvv8w+++yzbLtuttlmdO3albfffrvJZW+77bbMnTuX2bNn89hjj7HLLrsst95S+/zggw9m4sSJLFq0iJtuuomhQ4cyYMCAVn7S2rPZZpsxduxYXnvtNZ5//nneeOMNjj32WF599VU22qjpoVyuvPJKpk2bxmc+8xm22WabshpqF36vPvroIwYMGLBs///gBz9g5syZy6Yfc8wxy6atvvrqRMRy7abef/99hg0bxrbbbsvJJ5/cgq1g9cpBjdWs1VZbjffee2/Z+7feeqvseUuNzZOfvt5663HIIYcsF1QtWrSIk046CYDdd9+d+++/nzfffJPPfOYzHH744QCsvfbaXH755bzxxhtcdtllHHnkkcu1/ekMtttuO7p378748eOLTl9vvfWYNGnSctv2gw8+YJ111mm3Mq2zzjpst9123HbbbVx33XUccsgh7bauWvGZz3yGkSNH8vzzz7Peeust1zaplE022YRx48Yxc+ZMTjzxRPbbbz8WLVrU6DyF36vu3bsze/bsZft+/vz5TJ06ddn0yy67bLlj4/3332f77bcH4MMPP2TEiBGsu+66XHbZZa349FaPHNRYzdp666258cYb+fjjj7nnnnt46KGHyp53rbXW4j//+U+jeXJX9vfeey8ff/wxH3zwAQ8++CCvvfYab7/9NnfccQeLFi2ie/fuNDQ0kHvUxM0338xrr70GQN++fZFEZ3sMRe/evTnjjDM46qijGD9+PO+99x4fffQRkyZN4oQTTuCII47glFNOWXZLYdasWdxxxx3tXq5DDz2Uc845h+eee45vfOMb7b6+juZf//oX55577rLj89VXX2XcuHFsu+22fP/73+c3v/kNTz31FBHB9OnTV7jlA3D99dcza9YsunTpQp8+fQCadXwPGDCA3XbbjeOPP5758+ezdOlSXnrppWXf3yOOOIL/+7//WxbkzJs3j5tvvhmAjz76iP32249VV12Va665ptN9r6xpPiKsZp1//vlMnDiRPn36cMMNNzBixIiy5/3e977HCy+8QJ8+fUrOt9566y1rtNi/f3/WW2+9ZT13li5dynnnncenPvUpVl99dR566CEuueQSIPW0+dKXvkRDQwPDhw/n/PPPZ8MNN2yDT1xbjj/+eM477zzOPPPMZdvvoosuYsSIERxzzDEMHz6c3XbbjZ49e7LttttWZLyRffbZZ9mtr9VWW63d19fR9OzZk8cff5wvfelL9OjRg2233ZYtt9ySc889l29+85uccsopHHjggfTs2ZMRI0YU7ZF2zz33sMUWW9DQ0MAxxxzDH/7wB1ZdddVmlePaa69l8eLFbL755vTt25f99tuPN998E0j76MQTT+Tb3/42vXr1Ysstt2TSpEkAPProo9x5553cd9999OnTZ9k4OJMnT279xrG6oIiodhmsip566qn4whe+UO1imFXMRhttxGWXXbZcWxyrX0899RS/+MUvxgJXTZgwwdFPnXNNjZl1GrfeeiuS+OpX/dBms3pUe0Ommpm1wM4778wLL7zAdddd57YYZnXKQY2ZdQp18EBGM2uCL1fMzMysLjiosXBjcTOrRxGBz2+di4OaTk7SgsWLF1e7GGZmbW7x4sVExIfVLodVjoOaTk7SZTNmzFi8dOnSahfFzKzNLF26lP/+978fvfTSS09mSR83OoPVBTcU7uSWLl162jvvvPO1BQsWfLbUowPMzGpNRMT8+fPfvOWWWx4B1gKafrCY1TwPvmcMHz68K3AwsDO+mjGz+tIFuH7ChAl/rnZBrP05qDEAhg8f3gVYG+hZ7bKYmbWRAOZMmDBhdrULYpXhoMbMzMzqghsKm5mZWV1wUGNmZmZ1wUGNmZmZ1YX/D8XNBpxghK0vAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# time plot\n",
    "font = {\n",
    "'size' : 22}\n",
    "\n",
    "# matplotlib.rc('font', **font)\n",
    "tab10 = sns.color_palette('tab10')\n",
    "palette = [tab10[1], tab10[2], tab10[4]]\n",
    "\n",
    "fig, axes = plt.subplots(figsize=(8,4), nrows=1, ncols=2)\n",
    "\n",
    "g1 = sns.barplot(ax=axes[0], data=df_result[df_result['Measure'].isin(['Ancestor-Descendent F1-Score', 'Different-Lineage F1-Score'])], x='Measure', y='Value', hue='Method', width=0.4, palette=palette, hue_order=['HUNTRESS', 'CellPhy', 'ScisTree'])\n",
    "axes[0].set_title('Performance',  fontsize=16)\n",
    "axes[0].set_xlabel('')\n",
    "axes[0].set_xticks([0, 1])\n",
    "axes[0].set_xticklabels(['AD', 'DL'],  fontsize=16)\n",
    "axes[0].set_ylabel('F1 Score',  fontsize=16)\n",
    "axes[0].set_yticklabels(axes[0].get_yticklabels(), fontsize=12)\n",
    "axes[0].legend_ = None\n",
    "# axes[0].legend_.set_title('')\n",
    "# axes[0].legend(loc='lower right')\n",
    "handles, labels = axes[0].get_legend_handles_labels()\n",
    "g2 = sns.barplot(ax=axes[1], data=df_time, x='Measure', y='Value', hue='Method', width=0.4, palette=palette, hue_order=['HUNTRESS', 'CellPhy', 'ScisTree'])\n",
    "axes[1].set_title('Running Time',  fontsize=16)\n",
    "axes[1].set_xlabel('')\n",
    "axes[1].set_xticklabels(['User Time', 'Elapsed Time'], fontsize=16)\n",
    "axes[1].set_ylabel('Seconds',  fontsize=16)\n",
    "axes[1].set_yscale('log')\n",
    "axes[1].set_yticklabels(axes[1].get_yticklabels(), fontsize=12)\n",
    "# axes[1].legend_.set_title('')\n",
    "# axes[1].legend(loc='upper right')\n",
    "axes[1].legend_.remove()\n",
    "fig.legend(loc='lower center', bbox_to_anchor=(0.5, -0.1),\n",
    "          fancybox=True, shadow=True, ncol=5, handles=handles, labels=['Huntress', 'CellPhy', 'ScisTree2'], fontsize=12)\n",
    "plt.suptitle('Comparion of performance and efficiency on HGSOC', fontsize=18)\n",
    "plt.tight_layout()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "work",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}