{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mito Clones \n",
    "\n",
    "Using vireo for clonal reconstruction - mitochondrial mutations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Date: 15/03/2022"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mitochondrial mutations data set is extracted from [Ludwig et al, Cell, 2019](https://doi.org/10.1016/j.cell.2019.01.022), the 9 variants used here are from Supp Fig. 2F (and main Fig. 2F).\n",
    "\n",
    "Generally, you can use [cellSNP-lite](https://github.com/single-cell-genetics/cellsnp-lite) to genotype mitochondrial genomes and call clonal informed mtDNA variants with [MQuad](https://github.com/single-cell-genetics/MQuad)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the filtered variants at hand, we can use the `vireoSNP.BinomMixtureVB` class to reconstruct the clonality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5.6\n"
     ]
    }
   ],
   "source": [
    "import vireoSNP\n",
    "import numpy as np\n",
    "from scipy import sparse\n",
    "from scipy.io import mmread\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(vireoSNP.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.set_printoptions(formatter={'float': lambda x: format(x, '.5f')})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "AD = mmread(\"../data/mitoDNA/cellSNP.tag.AD.mtx\").tocsc()\n",
    "DP = mmread(\"../data/mitoDNA/cellSNP.tag.DP.mtx\").tocsc()\n",
    "# mtSNP_ids = np.genfromtxt('../data/mitoDNA/passed_variant_names.txt', dtype='str')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identify clones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vireoSNP import BinomMixtureVB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-190779.74335041404\n"
     ]
    }
   ],
   "source": [
    "_model = BinomMixtureVB(n_var=AD.shape[0], n_cell=AD.shape[1], n_donor=3)\n",
    "_model.fit(AD, DP, min_iter=30, n_init=50)\n",
    "print(_model.ELBO_iters[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Check the model fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAEYCAYAAADMNRC5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2UUlEQVR4nO3deZhcZZ33//cn6STshiVihoCJjygDiCwt4iOOyGZUNh8RYRwJwhgRcBwdgSD+BLcRV0YHtyjB6KCAKBJZhIAobgEChECAQIAwJiKEnYCmq6u/vz/OXclJUZVUd7r6dJ3+vK7rXH3OfbZvna6+q++6N0UEZmZmZmZmrRhVdABmZmZmZtY5XIAwMzMzM7OWuQBhZmZmZmYtcwHCzMzMzMxa5gKEmZmZmZm1rKvoAFqxzTbbxOTJk4sOw8ysMLfeeuvjETGh6DhGEn/2mNlI1+yzpyMKEJMnT2b+/PlFh2FmVhhJDxcdw0jjzx4zG+maffa4CZOZmZmZmbXMBQgzMzMzM2uZCxBmZmZmZtYyFyDMzMzMzKxlLkCYmZmZmVnLXIAwMzMzM7OWta0AIWkjSTdLukPSIkmfTuk/kPSQpAVp2b1dMZiZmZmZ2eBq5zwQq4D9I2KlpDHA7yVdnfadGhGXtvHeZmZmZmbWBm0rQEREACvT5pi0RLvuZ2ZmNhy90NPLLUuf4omVq3jy+R4eX9nDEytXsaq3r+jQzGwE2Gni5py03ysH9ZptnYla0mjgVuCVwDcj4iZJHwI+L+lTwPXAjIhY1c44zMzMivKlXy3mB39cunp7zGix1aZj2WRsWz+CzcwA2GjM4PdYaGvuFRFVYHdJ44HLJO0KnAH8FRgLzAROBz5Tf66k6cB0gB122KGdYZrZCDR5xpVDfs+l57xjyO9pxXv2bxW23WIcF01/A1tvNpbNx3UhqeiwzMwGbEhGYYqIp4EbgKkR8UhkVgEXAHs3OWdmRHRHRPeECROGIkwzM7NB11PtY9NxXUzZZlO22GiMCw9m1vHaOQrThFTzgKSNgYOAeyVNTGkCjgDualcMZmZmRatU+xg72qOmm1l5tLMJ00RgduoHMQq4JCKukPRrSRMAAQuAE9sYg5mZWaEq1aBrtGsdzKw82jkK00Jgjwbp+7frnmZmZsNNpdrHGNdAmFmJOEczMzNrIxcgzKxsnKOZmZm1UaUa7gNhZqXiHM3MzKyNKtU+94Ews1JxAcLMzKyNKtVwEyYzKxXnaGZmZm3kYVzNrGyco5mZmbVR1onaTZjMrDxcgDAzM2ujSm8fXa6BMLMScY5mZmbWRpU+94Ews3JxjmZmZtZGWR8IN2Eys/JwAcLMzKyNKr2eSM7MysU5mpmZWRtVquE+EGZWKs7RzMzM2iQiqPS5CZOZlYsLEGZm1naSvizpXkkLJV0maXxu3xmSlkhaLOmtufSpKW2JpBm59CmSbkrpF0sam9LHpe0laf/k9d2j3ap9QQRuwmRmpeIczczMhsJcYNeI2A24DzgDQNLOwNHALsBU4FuSRksaDXwTeBuwM3BMOhbgi8C5EfFK4CnghJR+AvBUSj83Hdf0Hm1+vUDWfAlgTJc/bs2sPJyjmZlZ20XEtRHRmzbnAZPS+uHARRGxKiIeApYAe6dlSUQ8GBE9wEXA4ZIE7A9cms6fDRyRu9bstH4pcEA6vtk92q6n2gdA1yg3YTKz8nABwszMhtrxwNVpfTvgz7l9y1Jas/StgadzhZFa+lrXSvufScc3u9aLSJouab6k+StWrBjQi8vrTQWIsa6BMLMS6So6ADMzKwdJ1wEva7DrzIi4PB1zJtALXDiUsbUqImYCMwG6u7tjQ6+3ugmT+0CYWYm4AGFmZoMiIg5c135JxwGHAAdERO2f8+XA9rnDJqU0mqQ/AYyX1JVqGfLH1661TFIX8JJ0/Lru0VaVVAPhAoSZlYlzNDMzaztJU4HTgMMi4oXcrjnA0WkEpSnAjsDNwC3AjmnEpbFknaDnpILHDcCR6fxpwOW5a01L60cCv07HN7tH2/WsLkC4D4SZlYdrIMzMbCicB4wD5mb9mpkXESdGxCJJlwB3kzVtOjkiqgCSTgGuAUYDsyJiUbrW6cBFkj4H3A6cn9LPB34kaQnwJFmhg3Xdo9163YTJzErIBQgzM2u7NLRqs32fBz7fIP0q4KoG6Q/SYBSliPg78O7+3KPd3ITJzMrIOZqZmVmbuAmTmZWRCxBmZmZtUul1DYSZlU/bcjRJG0m6WdIdkhZJ+nRKnyLpJklLJF2cOseZmZmVTm+f+0CYWfm0M0dbBewfEa8FdgemStoH+CJwbmoP+xRwQhtjMDMzK4ybMJlZGbWtABGZlWlzTFoC2B+4NKXPBo5oVwxmZmZFchMmMyujtuZokkZLWgA8BswFHgCeTpP/ACwDtmtnDGZmZkXxTNRmVkZtzdEiohoRu5PN+rk3sFOr50qaLmm+pPkrVqxoV4hmZmZt09vnJkxmVj5D8pVIRDxNNnPoG4DxkmrzT0wCljc5Z2ZEdEdE94QJE4YiTDMzs0HV4yZMZlZC7RyFaYKk8Wl9Y+Ag4B6ygsSR6bBpwOXtisHMzKxItSZMY7tcgDCz8mjnTNQTgdmSRpMVVC6JiCsk3Q1cJOlzwO3A+W2MwczMBkDSdsDLyX1ORMSNxUXUmWozUXeNchMmMyuPthUgImIhsEeD9AfJ+kOYmdkwJOmLwHuAu4FqSg7ABYh+qhUgxrgGwsxKpJ01EGZm1pmOAF4dEauKDqTTrW7C5D4QZlYiztHMzKzeg2Rz99gGchMmMysj10CYmVm9F4AFkq4HVtdCRMS/FRdSZ6pU+5BgtAsQZlYiLkCYmVm9OWmxDVSpBmNGj0JyAcLMysMFCDMzW0tEzJY0FnhVSlocEZUiY+pUlWqf+z+YWem4AGFmZmuRtB8wG1gKCNhe0jQP49p/lWofXZ6F2sxKxgUIMzOr91Xg4IhYDCDpVcBPgL0KjaoDVap9noXazErHuZqZmdUbUys8AETEfXhUpgGpVMNNmMysdFwDYWZm9eZL+j7wP2n7vcD8AuPpWFkNhJswmVm5uABhZmb1PgScDNSGbf0d8K3iwulcWR8I10CYWbm4AGFmZmtJM1B/LS22AXp6w30gzKx0XIAwMzMAJF0SEUdJuhOI+v0RsVsBYXW03r4+xroJk5mVjAsQZmZW85H085BCoygRj8JkZmXkXM3MzACIiEfS6kkR8XB+AU4qMrZOVekNzwNhZqXjAoSZmdU7qEHa24Y8ihLocQ2EmZWQmzCZmRkAkj5EVtPwCkkLc7s2B/5QTFSdLesD4QKEmZWLCxBmZlbzY+Bq4AvAjFz6cxHxZDEhdbaKR2EysxJyAcLMzACIiGeAZ4BjACS9FNgI2EzSZhHxv0XG14myeSDcB8LMysVfi5iZ2VokHSrpfuAh4LfAUrKaCeunnqqbMJlZ+ThXMzOzep8D9gHui4gpwAHAvGJD6ky9VTdhMrPyca5mZmb1KhHxBDBK0qiIuAHoLjqoTlSp9jGmy02YzKxc3AfCzMzqPS1pM+BG4EJJjwHPFxxTR+qp9tE1yt/VmVm5OFczM7N6hwMvAB8FfgU8ABxaaEQdqlLtY2yXP2rNrFzalqtJ2l7SDZLulrRI0kdS+tmSlktakJa3tysGMzMbkI8B20VEb0TMjohvAO8qOqhOlPWBcBMmMyuXdn4t0gv8R0TsTNYZ72RJO6d950bE7mm5qo0xmJlZ/30Y+JWkt+TSTiwqmE7V1xf09rkTtZmVT9tytYh4JCJuS+vPAfcA27XrfmZmNmiWA28DzpF0akrz1+j9VOnrA3ABwsxKZ0hyNUmTgT2Am1LSKZIWSpolacuhiMHMzFqXJo17M7CzpJ8CG2/I9SR9WdK9Ke+/TNL4lH6QpFsl3Zl+7p87Z6+UvkTSNyQppW8laa6k+9PPLVO60nFL0n32zF1rWjr+fknTNuS1tKpSDQA3YTKz0ml7ASKN5PEz4N8j4lng28D/AXYHHgG+2uS86ZLmS5q/YsWKdodpZmZrzAeIiL9HxPuB3wBjN/Cac4FdI2I34D7gjJT+OHBoRLwGmAb8KHfOt4EPADumZWpKnwFcHxE7AtenbchqTWrHTk/nI2kr4Czg9cDewFlD8eVVb9U1EGZWTm3N1SSNISs8XBgRPweIiEcjohoRfcD3yDLzF4mImRHRHRHdEyZMaGeYZmaWExEfqNv+ZkS8YgOveW1E9KbNecCklH57RPwlpS8CNpY0TtJEYIuImBcRAfwQOCIddzgwO63Prkv/YWTmAePTdd4KzI2IJyPiKbLCTK0w0jY9LkCYWUm1bR6IVNV8PnBPRHwtlz4xIh5Jm+8E7mpXDGZm1jpJl0TEUZLuBKJ+f6o9GAzHAxc3SH8XcFtErJK0HbAst28Za/rRbZv7HPkrsG1a3w74c4NzmqW3lZswmVlZtXMiuTcC7wPulLQgpX0COEbS7mQfTkuBD7YxBjMza91H0s9DBnKypOuAlzXYdWZEXJ6OOZNslL4L687dBfgicHB/7hkRIelFhZ2BkjSdrPkTO+ywwwZdq9LrGggzK6e2FSAi4vc0HrXDw7aamQ1DtW/1I+LhAZ5/4Lr2SzqOrHByQGqWVEufBFwGHBsRD6Tk5aRmTsmklAbwaK02OzVReix3zvYNzlkO7FeX/psmr2EmMBOgu7t7gwomvR6FycxKyrmamZkBIOk5Sc82WJ6T9OwGXnsqcBpwWES8kEsfD1wJzIiIP9TSU2HmWUn7pCaxxwKXp91zyDpck37m049NozHtAzyTrnMNcLCkLVPn6YNTWlv19NaaMPmj1szKpZ1NmMzMrINExOZtvPx5wDhgbhqNdV5EnAicArwS+JSkT6VjD46Ix4CTgB+QDSF7dVoAzgEukXQC8DBwVEq/Cng7sAR4AXh/el1PSvoscEs67jMR8WSbXudqldWdqN0HwszKxQUIMzNrSNJLgY1q22luiAGJiFc2Sf8c8Lkm++YDuzZIfwI4oEF6ACc3udYsYFY/Qt5gFY/CZGYl5VzNzMzWIukwSfcDDwG/JRvw4up1nmQvsmYUJn/Umlm5OFczM7N6nwX2Ae6LiClk3/bPKzakzlOrgRjb5SZMZlYuLkCYmVm9SmomNErSqIi4AeguOqhOUytAdI3yR62ZlYv7QJiZWb2nJW0G3AhcKOkx4PmCY+o47gNhZmXlXM3MzOodDvwN+CjwK+AB4NBCI+pAtT4QbsJkZmXjGggzM1tLRORrG2YXFkiHcw2EmZWVCxBmZgaApN9HxL6SngPyszCLbJTULQoKrSOt7gPhAoSZlYwLEGZmBkBE7Jt+tnNCuRGjZ/Uwrm7CZGbl4q9FzMxsLZJ+1EqarVtvbRhX10CYWck4VzMzs3q75DckdQF7FRRLx3IfCDMrK+dqZmYGgKQzUv+H3SQ9m5bngEeBywsOr+PURmHqchMmMyuZlgoQkl7T7kDMzKxYEfGF1P/hyxGxRVo2j4itI+KMouPrNKtrIDyRnJmVTKudqL8laRzwA+DCiHimfSGZmVkRJO0UEfcCP5W0Z/3+iLitgLA6VqXaR9coMWqUayDMrFxaKkBExJsk7QgcD9wq6WbggoiY29bozMxsKH0MmA58tcG+APYf2nA6W6Ua7v9gZqXU8jCuEXG/pE8C84FvAHtIEvCJiPh5uwI0M7OhERHT08+3FB1LGfT09rn/g5mVUksFCEm7Ae8H3gHMBQ6NiNsk/QPwJ8AFCDOzEpH0f4HJ5D4nIuKHhQXUgXr7+jyEq5mVUqs1EP8NfJ+stuFvtcSI+EuqlTAzs5JIcz78H2ABUE3JAbgA0Q+VXjdhMrNyarUA8Q7gbxFRBZA0CtgoIl6ICE8uZGZWLt3AzhERRQfSySrVPsZ0uQmTmZVPq1+NXAdsnNveJKWZmVn53AW8rOggOl1Ptc9DuJpZKbVaA7FRRKysbUTESkmbtCkmMzMr1jbA3WnEvVW1xIg4rLiQOk+vR2Eys5JqtQDxvKQ9a2OAS9oL+Nu6TpC0PVl72W3J2s7OjIivS9oKuJisc95S4KiIeGpg4ZuZWRucXXQAZeAmTGZWVq0WIP6dbGKhvwAiq9p+z3rO6QX+I43WtDnZ/BFzgeOA6yPiHEkzgBnA6QMJ3szMBl9E/LboGMqgp9rnGggzK6VWJ5K7RdJOwKtT0uKIqKznnEeAR9L6c5LuAbYDDgf2S4fNBn6DCxBmZoWT9PuI2FfSc2Q1x6t3ARERWxQUWkequA+EmZVUyxPJAa9jzZjge0pqeUxwSZOBPYCbgG1T4QLgr2RNnMzMrGARsW/6uXnRsZRBbzUYN8YFCDMrn1YnkhvwmOCSNgN+Bvx7RDybTV6dLhARkhoOEyhpOjAdYIcddmglTDMzs2GjUu1js4368z2dmVlnaDVnG9CY4JLGkBUeLoyI2mzVj0qaGBGPSJoIPNbo3IiYCcwE6O7u9ljkZmbWUXo8CpOZlVSrOVu/xwRXVtVwPnBPRHwtt2sOMC2tTwMu7891zczMOkGl2seY0R6FyczKp9UaiIGMCf5G4H3AnZIWpLRPAOcAl0g6AXgYOKq/QZuZWXtJejmwY0RcJ2ljoCsinis6rk7S61GYzKykWi1AnN3fC0fE78lG7mjkgP5ez8zMhoakD5D1QduKrP/bJOA7OO/ul4qbMJlZSbU6jOtv676N2gQY3d7QzMysICcDe5ONnEdE3C/ppcWG1Hk8D4SZlVVLOVv6NupS4LspaTvgF22KyczMirUqInpqG5K6WHteCGuB+0CYWVm1+tXIyWR9Gp6F7NsowN9GmZmV028lfQLYWNJBwE+BXxYcU8fpdRMmMyupVnM2fxtlZjZyzABWAHcCHwSuAj5ZaEQdyE2YzKysWu1EXf9t1En42ygzs1KKiD7ge2mxAYgIKtU+xroJk5mVUKsFiBnACaz9bdT32xWUmZkNPUl3so7a5YjYbQjD6WjVviACulwDYWYl1OooTP42ysys/A4pOoCy6O3LymFuwmRmZdRSAULSQzT4VioiXjHoEZmZWSEi4uGiYyiLnmofgEdhMrNSarUJU3dufSPg3WQTDJmZWclIeo4Xf2n0DDAf+I+IeHDoo+osld6sADG2yzUQZlY+LeVsEfFEblkeEf8FvKO9oZmZWUH+CziVbM6fScDHgR8DFwGzBnJBSV+WdK+khZIukzS+bv8OklZK+ngubaqkxZKWSJqRS58i6aaUfrGksSl9XNpekvZPzp1zRkpfLOmtA3kN/VGpZuWvrlEuQJhZ+bQ6kdyeuaVb0om0XnthZmad5bCI+G5EPBcRz0bETOCtEXExsOUArzkX2DV1xL4POKNu/9eAq2sbkkYD3wTeBuwMHCNp57T7i8C5EfFK4CmyQT5IP59K6eem40jnHQ3sAkwFvpWu3zYVN2EysxJr9auRr+aWLwB7AUe1KygzMyvUC5KOkjQqLUcBf0/7BjQHUERcGxG9aXMeWc0GAJKOAB4CFuVO2RtYEhEPpnmILgIOlyRgf+DSdNxs4Ii0fnjaJu0/IB1/OHBRRKyKiIeAJen6bVMrQLgJk5mVUaujML2l3YGYmdmw8V7g68C3yAoM84B/kbQxcMogXP944GIASZsBpwMHkTWVqtkO+HNuexnwemBr4OlcYWRZOnatcyKiV9Iz6fjt0mugwTlrkTQdmA6www47DOzVsaYJk0dhMrMyanUUpo+ta39EfG1wwjEzs6KlTtKHNtn9+2bnSboOeFmDXWdGxOXpmDOBXuDCtO9ssuZIK7PKgmKl5lozAbq7uwdU2wJraiC6RhX/mszMBlt/RmF6HTAnbR8K3Azc346gzMysOJImAB8AJpP7nIiI49d1XkQcuJ7rHkc218QBEVH75/z1wJGSvgSMB/ok/R24Fdg+d/okYDnwBDBeUleqhailk35uDyyT1AW8JB2/vMm12mZ1Hwg3YTKzEmq1ADEJ2DMingOQdDZwZUT8S7sCMzOzwlwO/A64DqgOxgUlTQVOA94cES/U0iPiTbljzgZWRsR5qQCwo6QpZP/sHw38c0SEpBuAI8n6RUxL8UL2Jdc04E9p/6/T8XOAH0v6GvAPwI5kX4K1Ta0J01g3YTKzEmq1ALEt0JPb7klpZmZWPptExOmDfM3zgHHA3NRUaV5EnNjs4NSH4RTgGmA0MCsiap2sTwcukvQ54Hbg/JR+PvAjSUuAJ8kKHUTEIkmXAHeTNZ86OSIGpWDUzJpRmFyAMLPyabUA8UPgZkmXpe0jWDPShZmZlcsVkt4eEVcN1gXT0KrrO+bsuu2rgBfFkPpovGgUpYj4O9lEp42u/Xng8y2Gu8FqM1F3eRhXMyuhVkdh+rykq4FaVfP7I+L29oVlZmYF+gjwCUmrgAogICJii2LD6hy9bsJkZiXWn8ngNgGejYgLJE2QNCWNp21mZiUSEZsXHUOncxMmMyuzVodxPYtsJKZXAxcAY4D/Ad7YvtDMzGwoSdopIu6VtGej/RFx21DH1Kk8E7WZlVmrNRDvBPYAbgOIiL9I8jdUZmbl8jGySdS+2mBfkM0AbS3o6XUNhJmVV6sFiJ40FF4ASNq0jTGZmVkBImJ6+vmWomPpdL19nonazMqr1ZztEknfJZu85wNkY4N/r31hmZlZUSS9u1bLLOmTkn4uaY+i4+okbsJkZmW23gKEsgG7LwYuBX5G1g/iUxHx3+s5b5akxyTdlUs7W9JySQvS8vYNjN/MzAbf/xcRz0naFziQbH6F7xQcU0dZ3YTJM1GbWQmttwlTarp0VUS8Bpjbj2v/gGzioB/WpZ8bEV/px3XMzGxo1SZZewcwMyKuTJO2WYtqM1GPGeUChJmVT6s5222SXtefC0fEjWQzgZqZWWdZnpqtvge4StI4Wv+8MKDXTZjMrMRa/UB4PTBP0gOSFkq6U9LCAd7zlHSNWZK2bHaQpOmS5kuav2LFigHeyszMBuAo4BrgrRHxNLAVcGqhEXWYSrUPCUaPcgHCzMpnnU2YJO0QEf8LvHWQ7vdt4LNkwwF+lmyowOMbHRgRM4GZAN3d3TFI9zczs/WIiBeAn+e2HwEeKS6iztNTDcaMGkXWjdDMrFzW1wfiF8CeEfGwpJ9FxLs25GYR8WhtXdL3gCs25HpmZmbDUaXa5+ZLZlZa62vClM/9XrGhN5M0Mbf5TuCuZseamZl1qt5qn0dgMrPSWl8NRDRZXy9JPwH2A7aRtAw4C9hP0u7pWkuBD/bnmmZmZp2gpxqeRM7MSmt9BYjXSnqWrCZi47RO2o6I2KLZiRFxTIPk8wcWppmZDRVJ+wD/DfwjMBYYDTy/rjzf1lap9jHGHajNrKTWWYCIiNFDFYiZmQ0b5wFHAz8FuoFjgVcVGlGHqbgJk5mVmHM3MzN7kYhYAoyOiGpEXABMLTqmTtLrJkxmVmLrnYnazMxGnBckjQUWSPoS2RCu/m+4H3qqfS5AmFlpOXczM7N67yP7fDgFeB7YHtigYbxHGg/jamZl5hoIMzNbS0Q8nFb/Dny6yFg6VcU1EGZWYs7dzMzMBlmlGq6BMLPScgHCzMxskLkGwszKzLmbmZnZIHMBwszKzH0gzMxsLZJeBZwKvJzc50RE7F9YUB2m0usmTGZWXi5AmJlZvZ8C3wG+B1QLjqUjVfpcA2Fm5eUChJmZ1euNiG8XHUQnq1T7GOsChJmVlHM3MzOr90tJJ0maKGmr2lJ0UJ2k0ht0uQmTmZWUayDMzKzetPTz1FxaAK8oIJaO5E7UZlZmLkCYmdlaImJK0TF0OhcgzKzMXIAwMzMAJO0fEb+W9P8a7Y+Inw91TJ2qUg3GdrkAYWbl5AKEmZnVvBn4NXBog30BuADRokq1j65R7gNhZuXkAoSZmQEQEWeln+8vOpZOFhH09oWbMJlZaTl3MzMzG0SVagC4CZOZlZZzNzMzs0FUqfYBeCZqMystFyDMzMwGUa0A0TXKH7FmVk7uA2FmZqtJeilwMrBLSloEfCsiHi0uqs7SU6uBcBMmMysp525mZgaApDcCt6TNH6YF4Ka0z1rQW+sD4SZMZlZSLkCYmVnNV4EjIuKsiJiTlrOAI4CvbciFJX1Z0r2SFkq6TNL43L7dJP1J0iJJd0raKKXvlbaXSPqGJKX0rSTNlXR/+rllSlc6bkm6z565e0xLx98vaRpttKYPhD9izayc2pa7SZol6TFJd+XSGmb6ZmY2LGwREbfXJ0bEAmDzDbz2XGDXiNgNuA84A0BSF/A/wIkRsQuwH1BJ53wb+ACwY1qmpvQZwPURsSNwfdoGeFvu2OnpfCRtBZwFvB7YGzirnZ8/q/tAuABhZiXVztztB6zJ7GuaZfpmZlY8NfrHOv0DvkGfFxFxbUT0ps15wKS0fjCwMCLuSMc9ERFVSRPJCjTzIiLImlMdkc45HJid1mfXpf8wMvOA8ek6bwXmRsSTEfEUWWGm/vNp0PT0ugmTmZVb2woQEXEj8GRdcrNM38zMincucK2kN0vaPC37AVenfYPl+HRNgFcBIekaSbdJOi2lbwcsy52zLKUBbBsRj6T1vwLb5s75c4NzmqW/iKTpkuZLmr9ixYr+vzKgt89NmMys3IZ6FKZmmf6LSJpOVgXNDjvsMAShmZmNbBExU9JfgM+y9ihMn4uIX67vfEnXAS9rsOvMiLg8HXMm0AtcmPZ1AfsCrwNeAK6XdCvwTIsxh6Ro5dgWrzcTmAnQ3d09oOu6D4SZlV1hw7iuL9MfjEzczMz6JyKuAK4Y4LkHrmu/pOOAQ4ADUrMkyGoDboyIx9MxVwF7kvWLmJQ7fRKwPK0/KmliRDySmig9ltKXA9s3OGc5Wd+KfPpv+vPa+qPWhKnLTZjMrKSG+uuRR1NmT12mb2ZmBZO0jaSzJH1Y0maSviXpLkmXS3rlBl57KnAacFhEvJDbdQ3wGkmbpA7VbwbuTrXVz0raJ42+dCxweTpnDlAbSWlaXfqxaTSmfYBn0nWuAQ6WtGXq43FwSmuLWg3EWNdAmFlJDXXu1izTNzOz4v0YGEfWL+FmYClwJFmNxPc38NrnkY3kNFfSAknfAUidmr9GNv/EAuC2iLgynXNSuu8S4AHW9Js4BzhI0v3AgWkb4CrgwXT899L5RMSTZM2ybknLZ1JaW7gPhJmVXduaMEn6CVmV8TaSlpENoXcOcImkE4CHgaPadX8zM+u3bSPiE+kb/4cj4ksp/V5JJ2/IhSOiaQ1GRPwPWZOl+vT5wK4N0p8ADmiQHmSzaDe6xyxgVj9CHrBaEyYXIMysrNpWgIiIY5rselGmb2Zmw0IVVvdRe7xuX18B8XSkNZ2o3QfCzMqpsE7U1vkmz7hy/Qe1wdJz3lHIfc1GgFdImgMot07anlJcWJ3FozCZWdm5AGFmZjWH59a/Urevftua6K2mJkxdLkCYWTm5AGFmZgBExG+b7ZN0MdB0v63R4yZMZlZy/nrEzMxa8YaiA+gUq5swjfJHrJmVk3M3MzOzQbS6AOEmTGZWUm7CZGZmAEjas9kuYMxQxtLJKrU+EG7CZGYl5QKEmZnVfHUd++4dsig6nJswmVnZuQBhZmYARMRbio6hDCrVPkaPEqNGuQbCzMrJX4+YmRkAkk7Lrb+7bt9/Dn1EnalSDTdfMrNScwHCzMxqjs6tn1G3b+pQBtLJKtU+TyJnZqXmHM7MzGrUZL3RtjVRqfYx1gUIMysx53BmZlYTTdYbbVsTld6gy02YzKzE3InazMxqXivpWbLaho3TOml7o+LC6ixuwmRmZecChJmZARARo4uOoQwqfeEmTGZWas7hzMzMBlGl1zUQZlZuzuHMzMwGUaXa5z4QZlZqLkCYmZkNoh73gTCzknMOZ2ZmNoh6q+4DYWbl5hzOzMxsEFWqfYzpchMmMysvFyDMzMwGUaXaR9cof7yaWXk5hzMzMxtEPdVwHwgzKzXncGZmZoOot9rHWDdhMrMScwHCzMxsEHkmajMru0Jmopa0FHgOqAK9EdFdRBxmZmaDrVIN94Ews1IrpACRvCUiHi/w/mZmZoOux02YzKzk/BWJmZnZIOp1EyYzK7micrgArpV0q6TpjQ6QNF3SfEnzV6xYMcThmZmZDUzFozCZWckVlcPtGxF7Am8DTpb0T/UHRMTMiOiOiO4JEyYMfYRmZmYD0FPto2u0mzCZWXkVUoCIiOXp52PAZcDeRcRhZmY22CrVPsa6BsLMSmzIczhJm0ravLYOHAzcNdRxmJmZDbZqXxCBmzCZWakVMQrTtsBlkmr3/3FE/KqAOMzMzAZVpdoHuABhZuU25AWIiHgQeO1Q39fMzKzdelYXINwHwszKq8h5IMwGZPKMK4f8nkvPeceQ39PMOk9vNQDXQJhZuTmHMzMzGyRuwmRmI4FzODMzs0HS0+smTGZWfi5AmJmZDRLXQJjZSOAczszM2k7SlyXdK2mhpMskjU/pYyTNlnSnpHsknZE7Z6qkxZKWSJqRS58i6aaUfrGksSl9XNpekvZPzp1zRkpfLOmt7XqdvX3uA2Fm5ecczszMhsJcYNeI2A24D6gVFN4NjIuI1wB7AR+UNFnSaOCbwNuAnYFjJO2czvkicG5EvBJ4CjghpZ8APJXSz03Hkc47GtgFmAp8K11/0LkJk5mNBC5AmJlZ20XEtRHRmzbnAZNqu4BNJXUBGwM9wLPA3sCSiHgwInqAi4DDlU0itD9waTp/NnBEWj88bZP2H5COPxy4KCJWRcRDwJJ0/UG3uglTlz9ezay8nMOZmdlQOx64Oq1fCjwPPAL8L/CViHgS2A74c+6cZSlta+DpXGGklk7+nLT/mXR8s2u9iKTpkuZLmr9ixYp+v7BKbRjXUf54NbPy8jwQZmY2KCRdB7yswa4zI+LydMyZQC9wYdq3N1AF/gHYEvhduk4hImImMBOgu7s7+nt+ryeSM7MRwAUIMzMbFBFx4Lr2SzoOOAQ4ICJq/5z/M/CriKgAj0n6A9BNVmOwfe70ScBy4AlgvKSuVMtQSyf93B5YlppEvSQdv7zJtQZdj5swmdkI4BzOzMzaTtJU4DTgsIh4Ibfrf8n6NCBpU2Af4F7gFmDHNOLSWLJO0HNSweMG4Mh0/jTg8rQ+J22T9v86HT8HODqN0jQF2BG4uR2vs9aEaaxHYTKzEnMNhJmZDYXzgHHA3KxfM/Mi4kSykZYukLQIEHBBRCwEkHQKcA0wGpgVEYvStU4HLpL0OeB24PyUfj7wI0lLgCfJCh1ExCJJlwB3kzWfOjkiqu14kbVO1F1uwmRmJeYChJmZtV0aWrVR+kqyoVwb7bsKuKpB+oM0GEUpIv6+jmt9Hvh8P0IeEE8kZ2YjQekLEJNnXDnk91x6zjuG/J5mZlY8N2Eys5HAOZyZmdkgcQ2EmY0EzuHMzMwGySZjRzNlm00Z61GYzKzESt+EyczMbKgcvvt2HL57wznqzMxKw1+RmJmZmZlZy1yAMDMzMzOzlrkAYWZmZmZmLXMBwszMzMzMWuYChJmZmZmZtcwFCDMzMzMza1khBQhJUyUtlrRE0owiYjAzMzMzs/4b8gKEpNHAN4G3ATsDx0jaeajjMDMzMzOz/iuiBmJvYElEPBgRPcBFwOEFxGFmZmZmZv2kiBjaG0pHAlMj4l/T9vuA10fEKXXHTQemp81XA4uHNNCB2QZ4vOgghhE/jzX8LNbm57FGq8/i5RExod3B2BqSVgAPD+DUTn9/O/5iOf5iOf61Nfzs6RrEGwyqiJgJzCw6jv6QND8iuouOY7jw81jDz2Jtfh5r+FkMXwMtsHX679TxF8vxF8vxt6aIJkzLge1z25NSmpmZmZmZDXNFFCBuAXaUNEXSWOBoYE4BcZiZmZmZWT8NeROmiOiVdApwDTAamBURi4Y6jjbpqCZXQ8DPYw0/i7X5eazhZ1E+nf47dfzFcvzFcvwtGPJO1GZmZmZm1rk8E7WZmZmZmbXMBQgzMzMzM2uZCxB1JH1Z0r2SFkq6TNL4lL63pAVpuUPSO3PnTJW0WNISSTNy6VMk3ZTSL06dxpE0Lm0vSfsn5845I6UvlvTWoXvlja3jeRwk6VZJd6af++fO2SulL5H0DUlK6VtJmivp/vRzy5SudNySdJ89c9ealo6/X9K0IX75a1nHs9ha0g2SVko6r+6cUj6LFE/D55H2NXwfl/xv5d2SFknqk9SdSx8r6YL0PrhD0n65faV9f4xUzd7jw5WkWZIek3RXLq3h+284krR9yn/vTn9/H0npHfEaJG0k6eaUNyyS9OmU3jBPHI4kjZZ0u6Qr0nbHxA4gaWnKhxdImp/SOuL9AyBpvKRL0+fxPZLeMCTxR4SX3AIcDHSl9S8CX0zrm+TSJwKPkXVCHw08ALwCGAvcAeycjrsEODqtfwf4UFo/CfhOWj8auDit75zOHwdMSdcdPUyfxx7AP6T1XYHluXNuBvYBBFwNvC2lfwmYkdZn5K719nSc0nk3pfStgAfTzy3T+pbD8FlsCuwLnAicV3dOKZ/Fep5Hw/fxCPhb+UeySS9/A3Tn0k8GLkjrLwVuBUaV/f0xEpd1vceH6wL8E7AncFcureH7bzguZJ/He6b1zYH7Uv7QEa8h/S1vltbHADelv+2GeeJwXICPAT8GrkjbHRN7inEpsE1dWke8f1J8s4F/TetjgfFDEX/hL3w4L8A7gQsbpE8BHiUrQLwBuCa374y0iGwmwNo/WKuPIxuB6g1pvSsdp9q5uWutPm44LOt4HgKeJPtnbiJwb27fMcB30/piYGJanwgsTuvfBY7JnbM47V99bqPjhtuzAI4jV4AYKc+i/nk0ex+PlL8VXlyA+Cbwvtz29cDeI+n9MVKWZu/xouNqIe7JrF2AaPj+64QFuBw4qBNfA9kXlbcBr2+WJw63hWwur+uB/YEr1pWfD9eFxgWIjnj/AC8BHiINijSU8bsJ07odT/ZtHwCSXi9pEXAncGJE9ALbAX/OnbMspW0NPJ2OyaeTPyftfyYd3+xaw8VazyPnXcBtEbGKLN5luX3517BtRDyS1v8KbJvWm73u4fw8mj2LvJHyLGDt59Hf11DGv5W8O4DDJHVJmgLsRTaZ5kh6f4wUZfk9NHv/DWupieMeZN/id8xrSE2AFpC1bJhLVovVLE8cbv4LOA3oS9vrys+HqwCuVdYce3pK65T3zxRgBXBBakb2fUmbMgTxD/k8EMOBpOuAlzXYdWZEXJ6OORPoBS6s7YyIm4BdJP0jMFvS+v6B7AgDfR4pfRey5isH9+eeERGSht0YwhvyLAZquD4LKOZ5DGetPI8GZpE1b5oPPAz8Eai2es/h/P6w8uuU95+kzYCfAf8eEc+m7kTA8H8NEVEFdlfWj+wyYKdiI2qNpEOAxyLi1nzfrg60b0Qsl/RSYK6ke/M7h/n7p4usCeKHI+ImSV8na7K0WrviH5EFiIg4cF37JR0HHAIcEKn+p+78eyStJLX9J/s2sWZSSnsCGC+pK5XEa+nkzlkmqYusCuqJdVyrrQb6PCRNIsvsjo2IB1LycrK4a/Kv4VFJEyPiEUm1fiS1cxq97uXAfnXpv2n5hQ3Ahr436nT0s4ABP491vY9L/bfS5Jxe4KO1bUl/JGun/RQd/v6wFynkfdkGzd5/w5KkMWSFhwsj4ucpuaNeA0BEPC3pBrJmP83yxOHkjWS1q28HNgK2AL5OZ8S+WkQsTz8fk3QZWRPTTnn/LAOWpS+4AS4lK0C0PX43YaojaSpZddxhEfFCLn1K+gcGSS8n+4ZgKXALsGPaP5aso+ec9M/UDcCR6RLTyNpmAsxJ26T9v07HzwGOVjbyzBRgR7JOloVZx/MYD1xJ1knnD7X0VGX2rKR9lH0FdCyNX3f98zhWmX2AZ9J1rgEOlrRlGkHg4JRWiGbPopkyPwtY5/No9j4u9d9KM5I2SVXKSDoI6I2Iu8v+/hihGr7HC45pIJq9/4ad9LdzPnBPRHwtt6sjXoOkCVozot/GZP037qF5njhsRMQZETEpIiaTvdd/HRHvpQNir5G0qaTNa+tkeedddMj7JyL+CvxZ0qtT0gHA3QxF/IPdqaLTF2AJWRvWBWmpjQDzPmBRSrsNOCJ3ztvJvlF8gKwpQy39FWT/1CwBfgqMS+kbpe0laf8rcuecma6zmDQiyzB9Hp8Ens+lLwBemvZ1k/0BPgCcx5oZz7cm62x1P3AdsFVKF1lH0wfI+pfkO6Aen2JYArx/OD6LtG8pWUfylWTfCNRGFyrls2jheTR8H5f8b+Wd6Xe/imyQhVpH8MkpxnvS7/rluXNK+/4YqUuz9/hwXYCfAI8AlfT+PaHZ+284LmQj4AWwMJcXvb1TXgOwG3B7iv8u4FMpvWGeOFwXshrQ2ihMHRN7ivWOtCyq/c12yvsnxbo7WRPZhcAvyEbia3v8tQ8rMzMzMzOz9XITJjMzMzMza5kLEGZmZmZm1jIXIMzMzMzMrGUuQJiZmZmZWctcgDAzMzMzs5a5ADGCSapKWpBbZqT030jqrjt2P0nPpOMWSrouzdpY2z9d0r1puVnSvk3u+RlJ/Z6Ma6BS3Ffk1v9vbt+Jko5dz/nHSTqvH/frlvSNFo77Y/o5WdI/9+f8utd0WO331h+Sdk+T/9S2B3QdM7ORTtnEsi/Kzwfp2p+o2/7jYF7fbKBG5EzUttrfImL3fhz/u4g4BEDSF4CTgbOUTWf/QbLp4B+XtCfwC0l7RzbJyWoR8alBin0g9iObp+GPKZbvDPYNImI+2XjM6zuuVpCZDPwz8OP+nJ+7zhwGNlHV7mRzEFy1gdcxM7PMZHL5eStyMzY38wngP2sbuc8Os0K5BsL6Lc38uTnwVEo6HTg1Ih4HiIjbgNlkBYz6c38g6ci0vlTSpyXdJulOSTs1OP44Sb+QNDcdf4qkj0m6XdI8SVul41bXmkjaRtLSuutMBk4EPppqUd4k6WxJH8+d//W07y5JezeIZYKkn0m6JS1vbHBMvnbgbEmz0rUflPRvueNWptVzgDel+3607vy9Jf0pvdY/5maarH8+56X1fG3S3yS9udE1lM2Q+xngPenY99RdZ7KkX6eapusl7ZD73X0jXefB3O9xoqQbc8/uTfVxmpmNAPX5+WhJX06fFwslfRBWf078TtIcslmDSZ9zt0paJGl6SjsH2Dhd78KUVqvtULr2Xenz8z25a/9G0qXKWgRcmD6zkXSOpLtTLF8Z8qdjpeIaiJFtY0kLcttfiIiL13H8m9LxW5PNQl2rWt0FuLXu2PmsmUZ9XR6PiD0lnQR8HPjXBsfsCuxBNivxEuD0iNhD0rnAscB/re8mEbFU0neAlRHxFQBJB9QdtklE7C7pn4BZ6b55XwfOjYjfp3+qrwH+cT233gl4C1mBa7Gkb0dEJbd/BvDxXM3Ofrl99wJvioheZc2+/hN41zpe4+7pGocCp5HVtGxcf42IeJekT5HNYnxKOue43KX+G5gdEbMlHQ98Azgi7ZtINvPrTmQ1FpeSfeN2TUR8XtJoYJP1PBMzszKqz8+nA89ExOskjQP+IOnadOyewK4R8VDaPj4inpS0MXCLpJ9FxAxJpzRpKfD/yGqSXwtsk865Me3bg+xz+S/AH4A3SroHeCewU0SEpPGD+9JtpHEBYmTbkCZMpwNfIvtWf0P8PP28lSxDbOSGiHgOeE7SM8AvU/qdwG4beP+8nwBExI2StmiQwR4I7Jy+zAHYQtJmEbGS5q6MiFXAKkmPAdsCy1qM5yXAbEk7AgGMWd8J6dgvA2+JiIqkl/X3GsAbWPO7+BHZ77nmFxHRB9wtaduUdgswS9KYtH9BC/cwMyu7g4HdarW1ZHn6jkAPcHOu8ADwb5Lemda3T8c9sY5r7wv8JCKqwKOSfgu8Dng2XXsZZDXTZE2r5gF/B85PtdxXbPjLs5HMTZhsoOYA/5TW7wb2qtu/F7CoheusSj+rNC/Qrsqt9+W2+3Ln9LLm/bxRC/dtJNazPQrYJyJ2T8t26yk8wNqxr+s1NvJZssLTrsChrOd1SdoMuAT4QEQ8MpBrtCD/egRZgYvsvbAc+IHW0zHdzGyEEPDh3GfGlIio1UA8v/qgrOb5QOANEfFa4HY2LK9+0edO6mexN1mt8SHArzbg+mYuQNiA7Qs8kNa/BHxR0taQjfADHAd8awjjWcqaQsyRTY55jqwpUTO1NqT7klU7P1O3/1rgw7WN9Do31LpiegnZP+WQPc/1mQVcEBG/a+Ea67rvH4Gj0/p7gd81OQ4ASS8HHo2I7wHfJ6uaNzMbaerz1WuAD6XaWSS9StKmDc57CfBURLygrC/gPrl9ldr5dX5H1o9ttKQJZF/i3NwssPQF00si4irgo2RNn8wGzE2YRrb6PhC/iojaUJ5XSqq11f8T8E3W9IEQ8Aypv0JEzJG0HfBHSUGWif5L7lvwofAV4JLU5vTKJsf8ErhU0uHkCgI5f5d0O1kzn+Mb7P834JuSFpL97dzIhjfhWghUJd0B/IDsm6eaL5E1P/okzV8TsPqf+COBV6V+C5D9fppd4wZgRvp9fqHuch8GLpB0KrACeP96XsN+wKnp/bKSrF+KmdlIU5+ff52s+dBtqSPzCtb0J8v7FXBi6qewmKy5Uc1MYKGk2yLivbn0y8iam95BVlt+WkT8VQ0GI0k2By6XtBHZZ/jHBvQKzRJF1LfSMBt5JP2GrPNby0OompmZmY1EbsJkZmZmZmYtcw2EmZmZmZm1zDUQZmZmZmbWMhcgzMzMzMysZS5AmJmZmZlZy1yAMDMzMzOzlrkAYWZmZmZmLfv/ATC0aBCtmnZ/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 792x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(11, 4))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.hist(_model.ELBO_inits)\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.xlabel(\"ELBO in multiple initializations\")\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(_model.ELBO_iters)\n",
    "plt.xlabel(\"Iterations\")\n",
    "plt.ylabel(\"ELBO in a single initialization\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### Visualize assignment probability and allele frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In mitochondrial, allele frequency is highly informative between 0.01 to 0.1, \n",
    "# so we rescale the colour to give more spectrum for this region.\n",
    "# You can design/choose your own colors from here: \n",
    "# https://matplotlib.org/stable/tutorials/colors/colormaps.html\n",
    "\n",
    "from matplotlib import cm\n",
    "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n",
    "\n",
    "raw_col = cm.get_cmap('pink_r', 200)\n",
    "new_col = np.vstack((raw_col(np.linspace(0, 0.7, 10)),\n",
    "                     raw_col(np.linspace(0.7, 1, 90))))\n",
    "segpink = ListedColormap(new_col, name='segpink')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArIAAAGGCAYAAACHemKmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABO3UlEQVR4nO3dd7wcdb3/8df7HJIAIQGlJCqiKIiKepEggjUaC1e9XiwIglKsYFADKhoEKXrF9guWgAWV4rWgonDBSzEiFsRgcqUoorRITQCBFEiB5PP74zsLk82es7N7ds/M7nk/ecwjmdnvzHx2w/mez37nWxQRmJmZmZn1moGyAzAzMzMza4cTWTMzMzPrSU5kzczMzKwnOZE1MzMzs57kRNbMzMzMepITWTMzMzPrSU5kzczMzKwnOZE1MzMzs57kRNbMzMzMepIT2R4g6XhJXoKtQiQ9VVJI+mgHr3l8ds2tCpRdJOmM3P707NzpuWNnSFrUqfjMrDokXSbpstx+rU46uI1r9Wz9MZL3bf3BiewISPpA9gM0v+xY+pWkoyXtXXYc/UDSplmyPL3sWMw6SdLBWV0ckl7S4HVJui17/YIyYrSRkbS/pFllx2HV40R2ZA4AFgG7S9qhi/f5DLBJF69fZUcDe5cdRAXtBLy3SZn3ZuVqNgWOA6Z3KSazsq0C9m9w/OXAtsDq0Q2np9XXH2XbH5jV4Pg/Sb8fvzeq0VhlOJFtk6TtgRcBRwL3kJLaroiIRyJiVbeub4mkTcuOoaiIWB0RDzcp83BE+Be3jSX/C+wjaaO64/sDC4HFox9Sb+p2/dGp+jaSVRGxthPXs97jRLZ9BwD3A78AfsoQiayk/SQtlLRc0jJJ10r6cO71cZKOk3SDpFWS/iXp95JenSuzQR9ZSZtI+qqke7Nr/4+kJ2WPzo6vP1fSDlmfpwckLZV0en1FkpWbK2kfSddJWinpCknPzV5/v6Qbszgvk/TUBu/3hZIuyu7xkKTfSHpxXZlCMWXveSJwUO6x4RlD/YPk+nntK+mzkhZLejD7bJ5cV/YySX+RNE3SbyU9BHw2e20bSd+RtCR7r1dLOmiY+x4h6Z/Z5/UbSc+pe/152fu8ObveYknflbTlEJfcStKPs/9f/iXpK5I2rrvmen1kh4jr0T5u2b/VPdlLx+U+z+MlHZL9/fkNrnG0pLWSnjTcvcwq4ofAlkC+/hwPvBX4QaMTJA1ImiXpr9nP5xJJ35T0uLpy/ynpF5LulLRa0k2SjpU0WFeuVrc8W9Kvs3rwDklHFXkD2c/jpZLuzu5znaTDWv0gctd7pqSfSrove38LJL2xwHkb9JHNPqsPK/0eWyXpnqy+363JtYarb5t+rkp9gV8PPCVXdy3KXmvYR1bSKyX9Lvsd8ICk8yQ9q9CHZj2l/lurFXcA8LOIWCPph8Bhkl4QEX+qFVBKRn8I/Ar4eHb4WcCLga9k+8cDs4FvA1cCk4HdgF2BXw5z/zOAt5Eep/yR9OjsF8OU/zFwS3avXYH3AHfn4qp5KfBG4JRsfzZwgaQvAB8ATgUeBxwFfBd4Ze79vhK4kNTycQKwDjgEuFTSSyPiyhZjeiePfS7fyo7dNMx7rPkkEMDngW1Ij6PmSdolIlbmym2Zxfsj4L+BJZI2AS4DdgDmZvHtA5whaYuI+ArrOxCYRPq8NgY+nL3f50bEkqzMq4GnAaeTWoR2Bt4H7Cxpj4ioH8j3Y1KXldnAHsCHSJ/5gQXe+1DuAQ4Dvg78HPhZdvya7D2eQvp/+s915x0AXBYRd4zg3majZRFwBfB20s82wL8Dm5N+zj/U4JxvAgeTfj6/CmwPHA48X9KLc08+DgZWAHOyP18JnEiqsz9Wd83HAReRfs5+TEqkPy/p2oi4kOEdBvwV+B/gEeA/gFMlDUTEKcOeWUfSzsDlwB3A54AHSb83zpX0loj4eSvXA75D+hwuJNXNG5F+Z+wBLGhy7gb1bXb8YJp/rv9F+jfcFjgiO7ZiqBtJelV2r5tJv2M3AT4IXC5p14hY1PSdWu+ICG8tbsA0UqL0qmxfwG3Al+vKfRlYCgwOc62rgAua3O/49E/16P6u2f1Prit3enb8+Ppzge/Ulf0ZcG/dsSD1MXtq7tj7suN3AZNyxz+bHX9q7jP4B6nyVq7cJqTK5JI2Y1oBnFHw32V6dt3b62LdJzv+odyxy7Jj76+7xoez4wfkjo0D/gAsr10XeGpW7iHgSbmyu2fH5+Q/gwax7peVe2mDz+W8urKnZMeflzu2KP+55N779NyxM4BFuf2t6v//yL32A9Ivu4Hcsedn5Q8u+2fOm7fhNlIyFKRGgJnAstrPHSmRvDT7+yJy9S3wkuy8/euu99r640P8HH+DlBxOyB2r1S3vzB0bT6pDf1rgvTS6z0XATXXHLiN9yazt1+qkg3PH5pG+rObjEym5/UfuWJH64xVZma80iE9N3lPtM3l/g9eKfq4X5ONp8r7/TEqUH5879jxgLXBm2f+/euvs5q4F7TmA9EPya8gyTDgb2K/uMdMDpEfjr66/QF2ZnSXt2ML998r+PLXu+NeGOecbdfu/A7aUNLnu+K9i/W+rtRkZzomI5Q2OPy37cxdgR1JCtKWkrZSmkZpIapF+maT6/9+KxtSqs+pi/Snpl8jr6sqtJiX/ea8jtZr+sHYgUovMV4HNSC3feedGrrUyUqvz/Py9ItcKLGnj7HP5Y3Zo1wbx17e61P5d6+PvpLOAJ5J+WdUcAKwEzunifc067cekL9BvkDQJeANDdCsgfcldCvyyVmdlP58LSV+iH/15qPs5npSV+x1pEOUz6667gtTqWDt3DenJ0tNoou4+m2f3+Q3wNEmbNzs/d+7jSa2bPwYm5d7blsDFwI4tdhl6CylhPKFBzEWmh2xU37b6uTYl6Qmk30dnRMR9uftcQ3rK2c161ErgRLZFWaK6HymJ3V6pn+cOpORlCjAjV/xUUivlhZJuV+oXuVfdJT8FbAH8I+t39EVJz2sSxlNIj+1vqTt+4zDn3Fq3f3/25+OalFua/XnbEMdr59cS8TNJj7Hz23uACaRHQ+3E1Kob8jtZJXsj6Zt73h3ZL5i8pwA3RMS6uuN/y70+5L0y/8jfS9Ljlfq5LiElhvfw2L9do19M9de8ifTvXR9/J/2SlOwfAKkvHOnx7Hl1XwrMKi0i7iG1RO4PvBkYJH2ZbWRH0s/g3WxYb21G6poEpMf0kn4uaSmpxfceHktW63+Ob2+Q3N1PgbpN0oslzZP0IKmh4x6y/qQN7jOcHUitr59mw/dWS0a3aXxqQ08H7swnhy1qVN+2+rkWUauj/97gtb+RxiBMbOO6VlHuI9u6VwJPICWz+zV4/QDgEoCIuFvSLqTHVP+ebYdIOisiDsrK/FbS04H/BF5DSvqOkHRoRHy7g3EPNaJTBcs1O7/2pehjpO4SjdT3aSoaU7esbF6kI35MmuHii6TPZgXp87qIYl8mi7R2jEhErJX0A+C9kj5A6sf9RHKtSmY95AfAacBU4MKIeGCIcgOkJHaoWWfuAZC0BalVdBmp8eEmUjesXUl98et/jtuq27LfBb8CrifNiHMbsIbUinhEg/sMp1b2S6QW2EaGa/zotA3q2zY+V7MNOJFt3QGkim9mg9feDLwpS0JXwqOPlM4Hzs9auU4F3i/p0xFxY1bmPtIjl9MlbQb8ltRfcqhE9p+kH/DtWb/1rptz2TZTG4S1LCLmdfC67SRx63XTkCTSZ3NNgXP/CTwvG1iRb5V9Zu71Ie+VeQapLx5KI59nAMdFxIm5mIbrSrIj67e270D6915UIP7hNPsszwI+Qhpc8u+kX+JD/QI0q7KfkwZx7QHsO0y5m4BXAZfH+gNB600nPZJ/c0T8tnZQaRrGTvoP0tOrN0bEo0+sJL1i6FOGdHP258MdqpNvAl4r6fEjaJWtN53in2vR3wW1OrrRHLjPJI3DeLCVIK3a/G2nBdmI9jeTBgv8tH4jjXKfRBr1j+qmV8oSo1oyNWGIMitI35InDBNKLbn4QN3xD7b+rjpmIami+2iWjK9H0tZtXvdBUteLVhyY9Y2reSupFb3ZaGFI81BOJffLT2lOyg+SWlJ/U1d+73w/M0m7Ay/M3avWMlPfEjNrmBjqvyTV/l2LxD+ch7I/t2j0YtaH7BrSU4G3AD+KiEdGeE+zUZfVo4eRGgTOH6boj0ldD46tf0HSRlmLITT4OVaa1qu+Dh6pRvfZnDT7S0si4m7SIKv3Z/1G19NGnXxOFtdxDa7V7lO0Vj7XBynQ1SAi7iI9+Too9++H0rSIryHV8dZH3CLbmjeSEtX/GeL1P/LY4ghnA9/OOtxfShpJ/xRSUnIVj/W5vE5pjryFwH2kUbdvJSXFDUXEQknnALOyRLg2/dYzakXae3vti4h1kt5DSrb+Kul00ij4J5EGTCwjtTa0aiHwKklHAncCt0REsyWB7wN+n8UwhZQ03kh61NjMt4D3k6bbmkZqBX0r6VH7rAb9RW/M7vV10pePWcC/gC8ARMQySb8FjpI0jvSZvIbUmj6U7SX9D6nrwZ7AO4AfRMTVBeIfUkSslHQdsK+kf5A+p79ExF9yxc4iPYoEdyuwHhYRZxYo8xtJ3wRmZ93ALgEeJj0V2Yc0i8lPSbOW3A+cKemrZLMS0PluUJeQuhKcn8W1GWmFrbtJX8ZbNRP4PXCtpNNIrbRTSPXKtsC/Fb1QRPxa0veAD2VPlGpdo15KGjMy5O+sYbTyuS4k1V1zgD8BKyJiqC8pHyP9LrpC0nd4bPqtpaQvN9ZHnMi25gBS/52G87tmydwvgAOyBPO/SdNXfYDUCraYlOAen3ts/VVSgvwaUiL0T+AYUn/K4RyYXe/twJtIgxv2JXVwL2UVsIi4TNKepNaNw0mV8GLSQLhvtnnZI0nJZW2Z3jN5bMaEoXyWNNXKbNIXj18BH4iIh4Y9i0eTvemkORcPIs1l+HfgkIg4o8EpZ5EGYs0iDZy4Ejg8axWo2Z8088BMUgV9CenR/Z1DhLEvaR7Fz5HmkZzLhvNUtus9WSwnk6YEOgHIJ7LfJ/VNuyk2nPfXrO9ExKGSFpK+wH6W9DO3iFR/X56V+ZekNwD/j1QX3Z+9/is62P0mIv4u6a3ZPb5Eqj+/Tmog+W4b17tOabGC40hTlG1JSor/TKpjWnUI6anNu0m/o5aS5o/9QxvXavVzPZU0G8EhpP7C/2SI1vaImJcNrD6B9D4fJj1N+3hE1A+Sth6nYrNmWC/IWhT+DLwjIr5fcjijLktAfw3sk3X1sBZlU9/cBZwYEZ8uOx4zM7PhuI9sj8r669abRWod/G2D18yKOJjUZ/B7JcdhZmbWlLsW9K6jsj6cvyY9CqtN7/WtiKif89VsWErLCz+btLzvueElHM3MrAc4ke1dfyCtGHYsqS/qraRO7P9VYkzWuz5Fmuv2csqd/cLMzKww95E1MzMzs57kPrJmZmZm1pN6IpGVNFPSIkmrJM3PJp03MzMzszGs8l0LJO1LmqvzUNL8obNIE1XvlK1c0ux8kdaMr5/I3sxG3yTgzqh6xTPGud40qxTXm8PohUR2PvCniDg82x8AbgO+FhGfK3D+k0iraplZNWwbEXeUHYQNzfWmWeW43hxCpWctyNZcngacVDuWrZ41j7TEXqNzJpBWyFrPtP/8BIPjNu5WqH3lgm8eUXYI1oeWL1vG07d/MriVrxcsB7jtttuYPHly2bGYjVnLli3jyU92vTmcSieywFakydmX1B1fAjxziHNmk5bjW8/guI3ZyIlsIf7FZWaQ6gLXB2ZWZT0x2KtFJwGb57Ztyw3HzMzMzLqh6i2y9wJrgSl1x6cAixudEBGrgdW1/TRmwczMzMz6TaVbZCNiDbAQmFE7lg32mgFcUVZcZmZmZla+qrfIAswBzpS0ALiSNP3WROD0MoMyMzMzs3JVPpGNiLMlbQ2cCEwFrgL2ioj6AWBmZmZmNoZUPpEFiIi5wNyy4zAzMzOz6qh0H1kzMzMzs6E4kTUzMzOznuRE1szMzMx6khNZMzMzM+tJTmTNzPqQpJmSFklaJWm+pN3LjsnMrNOcyJqZ9RlJ+5Lm4D4B2BW4GrhY0jalBmZm1mFOZM3M+s+RwGkRcXpEXAccCjwEvKvcsMzMOsuJrJlZH5E0HpgGzKsdi4h12f6eQ5wzQdLk2gZMGpVgzcxGyImsmVl/2QoYBOpXP1xCWh2xkdnA0tx2e9eiMzPrICeyZmZ2ErB5btu23HDMzIrpiSVqzcyssHuBtcCUuuNTgMWNToiI1cDq2r6krgVnZtZJbpE1M+sjEbEGWAjMqB2TNJDtX1FWXGZm3eAWWTOz/jMHOFPSAuBKYBYwETi9zKDMzDrNiayZWZ+JiLMlbQ2cSBrgdRWwV0TUDwAzM+tpTmTNzPpQRMwF5pYdh5lZN7mPrJmZmZn1JCeyZmZmZtaTnMiamZmZWU9yImtmZmZmPcmJrJmZmZn1JCeyZmZmZtaTnMiamZmZWU9yImtmZmZmPcmJrJmZmZn1JCeyZmZmZtaTvEStmY0pkl4GfAyYBjwBeFNEnNvknOnAHGBn4DbgMxFxRjfjrILr//dsNtt0k7LD6AnPfuM7yg7BbEwqtUVW0ssknS/pTkkhae+61yXpREl3SVopaZ6kHUsK18z6w0TgamBmkcKStgd+Afwa2AX4MvBtSa/tUnxmZlZQ2V0Lmv1COQr4EHAo8ELgQeBiSRuPTnhm1m8i4sKIOCYifl7wlEOBWyLiIxHxt4iYC/wUOKJ7UZqZWRGldi2IiAuBCwEkrfea0oFZpEd452XHDgSWAHsDPxrFUM1sFGVfVseP8DKrI2J1B8LZE5hXd+xiUsusmVkltFFvromIVd2KZ7RUuY/s9sBUcr9AImKppPmkXyxOZM36kKSNx208aeXDq5a3ctoKYLO6YycAx3cgpKmkL9B5S4DJkjaJiJUduIeZWdskbTxho41Wrn7kkVZOWyxp+15PZqucyE7N/mz0C2QqQ5A0AZiQOzSpw3GZWXeNf3jVcqbtfTSD45r3Ilr78CoWnvvZzYBtgXz224nWWDOzXjB+9SOP8MrnPpeNBgebFn5k7VouvfbaqaQWXCeyFTMbOK7sIMxsZAbHb8JGBRJZHuuWtDwilnUhlMXAlLpjU4Blbo01syrZaHCQcQUS2X5S9mCv4SzO/mz0C2QxQzsJ2Dy3bdv50Mys+5SS1GYbanqlEboCmFF37NXZcTMzK1GVE9lbSAnro79AJE0mzV4w5C+QiFgdEctqG+s/ajSzXlEkiX00mW3lstpM0i6SdskObZ/tb5e9fpKks3KnfAN4mqQvSHqmpA8AbwNO7sC7NDOzESi1a4GkzYAdcoe2z3653BcRt0r6MnCMpBtIie2ngTuBc0c5VDMbbRpMW5FyrdmNNCdszZzszzOBg0mLJGxXezEibpH0elLi+mHgduA9EXFxqzc2M+umcePHM26jAqlda4PCKq3sPrLNfqF8gTTX7LeALYDfA3v1+gg7MyugaGtriy2yEXEZw/RHiIiDhzjn+S3dyMzMuq7seWQvY/hfKAF8KtvMbCwRBRPZrkdiZtYTxk2YUKxFto8GhJXdImtm1pgG0laknJmZjUn+DWBm1mckvUzS+ZLulBSS9i47JjOzbnAia2bVVGuRLbJZvYnA1cDMsgMxM+smdy0ws2rq0mCvsSAiLgQuBJA/H7Mxo/CsBQP90wDgRNbMqsmJ7Kjx0t5m1qv6JyU3s/4yMFh8s5GaDSzNbbeXG46ZWTFOZM2smrq0spc15KW9zfrA4MAAGw0ONt0G2+haIGmmpEWSVkmaL2n3Ycq+WdICSQ9IelDSVZLeWVfmjGwwan67qNW43LXAzKrJXQtGTUSsBlbX9t2v1szyJO1LWrTqUGA+MAu4WNJOEXF3g1PuA/4LuB5YA7wBOF3S3XWrIl4EHJLbX02LnMiaWTV5Hlkzs6o4EjgtIk4HkHQo8HrgXcDn6gtnC17lfUXSQcBLgHwiuzoiFo8kMP8GMLOKUgub5UnaTNIuknbJDm2f7W9XZlxmVjmTJE3ObRPqC0gaD0wD5tWORcS6bH/PZjdQMgPYCfht3cvTJd0t6e+Svi5py1bfgFtkzaya3LVgJHYDfp3bn5P9eSZw8KhHY2ajYnBwkMECy88OrltX+2v9wM4TgOPrjm0FDAJL6o4vAZ451D0kbQ7cQZoRZS3wgYj4Za7IRcDPgFuApwOfBS6UtGdErG36JjJOZM2smpzIti17rOcPxsya2RZYnttvuY/qMJYDuwCbATOAOZJurnU7iIgf5cpeK+ka4CZgOvCrojdxImtm1SQV7CPrfM3MrE3LI2JZkzL3klpUp9QdnwIM2b81635wY7Z7laRnkab6u2yI8jdLuhfYgRYSWfeRNbNqchdZM7OWFJl6q7YVFRFrgIWkVlUAJA1k+1e0EN4A6y+8sh5J2wJbAne1cE23yJpZVQ2CilS2XhDBzKzL5gBnSloAXEmafmsiUJvF4CzgjoiYne3PBhaQugpMAF4HvBM4LHt9M+A44BxSq+7TgS+QWnDzsxo05UTWzKrJfWTNzCohIs6WtDVwIjAVuArYKyJqA8C2A9blTpkInErqg7uSNJ/sOyLi7Oz1tcDzgIOALYA7gUuAY7N5rQtzImtm1eR5ZM3MKiMi5gJzh3htet3+McAxw1xrJfDaTsTlRNbMqsktsmZmLRkcGCi0/Gw7S9RWlRNZM6uooiO5nMiamY1VTmTNrJrcImtmZk04kTWzanIfWTOzlrSxslfPcyJrZtXkFlkzM2ti7CSyAxvB4Nh5uzZ61kWUHUJPaPVzkoQKJKlFylh7xj1uY8ZP3KTsMMzMhuTMzsyqyWO9zMysCSeyZlZJbpE1M2vNgMRggTpxoI/qTSeyZlZJTmTNzKyZMZPIXvD1DzJ58uSyw+gJrzzw82WH0FMuPevjZYfQE1ptAXAia2bWmrE4a0Gp89ZImi3pT5KWS7pb0rmSdqors7GkUyT9S9IKSedImlJWzGZmZmZWDWVPwPhy4BRgD+DVwDjgEkkTc2VOBv4D2Ccr/0TgZ6Mcp5mNNrWwmZnZmFRq14KI2Cu/L+lg4G5gGvBbSZsD7wb2j4hLszKHAH+TtEdE/HGUQzazUeKuBWZm1kzV+shunv15X/bnNFIr7bxagYi4XtKtwJ7ABomspAnAhNyhSd0J1cy6yYmsmZk1U5lEVtIA8GXg8oj4S3Z4KrAmIh6oK74ke62R2cBx3YjRzEZPSmSb935yImtmlgxIhQbW9tP0W2X3kc07BXgOsN8Ir3MSqWW3tm07wuuZWRncR9bMzJqoRCIraS7wBuAVEXF77qXFwHhJW9SdMiV7bQMRsToiltU2YHk3Yjaz7pKKb61fWzMlLZK0StJ8Sbs3KT9L0t8lrZR0m6STJW3c7nszM+uGAYmBgYHmm1tkO0PJXOBNwCsj4pa6IguBh4EZuXN2ArYDrhi1QM1s1Ak92k922K3FJllJ+wJzgBOAXYGrgYslbTNE+f2Bz2Xln0UagLov8Nn23133FJnW0MysX5TdInsK8A5gf2C5pKnZtglARCwFvgPMkfQKSdOA04ErPGOBWX9TC/+16EjgtIg4PSKuAw4FHgLeNUT5F5H67v8gIhZFxCXAD4FhW3FLVGRaQzOzvlB2InsYqR/rZcBduW3fXJkjgAuAc4DfkroUvHlUozSz0deFPrKSxpNmQ8nPhLIu299ziNP+AEyrdT+Q9DTgdcD/tvJ2RktE7BURZ0TEXyPiauBg0lOsaeVGZma9rJUuWZLeLGmBpAckPSjpKknvrCsjSSdKuivrtjVP0o6txlX2PLJNfwVFxCpgZraZ2RjRxvRbk+rKr46I1XXFtwIGSTOf5C0Bntno+hHxA0lbAb9XusFGwDciopJdCxqon9bQzKwluS5ZhwLzgVmkLlk7RcTdDU65D/gv4HpgDWkc1OmS7o6Ii7MyRwEfAg4CbgE+nV3z2VnuV0jZLbJmZo21PtrrdmBpbpvdmTA0HTga+ACpT+2bgddLOrYT1++mIaY1bFRugqTJtQ3Pv23WkwalwluLWuqSFRGXRcTPI+JvEXFTRHwFuAZ4CaTWWFIy/JmIOC8irgEOJK3euncrgVVmHlkzs7w25kPclvVnKalvjQW4F1hLmvkkb8iZUEitBN+LiG9n+9dm/U2/Jem/sq4JVVWb1vAlTcp5/m2zsanpk6xcl6yTasciYp2k4bpk5c8X8EpgJ+Dj2eHtSesB5Lt5LZU0P7vmj4q+AbfImlk1td5Hdnl+6r0G3QqIiDWk2VDyM6EMZPtDzYSyKVCfrK7NRVlJw0xr2Ijn3zbrA4Wm3sq2TJEnWcN1yRpqcSokbS5pBalrwS+AD0bEL7OXa+e1dM1G3CJrZpXUxSVq5wBnSloAXEl6vDWRNCMKks4C7oiIWoV+PnCkpD+T+obtQGqlPT8i1lIxWevH10jTGk5vMK3hBrKk/9HE36ulmY0ZRZ5ktWs5sAuwGamxYI6kmyPisg7ew4msmVVUwUS21RURIuJsSVsDJ5K++V8F7BURtZaB7Vi/BfYzQGR/Pgm4h5TcfrKlG4+eU0hTGv4n2bSG2fGlEbGyvLDMrIKWZ4tHDaedLlm1GWFuzHavkvQsUovvZbnzppBmq8pf86oigdc4kTWzaio6tVYbjYcRMReYO8Rr0+v2HyEthnBC63cqxWHZn5fVHT8EOGNUIzGznhcRayTVumSdC+t1yWpYjw5hAJiQ/f0WUjI7gyxxzQaavhD4eivxOZE1M+sjRaY1NDNrUUtdsiTNBhYAN5GS19cB7yT7oh0RIenLwDGSbuCx6bfuJEuWi3Iia2aV1MU+smZmfWlwYIDBgebj+IuUyWujS9ZE4FRSH9yVpPlk3xERZ+fKfCEr9y1gC+D32TULzyELTmTNrKKcyJqZVUeLXbKOAY5pcr0APpVtbXMia2aVVDdFzLDlzMysrfm3e55/A5iZmZlZT3KLrJlVUxdnLTAzs/7gRNbMKsl9ZM3MrBknsmZWSVKxtQ6cx5qZjV1OZM2sktwia2bWooEBVGQAbB8NknUia2aV5C6yZmbWjBNZM6sm9y0wM2vJoMRggTqxSJle4UTWzCpJFMxjux6JmZlVlRNZM6ukwQExONA8TV1XoIy158nTXs7kyZPLDqMnrFhxU9khWB9asWJ52SFUnhNZM6smd5I1M7MmnMiaWSW5i6yZWWvG4mwvTmTNrJLGYoVsZmatcSJrZpWUWmSLJLKjEIyZmVWSE1kzqyR3kTUza83gwACDBRY7KFKmVziRNbNqcidZMzNrwomsmVWS81gzM2um1LZlSYdJukbSsmy7QtK/517fWNIpkv4laYWkcyRNKTNmMxsdGhADBTZ5HlkzszGr7E4StwOfAKYBuwGXAudJ2jl7/WTgP4B9gJcDTwR+VkKcZjbK1MJmZmaPzfZSZOsXpSayEXF+RPxvRNwQEf+IiE8CK4A9JG0OvBs4MiIujYiFwCHAiyTtUWbcZjYKnMmamVWGpJmSFklaJWm+pN2HKfteSb+TdH+2zasvL+kMSVG3XdRqXGW3yD5K0qCk/YCJwBWkVtpxwLxamYi4HrgV2LOUIM1s1IzFlgUzsyqStC8wBzgB2BW4GrhY0jZDnDId+CHwClLOdhtwiaQn1ZW7CHhCbnt7q7GVnshKeq6kFcBq4BvAmyLiOmAqsCYiHqg7ZUn22lDXmyBpcm0DJnUpdDProtpgryKbPabZ2AMz618DUuGtRUcCp0XE6VmOdijwEPCuRoUj4oCIODUirsoaId9Dyjln1BVdHRGLc9v9rQZWeiIL/B3YBXgh8HXgTEnPHsH1ZgNLc9vtIw3QzEafW2Tb1mzsgZlZYZLGk+qT/BPyddl+0Sfkm5Kest9Xd3y6pLsl/V3S1yVt2Wp8pU+/FRFrgBuz3YWSXgB8GDgbGC9pi7pW2SnA4mEueRKp+btmEk5mzXqOF0RoT0ScX3fok5IOA/YA/lpCSGZWXZPqGgNWR8TqujJbAYOkJ+J5S4BnFrzP54E7ySXDpG4FPwNuAZ4OfBa4UNKeEbG24HXLT2QbGAAmAAuBh0nN0OcASNoJ2I7Uh7ah7B/g0X8Et9aY9ShPJDtikgZJs77Uxh6YmeXVN/SdABzfyRtI+gSwHzA9IlbVjkfEj3LFrpV0DXATqX/tr4pev9REVtJJwIWkAVyTgP1Jb+C1EbFU0neAOZLuA5YBXwOuiIg/lhSymY0S57Htk/RcUuK6MWkmmNrYg6HKTyA1INR4bIFZDxoYGGCgwPKzuTLbAstzL9W3xgLcC6wlPRHPa/aEHEkfJXV1elVEXDNc2Yi4WdK9wA70SiILbAOcRRqpthS4hpTE/jJ7/QhgHalFdgJwMfCBEuI0s1E2oLTgQZFytoHa2IPNgbeSxh68fJhkdjZw3CjFZmbVsTwilg1XICLWSFpIekJ+LoCk2sCtuUOdJ+ko4JOkvG5Bs0AkbQtsCdxVOHpKTmQj4t1NXl8FzMw2MxtD3Ee2fcOMPXj/EKd4bIFZHyg6I0EbDQBzSF+IFwBXArNIXZZOB5B0FnBHRMzO9j8OnEh60r5IUm22qRURsULSZqQvz+eQWnWfDnyBVG9d3EpgZbfImpk15r4FnVQbe9CQxxaY2XAi4mxJW5OS06nAVcBeEVEbALYd6Ql6zWHAeOCndZeq9cFdCzwPOAjYgjQQ7BLg2AaDzYblRNbMKsl5bHuGG3tQYlhm1uMiYi5DdCWIiOl1+09tcq2VdKhOciJrZpVUdI5Ytx5uoNnYAzOzvuFE1swqyS2y7Wk29sDMrJ9UYWUvM7MNqIX/Wr62NFPSIkmrJM2XtHuT8ltIOkXSXZJWS/qHpNe1/ebMzLpggILL1JYdaAe5RdbMKqlbLbKS9iWNwD0UmE8afXuxpJ0i4u4G5ccDvwTuJk1ldQfwFOCB1u5sZmad1nIiK2kTQBHxULb/FOBNwHURcUmH4zOzMWpwQAwWmEe2SJk6RwKnRURt2phDgdcD7wI+16D8u4DHAy+KiIezY4tavamZWbeNxbEF7bQunwccCOlxG6lF4yPAedl63mZmI6cWtqKXTK2r08it9x0R67L9PYc47Y2kVbJOkbRE0l8kHZ0t/2pmZiVqJ5HdFfhd9ve3AktIj9kOBD7UobjMzFrNYSdJmpzbGs2buhUwSKq38paQ5kZs5Gmkum4QeB3wadKX92Naf0dmZtZJ7SSym/LYuryvAX6WtWj8kZTQmpmNWBuDvW4nTTdV22Z3KJQBUv/Y90XEwog4G/gvUh/bjpF0kKTX5/a/IOkBSX/IunCZmVmddhLZG4G9JT2ZNJltrV/sNsCw6/WamRU1oOJbZltg89x2UoPL3ktaUWZK3fEppGUSG7kL+EdErM0d+xswNeuq0ClHAysBJO1JWpr7qCzmkzt4HzOzvtFOInsi8CXSYIf5EXFFdvw1wJ87FJeZjXW1aQuKbMnyiFiW2zZY5jAi1gALgRmP3UYD2f4V9eUzlwM7ZOVqngHclV2vU55MaigA2Bs4JyK+RWpZfmkH72NmfUoDAwwU2DTQPxNwtfxOIuKnpDV1dwP2yr30K+CIDsVlZmNc63lsYXOA92aP8p8FfB2YCNRmMTgrW+a15uukWQu+IukZ2eP/o4FTRvoe66wAtsz+/hrSlF8Aq4BNOnwvM7O+0NY8shGxmLrHcBFxZUciMjPjsUm9i5RrRUScLWlr0tOlqcBVwF4RURsAth2wLlf+NkmvJT3ev4Y0j+xXgM+3dOPmfgl8W9KfSS2+/5sd3xlP92VmBRSdyKV/Jt8qmMhK+lnRC0bEm9sPx8ws6WaFHBFzgblDvDa9wbErgD3auFUrZgKfIXUxeEtE/Cs7Pg34YZfv3dCaNfexZs3DzQuamZWkaIvs0q5GYWZWRwNCBRY7KFKm6pRmJ9+a1I3h7xHxSO21iDiutMDMzCquUCIbEYd0OxAzs7y6GQmGLdfLJG0P/A/w7OzQ7ZLeEhELSgzLzKwn9M+wNTPrK11Y2KuqvkhqVHgHaeGF24FvlhqRmVmPKNpH9s9AFCkbEbuOKCIzMyg+JUHvrxn+EuCtEfF7AEl/JLXKToyIB8sNzcx6yaDEYIE6sUiZXlG0j+y53QzCzKzeAMUeGfXBY6VtgBtqOxFxl6SV2fFbSovKzKwHFO0je0K3AzEzy5OECrQaFClTcQFsliWvNeuASZImP1oowisnmtmwVHCxgzG9IAKApC0kvUfSSZIenx3bVdKTOhuemY1VXVwQoWoE/AO4P7dtRlop8X7ggexPM7PSSJopaZGkVZLmS9p9mLLvlfQ7Sfdn27z68kpOlHSXpJVZmR1bjavlBREkPQ+YR5qS66nAacB9wJtJE4kf2Oo1zcw2VDRL7flM9hVlB2Bm/aFbC8lI2pe0KuKhwHxgFnCxpJ0i4u4Gp0wnzX/9B9LqhB8HLpG0c0TckZU5CvgQcBCpG9Wns2s+OyJWFY2tnZW95gBnRMRRkpbnjv8v8IM2rmdmtoEBFRuQ0OvTb0XEb8qOwcysiSOB0yKitpT3ocDrgXcBn6svHBEH5PclvQd4CzADOCubO3sW8JmIOC8rcyCwBNgb+FHRwNrpWvACGk8NcwdpuUczMzMz6wOSxpNWGJxXOxYR67L9PQteZlNgHOkJPsD2pJwxf82lpNbeotcE2muRXQ1MbnD8GcA9bVzPzGwDY2hBhHU0n94wIqKd+trMbDiT6gbMro6I1XVltgIGSa2leUuAZxa8z+eBO3ksca01fDa6ZkuNou1UjP8DfErS27L9kLRdFuQ5bVzPzGwDY2caWd40zGt7kvqQ9c8QYzOrktvr9k8Aju/kDSR9AtgPmN5K39ei2qkcP0IaUXs3sAnwG+BGYDnwyXYDkfQJSSHpy7ljG0s6RdK/JK2QdI6kKe3ew8x6R236rSJbL4uI8+o34HrgYOCjwE+Andq9fqO61cz6UxsrIm4LbJ7bTmpw2XuBtUB9/jUFWDxsPNJHgU8Ar4mIa3Iv1c5r+Zr1Wk5kI2JpRLwaeAOppWAu8LqIeHm7q9BIegHwfuCaupdOBv4D2Ad4OfBE4Gft3MPMeswYWqO2RtITJZ0GXEt6YrZLRBwUEf9s83pD1a1mZgDLI2JZbqvvVkBErAEWkgZqASBpINu/YqgLSzoKOBbYKyIW1L18CylhzV9zMvDC4a7ZSNt9riLicuDyds+vkbQZ8H3gvcAxueObA+8G9o+IS7NjhwB/k7RHRPxxpPc2s+oaQAwUyFKLlKm6rL47GvggcBUwIyJ+N8JrNqxbzax/dXFBhDnAmZIWAFeSZhyYCNRmMTgLuCMiZmf7HwdOBPYHFkmq9XtdERErIqL2lOgYSTfw2PRbd9LiarItvxNJX5X0oQbHD2/z0dUpwC8iYl7d8WmkEW75EW3XA7fS4og2M+s9Y2VBhKzV4mbSU663R8SLRprEZoaqW83MWhIRZ5O6Op1I+rK9C6mltTZYazvgCblTDgPGAz8F7sptH82V+QLwNeBbwJ9I3Vb3arUfbTstsm8B3tjg+B9I/SBmFb2QpP2AXUlTetWbCqyJiAfqjg87ok3SBGBC7tCkovGYWXWMoSVqPwesJI01OEjSQY0KRcSbi16wSd3aqLzrTTMbVkTMJXUnbfTa9Lr9pxa4XgCfyra2tZPIbkla1aveMtIUDYVIejLwFeDVHR7FNhs4roPXM7MSdGuFmgo6i+bTbxXWZt3qetPMelI7ieyNwF5smJX/O+nxWFHTgG2A/8u1qAwCL5N0OPBaYLykLepaZZuNaDuJ1JejZhIbTi9hZhU3QLG+T70+L1VEHNzhSzarWydExNq6c1xvmllPaneJ2rmStgYuzY7NIE3LNauF6/wKeG7dsdNJ0858HrgNeDi79jkAknYi9cMYckRbNuLu0VF3ffDY0WxMGkPzyDYk6SmkwRTXZ6voFDVs3dogiXW9adYnxtCTrEe1nMhGxHez/lSfJE2rALAIOCwizmrhOsuBv+SPSXoQ+FdE/CXb/w4wR9J9pK4LXwOu8IwFZv1vrCSykt4FbBERc3LHvkWatQXg75JeGxG3FblekbrVzKxftPVULiK+HhHbkh7zT46Ip7WSxLbgCOACUovsb0ldCgoPeDCz3lVborbI1uPeB9xf25G0F3AIcCBpsNYDuP+qmRVQa5EtsvWLEa3dHRH3dCqQ7HrT6/ZXATOzzczGkqKLHfR+fbwjkJ8s/D+B8yLi+wCSjiabq7Fd9XWrmVm/6PVxEmbWp5QtiNBsU+9nspuQuk7VvIj0BKrmZoaZctDMbCxzImtmlTTQwtbj/kmaaQBJWwE7s/6qiVNpPOWhmdmYN6KuBWZm3TIwIAYKdIAtUqbizgROkbQz8ErSLAULc6+/iLrBW2ZmljiRNbNKGiuzFpCWadyUNJB1MbBP3esvBn442kGZWe8ZQysiPqpjiaykKcD7I+LETl3TzMausZLIZnPEDrlMY0TUJ7ZmZpbpZPeyqXiKGDPrkCIDvWqbmZl5+q1hSXpekyI7jTAWM7NHjZUWWTMza18rXQuuAoLGszbWjkcHYjIzS9PIFklkux6JmZlVVSuJ7H3AUaR1vBvZGTh/xBGZmUE2Q2yBQQtOZc3MxqxWEtmFwBMj4p+NXpS0BW4cMbMOKbr8bO/PvmVmZu1qJZH9BjBxmNdvJa0PXklX3HIvEzdbXXYYPeHE495Vdgg95fc3dXSl5r714IrlLZUfS31kJT0BOAx4CfAEYB1pRa9zgTMiYm0ZcT1w642s3WzTMm7dc8ZtPqHsEKwPPbj8wdZOKDj9Vl9UnJnCsxZExM8j4r+Hef3+iDizM2GZ2Vg3qOJbL5O0G/A34HXAOGBH0hOwB4EvAb+VNKm8CM3MqqsPVnc0s76kx1plh9va6dAkaaakRZJWSZovafeC5+0nKSSd2/pdh/Rl4OSI2C0iXgocDDwjIvYDnkZaLOEzHbyfmfWp2oIIRbZ+0VIiK+lwSWdJ2i/bf6ek6yRdL+mzkrxSmJl1RLcqZEn7AnOAE4BdgauBiyVt0+S8p5JaSH/X1hsa2q7A93L7PwB2lTQlIu4nDbJ9a4fvaWbWklYaACTtLOmcrHxImtWgzPHZa/nt+lbjamUe2WNIFeolwMmSngJ8DDiZ1J/rCOBhKrooQr99A+mmdeFZ1KzzWv2/qouDvY4ETouI0wEkHQq8HngX8LlGJ0gaBL5Pqt9eCmzR8l2HdjepX+zN2f4UUt28LNu/AXh8B+9nZtaSXAPAocB8YBapAWCniLi7wSmbkuq0n5DyxKH8FXhVbv+RVmNrpQX1YODgiPiZpH8j9eE6KCK+D5Bl0V+goomsmfWWor0GcmUm1X1ZXR0R643wlDQemAacVDsWEeskzQP2HOY2nwLujojvSHppgbBacS7wDUkfA1YDxwK/iYiV2es7AXd0+J5mZq1oqQEgIv4E/Ckr27CBIPNIRCweSWCtdC14IrAAICKuJrXCXpV7/f+yMmZmI6YW/svcDizNbbMbXHYrYBBYUnd8CWmZ7Q3jkF4CvBt4bwfeViPHANeR5uH+FTCB9MuhJmj8XszM1jPQwlZUrgFgXu1YRKzL9odrAChiR0l3SrpZ0vclbdfqBVppkV0MPBu4VdKOpF8GzyY1C0NaEKFR87KZWcvamH5rWyA/x9eI59vLZgv4HvDeiLh3pNdrJCJWAPtK2hjYKNvPv35JN+5rZkaBJ1kM3wDwzBHcez7paf/fSd2rjgN+J+k5EVF4vsZWEtnvA2dJOg+YQepG8CVJW5JaDD4J/LSF65mZDamNPrLLI2LZMEUB7gXWkvqh5k0hfVmv93TgqcD5ucp+AEDSI8BOEXFT8yibi4hVnbiOmVkLbq/bPwE4fjRuHBEX5navkTQf+CfwNuA7Ra/TSiJ7HLCS1Ix8GqlPxNWkhHZT0mOxY1u4npnZkAYGxECBTLZImZqIWCNpIenL+LkAkgay/bkNTrkeeG7dsc8Ak4APA7cVvrmZWZdpYAANNO84kCtT5ElWqw0AbYmIByT9A9ihlfMKJ7JZf4jP1h3+UbaZmXXUAGKgwHCvImXqzAHOlLQAuJI0+nYiUBvEcBZwR0TMzlpJ/5I/WdIDABGx3nEzsx7U9ElWGw0AbZG0Gekp2Pealc3zvK9mVk1tTFtQREScLWlr4ETSAK+rgL0iotb/azvSYFYzM0sKNwBk++NJ46gAxgNPkrQLsCIibszKfIn0NP+fpMkCTiC1/P6wlcCcyJpZJRUdWdvO8oQRMZchWhIiYnqTcw9u45ZmZj2rjQaAJwJ/zu1/NNt+A0zPjm1LSlq3BO4Bfg/sERH3tBKbE1kzq6Q2Zi0w0mo5bDif998jYiSji82sBwwAAwUqxW43AETEIpo8L8uW4R4xJ7JmVkld6lkwVox4tRwzs17gRNbMKqnostJeerqhEa+WY2bWC5zImlkluUV2RHaUdCewCrgCmB0Rtw5VWNIE0opiNZO6HJ+ZdcFYrDfb6SbRMZKOlxR12/W51zeWdIqkf0laIekcSfXzmJlZP9Jj/WSH2/qqRu6M2mo5ewGHAduTVssZLjmdzfrL+9ZPkm5mVkmlJrKZv5KWJqttL8m9djLwH8A+wMtJo+B+NtoBmtnoq80jW2Szx0TEhRHxk4i4JiIuBl4HbEFaLWcoJwGb57Ztux6omVkHVKFrQcO+XJI2B94N7B8Rl2bHDgH+JmmPiPhjKzcputylwbooOwIz95HtlCKr5WRrqz+6oo8/UzPrFVVokd1R0p2Sbpb0fUnbZcenAeOAebWCEXE9cCtpmdyGJE2QNLm24b5eZj1JLWw2tNxqOXeVHYuZdVetAaDI1i/KTmSH68s1FVgTEQ/UnbMke20o7utl1gfGYoXcCZK+JOnlkp4q6UXAz2ljtRwzs15QateCiLgwt3uNpPmkpcreBqxs87InkZZSq5mEk1mznjMWR992SEdWyzGz3jMgFVsQoY8aAKrQR/ZRdX25fgmMl7RFXavsFGDI+RHd18usP3hlr/Z0arUcM7NeUHbXgvXU9eVaCDwMzMi9vhNpPd8rSgnQzEaN+8iamVkzpbbISvoScD6pO8ETgRPI+nJFxFJJ3wHmSLoPWAZ8Dbii1RkLzKxHOUs1M7NhlN21oFlfriOAdcA5pFVnLgY+UEKcZjbKBhGDBTLZImXMzKw/lT3Ya9i+XBGxCpiZbWY2hqRuAwXmke1+KGZmPWEszr9ddousmVljnrbAzMyacCJrZpVUdPlZL1FrZpa4RdbMrCI8/ZaZmTUzZhLZAYpNEmywNtaVHYL1obURLZV3ImtmZs2MmUTWzHqNO8mamdnwnMiaWSVJMOAWWTMzG8aYSWR3e8rjmTx5ctlhWB965YGfLzuEnvDIw6taKu/22PLNP3cBm06YUHYYPWHfE04oOwTrQ5suW9ZS+W7Wm5JmAh8DpgJXAx+MiCuHKLszcCIwDXgKcEREfHkk1xxKpZaoNTOrUQv/mZlZ90jaF5hDWoF1V1LSebGkbYY4ZVPgZuATwOIOXbMhJ7JmVkkDKr6ZmRkMSIW3Fh0JnBYRp0fEdcChwEPAuxoVjog/RcTHIuJHwOpOXHMoTmTNrJLcImtm1nWTJE3ObRv0JZI0ntRFYF7tWESsy/b3bOemnbymE1kzq6Ta9FtFNjMza8vtwNLcNrtBma2AQWBJ3fElpL6t7ejYNcfMYC8zMzMzW8+2wPLc/lDdACrLiayZVVLRflxe6MTMrG3LI6LZ1Aj3AmuBKXXHpzDEQK4COnZNdy0ws0pSC5uZmYGkwltREbEGWAjMyN1nINu/op04O3lNt8iaWTV5Ilkzs6qYA5wpaQFwJTALmAicDiDpLOCOiJid7Y8Hnp2dOx54kqRdgBURcWORaxblRNbMKqnojASetcDMLCna2tpKiyxARJwtaWvSIgdTgauAvSKiNlhrO2Bd7pQnAn/O7X80234DTC94zUKcyJpZJblB1sysOiJiLjB3iNem1+0vokD1PNw1i3Iia2bV5EzWzMyacCJrZpXkPNbMzJrxrAVmVkleEKF9kp4k6b8l/UvSSknXStqt7LjMzDrNiayZVVI3p9+SNFPSIkmrJM2XtPswZd8r6XeS7s+2ecOVL5ukxwGXAw8D/04aOfwR4P4y4zKz7lM2/3azrdXBXlXmrgVmVk1d6lsgaV/StC+HAvNJU75cLGmniLi7wSnTgR8CfwBWAR8HLpG0c0Tc0drdR8XHgdsi4pDcsVvKCsbMrJvcImtmFdW1NtkjgdMi4vSIuI6U0D4EvKtR4Yg4ICJOjYirIuJ64D2kunNGo/IV8EZggaSfSLpb0p8lvXe4EyRNkDS5tgGTRidUM+uksbiQjBNZM6ukNvrITsonY5ImbHhNjQemAfNqxyJiXba/Z8HQNgXGAfeN5P110dOAw4AbgNcCXwe+KumgYc6ZDSzNbbd3O0gzs05wImtmldRGy8LtrJ+MzW5w2a2AQaB+wu0lpAm5i/g8cCe5ZLhiBoD/i4ijI+LPEfEt4DRSy/NQTgI2z23bdj9MM7ORcx9ZM6ukNlb22hZYnntpdcdjkj4B7AdMj4hVnb5+h9wFXFd37G/AW4Y6ISJWk/u8+mkgiJn1t9JbZJtNE6PkREl3Za/Pk7RjmTGb2Sgo2q3gsZxreUQsy22NEtl7gbXAlLrjU4DFw4YjfRT4BPCaiLhmRO+tuy4Hdqo79gzgnyXEYmbWVaUmsgWniTkK+BDpsdgLgQdJI4w3Ht1ozazXRcQaYCG5gVqSagO3rhjqPElHAceS1gFf0O04R+hkYA9JR0vaQdL+wPuAU0qOy8y6TNnUWkW2flF214Jhp4lR+qRnAZ+JiPOyYweS+rPtDfxo1CI1s9HVvaW95gBnSloAXEmqYyYCpwNIOgu4IyJmZ/sfB04E9gcWSar1pV0REStavnuXRcSfJL2J1O/1U6Q6dVZEfL/cyMzMOq/sRPaNpNbVnwAvB+4ATo2I07LXtycNwMiPMF4qaT5phLETWbM+1UYf2UIi4mxJW5OS06nAVaSW1toAsO2AdblTDgPGAz+tu9QJwPEt3XyURMQFwAVlx2Fmo6toa6tbZDunNk3MHOCzwAtI08SsiYgzeWwUceERxtmUO/lpdzwfolkPKrr8bDv1cUTMBeYO8dr0uv2ntn4HMzMbDWUnsgPAgog4Otv/s6TnkPrDntnmNWcDx3UiODMrT/d6FpiZWb8oe9aCoaaJ2S77e20UcSsjjD0folkfGIuDFszMrDVlJ7LNpom5hZSw5kcYTybNXtBwhHFErM5PwcP680qamZmZ9aUBqfDWL8ruWnAy8AdJRwM/BnYnTRPzPoCICElfBo6RdAMpsf00aVWdc8sI2MxGh7sWmJlZM6UmsgWnifkCaWqcbwFbAL8njTCu6qo6ZtYJzmTNzKyJsrsWEBEXRMRzI2LjiHhWbuqt2usREZ+KiKlZmVdFxD/KitfMRoda2MzMjGLLIRadEmaDS2umpEWSVkmaL2n3JuX3kXR9Vv5aSa+re/0MSVG3XdRqXKUnsmZmjYiC9XHZgZqZ9TlJ+5KmSj0B2BW4mrQOwDZDlH8R8EPgO8DzSd1Bz81mpsq7CHhCbnt7q7E5kTWzSnKLrJlZZRwJnBYRp0fEdaRpUh8C3jVE+Q8DF0XEFyPibxFxLPB/wOF15VZHxOLcdn+rgTmRNbNqciZrZtZtkyRNzm0T6gtIGg9MY/1VVtdl+3sOcd098+UzFzcoP13S3ZL+LunrkrZs9Q2UPWuBmVlDY3Gpxap5y9FHM3ny5LLD6Alv2G23skPoKRcsWFB2CH1pgGItlLkyt9e91Gjp7a2AQRqvsvrMIW4xdYjy+VVZLwJ+Rhro/3TSCq8XStozItYOF3+eE1kzqyRPWmBm1nXbsv58+6tH68YR8aPc7rWSrgFuAqYDvyp6HXctMDMzMxublucXkYqIRonsvcBaWltldXGL5YmIm7N77VAo8owTWTOrJLXwn5mZdWdp74hYAyxk/VVWB7L9hqusZsdn1B179TDlkbQtsCVwV+HgcNcCM6uoolMduousmVnXzQHOlLQAuBKYRVqs6nQASWcBd0TE7Kz8V4DfSPoI8AtgP2A3spVbJW0GHAecQ2qlfTppAawbSYPCCnMia2ZmZmZDioizJW0NnEgasHUVaZXV2oCu7YB1ufJ/kLQ/8BnSIK4bgL0j4i9ZkbXA84CDSKu23glcAhw7RPeGITmRNbNKcousmVl1RMRcYO4Qr01vcOwnwE+GKL8SeG0n4nIia2YV5XkLzMxseE5kzayS3CJrZmbNOJE1s0oaUNqKlDMzM9h80kQ222STpuU2Gjc4CtGMDieyZlZR7lpgZmbD8zyyZmZmZtaT3CJrZpXk9tj2SFoEPKXBS6dGxMxRDsfMRtGErTdlwqabNi338EP9U3M6kTWzaio42MuZ7AZeAOQ7wD0H+CVDTINjZtbLnMiaWSW5RbY9EXFPfl/SJ4CbgN+UE5GZjZZxj9uY8RM3bl5uwrqmZXqF+8iaWTWphc0akjQeeAfw3YiIsuMxM+s0t8iaWSW5RbYj9iYt/3jGcIUkTQAm5A5N6lpEZtY1E7bYhAmbNZ9+a824/vle6xZZM6skSYU3G9K7gQsj4s4m5WYDS3Pb7d0OzMysE5zImlkl1Vb2KrLZhiQ9BXgV8O0CxU8CNs9t23YxNDOzjnHXAjOz/nQIcDfwi2YFI2I1sLq271Zus940bvIExk0qMNhrYO0oRDM6nMiaWSW5j2z7JA2QEtkzI+KRsuMxM+sWJ7JmVknK/itSzjbwKmA74LtlB2Jmo2f8xIlM2GyzpuVW98/sW05kzayi3CTbtoi4BH8yZjYGOJE1s0pyHmtm1pqNNno8G23UfPa8jTYaPwrRjI5SZy2QtEhSNNhOyV7fWNIpkv4laYWkcyRNKTNmMxslXhDBzMyaKLtFttma4CcDrwf2Ic1tOBf4GfDiUYzRzErgPrJmZq2ZMOHxTJgwuUC5caMQzegotUU2Iu6JiMW1DXgD2ZrgkjYnTeZ9ZERcGhELSaNwXyRpjxLDNrNRMKDiW6skzcyeCK2SNF/S7k3K7yPp+qz8tZJe1+77MjOzzqnMgggN1gSfBowD5tXKRMT1wK3AnsNcZ4KkybUNL7Vo1pu61LVA0r7AHOAEYFfgauBiSdsMUf5FwA+B7wDPB84FzpX0nBbfkZmZdVhlElk2XBN8KrAmIh6oK7cke20oXmrRrC90rZPskcBpEXF6RFwHHAo8BLxriPIfBi6KiC9GxN8i4ljg/4DDW72xmVk3jRu3ReGtX5TdRzav6JrgzZxEam2pmYSTWbOes3zZskIp6vJly2p/nVS3ItXqbMWqR2VPfqaR6gkAImKdpHkM/aRnT9avUwAuJn35NjOrjGWP1YcdKdcLKpHI5tYEf3Pu8GJgvKQt6lplp2SvNeSlFs163hpg8Q7bP3m4Jy/1VrDhF9YTgOPrjm1FGmC6pO74EuCZQ1x76hDlW4nPzKyb1gCLn/zklurNxdl5Pa0SiSyN1wRfCDwMzADOAZC0E2m1mitGO0AzGx0RsUrS9sBIJzpc3byImVnva7PeXBMRq7oV02gpPZEdak3wiFgq6TvAHEn3AcuArwFXRMQfy4nWzEZDVrl2o4K9F1hLerKTN9yTnsUtljczG3VdrDcrrQqDvYZbE/wI4AJSi+xvSb843tygnJlZUxGxhvS0Z0btWPZlegZDP+m5Il8+8+phypuZ2SgpvUV2uDXBs28XM7PNzKwT5gBnSloAXAnMAiYCpwNIOgu4IyJmZ+W/Qprb+iOk7k/7AbsB7xvluK3CLliwoOwQesopRx5Zdgg9YeVq95BqpvRE1sxsNEXE2ZK2Bk4kDdi6CtgrImoDurYD1uXK/0HS/sBngM8CNwB7R8RfRjVwMzPbgBNZMxtzImIuacnrRq9Nb3DsJzy2dLaZmVVEFfrImpmZmZm1zImsmZmZmfUkJ7JmZmZm1pOcyJqZmZlZT3Iia2ZmZmY9yYmsmZmZmfUkJ7JmZmZm1pOcyJqZmZlZT3Iia2ZmZmY9yYmsmVkfkTQo6dOSbpG0UtJNko6VpLJjMzPrNC9Ra2bWXz4OHAYcBPwV2A04HVgKfLXEuMzMOs6JrJlZf3kRcF5E/CLbXyTp7cDuJcZkZtYV7lpgZtZf/gDMkPQMAEn/BrwEuLDUqMzMusAtsmZm/eVzwGTgeklrgUHgkxHx/aFOkDQBmJA7NKm7IZqZdYZbZM3M+svbgAOA/YFdSX1lPyrpoGHOmU3qQ1vbbu92kGZmneAWWTOz/vJF4HMR8aNs/1pJTyElq2cOcc5JwJzc/iSczJpZD3Aia2bWXzYF1tUdW8swT+AiYjWwurbvmbrMrFc4kTUz6y/nA5+UdCtp+q3nA0cC3y01KjOzLnAia2bWXz4IfBo4FdgGuBP4JnBimUGZmXWDE1kzsz4SEcuBWdlmZtbXPGuBmZmZmfUkJ7JmZmZm1pOcyJqZmZlZT3Iia2ZmZmY9qdREVtKgpE9LukXSSkk3STpWuUkMlZwo6a6szDxJO5YZt5mZmZmVr+wW2Y8DhwGHA8/K9o8iTR9TcxTwIeBQ4IXAg8DFkjYe3VDNzMzMrErKnn7rRcB5EfGLbH+RpLcDu0NqjSVNIfOZiDgvO3YgsATYG/hR/QXNzMzMbGwou0X2D8AMSc8AkPRvwEuAC7PXtwemAvNqJ0TEUmA+sOfohmpmZmZmVVJ2i+zngMnA9ZLWAoPAJyPi+9nrU7M/l9SdtyT32nokTQAm5A5N6ly4ZmZmZlYVZbfIvg04ANgf2BU4CPiopINGcM3ZwNLcdvtIgzQzMzOz6im7RfaLwOciotbX9VpJTyElo2cCi7PjU4C7cudNAa4a4ponAXNy+5OA25cvW9apmM3W88jDq8oOoSes9efUc5a53rQuWbl6ddkh9IRVa9aUHULllZ3Ibgqsqzu2lsdaim8hJbMzyBJXSZNJsxd8vdEFI2I18OhPiKRJAE/f/skdDNvMRmAS4Ayp2iYBPPnJrjfNKsL15hDKTmTPBz4p6Vbgr8DzgSOB7wJEREj6MnCMpBtIie2ngTuBcwve405gW2B5RyMfuUmkbg9VjK2K/HkVV+XPahLpZ9KqzfVm7/Nn1Zoqf16uN4dRdiL7QVJieiqwDekf6pvAibkyXwAmAt8CtgB+D+wVEYWeU0ZEAHd0LuTOyK35sDwi/C2rCX9exVX8s6paPNaA683e58+qNRX/vKoWT6Uo1Vc22rIuEkuBzSv4Q1M5/ryK82dl/cr/bxfnz6o1/rx6V9mzFpiZmZmZtcWJbHlWAyeQG5hmw/LnVZw/K+tX/n+7OH9WrfHn1aPctcDMzMzMepJbZM3MzMysJzmRNTMzM7Oe5ETWzMzMzHqSE9kSSJopaZGkVZLmS9q97JiqStLLJJ0v6U5JIWnvsmOqKkmzJf1J0nJJd0s6V9JOZcdl1gmuN4tzvVmc683e50R2lEnaF5hDGh25K3A1cLGkbUoNrLomkj6jmWUH0gNeDpwC7AG8GhgHXCJpYqlRmY2Q682Wud4szvVmj/OsBaNM0nzgTxFxeLY/ANwGfC0iPldqcBUnKYA3RcS5ZcfSCyRtDdwNvDwiflt2PGbtcr3ZPtebrXG92XvcIjuKJI0HpgHzasciYl22v2dZcVnf2jz7875SozAbAdebNspcb/YYJ7KjaytgEFhSd3wJMHX0w7F+lbVYfRm4PCL+UnI4ZiPhetNGhevN3rRR2QGYWVecAjwHeEnZgZiZ9QjXmz3IiezouhdYC0ypOz4FWDz64Vg/kjQXeAPwsoi4vex4zEbI9aZ1nevN3uWuBaMoItYAC4EZtWPZo4wZwBVlxWX9Qclc4E3AKyPilrJjMhsp15vWTa43e59bZEffHOBMSQuAK4FZpKlSTi8zqKqStBmwQ+7Q9pJ2Ae6LiFvLiaqyTgH2B/4TWC6p1n9waUSsLC8ssxFzvdkC15stcb3Z4zz9VgkkHQ58jDRQ4SrgQxExv9SgKkrSdODXDV46MyIOHtVgKi6bZqeRQyLijNGMxazTXG8W53qzONebvc+JrJmZmZn1JPeRNTMzM7Oe5ETWzMzMzHqSE1kzMzMz60lOZM3MzMysJzmRNTMzM7Oe5ETWzMzMzHqSE1kzMzMz60lOZM3MzMysJzmRtUqQFJL2LjsOM7Ne4DrTLHEia6NC0lRJX5N0s6TVkm6TdL6kGWXHZmZWNa4zzYrZqOwArP9JeipwOfAAaa30a4FxwGuBU4BnlhWbmVnVuM40K84tsjYaTgUC2D0izomIf0TEXyNiDrBHoxMkPVfSpZJWSvqXpG9J2iz3+hmSzpX0UUl3ZWVOkTQuV2aCpC9JukPSg5LmS5re5fdqZjZSrjPNCnIia10l6fHAXsApEfFg/esR8UCDcyYCFwP3Ay8A9gFeBcytK/oK4OnZnwcBB2dbzVxgT2A/4HnAT4CLJO04grdkZtY1rjPNWuNE1rptB0DA9S2csz+wMXBgRPwlIi4FDgfeKWlKrtz9wOERcX1EXAD8ApgBIGk74BBgn4j4XUTcFBFfAn6fHTczqyLXmWYtcB9Z6za1cc6zgKvrWiMuJ33x2glYkh37a0SszZW5C3hu9vfnAoPAP6T1QpgA/KuNmMzMRoPrTLMWOJG1bruB1NerG4MTHq7bDx57yrAZsBaYlv2Zt6ILsZiZdYLrTLMWuGuBdVVE3EfquzUz68e1HklbNDjtb8C/1ZV/MbAO+HvBW/+Z1LqwTUTcWLctbulNmJmNEteZZq1xImujYSapgrxS0lsk7SjpWZI+BFzRoPz3gVXAmZKeI+kVwNeA70XEkgblNxAR/8iuc5akN0vaXtLukmZLen1n3paZWVe4zjQryImsdV1E3AzsCvwa+H/AX4BfkgYZHNag/EOk+RIfD/wJ+CnwK9LghVYcApyV3fPvwLmkEb23tvE2zMxGhetMs+IUEWXHYGZmZmbWMrfImpmZmVlPciJrZmZmZj3JiayZmZmZ9SQnsmZmZmbWk5zImpmZmVlPciJrZmZmZj3JiayZmZmZ9SQnsmZmZmbWk5zImpmZmVlPciJrZmZmZj3JiayZmZmZ9SQnsmZmZmbWk/4/ZU/n5qzfyrEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 700x400 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from vireoSNP.plot import heat_matrix\n",
    "\n",
    "fig = plt.figure(figsize=(7, 4), dpi=100)\n",
    "plt.subplot(1, 2, 1)\n",
    "im = heat_matrix(_model.ID_prob, cmap=\"Blues\", alpha=0.8,\n",
    "                 display_value=False, row_sort=True)\n",
    "plt.colorbar(im, fraction=0.046, pad=0.04)\n",
    "plt.title(\"Assignment probability\")\n",
    "plt.xlabel(\"Clone\")\n",
    "plt.ylabel(\"%d cells\" %(_model.n_cell))\n",
    "plt.xticks(range(_model.n_donor))\n",
    "\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "im = heat_matrix(_model.beta_mu, cmap=segpink, alpha=0.8, \n",
    "                 display_value=False, row_sort=True)\n",
    "plt.colorbar(im, fraction=0.046, pad=0.04)\n",
    "plt.title(\"Mean allelic ratio\")\n",
    "plt.xlabel(\"Clone\")\n",
    "plt.ylabel(\"%d SNPs\" %(_model.n_var))\n",
    "plt.xticks(range(_model.n_donor))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "# plt.savefig(\"you_favorate_path with png or pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagnosis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Double check if global optima is found\n",
    "We are using multiple initializations via `n_init` and choose the one with highest ELBO. However, this doesn't guarantee to be the global optima. To double check it, we can run the same scripts multiple time (without fixed random seed), and check if the same (best) ELBO is found.\n",
    "\n",
    "If yes, it is likely to be the global optima, otherwise, we need increase `n_init`, e.g., 300 for more search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rerun 0: -190779.74335041404\n",
      "rerun 1: -190779.74335041404\n",
      "rerun 2: -190779.74335041404\n"
     ]
    }
   ],
   "source": [
    "n_init = 50\n",
    "for i in range(3):\n",
    "    _model = BinomMixtureVB(n_var=AD.shape[0], n_cell=AD.shape[1], n_donor=3)\n",
    "    _model.fit(AD, DP, min_iter=30, n_init=n_init)\n",
    "    print(\"rerun %d:\" %i, _model.ELBO_iters[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Identify the suitable number of clones\n",
    "\n",
    "It is generally difficult to identify the number of clones, which is a balance between subclone resolution and analysis reliability. More clones maybe preferred, but there could be higher risk that the subclones are not genuine but rather technical noise.\n",
    "\n",
    "Here, we could use ELBO for different number of clones as an indictor for model selection. However, this is still imperfect. One empirical suggestion is to choose the `n_clones` when ELBO stops increasing dramatically, for example in the case below, we will pick 3 clones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_init = 50\n",
    "n_clone_list = np.arange(2, 6)\n",
    "\n",
    "_ELBO_mat = []\n",
    "for k in n_clone_list:\n",
    "    _model = BinomMixtureVB(n_var=AD.shape[0], n_cell=AD.shape[1], n_donor=k)\n",
    "    _model.fit(AD, DP, min_iter=30, n_init=n_init)\n",
    "    _ELBO_mat.append(_model.ELBO_inits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEHCAYAAADoL5IPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlZ0lEQVR4nO3de3xU9Z3/8dfbhCQSREARKJFCq7YI2lYisK27673Y9ad267ZQW0ApbLfqdtvdn5elv9VeH623XtTWoqhoAXW7dWWtVlF0W9sGDEq9gRpvFURFuQkpl4TP74/5BoaYEIQkZyZ5Px+Pecw5n+/3nPnO6PDOOec7M4oIzMzMOts+WQ/AzMy6JweQmZllwgFkZmaZcACZmVkmHEBmZpaJ0qwHUCwOPPDAGDp0aNbDMDMrKosXL34rIvq31OYA2k1Dhw6ltrY262GYmRUVSa+01uZTcGZmlgkHkJmZZcIBZGZmmXAAmZlZJhxAZmaWCQeQmZllwgFkZmaZ8OeAzMwKUOO2YMPmBjam2zvpfsOmhu31SR8fiqSsh7rHHEBmZu2koXEbGzc3smHLzkGxoem2Ka1vyVve3NSvcfty7TdO2q3HO3s3x1Wov/vmADKzbm1r47bcEcamBjZuadixvLmRDZu3smFzY8shsvndAbNp67bdeszy0n3oVV5Kr4pSKstK6VVeSv/9yhl6YCW9yks54+6nqSzP1bf3y18vL2XIAZVs27bNR0BmZp1pS8O2dx9d7DIcGnc6fZXftrlh90Kjosc+2//xbwqDgb0rcssVpTu17ZfuK8tL2K9ZeFSWl9KjpH0uvxdz+IADyMw6yeaGxhQQjbyzeSsbUyi0dG1jQ7MQ2XF0ktvHlsbdC419e5Q0C4cS3tenYkeIVJTSq6y0xRDptVN7CaXtFBq2gwPIzN6zzQ2NrK3fyuqNW1izcQtr6reyuj63vHrjFtbWb2F1/dad1jduadytfVeWlewUAJVlpRzcr+f2AOlV3oNe5SU7jjYqmo42dhx5NG1Xsk9xHyF0dQ4gs25uS8O2FBhbUqBsZU1TmNTvCJg129vfHSav/ODUdh1ToV4031v9+vVjzZo17ba/9joF17dvX1avXt0u+3ovMgkgSf8AXAoMB0ZHRG2qlwE/B6qBbcBXI+Lh1DYKuBnYF7gntYWkfsDtwFDgZeCzEbFGuf8yPwY+BdQDkyPisbSvScA30nC+ExGzOvYZm3WO/DBpCpLtRyTN11OfDZsbWt3ffuWl9K0so2/PHvSrLOOQ/r22r/etLKNfz7Lc+r+sp29lD/rsW0ZZaeunqiR12XDZHWvWrCnI55/VtaSsjoCeAv6eXNjkmwoQEUdIOgi4V9LREbEN+FlqX0gugMYB9wIXAQ9GxPclXZTWLwROAQ5NtzFp+zEpsC4hF3IBLJY0LyLa788Ss3awtXFbOhLZ2iw0cmGy03r9FtZu3Mo7uwiTXuWl9K3sQb+eZfTpWcYH+veib88y+lX2oE/PMvpVlqX1XMD06ZkLk0L8qz2rv9j3VlzSGy7dP+thvEtc0juTx80kgCJiKbT4P+LhwILU501Ja4FqSa8CvSOiJm13C3AGuQA6HTg2bT8LeJhcAJ0O3BK5PzdqJPWRNCj1nR8Rq9O+5pMLs7nt/0zNcrY2bmNts9NYzU9rrcm7brJm45ZdhkllWUnuCCSFxrADK3c+IulZlgubVNu/Zw/KS0v2aOyF+Fd7sc7+0jfXF9xrCenI9NLOf9xCuwb0J+A0SXOBg4FR6X4bsDyv33JgcFoeEBEr0/LrwIC0PBh4tYVtWqu/i6RpwDSAIUOG7Nkzsi6noXEba+rTEUhTcLRx3eSdTbsOk+1HIJVlDDug507r/XrmnfKqLKPPXoSJWSHpsACS9AAwsIWm6RFxVyub3UjuulAt8ArwB2D3ps4A6ZpQu/15EREzgBkA1dXVhfdni3WYiOCO2ldZ9NKaZtdRtrB+F2HSs6xk+2msPj178P4DeuaORtKprqZA6ZPXp6KHw8S6pw4LoIg4cQ+2aQC+1rQu6Q/Ac8AaoCqvaxWwIi2/IWlQRKxMp9jeTPUV5I6emm+zgh2n7JrqD7/XsVoBa4dz7AI+l27vUtHGxpvS7e0W2i5dt3cDy0AhXrfI6ppFeyjE04d9+/bN5HEL6hScpJ6AImKjpJOAhoh4JrWtlzSW3CSEicDVabN5wCTg++n+rrz6eZJuIzcJYV0KqfuA70lqesVPBi7uhKdnnWUv/5H/0QPP8aMHnuez1VV8/++PZJ/u/lmSIgzNQlWI13+ylNU07E+TC5D+wK8lLYmITwIHAfdJ2kbuSOWLeZt9hR3TsO9NN8gFzx2SppA7bffZVL+H3BTsOnLTsM8GiIjVkr4NPJr6fatpQoLZTx58nh898DxnjnL4mHU0OZF3T3V1ddTW1mY9DOtA1yx4nivuf47PHFXFZWce6U/Rm7UDSYsjorqlNn+5kRlw7UN1XHH/c/z9xwY7fMw6iQPIur2fPfwCl9/3LGd89H1c/g8fcfiYdRIHkHVrM377Aj/4zTJO+8j7uMLhY9apHEDWbd3wuxf53j3LOPXIQVz12Y/46/bNOpnfcdYtzXzkJb7z66X83RGD+NHnPurwMcuA33XW7dz0+5f49t3P8KkjBvKj8Q4fs6z4nWfdyqw/vMw3/+cZxo0YyI/Hf6zdfhrZzN47v/us27jljy9zybynOfnwAfxkgsPHLGt+B1q38IuaV/iPu57mxOEDuObzR+3yR9PMrHP4XWhd3pyFf+Yb//0UJw4/iJ+e5fAxKxR+J1qXdtuiP/Pvdz7J8R8+iGsdPmYFxe9G67LuePRVLvrVkxz7of787AtH+UfczAqMA8i6pDtqX+XCXz3B3xzWn+u+MMrhY1aAHEDW5fxy8XIu/K8nOOaQA5nxxVH+xVGzAuUAsi7lV48t5//+8k984oMHcv3EaoePWQFzAFmX8d+Pr+Df/vNPfPyDBzh8zIqAA8i6hLuWrODrdyxhzLADuGHi0exb5vAxK3QOICt6//On1/ja7Us4emg/Zk6udviYFQkHkBW1u594jX+5fQnVQ/tx09lH07OsNOshmdlucgBZ0brnyZV89bYlHDWkDzdNdviYFRsHkBWl3zy1kn+e+zgfO7gPN509mspyh49ZsckkgCRdLmmZpCck3SmpT17bxZLqJD0r6ZN59XGpVifporz6MEkLU/12SWWpXp7W61L70LYew4rDfU+/znlzHufIqv25+ZzR9HL4mBWlrI6A5gMjI+JI4DngYgBJhwPjgRHAOOCnkkoklQDXAqcAhwMTUl+AHwA/jIhDgDXAlFSfAqxJ9R+mfq0+Rgc/X2sn9z/9OufOfowjqvZnlsPHrKhlEkARcX9ENKTVGqAqLZ8O3BYRmyPiJaAOGJ1udRHxYkRsAW4DTpck4Hjgl2n7WcAZefualZZ/CZyQ+rf2GFbgHnjmDc6d8xgjBufCZ7+KHlkPycz2QiFcAzoHuDctDwZezWtbnmqt1Q8A1uaFWVN9p32l9nWpf2v7sgL24NI3+KfZizl8UG9uOWc0vR0+ZkWvw85fSHoAGNhC0/SIuCv1mQ40ALM7ahx7Q9I0YBrAkCFDMh5N9/XQsjf5p188xvBBvbllyhj239fhY9YVdFgARcSJu2qXNBk4FTghIiKVVwAH53WrSjVaqb8N9JFUmo5y8vs37Wu5pFJg/9R/V4/R/DnMAGYAVFdXR0t9rGM9/Oyb/OMvFnPYwF7ceo7Dx6wryWoW3DjgAuC0iKjPa5oHjE8z2IYBhwKLgEeBQ9OMtzJykwjmpeB6CDgzbT8JuCtvX5PS8pnAgtS/tcewAvPb51Yx7dbFHNK/F7+YMob9ezp8zLqSrKYQXQOUA/Nz8wKoiYgvR8TTku4AniF3au7ciGgEkHQecB9QAtwYEU+nfV0I3CbpO8DjwMxUnwncKqkOWE0utNjVY1jh+N3zq5h6Sy0f7N+L2V8aQ5+eZVkPyczamXac/bJdqa6ujtra2qyH0S38vu4tzrn5UYYdWMmcqWPpV+nwMStWkhZHRHVLbYUwC85suz/UvcWUWQ4fs+7AAWQF448vvM05sx5lSL+ezP7SGIePWRfnALKCsPDFtznn5kc5uG9P5kwdywG9yrMekpl1MAeQZW7RS6s5++ZHGdx3X+ZMHcuBDh+zbsEBZJl69OXVTL5pEYP2r2DO1DH038/hY9ZdOIAsM4tfWc3kGxcxsHcFc6eO5aD9KrIekpl1IgeQZeKxP69h0o2PclDvCuZOG8tBvR0+Zt2NA8g63eN/XsOkmYs4sFcZc6eOZYDDx6xbcgBZp1ry6lomzlxEv15lzJ02loH7O3zMuisHkHWaJ5av5YszF9K3MnfkM2j/fbMekpllyAFkneLJ5ev4wg0L6dOzB3OnjeV9fRw+Zt2dA8g63FMr1vGFmQvpvW8P5k4dy2CHj5nhALIO9vRr6zjrhoX0Ki9l7tSxVPXtmfWQzKxAOICswzzz2nrOumEhlWUlzJ06loP7OXzMbAcHkHWIpSvXc9YNNezbo4S508Yy5ACHj5ntzAFk7e7Z19/hrBsWUl6aO/J5/wGVWQ/JzAqQA8ja1XNvvMPnr6+hR4mYO20sQw90+JhZyxxA1m6eT+FTso+YO3Uswxw+ZrYLDiBrF3VvbmDC9QuRxJypY/lA/15ZD8nMCpwDyPbaC6s2MOH6GgDmTh3DIQc5fMysbQ4g2ysvrtrAhBk1REQKn/2yHpKZFQkHkO2xl97ayITra2jcFsyZOpZDBzh8zGz3ZRJAki6XtEzSE5LulNQn1Q+Q9JCkDZKuabbNKElPSqqT9BNJSvV+kuZLej7d9011pX516XGOytvXpNT/eUmTOvGpdxkvv7WRCTNq2NqYC5/DHD5m9h5ldQQ0HxgZEUcCzwEXp/om4P8B/9bCNj8DpgKHptu4VL8IeDAiDgUeTOsAp+T1nZa2R1I/4BJgDDAauKQptGz3vPJ27shnc0Mjs780hg8NdPiY2XuXSQBFxP0R0ZBWa4CqVN8YEY+QC6LtJA0CekdETUQEcAtwRmo+HZiVlmc1q98SOTVAn7SfTwLzI2J1RKwhF4ZNYWZt+PPb9UyYUcNftjYy+0tjGT6od9ZDMrMiVQjXgM4B7m2jz2Bged768lQDGBARK9Py68CAvG1ebWGb1urWhldX1zPh+ho2bskd+Rz+PoePme250o7asaQHgIEtNE2PiLtSn+lAAzC7PR4zIkJStMe+ACRNI3f6jiFDhrTXbovS8jX1jJ9Rw4bNDcz+0hhGvG//rIdkZkWuwwIoIk7cVbukycCpwAnptNqurCCdpkuqUg3gDUmDImJlOsX2Zt42B7ewzQrg2Gb1h1t5DjOAGQDV1dXtFmzFZsXavzB+Rg3vbNrK7C+NZeRgh4+Z7b2sZsGNAy4ATouI+rb6p1Ns6yWNTbPfJgJ3peZ5QNNMtknN6hPTbLixwLq0n/uAkyX1TZMPTk41a8Fra//C+Bl/ZN1ftnLrlDEcUeXwMbP20WFHQG24BigH5qfZ1DUR8WUASS8DvYEySWcAJ0fEM8BXgJuBfcldM2q6bvR94A5JU4BXgM+m+j3Ap4A6oB44GyAiVkv6NvBo6vetiFjdUU+0mK1clzvyWbtxK7d+aQwfObhP1kMysy5EbZ/9Msidgqutrc16GJ3m9XWbGD/jj7y1YQu3ThnNx4Z4prqZvXeSFkdEdUtthTALzgrMG+s3MeH6Gt7asIVZ5zh8zKxjOIBsJ2+u38SEGTW8uX4Ts845mlHvd/iYWcfI6hqQFaA339nE+OtreH39JmadM5pR7++X9ZDMrAvzEZABsOqdzUyYUcPr6zZx89mjOXqow8fMOpYDyHhrw2Y+f30Nr63dxI2Tj2b0MIePmXU8B1A393YKn1fX1HPj5KMZ+4EDsh6SmXUTDqBu7O0NmznrhoW88nY9N046mr/6oMPHzDqPA6ibWr1xC2fdsJCX3trIzElH8/FDDsx6SGbWzTiAuqE1KXxefGsjN0yq5phDHT7WsebOncvIkSMpKSlh5MiRzJ07N+shWQHwNOxuZm19LnxeWLWB6ydW89eH9s96SNbFzZ07l+nTpzNz5kyOOeYYHnnkEaZMmQLAhAkTMh6dZclHQN3IuvqtfGHmQure3MCML47ibw9z+FjH++53v8vMmTM57rjj6NGjB8cddxwzZ87ku9/9btZDK0pd6WjSAdRNrPtLLnyee30DP//iKI790EFZD6mgdaU3edaWLl3K8uXLd3o9ly9fztKlS7MeWtFpOpq8+uqr2bRpE1dffTXTp08v3v8/I8K33biNGjUqitXa+i1x2tW/i0P+/dfxwDOvZz2cgjdnzpwYNmxYLFiwILZs2RILFiyIYcOGxZw5c7IeWlGqqqqKQYMG7fR6Dho0KKqqqrIeWtEZMWJELFiwYKfaggULYsSIERmNqG1AbbTy72qbR0CSDpL0TUm/TLdvShrQ1nZWGNZv2srEGxfxzMr1/OysUZww3P/p2uJTRu0vmn3rfvN12z1Lly7lmGOO2al2zDHHFO3R5C4DSNIn2PG7ObekG8DC1GYF7J1NW5k4cxFPr1jHtZ8/ihMPd/jsjq72Js/aa6+9xmWXXcb5559PRUUF559/PpdddhmvvfZa1kMrOsOHD+eRRx7ZqfbII48wfPjwjEa0d9o6AroSOCMiLomIeel2CXAGcFWHj8722IbNDUy6cRFPrVjHNZ8/ipNHDMx6SEWjq73JszZ8+HCqqqp46qmnaGxs5KmnnqKqqsqv5x6YPn06U6ZM4aGHHmLr1q089NBDTJkyhenTp2c9tD3T2rm5dIj8zJ60dcVbMV0DemfT1vjMT38fH7j413Hvk69lPZyi42tA7cuvZ/uaM2dOjBgxIvbZZ58YMWJEwb+O7OIaUFufA5KkvhGxplmxH55BV5A2bm7g7JsW8fira7l6wscYN3JQ1kMqOhMmTGDixIkcf/zx22ulpaX+zMoemjBhAjfffDMnnHACEYEkTjrpJL+e1maI/BC4X9LfStov3Y4F7k1tVkBy4fMoj/15LT8e/1E+dYTDZ09UVFTQ0NDAgAEDWLp0KQMGDKChoYGKioqsh1aUzj//fBYsWMAVV1zBxo0bueKKK1iwYAHnn39+1kMrOt1uGjZwKvBb4O10+y3wf9rarqvdiuEU3NdufzyGXXR33LVkRdZDKWpADBgwYKfagAEDIvd2sfeqvLw8rrzyyp1qV155ZZSXl2c0ouLV1aZhK9dubamuro7a2tqO2fml+3fMftvDpeuyHkGnk8TSpUv58Ic/vL22bNkyhg8fjt8v750kNm7cSM+ePbfX6uvrqays9Ov5HpWUlLBp0yZ69OixvbZ161YqKipobGzMcGStk7Q4IqpbatvlNSBJBwLnAquBm4DLgL8BXgD+NSLq2nms3ZK+ub4g34iSiEuzHkU2jj32WF5//fWd1m3PlJeXc9111/H1r399e+26666jvLw8w1EVp6YZmscdd9z2WjHP0GzrGtAcoBw4DFgEvAycCdwN3LCnDyrpcknLJD0h6U5JfVL9JEmLJT2Z7o/P22ZUqtdJ+okkpXo/SfMlPZ/u+6a6Ur+69DhH5e1rUur/vKRJe/o8rGsqLy/njTfeYODAgSxbtoyBAwfyxhtv+B/MPTR16lQuvPBCrrrqKurr67nqqqu48MILmTp1atZDKzrdbRr2n9K9gD83a1uyq23b2O/JQGla/gHwg7T8MeB9aXkksCJvm0XA2DSWe4FTUv0y4KK0fFHevj6V+ilttzDV+wEvpvu+ablvW2PuyGtAFOi1hUIdV2coLy8PYPvN1yv2znnnnbf9NS0vL4/zzjsv6yEVra40DXuX14AkPRYRRzVfbml9T0n6NHBmRJzVrC5ykx4GpbB4KCI+nNomAMdGxD9KejYtr5Q0CHg4Ij4k6edpeW7a5lng2KZbRPxjqu/UrzUdeQ1IUuGegivAcZlZ8djja0DAByTNI3cU0bRMWh/WTuM7B7i9hfpngMciYrOkwcDyvLblwOC0PCAiVqbl14Gm75sZDLzawjat1d9F0jRgGsCQIUN29/mYmdluaCuATs9bvqJZW/P1nUh6AGjp+1+mR8Rdqc90oAGY3WzbEeROzZ3cxvh2EhEhqd3+ZI+IGcAMyB0Btdd+zcysjQCKiP9trU3S7UCr7RFx4q72LWkyuc8YnRB553kkVQF3AhMj4oVUXgFU5W1elWoAb0galHcK7s28bQ5uYZsV5E7D5dcf3tVYzcys/e3N1+n81Z5uKGkccAFwWkTU59X7AL8mN6ng9031dIptvaSx6drQROCu1DwPaJrJNqlZfWKaDTcWWJf2cx9wsqS+acbcyalmZmadKKvvc7sG2A+YL2mJpOtS/TzgEOA/Un2JpKaf7vwKuanfdeQ+h3Rvqn8fOEnS88CJaR3gHnIz3OqA69P2RMRq4NvkfmbiUeBbqWZmZp2orVlwrc1yE3B3RHSbLxvzLDgzs/dub2bBXbmLtmV7PiQzM+vu2pqEcNyu2s3MzPZUWz/JfUHe8j80a/teRw3KzMy6vrYmIYzPW764Wdu4dh6LmZl1I20FkFpZbmndzMxst7UVQNHKckvrZmZmu62tWXAfkbSe3NHOvmmZtO7fJzYzsz3W1iy4ks4aiJmZdS9ZfROCmZl1cw4gMzPLhAPIzMwy4QAyM7NMOIDMzCwTDiAzM8uEA8jMzDLhADIzs0w4gMzMLBMOIDMzy4QDyMzMMuEAMjOzTDiAzMwsEw4gMzPLRCYBJOlyScskPSHpTkl9Un20pCXp9idJn87bZpykZyXVSboorz5M0sJUv11SWaqXp/W61D40b5uLU/1ZSZ/svGduZmZNsjoCmg+MjIgjgeeAi1P9KaA6Ij4KjAN+LqlUUglwLXAKcDgwQdLhaZsfAD+MiEOANcCUVJ8CrEn1H6Z+pO3GAyPSY/w07d/MzDpRJgEUEfdHRENarQGqUr0+r17Bjp/9Hg3URcSLEbEFuA04XZKA44Ffpn6zgDPS8ulpndR+Qup/OnBbRGyOiJeAurR/MzPrRIVwDegc4N6mFUljJD0NPAl8OQXSYODVvG2Wp9oBwNq80Gqqk79Nal+X+re2r3eRNE1SraTaVatW7dWTNDOznXVYAEl6QNJTLdxOz+szHWgAZjfVImJhRIwAjgYullTRUWNsS0TMiIjqiKju379/VsMwM+uSSjtqxxFx4q7aJU0GTgVOiIho3h4RSyVtAEYCK4CD85qrUu1toI+k0nSU01Qnb5vlkkqB/VP/1vZlZmadKKtZcOOAC4DTIqI+rz4shQWS3g98GHgZeBQ4NLWXkZtEMC8F10PAmWkXk4C70vK8tE5qX5D6zwPGp1lyw4BDgUUd9mTNzKxFHXYE1IZrgHJgfm5eADUR8WXgGOAiSVuBbcBXIuItAEnnAfcBJcCNEfF02teFwG2SvgM8DsxM9ZnArZLqgNXkQouIeFrSHcAz5E7/nRsRjR39hM3MbGdq4eyXtaC6ujpqa2s7ZN+SKMT/DoU6LjMrHpIWR0R1S22FMAvOzMy6IQeQmZllwgFkZmaZcACZmVkmHEBmZpYJB5CZmWXCAWRmZpnI6oOo1kz6QG5B6du3b9ZDMLMuzAFUANrzw57+8KiZFQufgjMzs0w4gMzMLBMOIDMzy4QDyMzMMuEAMjOzTDiAzMwsEw4gMzPLhAPIzMwy4QAyM7NMOIDMzCwTDiAzM8uEA8jMzDLhADIzs0xkEkCSLpe0TNITku6U1KdZ+xBJGyT9W15tnKRnJdVJuiivPkzSwlS/XVJZqpen9brUPjRvm4tT/VlJn+z4Z2xmZs1ldQQ0HxgZEUcCzwEXN2u/Cri3aUVSCXAtcApwODBB0uGp+QfADyPiEGANMCXVpwBrUv2HqR9pu/HACGAc8NO0fzMz60SZBFBE3B8RDWm1BqhqapN0BvAS8HTeJqOBuoh4MSK2ALcBpyv3K27HA79M/WYBZ6Tl09M6qf2E1P904LaI2BwRLwF1af9mZtaJCuEa0Dmkox1JvYALgW826zMYeDVvfXmqHQCszQuzpvpO26T2dal/a/t6F0nTJNVKql21atUePTkzM2tZhwWQpAckPdXC7fS8PtOBBmB2Kl1K7nTaho4a13sRETMiojoiqvv375/1cMzMupQO+0nuiDhxV+2SJgOnAifEjt+QHgOcKekyoA+wTdImYDFwcN7mVcAK4G2gj6TSdJTTVCfdHwwsl1QK7J/6r2hlX2Zm1omymgU3DrgAOC0i6pvqEfHXETE0IoYCPwK+FxHXAI8Ch6YZb2XkJhHMS8H1EHBm2sUk4K60PC+tk9oXpP7zgPFpltww4FBgUcc9WzMza0mHHQG14RqgHJifmxdATUR8ubXOEdEg6TzgPqAEuDEimiYpXAjcJuk7wOPAzFSfCdwqqQ5YTS60iIinJd0BPEPu9N+5EdHY3k/QzMx2TTvOftmuVFdXR21tbdbDaJMk/N/UzAqFpMURUd1SWyHMgjMzs27IAWRmZplwAJmZWSYcQGZmlgkHkJmZZcIBZGZmmXAAmZlZJhxAZmaWCQeQmZllwgFkZmaZcACZmVkmHEBmZpYJB5CZmWXCAWRmZplwAJmZWSYcQGZmlgkHkJmZZcIBZGZmmXAAmZlZJhxAZmaWCQeQmZllwgFkZmaZyCSAJF0uaZmkJyTdKalPqg+V9BdJS9LturxtRkl6UlKdpJ9IUqr3kzRf0vPpvm+qK/WrS49zVN6+JqX+z0ua1MlP38zMyO4IaD4wMiKOBJ4DLs5reyEiPppuX86r/wyYChyabuNS/SLgwYg4FHgwrQOcktd3WtoeSf2AS4AxwGjgkqbQMjOzzpNJAEXE/RHRkFZrgKpd9Zc0COgdETUREcAtwBmp+XRgVlqe1ax+S+TUAH3Sfj4JzI+I1RGxhlwYNoWZmZl1kkK4BnQOcG/e+jBJj0v6X0l/nWqDgeV5fZanGsCAiFiZll8HBuRt82oL27RWfxdJ0yTVSqpdtWrVe3xaZma2K6UdtWNJDwADW2iaHhF3pT7TgQZgdmpbCQyJiLcljQL+W9KI3X3MiAhJsZdDz9/fDGAGQHV1dbvt18zMOjCAIuLEXbVLmgycCpyQTqsREZuBzWl5saQXgMOAFex8mq4q1QDekDQoIlamU2xvpvoK4OAWtlkBHNus/vB7fHpmZraXspoFNw64ADgtIurz6v0llaTlD5CbQPBiOsW2XtLYNPttInBX2mwe0DSTbVKz+sQ0G24ssC7t5z7gZEl90+SDk1PNzMw6UYcdAbXhGqAcmJ9mU9ekGW9/A3xL0lZgG/DliFidtvkKcDOwL7lrRk3Xjb4P3CFpCvAK8NlUvwf4FFAH1ANnA0TEaknfBh5N/b6V9xgFLb1W7dYvHXiamWVC/kdo91RXV0dtbW3WwzAzKyqSFkdEdUtthTALzszMuiEHkJmZZcIBZGZmmXAAmZlZJhxAZmaWCQeQmZllwgFkZmaZcACZmVkm/EHU3SRpFblvWih0BwJvZT2ILsSvZ/vy69l+iuW1fH9E9G+pwQHUxUiqbe1Tx/be+fVsX349209XeC19Cs7MzDLhADIzs0w4gLqeGVkPoIvx69m+/Hq2n6J/LX0NyMzMMuEjIDMzy4QDyMzMMuEA6iIkHSzpIUnPSHpa0lezHlOxklQhaZGkP6XX8ptZj6krkFQi6XFJd2c9lmIn6WVJT0paIqlofykzq5/ktvbXAPxrRDwmaT9gsaT5EfFM1gMrQpuB4yNig6QewCOS7o2ImqwHVuS+CiwFemc9kC7iuIgohg+itspHQF1ERKyMiMfS8jvk3uiDsx1VcYqcDWm1R7p5ts5ekFQF/B1wQ9ZjscLhAOqCJA0FPgYszHgoRSudLloCvAnMjwi/lnvnR8AFwLaMx9FVBHC/pMWSpmU9mD3lAOpiJPUC/gv4l4hYn/V4ilVENEbER4EqYLSkkRkPqWhJOhV4MyIWZz2WLuSYiDgKOAU4V9LfZD2gPeEA6kLS9Yr/AmZHxK+yHk9XEBFrgYeAcRkPpZh9AjhN0svAbcDxkn6R7ZCKW0SsSPdvAncCo7Md0Z5xAHURkgTMBJZGxFVZj6eYSeovqU9a3hc4CViW6aCKWERcHBFVETEUGA8siIgvZDysoiWpMk00QlIlcDLwVLaj2jOeBdd1fAL4IvBkunYB8O8RcU92Qypag4BZkkrI/ZF2R0R46rAVigHAnbm/OSkF5kTEb7Id0p7xV/GYmVkmfArOzMwy4QAyM7NMOIDMzCwTDiAzM8uEA8jMzDLhADIzs0w4gMyKjKShkoryg4dm+RxAZmaWCQeQWSdLRzBLJV2ffvDu/vSVPy31PUTSA+nH8R6T9MFm7RWSbko/Tva4pONSfbKkX0n6jaTnJV2Wt83Jkv6Y9vef6QtskfT99IOGT0i6oiNfAzNwAJll5VDg2ogYAawFPtNKv9mp30eAjwMrm7WfS+4njI4AJpD7CqGK1PZR4HPAEcDn0q/mHgh8AzgxfZtyLfB1SQcAnwZGRMSRwHfa52matc7fBWeWjZciYklaXgwMbd4hfeHk4Ii4EyAiNqV6frdjgKtT+zJJrwCHpbYHI2Jd2uYZ4P1AH+Bw4PdpP2XAH4F1wCZgZvrJbH/3nXU4B5BZNjbnLTcCLZ6Ca+fHKAVE7gf2JjTvLGk0cAJwJnAecHwHjMlsO5+CMytQ6afVl0s6A0BSuaSezbr9DjgrtR8GDAGe3cVua4BPSDokbVMp6bB0HWj/9O3pXwM+0q5PxqwFDiCzwvZF4J8lPQH8ARjYrP2nwD6SngRuByZHxGZaERGrgMnA3LTPPwIfBvYD7k61R4Cvt/cTMWvOP8dgZmaZ8BGQmZllwpMQzAqApGvJ/aptvh9HxE1ZjMesM/gUnJmZZcKn4MzMLBMOIDMzy4QDyMzMMuEAMjOzTPx/lIjJ51LrckMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(1, len(n_clone_list)+1), np.max(_ELBO_mat, axis=1))\n",
    "plt.boxplot(_ELBO_mat)\n",
    "plt.xticks(np.arange(1, len(n_clone_list)+1), n_clone_list)\n",
    "plt.ylabel(\"ELBO\")\n",
    "plt.xlabel(\"n_clones\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to visualise the raw allele frequency with annotation of cells, you may consider [seaborn.clustermap](https://seaborn.pydata.org/generated/seaborn.clustermap.html). We also wrap this function here as `vireoSNP.plot.anno_heat` for quick use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mtSNP_ids = ['mt_variant%d' %x for x in range(AD.shape[0])]\n",
    "cell_label = np.array(['clone1'] * 27 + ['clone2'] * 27 + ['clone3'] * 27)\n",
    "id_uniq = ['clone1', 'clone2', 'clone3']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.matrix.ClusterGrid at 0x7fd475067eb0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwcAAALICAYAAAAqkYyaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8bUlEQVR4nO3dfbxdVX0n/s83CTEoJIhgdAwYCyJCAB2UWgZGWnXEqqhtUQq0xgeCMgytj2W0VrTaHxVtVdDQaDWt1opRZohKpYrYEfzJLzCGJ3mQJwUKCqhAUBKSrN8fd2NPb8/JfUxucvN+v17nxdlr773W2vueG+7nrLX3rtZaAAAAZkx1BwAAgK2DcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkSWZNdQeAbc/ll1/++FmzZn0yyaL4kgG2FhuTXL1+/frXH3zwwT+Z6s4A2ybhABizWbNmffIJT3jC03ffffefzZgxo011f4Bk48aNdffdd+931113fTLJUVPdH2Db5Bs/YDwW7b777vcLBrD1mDFjRtt9993vy9CIHsC4CAfAeMwQDGDr0/1e+n87MG7+AQEAAJK45gCYBEcuve7gyazva2/c9/Lx7PfmN7/5P+20004b3vve9/54Mvtz+OGHP3X16tWPedaznrXmoosuunEy637EQ58+elLP4ZzXrNhqzuF3vvOdHU866aQnr1mzZuaMGTPa2972tjtPOOGEn01W/b3e8fmZk3oe/+KYDWM+j5vjHN5www2zX/7yl++1cePGWr9+fS1ZsuQnb3/72++erPoBHiEcAIzgrW99610PPvjgjE984hO7T3VftkU77bTTxs985jO3HHDAAWtvvfXWHZ797Gc//RWveMX9u+2224ap7tu2Ys8993z48ssvv27HHXds991334z99ttv/1e+8pU/X7hw4cNT3TdgejGtCNhmnXXWWY/bZ5999nva056238tf/vKn9K77zne+s+NBBx207z777LPfC17wgr3uvvvumUlyyCGHPO2Nb3zjkw444ICnL1y4cNHXvva1nZJk/fr1OfHEExcsWrTo6fvss89+Z5xxxm6P1PWyl73sgblz527cske3ZWyJc3jggQeuPeCAA9YmycKFCx/edddd1995553T5supLXEO58yZ03bccceWJL/85S9r48Zp+XEEtgLCAbBNuuyyy+Z88IMffOK//Mu/3HD99dd//2/+5m9+1Lt+8eLFT/mLv/iL22+44Ybv77///r/8kz/5k//0yLr169fXVVddde1f/uVf3vbe9773PyXJhz/84d3mzZu34eqrr772iiuuuPbv/u7vdr/uuutmb+nj2pKm4hxedNFFj3744Ydrv/32W7tljnLz2pLn8MYbb9xhn3322e8pT3nKgaeccspdRg2AzUE4ALZJF1xwwdyXvvSlP3viE5+4Pknmz5//qykq995778wHHnhg5otf/OI1SXLCCSfc+93vfnenR9YfffTRP0uSQw899MHbb799dpJ84xvfmPuFL3zhcfvuu+9+z3zmM5/+s5/9bNb3v//9OVv2qLasLX0Of/jDH+7wmte85tc+8YlP3Dpz5swtdZib1ZY8h3vvvffDN9xww/evvfbaqz/3uc/tdtttt02b0Rdg6+EfFmC7M2fOnJYks2bNyoYNGypJWmv1oQ996Ee/+7u/e//U9m7bMNZz+NOf/nTGi170or3f/e533/G85z3vwS3d363ReD+HCxcufHjffff95Te+8Y2dX/Oa12yWC7uB7ZeRA2Cb9MIXvvD+L3/5y4+96667ZibJj3/84199Ff24xz1uw9y5czc8Mo/7b//2bx/3G7/xG2s2Vd8LXvCC+5YuXbr72rVrK0muvPLKR91///3T+t/ILXUOH3rooXrxi1+89zHHHHPvdPtjdkudw5tuummHNWvWVJLcfffdM1etWrXT/vvv/9DmOzJge2XkAJiw8d56dCKe9axnPfSWt7zlzsMPP3zfGTNmtEWLFv3iyU9+8rpH1n/605++5Y1vfOOTTznllBl77rnn2n/8x3+8dVP1velNb7rn1ltvfdQBBxzw9NZa7brrrg+ff/75NyXJwQcf/LSbb755zi9/+cuZ8+fPP/DjH//4rZM9wjDeW49OxJY6h5/97Gd3WbVq1U4/+9nPZn3uc5/bLUk+9alP3XLooYf+crKPaTy3Hp2ILXUOr7zyyh3/5E/+ZEFVpbWWk08++a5DDjlk0s8fQLXmIafA2FxxxRW3HnTQQfdMdT+A/+iKK67Y7aCDDlo41f0Atk3TesgcAAAYPeEAAABIIhwA47Nx48aNNdWdAP697vfSE9KAcRMOgPG4+u67754nIMDWY+PGjXX33XfPS3L1VPcF2Ha5WxEwZuvXr3/9XXfd9cm77rprUXzJAFuLjUmuXr9+/eunuiPAtsvdigAAgCS+8QMAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6s6a6A2ydjlx6XZvqPsAg/3vOu6a6CzDQe3c8d6q7AAP9xTEbaqr7wNbNyAEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrVWpvqPrB18sGAPtatu6dv+ezZu23hnmy7nMPR217P1WQe99Z4Dqe4T7UlGmHbZeQAAABIIhwAAAAd4QAAAEgiHAAAAB3hAAAASCIcAAAAHeEAAABIIhwAAAAd4QAAAEgiHAAAAB3hAAAASCIcAAAAHeEAAABIIhwAAAAd4QAAAEgiHAAAAB3hAAAASCIcAAAAHeEAAABIIhwAAAAd4QAAAEgiHAAAAB3hAAAASCIcAAAAHeEAAABIIhwAAAAd4QAAAEgiHAAAMAZVdURVHbqZ2ziqqk4d5767VNVJw8peXVU/6F6vnpxeTk/CAQAAY3FEks0WDqpqVmttZWvt9HFWsUuSX4WDqto1ybuT/HqSQ5K8u6oeO+GOTlPCAQDAdqKqllTVZT2vJcPWL6yq66pqeVXdUFX/UFXPr6pLum/dD0nyhiRvqqrVVXV4nzbmVdUPq2pGt/yYqrqtqnaoqhOqalVVXVFVX6qqR3fbLK+qs6vq0iQfqKrFVXVWt+6lVXVpVX2vqr5RVfO78tOq6lNV9a2qurmqTum6cHqSvbr+nZHkhUm+3lr7aWvtZ0m+nuTIzXKCpwHhAABgO9FaW9Zae1bPa1mfzfZO8qEk+3avY5McluStSd6R5Owkf91ae0Zr7dt92rgvyeokz+2KXpLkgtbaw0nOba09u7V2UJJrk7yuZ9cFSQ5trb15WJUXJ3lOa+2ZST6f5O096/bN0B//j4wI7JDk1CQ3df17W5InJbmtZ5/buzL6mDXVHQAAYKtyS2vtqiSpqmuSXNhaa1V1VZKFGfrDfyTnJHlVkouSHJPk4135oqp6X4am/uyU5IKefVa01jb0qWtBknOq6olJZie5pWfdV1tra5OsraqfJJk/qiNkICMHAAD0WtvzfmPP8saM/ovllUmO7Ob7H5zkm1358iQnt9YOSPKeJHN69nlwQF1nJjmr2+fEYfv09nXDgP7dkWSPnuUFXRl9CAcAAIzFA0l23tQGrbU1SVYl+UiSr/SMCOyc5M5u+s9xo2xvXv7tj/nR3GloeP8uSPLfquqx3YXI/y3/fsSCHsIBAABj8eUkrxh0QXKPc5Ic3/33Ee9KcmmSS5JcN8r2TkuyoqouT3LPSBu31u5NcklVXV1VZ7TWfprkzzMUVlYleW9XRh/VWpvqPrB18sGAPtat6///pdmzd9vCPdl2OYejt72eq8k87q3xHE5xn2pLNMK2y8gBAACQxN2KAAAYp6p6Z5KjhxWvaK29fyr6w8QJBwAAjEsXAgSBacS0IgAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdKq1NtV9YCt0//3X9P1g3HbxxX233/+3T+xbvm7dPZto44q+5bvt9ry+5ffcc2Hf8ofuu69v+YK9fmdg2zet+mzf8ju/dlPf8gP/6Pf6lg86H0nys8vv6lt+2LvePXCffm6/6dyB6+67/u6+5Xscdljf8rlz9+9bfvGfv6dv+cJjDxjY9qbObz+Dzvlezz5+4D6DPj+zZ+82prbZtt1//zUD1w36TAMD1VR3gK2bkQMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHRmTXUHAABgW/bkqjbZdf6wtZrsOkfDyAEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQxN2KAABgQmbPmTPVXZg0Rg4AABi1qjqiqg7dzG0cVVWnjnPfXarqpGFlX6uqn1fVVyanh9OXcAAAwFgckWSzhYOqmtVaW9laO32cVeyS5KRhZWck+YMJdWw7YVoRAMB2oqqWJFnSU7SstbasZ/3CJF9L8t0MBYBVST6d5D1JHp/kuCRvSLKhqo5P8j9aa98e1sa8JFcmeUprbWNVPSbJdUl+Lcnirv3ZSW5M8gettV9U1fIkDyV5ZpJLqurKJM9qrZ1cVS9N8qfdPvcmOa619uOqOi3Jnl29eyb5cGvto0lOT7JXVa1O8vXW2ttaaxdW1RETO3uDPWrHHTdX1VuccAAAsJ3ogsCyETbbO8nRSV6boXBwbJLDkhyV5B1Jzk6yprX2wQFt3Nf9Yf7cJBcleUmSC1prD1fVua21TyRJVb0vyeuSnNntuiDJoa21DVW1uKfKi5M8p7XWqur1Sd6e5C3dun2T/GaSnZNcX1VLk5yaZFFr7RkjnxGGEw4AAOh1S2vtqiSpqmuSXNj9YX5VkoVJVo+ijnOSvCpD4eCYJB/vyhd1oWCXJDsluaBnnxWttQ196lqQ5JyqemKGRg9u6Vn31dba2iRrq+onSeaP6ggZyDUHAAD0WtvzfmPP8saM/ovllUmOrKpdkxyc5Jtd+fIkJ7fWDsjQVKXe2/w8OKCuM5Oc1e1z4rB9evu6YQz9YwDhAACAsXggQ9N4BmqtrcnQlKSPJPlKz4jAzknurKodMnT9wmjMS3JH9/7Vk9E/BhMOAAAYiy8neUVVra6qwzex3TlJju/++4h3Jbk0ySUZukh5NE5LsqKqLk9yz0gbt9buzdBFzVdX1RlJUlXfTrIiyfOq6vaqeuEo296qVdWRVXV9Vd3Y79avVbW4qu7uflaru2s2NsnQCwAASZLW2q1JFvUsLx6w7sBR1PXFJDWsbGmSpX22XTxseXmGpiCltXZekvP67HPasOXefh87bN2mQsw2qapmJvlYkhckuT3Jqqpa2Vr7/rBNz2mtnTzaeoUDAACYgCl6QvIhSW5srd2cJFX1+SQvSzI8HIyJaUUAAIxLVb2zZ8rKI693TnW/poOqWlJVl/W8lgzb5ElJbutZvr0rG+53q+rKqvpiVe0xUrtGDgAAGJfW2vuTvH+q+zEdjfKZFCP5cpJ/bK2traoTk/xdkt/a1A7CAQAATMAOs6bkT+o7kvSOBCzIv93VKcmvLs5+xCeTfGCkSk0rAgCAbc+qJE+tqqdU1ewMPWxuZe8G3YPjHnFUkmtHqtTIAQAAbGNaa+ur6uQMPWV6ZpJPtdauqar3JrmstbYyySlVdVSS9Ul+mmTxSPUKBwAAsA1qrZ2f5PxhZX/W8/5/JvmfY6nTtCIAACCJcAAAAHSEAwAAIIlrDgAAYEJmTc2tTDcLIwcAAEAS4QAAAOhMnzEQAACYAjvMnDnVXZg0Rg4AAIAkwgEAANARDgAAgCTCAQAA0BEOAACAJO5WBAAAEzLT3YoAAIDpRjgAAACSmFYEAAATMmsaTSsSDujrB1+6oG/5nPmPGVM9//NlfzBw3eFPf3rf8pf/1fP6ll+39OK+5Tsu2Llv+YK9Bvdr96c9s2/5o3Z9dN/ydevu6lu+x2GHDWxjj8GrxmTOvHkD163d/Rd9ywf1N9m/b+kTj+x/subu/rRN9m0sBp3bTfnJbf+nb/mCvX5not1hG3Lhuz85cN0r/vqvt2BPAKY/04oAAIAkwgEAANARDgAAgCTCAQAA0BEOAACAJO5WBAAAEzJrxvT5vn36HAkAADAhwgEAAJDEtCIAAJiQmdPoCclGDgAAgCTCAQAA0BEOAACAJMIBAADQEQ4AAIAkwgEAANBxK1MAAJiAmZ6QDAAATDfCAQAAkMS0IgAAmBDTigAAgGlHOAAAAJIIBwAAQEc4AAAAkggHAABAx92KAABgAtytCAAAmHaEAwAARq2qjqiqQzdzG0dV1anj3HeXqjqpZ/kZVfX/VtU1VXVlVb1q8no6/QgHAACMxRFJNls4qKpZrbWVrbXTx1nFLklO6ln+RZI/bK3tn+TIJB+uql0m1svpyzUHAAAkSapqYZKvJfluhgLAqiSfTvKeJI9PclySNyTZUFXHJ/kfrbVvD6tjXpIrkzyltbaxqh6T5Lokv5ZkcZIlSWYnuTHJH7TWflFVy5M8lOSZSS6pqiuTPKu1dnJVvTTJn3b73JvkuNbaj6vqtCR7dvXumeTDrbWPJjk9yV5VtTrJ11trb3ukb621f62qnyTZPcnPJ+m0ueYAAIBtT1UtqarLel5L+my2d5IPJdm3ex2b5LAkb03yjiRnJ/nr1tozhgeDJGmt3ZdkdZLndkUvSXJBa+3hJOe21p7dWjsoybVJXtez64Ikh7bW3jysyouTPKe19swkn0/y9p51+yZ5YZJDkry7qnZIcmqSm7r+va23oqo6JEMh46bBZ2n7ZuQAAGA70VpblmTZCJvd0lq7Kkmq6pokF7bWWlVdlWRhhv7wH8k5SV6V5KIkxyT5eFe+qKrel6GpPzsluaBnnxWttQ196lqQ5JyqemKG/rC/pWfdV1tra5Os7UYE5g/qULf/Z5K8urW2cRTHsF0ycgAAQK+1Pe839ixvzOi/WF6Z5Miq2jXJwUm+2ZUvT3Jya+2ADE1VmtOzz4MD6jozyVndPicO26e3rxsG9a+q5ib5apJ3tta+O8pj2C4JBwAAjMUDSXbe1AattTUZul7hI0m+0jMisHOSO7vpP8eNsr15Se7o3r96rP2rqtlJ/leSv2+tfXGUbW63hAMAAMbiy0leUVWrq+rwTWx3TpLju/8+4l1JLk1ySYYuUh6N05KsqKrLk9wz0sattXszdFHz1VV1RpJXJvmvSRZ3fV5dVc8YZdvbHdccAACQJGmt3ZpkUc/y4gHrDhxFXV9MUsPKliZZ2mfbxcOWl2doClJaa+clOa/PPqcNW+7t97HDNv/sSP2dCHcrAgAAplRVHVlV11fVjZt6aFxV/W5Vtap61kh1GjkAAGBcquqdSY4eVryitfb+qejP9qSqZib5WJIXJLk9yaqqWtla+/6w7XZO8kcZms41IuEAAIBx6ULAdh8Epmha0SFJbmyt3ZwkVfX5JC9L8v1h2/15kr9M8raMgmlFAACw7XlSktt6lm/vyn6lqv5zkj1aa18dbaXCAQAAbGVG+TTrTe0/I8lfJXnLWPYzrQgAALYyo3ia9R1J9uhZXpB/ex5EMvSsh0VJvlVVSfKEJCur6qjW2mWDKjVyAAAA255VSZ5aVU/pHvR2TIaeTJ0kaa3d11rbrbW2sLW2MMl3k2wyGCTCAQAAbHNaa+uTnJzkgiTXJvlCa+2aqnpvVR013npNKwIAgG1Qa+38JOcPK/uzAdseMZo6hQMAAJiAGZ6QDAAATDfCAQAAkMS0IgAAmJApekLyZjF9jgQAAJgQ4QAAAEgiHAAAAB3hAAAASCIcAAAAHeEAAABI4lamAAAwITOqproLk8bIAQAAkEQ4AAAAOqYVAQDABEynJyRXa22q+8BWaN26e/p+MB566Md9t587d/8xt3H//deMqa516+7pWz6ePg1qe86c+QP3mSyzZ+82pu0HHfemjPWcjPXcbqquQQa1sanzMdZ+bYmf31jbvv/+KwbWNXv2E8ZU15Yw6PjWrbtrzHWN9fgm8+c6Wb8DmzLoZzt37kFj2n5T+4zVeM7hpn7Px1rXVBr0b8lk/Xu/qX+rBrUxyHjO4Vg/b4PMnr3b9JkcvxU56XnPm/Q/qD9+4YVT8rOaPjEHAACYEOEAAABIIhwAAAAd4QAAAEjibkUAADAhM6bR3Yqmz5EAAAATIhwAAABJhAMAAKDjmgMAAJiAGTV9ni1n5AAAAEgiHAAAAB3hAAAASCIcAAAAHeEAAABI4m5FAAAwITM9IRkAAJhuhAMAACCJaUUAADAhHoIGAABMO8IBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkMStTAEAYEJmeEIyAAAw3QgHAABAEuEAAIAxqKojqurQzdzGUVV16jj33aWqTupZfnJV/d+qWl1V11TVGyavp0NmVE36a6q45gAAgLE4IsmaJN/ZHJVX1azW2sokK8dZxS5JTkry8W75ziS/0VpbW1U7Jbm6qla21v514r2dfowcAABsJ6pqSVVd1vNaMmz9wqq6rqqWV9UNVfUPVfX8qrqkqn5QVYckeUOSN3XfxB/ep415VfXDqprRLT+mqm6rqh2q6oSqWlVVV1TVl6rq0d02y6vq7Kq6NMkHqmpxVZ3VrXtpVV1aVd+rqm9U1fyu/LSq+lRVfauqbq6qU7ounJ5kr65/Z7TW1rXW1nbrHhV//26SkQMAgO1Ea21ZkmUjbLZ3kqOTvDbJqiTHJjksyVFJ3pHk7CRrWmsfHNDGfVW1Oslzk1yU5CVJLmitPVxV57bWPpEkVfW+JK9Lcma364Ikh7bWNlTV4p4qL07ynNZaq6rXJ3l7krd06/ZN8ptJdk5yfVUtTXJqkkWttWc8UkFV7ZHkq92xvc2owWCSEwAAvW5prV3VWtuY5JokF7bWWpKrkiwcZR3nJHlV9/6YbjlJFlXVt6vqqiTHJdm/Z58VrbUNfepakOSCbp+3Ddvnq621ta21e5L8JMn8fp1prd3WWjswQ+Hg1Y+MPvAfCQcAAPRa2/N+Y8/yxox+1snKJEdW1a5JDk7yza58eZKTW2sHJHlPkjk9+zw4oK4zk5zV7XPisH16+7phpP51IwZXJ/kP06EYIhwAADAWD2RoGs9ArbU1GZqS9JEkX+kZEdg5yZ1VtUOGRg5GY16SO7r3rx5r/6pqQVXt2L1/bIamSF0/yrZHZcaMGZP+mirCAQAAY/HlJK8YdEFyj3OSHJ9/m1KUJO9KcmmSS5JcN8r2TkuyoqouT3LPSBu31u5NcklVXV1VZyR5epJLq+qKJP+S5IOttatG2fZ2xwXJAAAkSVprtyZZ1LO8eMC6A0dR1xeT1LCypUmW9tl28bDl5RmagpTW2nlJzuuzz2nDlnv7feywzUfsL0OMHAAAAEmMHAAAME5V9c4M3fa014rW2vunoj9TZSqfaDzZhAMAAMalCwHbVRDYmlTVkRm66Htmkk+21k4ftv4NSf57hu7ktCbJktba9zdVp2lFAACwjamqmUk+luRFSfZL8vtVtd+wzT7XWjugeyDcB5L81Uj1CgcAALDtOSTJja21m1tr65J8PsnLejdord3fs/iYJG2kSk0rAgCArUxVLUmypKdoWWttWc/yk5Lc1rN8e5Jf71PPf0/y5iSzk/zWSO0KBwAAsJXpgsCyETccuZ6PJflYVR2b5E8zwoPkhAMAAJiAmVNzt6I7kuzRs7wg//Yk6X4+nz7PmBjONQcAALDtWZXkqVX1lKqaneSYJCt7N6iqp/YsvjjJD0aq1MgBAABsY1pr66vq5CQXZOhWpp9qrV1TVe9NcllrbWWSk6vq+UkeTvKzjDClKBEOAABgQmbMmJrJOK2185OcP6zsz3re/9FY6zStCAAASCIcAAAAHeEAAABI4poDAACYkBlTcyvTzcLIAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJDE3YoAAGBCyt2KAACA6UY4AAAAkiTVWpvqPrB18sFgq/W8ffftW/70PfboW/6Df/3XMbfxxMc+tm/5/F126Vv+0zVr+pbvNGfOwDbe84UP9S3/25P+n77lP7jzzoF1DXLo057Wt/z+X/6yb/kXLr64b/mz99lnYBvX33HHmPr07L337lv+O699Qd/y/X/7xIF1/fqee/Ytv+u22/qW/5df//W+5T/++c8HtvHkxz++b/n9v/hF3/IvXnbZwLq2JSe/oP/P46yvf30L94QkOe13fqd/+bnnjrWq6TP/ZSvy/le+ctL/bnrnF74wJT8rIwcAAEAS4QAAAOgIBwAAQBLhAAAA6AgHAABAEuEAAADoeEIyAABMwAxPSAYAAKYb4QAAAEhiWhEAAEyIaUUAAMC0IxwAAABJhAMAAKAjHAAAAEmEAwAAoCMcAAAASdzKFAAAJmTGjOnzffv0ORIAAGBChAMAACCJaUUAADAhnpAMAABMO8IBAACQRDgAAAA6wgEAAJBEOAAAADruVgQAABPgbkUAAMC0IxwAADBqVXVEVR26mds4qqpOHee+u1TVSX3K51bV7VV11sR7OH0JBwAAjMURSTZbOKiqWa21la2108dZxS5J/kM4SPLnSf7PuDu2nRAOAAC2E1W1pKou63ktGbZ+YVVdV1XLq+qGqvqHqnp+VV1SVT+oqkOSvCHJm6pqdVUd3qeNeVX1w6qa0S0/pqpuq6odquqEqlpVVVdU1Zeq6tHdNsur6uyqujTJB6pq8SPf8FfVS6vq0qr6XlV9o6rmd+WnVdWnqupbVXVzVZ3SdeH0JHt1/Tuj2/bgJPOT/PNmOq+T/poqLkgGANhOtNaWJVk2wmZ7Jzk6yWuTrEpybJLDkhyV5B1Jzk6yprX2wQFt3FdVq5M8N8lFSV6S5ILW2sNVdW5r7RNJUlXvS/K6JGd2uy5IcmhrbUNVLe6p8uIkz2mttap6fZK3J3lLt27fJL+ZZOck11fV0iSnJlnUWntG186MJB9KcnyS549w7Ns94QAAgF63tNauSpKquibJhd0f5lclWZhk9SjqOCfJqzIUDo5J8vGufFEXCnZJslOSC3r2WdFa29CnrgVJzqmqJyaZneSWnnVfba2tTbK2qn6SodGB4U5Kcn5r7fap/EZ+W2FaEQAAvdb2vN/Ys7wxo/9ieWWSI6tq1yQHJ/lmV748ycmttQOSvCfJnJ59HhxQ15lJzur2OXHYPr193TCgf7+R5OSqujXJB5P8YVWN93qGac/IAQAAY/FAkrmb2qC1tqaqViX5SJKv9IwI7JzkzqraIclxSe4YRXvzerZ79Sj7t3NPX4575H03XelZrbVx3Qlpe2DkAACAsfhyklcMuiC5xzkZmud/Tk/Zu5JcmuSSJNeNsr3TkqyoqsuT3DPSxq21e5NcUlVXP3JBMqNn5AAAgCRJa+3WJIt6lhcPWHfgKOr6YpIaVrY0ydI+2y4etrw8Q1OQ0lo7L8l5ffY5bdhyb7+PHdCnX9U7mTwhGQAAmHaMHAAAMC5V9c4M3fa014rW2vunoj9MnHAAAMC4dCFguw8CM2ZMn8k40+dIAABgO1JVR1bV9VV1Y1X9hzswVdWbq+r7VXVlVV1YVU8eqU7hAAAAtjFVNTPJx5K8KMl+SX6/qvYbttn3MnTr1gOTfDHJB0aqVzgAAIBtzyFJbmyt3dxaW5fk80le1rtBa+2i1tovusXvZuhp05skHAAAwFamqpZU1WU9ryXDNnlSktt6lm/vygZ5XZJ/GqldFyQDAMBWprW2LMmyyairqo5P8qwkzx1pW+EAAAC2PXck2aNneUFX9u9U1fOTvDPJc1tra0eqVDgAAIAJmKInJK9K8tSqekqGQsExSf7dk6Gr6plJ/ibJka21n4ymUtccAADANqa1tj7JyUkuSHJtki+01q6pqvdW1VHdZmck2SnJiqpaXVUrR6rXyAEAAGyDWmvnJzl/WNmf9bx//ljrFA4AAGACamqmFW0WphUBAABJhAMAAKAjHAAAAEmEAwAAoCMcAAAASdytCAAAJmSKHoK2WRg5AAAAkggHAABARzgAAACSuOYAAAAmZDpdcyAc0Ne6dff0LZ89e7ct3JPNY7KO7/77rxm4bu7c/cdU11Ta1n7e/3TlxX3LH3rox33Lt9afxaDz/qbPfnZM22/KWH+Gr5/Ez8JY+zueNi790Y/GvM9kGXR8Jz3veX3LH96wYWBdBy1c2Ld85zlz+pbf/tOf9i3fuHFj3/K95s8f2PadP/953/K/+uo/9i3/6B/+Yd/yB9euHdjGILNn9f8z5F8HHN+9DzwwsK6Fj3983/Ld5s7tW75+wM/jUTvs0Lf84fXrB7a94+zZfcvXDthn0M/18ptvHtjGR//5n/uWL33ta/u3/fDDfcv/+DOfGdgGJKYVAQAAHeEAAABIIhwAAAAd4QAAAEjigmQAAJiQmkZ3KzJyAAAAJBEOAACAjmlFAAAwAdPpIWhGDgAAgCTCAQAA0BEOAACAJMIBAADQEQ4AAIAkwgEAANBxK1MAAJgAT0gGAACmHeEAAABIYloRAABMiCckAwAA045wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACRxK1MAAJgQtzIFAACmHeEAAIBRq6ojqurQzdzGUVV16jj33aWqThpWtqGqVnevlZPTy+nJtCIAAMbiiCRrknxnc1ReVbNaayuTjPeP+F2SnJTk4z1lv2ytPWOCXRuoTCsCAGC6qaqFVXVdVS2vqhuq6h+q6vlVdUlV/aCqDknyhiRv6r6FP7xPHfOq6odVNaNbfkxV3VZVO1TVCVW1qqquqKovVdWju22WV9XZVXVpkg9U1eKqOqtb99KqurSqvldV36iq+V35aVX1qar6VlXdXFWndF04PcleXf/O2AKnbVoRDgAAthNVtaSqLut5Lemz2d5JPpRk3+51bJLDkrw1yTuSnJ3kr1trz2itfXv4zq21+5KsTvLcruglSS5orT2c5NzW2rNbawcluTbJ63p2XZDk0Nbam4dVeXGS57TWnpnk80ne3rNu3yQvTHJIkndX1Q5JTk1yU9e/t3XbzemO97tV9fKRztP2zLQiAIDtRGttWZJlI2x2S2vtqiSpqmuSXNhaa1V1VZKFGfrDfyTnJHlVkouSHJN/m+KzqKrel6GpPzsluaBnnxWttQ196lqQ5JyqemKS2Ulu6Vn31dba2iRrq+onSeYP6M+TW2t3VNWvJflmVV3VWrtpFMex3TFyAABAr7U97zf2LG/M6L9YXpnkyKraNcnBSb7ZlS9PcnJr7YAk70kyp2efBwfUdWaSs7p9Thy2T29fNwzqX2vtju6/Nyf5VpJnjvI4tjvCAQAAY/FAkp03tUFrbU2SVUk+kuQrPSMCOye5s5v+c9wo25uX5I7u/avH2r+qemxVPap7v1uS/5Lk+6Nse7sjHAAAMBZfTvKKQRck9zgnyfHdfx/xriSXJrkkyXWjbO+0JCuq6vIk94y0cWvt3iSXVNXV3QXJT09yWVVdkaFpTqe31iY1HMyomvTXVHHNAQAASZLW2q1JFvUsLx6w7sBR1PXFJDWsbGmSpX22XTxseXmGpiCltXZekvP67HPasOXefh87bPMDRuovQ4wcAAAASYwcAAAwTlX1ziRHDyte0Vp7/1T0h4kTDgAAGJcuBGz3QcATkgEAgGlHOAAAAJIIBwAAsE2qqiOr6vqqurGqTu2z/r9W1f+tqvVV9XujqVM4AACAbUxVzUzysSQvSrJfkt+vqv2GbfajJIuTfG609bogGQAAtj2HJLmxtXZzklTV55O8LD1Pf+6eTZGq2jjaSoUDAACYgM3xROOqWpJkSU/Rstbasp7lJyW5rWf59iS/PtF2hQMAANjKdEFg2YgbTjLXHAAAwLbnjiR79Cwv6MomxMgBAABMwBQ9BG1VkqdW1VMyFAqOSXLsRCs1cgAAANuY1tr6JCcnuSDJtUm+0Fq7pqreW1VHJUlVPbuqbk9ydJK/qaprRqrXyAEAAGyDWmvnJzl/WNmf9bxflaHpRqNm5AAAAEgiHAAAAB3hAAAASCIcAAAAHRckAwDABGyOJyRPFSMHAABAEuEAAADomFYEAAATMEVPSN4sjBwAAABJhAMAAKBTrbWp7gNboU+ecELfD8ajZvWfifYHS5eOuY2lr31t3/I3fupTfcv/8eST+5YPukPAq848c2Dbg+pa89BDfctnzuifox/esGFgG79ct65v+R9/5jMD9+nnC6ecMnDdoP7uNGdO3/JXfvSjfcv/7g1v6Ft+3y9+MbDtU/7+7weu6+eTJ5zQt/z1n/jEwH0G9ev3P/q+vuWzZ+82pj5trdatu6dv+XQ5vrEa9PuaJL9/1lmT0sagc74pU/nz2BKfEZ/DbcM4fk7TZ/7LVuRLf/zHk/4H9e9++MNT8rMycgAAACQRDgAAgI67FQEAwAR4CBoAADDtCAcAAEAS4QAAAOi45gAAACbAE5IBAIBpRzgAAACSCAcAAEBHOAAAAJIIBwAAQEc4AAAAkggHAABARzgAAACSeAgaAABMyG6PnTfVXZg0Rg4AAIAkwgEAANARDgAAgCSuOQAAgAnZccHOU92FSWPkAAAASCIcAAAAHdOKAABgAubMf8xUd2HSGDkAAACSCAcAAEDHtCIAAJiAObubVgQAwHaoqo6oqkM3cxtHVdWp49x3l6o6aVjZnlX1z1V1bVV9v6oWTkpHpyHhAACAsTgiyWYLB1U1q7W2srV2+jir2CXJScPK/j7JGa21pyc5JMlPJtDFac20IgAAkiTdN+pfS/LdDAWAVUk+neQ9SR6f5Lgkb0iyoaqOT/I/WmvfHlbHvCRXJnlKa21jVT0myXVJfi3J4iRLksxOcmOSP2it/aKqlid5KMkzk1xSVVcmeVZr7eSqemmSP+32uTfJca21H1fVaUn27OrdM8mHW2sfTXJ6kr2qanWSr3f9n9Va+3qStNbWTOY5S5JH7froya5yyhg5AADYTlTVkqq6rOe1pM9meyf5UJJ9u9exSQ5L8tYk70hydpK/bq09Y3gwSJLW2n1JVid5blf0kiQXtNYeTnJua+3ZrbWDklyb5HU9uy5Icmhr7c3Dqrw4yXNaa89M8vkkb+9Zt2+SF2ZoNODdVbVDklOT3NT1721J9kny86o6t6q+V1VnVNXM0Zyv7ZGRAwCA7URrbVmSZSNsdktr7aokqaprklzYWmtVdVWShRn6w38k5yR5VZKLkhyT5ONd+aKqel+Gpv7slOSCnn1WtNY29KlrQZJzquqJGRo9uKVn3Vdba2uTrK2qnySZ32f/WUkOz9CoxI+6vi1O8rejOI7tjpEDAAB6re15v7FneWNG/8XyyiRHVtWuSQ5O8s2ufHmSk1trB2RoqtKcnn0eHFDXmUnO6vY5cdg+vX3dMKB/tydZ3Vq7ubW2Psn/TvKfR3kc2x0jBwAAjMUDSeZuaoPW2pqqWpXkI0m+0jMisHOSO7vpP8cluWMU7c3r2e7Vo+zfzj3Lq5LsUlW7t9buTvJbSS4bRT2jNmfevMmsbkoZOQAAYCy+nOQVVbW6qg7fxHbnJDm+++8j3pXk0iSXZOgi5dE4LcmKqro8yT0jbdxauzdDFzVfXVVndMHkrUku7KZGVZJPjLLt7Y6RAwAAkiSttVuTLOpZXjxg3YGjqOuLGfpDvLdsaZKlfbZdPGx5eYamIKW1dl6S8/rsc9qw5d5+Hzts3ddH02eEAwAAmJDZs58w1V2YNMIBAADjUlXvTHL0sOIVrbX3T0V/mDjhAACAcelCgCAwjQgHAAAwAXPm9Hu8wrbJ3YoAAIAkwgEAANAxrQgAACZg9uzdpqTdqjoyQw+am5nkk62104etf1SSv8/QU6rvTfKq7pa0Axk5AACAbUxVzUzysSQvSrJfkt+vqv2Gbfa6JD9rre2d5K+T/OVI9QoHAACw7TkkyY2ttZtba+uSfD7Jy4Zt87Ikf9e9/2KS51VVZROEAwAA2PY8KcltPcu3d2V9t2mtrU9yX5LHbapS4QAAACamJvtVVUuq6rKe15ItcSAuSAYAgK1Ma21ZkmWb2OSOJHv0LC/oyvptc3tVzUoyL0MXJg9k5AAAALY9q5I8taqeUlWzkxyTZOWwbVYmeXX3/veSfLO11jZVqZEDAADYxrTW1lfVyUkuyNCtTD/VWrumqt6b5LLW2sokf5vkM1V1Y5KfZihAbJJwAAAA26DW2vlJzh9W9mc97x9KcvRY6jStCAAASCIcAAAAHeEAAABIIhwAAAAd4QAAAEgiHAAAAB3hAAAASCIcAAAAnRrhCcpsv7a6D8b991/Tt/y2iy/uW77/b584aW3Mnbv/mLZPkjlz5vctnz17tzH1ad26ewauG2td42ljc7e9KYP69dBDP+5bPuicD9p+qI27+pbPnv2ESWlj0GdnU8b689jUz+L2m87tW/74Pf5r3/KxnttNtT+ZxzHIZLWxqd/l8fwMx9rGoM/h3LkHjamNLfHzG08bYzWez9RYfwcH1TWetqfyHI617dmzd6sxN8J2xcgBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAADrCAQAAkEQ4AAAAOsIBAACQRDgAAAA6wgEAAJBEOAAAYAyq6oiqOnQzt3FUVZ06zn13qaqTepZ/s6pW97weqqqXT1pnpxnhAACAsTgiyWYLB1U1q7W2srV2+jir2CXJr8JBa+2i1tozWmvPSPJbSX6R5J8n3NFpSjgAANhOVNWSqrqs57Vk2PqFVXVdVS2vqhuq6h+q6vlVdUlV/aCqDknyhiRv6r6FP7xPG/Oq6odVNaNbfkxV3VZVO1TVCVW1qqquqKovVdWju22WV9XZVXVpkg9U1eKqOqtb99KqurSqvldV36iq+V35aVX1qar6VlXdXFWndF04PcleXf/OGNa930vyT621X0ziaZ1WZk11BwAA2DJaa8uSLBths72THJ3ktUlWJTk2yWFJjkryjiRnJ1nTWvvggDbuq6rVSZ6b5KIkL0lyQWvt4ao6t7X2iSSpqvcleV2SM7tdFyQ5tLW2oaoW91R5cZLntNZaVb0+yduTvKVbt2+S30yyc5Lrq2ppklOTLOpGCoY7JslfjXD82zXhAACAXre01q5Kkqq6JsmF3R/mVyVZmGT1KOo4J8mrMhQOjkny8a58URcKdkmyU5ILevZZ0Vrb0KeuBUnOqaonJpmd5JaedV9tra1NsraqfpJk/qAOdfsfMKxNhjGtCACAXmt73m/sWd6Y0X+xvDLJkVW1a5KDk3yzK1+e5OTW2gFJ3pNkTs8+Dw6o68wkZ3X7nDhsn96+bhihf69M8r9aaw+P8hi2S8IBAABj8UCGpvEM1Fpbk6EpSR9J8pWeEYGdk9xZVTskOW6U7c1Lckf3/tUT6N/vJ/nHUba53RIOAAAYiy8necWgC5J7nJPk+O6/j3hXkkuTXJLkulG2d1qSFVV1eZJ7Rtq4tXZvkkuq6upHLkiuqoVJ9kjyL6Nsc7vlmgMAAJIkrbVbkyzqWV48YN2Bo6jri0lqWNnSJEv7bLt42PLyDE1BSmvtvCTn9dnntGHLvf0+dti6W5M8aaQ+Y+QAAADoGDkAAGBcquqdGbrtaa8VrbX3T0V/mDjhAACAcelCgCAwjZhWBAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAjnAAAAAkEQ4AAICOcAAAACQRDgAAgI5wAAAAJBEOAACAzqyp7gBbp3vuubBv+X233Nm3fK9nHz/mNtatu6dv+ezZu/Utv/v67/Utn7P7Y8bc9v33X9O3/LaLL+5b/tTnz+9fz93XD2xjzh799xmr+++/YuC62bOf0Ld83bq7+pbvttvzxtzGIIPqGuT2m87tW75gr98ZuM9tV3ytb/nuT3tm3/JBxzHoPCXJ3LkHjamuQed2kGvO7/+ZSpI9Djusb/mgz/oeBx3Zt3zQ72sy+PwO+v0b9Duw/2+fOLCNQb9PYz1XD913X9/yWz931cB9Dvyj3+tbPmdO/9+/yz/9V33Ln/zS/p+DJLnnnrEdx6DP26A+JYP/Lfnx//elvuVj/exsyqB/vx966MeT1sbA39kBxz1n3ry+5YN+X5PBn909+p+qgccxqK9z5+4/sO1Bv4ODPtNzd39a/+0HnPNk8Lla+9Nf9C1/1K6P7lu+qX9zITFyAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAACCJcAAAAHSEAwAAIIlwAAAAdIQDAAAgiXAAAAB0hAMAAEatqo6oqkM3cxtHVdWp49x3l6o6aVjZB6rqmqq6tqo+WlU1OT2dfoQDAADG4ogkmy0cVNWs1trK1trp46xilyS/CgddkPkvSQ5MsijJs5M8d6L9nK6EAwCA7URVLamqy3peS4atX1hV11XV8qq6oar+oaqeX1WXVNUPquqQJG9I8qaqWl1Vh/dpY15V/bCqZnTLj6mq26pqh6o6oapWVdUVVfWlqnp0t83yqjq7qi5N8oGqWlxVZ3XrXlpVl1bV96rqG1U1vys/rao+VVXfqqqbq+qUrgunJ9mr698ZSVqSOUlmJ3lUkh2S/HgznN5pYdZUdwAAgC2jtbYsybIRNts7ydFJXptkVZJjkxyW5Kgk70hydpI1rbUPDmjjvqpanaFv5y9K8pIkF7TWHq6qc1trn0iSqnpfktclObPbdUGSQ1trG6pqcU+VFyd5TmutVdXrk7w9yVu6dfsm+c0kOye5vqqWJjk1yaLW2jMeqaCqLkpyZ5JKclZr7doRzsF2y8gBAAC9bmmtXdVa25jkmiQXttZakquSLBxlHeckeVX3/phuOUkWVdW3q+qqJMcl2b9nnxWttQ196lqQ5IJun7cN2+errbW1rbV7kvwkyfzhO1fV3kme3tXzpCS/1W/EgyHCAQAAvdb2vN/Ys7wxo591sjLJkVW1a5KDk3yzK1+e5OTW2gFJ3pOh6T6PeHBAXWdm6Nv+A5KcOGyf3r5uGNC/VyT5bmttTWttTZJ/SvIbozyO7Y5wAADAWDyQoWk8A3V/hK9K8pEkX+kZEdg5yZ1VtUOGRg5GY16SO7r3rx5H/36U5LlVNatr97lJTCsaQDgAAGAsvpzkFYMuSO5xTpLj829TipLkXUkuTXJJkutG2d5pSVZU1eVJ7hlp49bavUkuqaqruwuSv5jkpgxNi7oiyRWttS+Psu3tjguSAQBIkrTWbs3Q7T4fWV48YN2Bo6jrixm6ALi3bGmSpX22XTxseXmGpiCltXZekvP67HPasOXefh87bPMTR+ovQ4wcAAAASYwcAAAwTlX1zgzd9rTXitba+6eiP0yccAAAwLh0IUAQmEZMKwIAAJIIBwAAQEc4AAAAkggHAABARzgAAACSCAcAAEBHOAAAAJIIBwAAQEc4AAAAkggHAABARzgAAACSCAcAAEBHOAAAAJIIBwAAQEc4AAAAkggHAABARzgAAACSCAcAAEBHOAAAAJIIBwAAQKdaa1PdBwAAYCtg5AAAAEgiHAAAAB3hAAAASCIcAAAAHeEAAABIIhwAAACd/x89GpSPUKxINAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vireoSNP.plot.anno_heat(AD/DP, col_anno=cell_label, col_order_ids=id_uniq, \n",
    "                        cmap=segpink, yticklabels=mtSNP_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PyEnv",
   "language": "python",
   "name": "pyenv"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}