{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Donor deconvolution\n",
    "\n",
    "Examples of using vireoSNP for donor deconvolution.\n",
    "\n",
    "**Please use the command line in [manual](https://vireosnp.readthedocs.io/en/latest/manual.html)**. This notebooks is mainly for development and diagnosis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vireoSNP\n",
    "import numpy as np\n",
    "\n",
    "from scipy import sparse\n",
    "from scipy.io import mmread"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load cellSNP data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Option 1 - VCF format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_vcf = vireoSNP.load_VCF(\"../data/cells.cellSNP.vcf.gz\", biallelic_only=True)\n",
    "cell_dat = vireoSNP.vcf.read_sparse_GeneINFO(cell_vcf['GenoINFO'], keys=['AD', 'DP'])\n",
    "\n",
    "AD = cell_dat['AD']\n",
    "DP = cell_dat['DP']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Option 2 - sparse matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "AD = mmread(\"../data/cellSNP_mat/cellSNP.tag.AD.mtx\").tocsc()\n",
    "DP = mmread(\"../data/cellSNP_mat/cellSNP.tag.DP.mtx\").tocsc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, `AD` and `DP` are sparse matrices, not `numpy.array`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Vireo\n",
    "For donor deconvolution, the `vireoSNP.vireo_wrap` function can be directly used, which contains multiple initialization and have options for different situtaions of donor genotype, fully missing, partially lack, or all available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[vireo] lower bound ranges [-29456.1, -25592.6, -19155.3]\n",
      "[vireo] allelic rate mean and concentrations:\n",
      "[[0.017 0.462 0.941]]\n",
      "[[28199.7 38020.2 21634.1]]\n",
      "[vireo] donor size before removing doublets:\n",
      "donor0\tdonor1\tdonor2\tdonor3\n",
      "241\t246\t234\t231\n"
     ]
    }
   ],
   "source": [
    "res = vireoSNP.vireo_wrap(AD, DP, n_donor=4, learn_GT=True,\n",
    "                          n_extra_donor=0, ASE_mode=False, fix_beta_sum=False,\n",
    "                          n_init=50, check_doublet=True, random_seed=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### The donor assignment probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from vireoSNP.plot.base_plot import heat_matrix\n",
    "\n",
    "fig = plt.figure(figsize=(2, 2.5), dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAHNCAYAAADys6skAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhdVZnv8e+vMpIQRMagMYgiiLaiggPYShwYpJ2HBptWGa4itNem6cYrjTIoCA4NqE3EKzQIoqKIAyiCMonCRYkDIGFQEyAMwQCSuapSee8faxec7BpS++Sc2jln/T7Ps5/i7PHdleLd66y9BkUEZmaWh566AzAzs/HjpG9mlhEnfTOzjDjpm5llxEnfzCwjTvpmZhlx0jczy4iTvplZRpz0zcwy4qTfASSFpBMbPs8p1s2pL6qhJJ0oKcsu3pKOlHRw3XEMR9LBxd/L7nXHYvVz0rdWOgfYo+4ganIkcHDdQZitz8S6A7DuERGLgEV1x9HtJG0CrA4PnGVNcEm/TSQ9X9K3JC2W1CvpPkkXSJrSsM9MSV+VtEhSn6QFkk6Q1JKHsaStJc2VdIek5ZIekXSNpNcMs+8Rkv5Q7LdM0p2SPtOwfZqkLxQxrpb0mKRbJL23YZ8h1TuSpkj6L0kPS1op6ReSdpO0UNL5DfsNVkG8TtJXJC2R9KikSyU9o3TOhZIul/RmSb+TtErSfElvbjjXfEkrJP16uGoNSbtL+lFxH6uL8/xjaZ8xxSRpIfBCYK9i/yjWjfZvE5L+W9Lhku4u/kbukHTgCDHsI+l/JP0VWAlMKbb/vaSri3+zlZJulPQPI1z26ZLOK+55haTLJD1ntDit+7ik3waSdgV+CSwBjgfuAbYD3gpMBnolzQR+DawFPgX8mVQ18gng2cAhLQhli+LnScDDwKbAO4DrJL0hIq4r4j0QmAt8GfiPIqYdgRc0nOt04H1FfL8DpgN/B2y5nhjOAw4APgdcU5zz+8BmI+x/DvBj4J+AZwGfB74BvL60367AqcApwBPACcClkk4F3gD8JxDAZ4HLJe0QEauK+30d8FPgZuDDxfEHAhdLmhYR51eM6R3AJcV5jizW9a7n9wLp7+F1pL+RFcWx35K0JiIuKe37P0UM7yP97vsl7QX8DLgVOKy45pHAZZLeGxEXl85xbrH/4H2cTPpbeHFE/G0M8Vo3iAgvLV6Aq4HHga1H2edsYBkwu7T+30nJ6gUN6wI4seHznGLdnIpxTSA96H8OXNqw/svA4+s59jbg++vZ58T0J/Xk5xcUcZ5W2u/AYv35DesOLtadVdr3mGL9zIZ1C0ml3Wc2rNu12O9BYFrD+rcV69/SsG4+8FtgYulalxXH9zQR0+3AdRX+LaK4h21L/z7zgXuG+b18fZhz3AQsBjYtneM24H5ApXNcWjp+z2L9cXX/P+Nl/BZX77SYpGnAXsB3IuKvo+z6ZuBa4EFJEwcX4Ipi+14tiufDkn4raTWwBugnlYR3adjt18DmRXXU2yRtNcypfg28SdJpSq2HNhnD5Qfv4Tul9ZcUsQznR6XPtxY/ty+t/31EPNDweX7x87qIWDnM+u0BJO0IPB+4qPjc+Lv/Cekb2c5NxlTV1RGxePBDRAwAFwM7SppV2vd7jR8kTQdeCVwSEctL57gQmMXQ+7io8UNE3AjcS/q2YZlw0m+9p5NKW+t7obkt8BZSEm5c/lhsHy7xViLpaOArpGqMdwGvAl5Oqtp4MmlHxIXAoaQk9j3gEUk3S9q74XQfJVWVvJ30sHpM0g8kPW+UEAarfhY3royINcCjIxxTXj9YTVJ+yDxWOmffcOuBwfVTi5/bFj+/wNDf/dxiW/l3P9aYqnp4lHXlarOHSp+fDmiY9ZC+rQx3jpGut74qOusirtNvvceAAVJJazRLSCXG40bY/uAI66v4Z1LJ94jGlZJmlHeMiPOA84oS5GtJ7wEul7RTRNwbEStI9eYnSNoWeBNwGqlK5PkjXH8wWW4LPFkqL0rVdSWaJcXPU4FLR9jnrnGKZeYo68oPmnJLncdJ7162G+Ycgy+Zl5TWj3S9P40So3UZl/RbLNLLwuuB94xQTTLoctKL0D9HxC3DLK1I+kHphaKkFzNKW/qIWBERV5BekE4mtUop77M40svObwE7F1Vaw/lF8fOA0vp3U1OBIyLuIr1Y33WE3/stEbGsiVP3Ur3k/4biAQqApAmk39WfIzV/HVHxEL4ZeGdjVZukHtLDfhFwd+mwgxo/SNqT9O3uuopxWwdzSb89jia13rlZ0mmkktS2pNYahxdJ5Xhgb+BGSV8ilS6nklru7A98eH3/44/B5cAnJZ1EehDtXFx3AQ3/9pK+BqwCfkWqLpgJHEtqjfKbYp+bi/PdSipl7kJqSXJTqQ79SRHxR0nfAv5d0gCp9c4LSS+rnyCVVOtwOHCFpCuB80nfQrYg3dPLIuI9TZzzNuBASQcAfyG1o79tPccsAa6R9Gmear3zfNKL7rE4ltQa51pJXyBVZR1JKky8NyLK3w52l3QO8F1S651TSPc+F8uGk34bRMQfJL2CVEVyKjCDVHd6DUUdc0Q8VLQf/ySpNcgsUmueBaQ698dbEMopwDRSc76PAXeQmii+g9QCaNANpBYe/0iqK15Cemi9v+Fl9DWkh9a/Fed8ALiguMZoDiE9SA4rjv19cZ2fArU0E4yIa4t/n+OAM0n3/Cjp91N+6TxWJ5CqWr5G+ve+l/QAH82PSO9wTgZmk5rtHhRDm1oOKyKul/R60t/Z+aRv7n8A3hoRlw9zyGGkB/W3Se38rwX+NSLK70Gsi2loYcCsvYpqhV+REtw3646nDkqd2M6KiI/UHYvlxSV9a6uiBdAewDxSFdKuwMdJ9eojvUg1szZx0rd2WwrsAxxFqvZYQuqLcGxErK4zMLMcuXrHzCwjbrJpZpYRJ30zs4w46ZuZZaSWF7mSROoq3kzPRzOzKmYADw7TWS1LdbXeeQaeYcnMxs8sGsZ/ylldSX8ZwKn7vYSpkybUFEL7HXqOm6Gb1WnpsmU8a8e/A9cqPKnWdvpTJ01gky5O+pttNtLkUGZm9fCLXDOzjDjpm5llxEnfzCwjTvpmZhlx0jczy4iTvplZRpz0zcwy4qRvZpYRJ30zs4w46ZuZZcRJ38wsI076ZmYZcdI3M8uIk76ZWUac9M3MMuKkb2aWESd9M7OMOOmbmWXESd/MLCNO+mZmGXHSNzPLiJO+mVlGnPTNzDLipG9mlhEnfTOzjDjpm5llxEnfzCwjTvpmZhlx0jczy4iTvplZRpz0zcwy4qRvZpYRJ30zs4w46ZuZZcRJ38wsI076ZmYZcdI3M8uIk76ZWUac9M3MMuKkb2ZWA0mvlXSZpAclhaS3j+GYvSTNk7Ra0l8kfbjqdZ30zczqMR34A/CRsewsaQfgJ8ANwEuBzwBfkvSuKhedWDFIMzNrgYi4ArgCQNJYDvkwcF9EHFV8ni9pd+A/gO+N9bpO+mZmDSRNBSa38JS9EdHbgvPsAVxVWnclcJikSRHRP5aTOOmbmRUkTd1syqRVS3vHlD+HsxzYtLTuJODEDYmrMBNYXFq3mJTHtwIeGstJnPTNzJ4yeWlvP8fssilTJoypyuVJvQPB5+cv3xSYBSxr3NTC+KL0WSOsH5GTvplZydQJYmrFpN+w97KIWNrikAAeJpX2G20DrAEeHetJnPTNzEqktFQ9ps1uAt5SWrcPcMtY6/PBTTbNzIboaXKpQtKmkl4i6SXFqh2Kz7OL7adKuqDhkLOB7SWdLmkXSYcChwFfqHJdl/TNzErGqaS/O3Btw+fTi59fBw4GtgNmD26MiAWS9gfOAP4FeBD4aESMubkmOOmbmQ2x5Q47scmkCZWOWdU/ALfOG/P+EXEd67wKGLL94GHWXQ+8rFJgJU76ZmYljy24u/KL3NUDY25AUysnfTOzko30RW5LOOmbmZWIUepdRjmmEzjpm5mVuKRvZpYRl/TNzDLikr6ZWUZc0jczy4gEPS7pm5nlwSV9M7OMbLHDTkyr2CN3Zf8A/HbsPXLr4qRvZlby+MK7WV2xR+4q98g1M+tMrt4xM8tID1H5RW7P2CevqpWTvplZiUv6ZmYZcdI3M8uIe+SamWXEJX0zs4z0NNEjt+r+dXHSNzMrcUnfzCwjm++wc+UeuVP6B+A37pFrZtZxnlh4F30TK/bIXeN2+mZmHcnVO2ZmGfGLXDOzjLikb2aWkyY6Z3VK1u+pOwAzMxs/LumbmZW4Tt/MLCOu0zczy8jTtt+Z6ZOrdc6a1DcAN7lzlplZx1l6312sqdg5a6U7Z5mZdSZX75iZZcQvcs3MMuKSvplZRlzSNzPLSA/Ve652Sk9XJ30zsxLPkWtmlhGX9M3MMuKSfps084bczForojM6FTWj2XubMXunyj1yJ/QNAL9t6nrjyUnfLHvdm/SbvbcV991FuEdu602dkBYzq9FAX90RtM9Af1OHqYkmm67eGYMDv3opm222WZ0htNVXD9q37hDa7vCLrqw7BNtAmji17hDaRhObe6C5c5aZWUbcOcvMLCMu6ZuZZcQlfTOzjHRzSb9TOpGZmY2bwc5ZVZfq19GRkhZIWi1pnqTXrGf/oyTdJWmVpPslnSGp0pt4l/TNzErGYxgGSQcAZwJHAr8CDgeukPSCiLhvmP0PAk4DDgVuBHYCzi82/9tYr+ukb2ZWMm129Tlyo3qP3KOBcyPinOLzUZL2BY4Ajh1m/z2AX0XEN4vPCyV9C3hFlYs66ZuZlay6/y56KvbIXfVUj9wZWreupzciehtXSJoM7EYquTe6CthzhEv8EvhnSa+IiF9Leg6wP/D1KnE66ZuZlWxg9c6i0qaTgBNL67YCJgCLS+sXAzOHO39EfFvS1sAvlZ4qE4GvRET5wTEqJ30zs5INbLI5C1jWsKl3yM5PKQ/Yo2HWpQ3SHOA40juAm4EdgS9KeigiPj3WOJ30zcxKNrDJ5rKIWLqe3ZcAAwwt1W/D0NL/oE8DFza8A7hN0nTg/0o6JSLWjiVON9k0MysRT5X2x7pUeUhERB8wD9i7tGlvUsuc4UwDyol9gIrPKJf0zcxKxqlz1unAhZJuAW4CPgTMBs4GkHQB8EBEDLbkuQw4WtLveKp659PAjyJiYKwXddI3MysZj2EYIuJiSVsCxwPbAbcD+0fEvcUus1m3ZH8yqb7/ZOCZwF9JD4LjqlzXSd/MrKRH0UTSrz6JSkTMBeaOsG1O6fMaUkugkypfqIGTvplZSTePveOkb2ZWssnsndmkYo/cgb4B4HftCaiFnPTNzEp677+LSZOqld17+z1HrplZR3L1jplZRjyJiplZRlzSNzPLSDOTojQziUodnPTNzEpcvWNmlhFX75iZ5URClet32hNKqznpm5mVTJm1M1OmVOuc1ddbebrEWjjpm5mV9D1wF30VO2f1uXOWmVmHaqr5TntCaTUnfTOzEjfZNDPLiJp4keukb2bWoZz0zcxy0kP1GcQ7ZMZxJ30zsxKX9M3MMuIXuWZmGXFJ38wsI5Oe8TwmTamWHif1rgF+056AWshJ38yspP+he+iv2CO33z1yzcw6k2iieqdNsbSak76ZWYlf5JqZZSSNp1+1pO/qHTOzztTFs6g46ZuZlTTXZLMzsr6TvplZiev0zcxy0sVZ30nfzKyki3O+k76ZWdnEZ+zExIo9cif2rgFuak9ALeSkb2ZWMvDQPQxMrjZW8kDf2jZF01pO+mZmJa7eMTPLiJtsmpnlxJ2zzMwy0sX1O076ZmYlXZzznfTNzIZook6/U7K+k76ZWYlL+mZmOenirO+kb2ZWMnHb5zFxasUeuavXAA+3J6AWctI3MysZeORP7pFrZpaLLq7dcdI3MytrbmL0zsj6TvpmZmXukWtmlg/19KCeanX6qrZ7bTokTDOzcTRYqV91qXwZHSlpgaTVkuZJes169t9c0lmSHiqOmS9p/yrXdEnfzKxsHN7kSjoAOBM4EvgVcDhwhaQXRMR9w+w/GfgZ8AjwbmAR8CxgWZXrOumbmZVIPahifU0TBf2jgXMj4pzi81GS9gWOAI4dZv9DgS2APSOiv1h3b9WLbnDSlzQBeBFwb0Q8vqHnMzOrW8+2O9JTsXNWz+o1wEKAGaWWP70R0du4oii17wacVjrNVcCeI1ziraT5GM+S9Dbgr8A3gc9GxMBY46yc9CWdCdwWEecWCf/6IsiVkt4cEddVPaeZ2cZk7eI/s3ZKtZL+2t4nO2ctKm06CTixtG4rYAKwuLR+MTBzhEs8B3g9cBGwP/A84CxSHv/UWONspqT/buAbxX+/BdgBeD7wfuAU4NVNnNPMbOOxYXX6s1i3nr136M5PivJZhlk3qIdUn/+homQ/T9IzgGNoc9LfiqcGmNgf+G5E3C3pXOCjTZzPzGyjknJ+1RezT/7nsohYup7dlwADDC3Vb8PQ0v+gh4D+UlXOfGCmpMkR0TeWOJtpsrkYeEFRtbMf8PNi/TTSTZiZdbaenuaWMSoS9Dxg79KmvYEbRzjsV8COWvcN807AQ2NN+NBc0j8P+A5wO+lryM+K9a8E7mzifGZmG5XBidGrLhWdDvwvSYdK2kXSGcBs4Owihgskndqw/1eALYEvStpJ0j8A/0mq1x+zytU7EXGipNtJ7UO/2/BWeoChb6LNzDrPOLTTj4iLJW0JHA9sRypI7x8Rg80wZwNrG/a/X9I+wBnArcADwBeBz1a5blNNNiPikmHWfb2Zc5mZbXTUU31cBY30/nVkETEXmDvCtjnDrLsJeFXlCzUYU9KXNOYXtBHxpebDMTOrXzPVNU1U79RirCX9fxvjfgE46ZtZZ+viAfXHlPQjYod2B2JmtrHQ1s9Bm0yqdsyqfuCO9gTUQh57x8ysJJYsIKZMqHZMb2e0WB9rnf7pYz1hRBzdfDhmZvVznT68dIz7VX99bWa20WlmfPwuSvoR8bp2B2JmtrFQj1BPxZJ+xf3r0nSdvqQdgecCv4iIVZIUES7pm1nnG6d2+nWoPAyDpC0lXQ3cDfyE1JMM4BxJ/9XK4MzMajFO0yXWoZmxd84A+kldhFc2rL+YNACbmVlHG6exd2rRTPXOPsC+EbGodJP3ANu3JCozszrl3jmrZDrrlvAHbcXokwWYmXUGqYk6/bXr32cj0EzS/wVplqxPFp+jGN/5GODaVgVmZlaXnq12oKdij9yeVf3ALe0JqIWaSfrHANdJ2h2YDHwOeCFplnZPlWhmHS8eXUhUnBg9Vq9pUzStVflFbkTcAbwY+DVpApXpwKXASyPiz60Nz8ysBj1qbukAzY6n/zBwQotjMTPbKEg9qGKdftX969JMO/1DJL1nmPXvkfSB1oRlZlYjt9Nfx8dJM7mXPUKar9HMrMOpyWXj10z1zvbAgmHW30vqsGWFzvgTsNx18+gpzd6bR9lc1yOkF7kLS+t3BR7d0IC6ycTOqOKz3K3tjHHgm7K2ybbzTY290xn/wzeT9L8NfEnSMlKbfYC9SLOyf7tVgXWDiZ3x4LfcRWc0NWxKs/fmHrnr+ASpiudqYPA32gNcgOv01/H+C66sO4S2++pB+9YdQlsdflH3/xtq4tS6Q2gbTexr7jhX7zwlIvqAAyR9AngJsAq4LSLubXVwZma12GJ7mDa52jErm3vAjLemx9OPiHtIg6yZmXWXx++H1RXT46rOqCbzxOhmZmWu0zczy4iTvplZRpoaWtlJ38ysM3VxSb+ZsXeGPUZSjyT3yDWzzuexd0DSZpK+A6yQtFjSSZImNOyyNcMPz2Bm1lkGe+RWXTpAleqdT5OGWngfsDmpk9Zukt5ZtN0HDzdjZt2gi6t3qiT9twMfiIjrACR9H/gxcJmktxb7dO/ITWaWj6fP7trOWVW+j2xFGkkTgIh4FNgbmAH8BJjW2tDMzGoimqjTrzvosalS0r8f2IWGevuIWCZpH+Aq4Pstjs3MrB5/ewB6q02Mzqr+9sTSYlVK+lcBh5RXRsRyYF9gdauCMjOrVRe33qlS0j8BeMZwG4oS/xuB3VoSlZlZnTyePkTE48Djo2xfDlzfiqDMzGrl1jtDSZoE/APwPOAh4PsRsaJVgZmZ1aaLS/pVOmfdKGnz4r+3BuYBFwMfBL4G3CHpmW2J0sxsPHVxnX6VR9OrgMGGq6cAA8D2EbETMAtYBHyqteGZmdWgi3vkNhvlXsAnIuJheLLN/nHA61sVmJlZbbq4pF+1Tn+wx+3mDB1nZwGw3QZHZGZWt6c9E6ZNqXbM5N72xNJiVZP++ZJ6gUmkydHvaNi2HfC3VgVmZlabpQ9Cf8VhGFZ1xjAMVZL+1xv++4fApqXt7wJ+v8ERmZnVrpnqmi6r3omIIb1xS04kvdw1M+twonoS77Kkvz5uo29mXcPt9BNJ20n6lKRrJM2XdLukyyQdVppQxcysc41T6x1JR0paIGm1pHmSXjPG4w6UFJJ+UPWaVTpn7Q7MB94CTAV2An4LrAC+ANwgaUbVAMzMNjrjkPQlHQCcSer39FLgBuCK9U07K2l7ipzbzK1VKemfCZwRES+NiD2BDwA7RcSBwHOATYCTmwnCzGyjMj6ds44Gzo2IcyJifkQcRRrC/ogRw0o1KheRBsD8SzO3ViXKlwEXNnz+JvAySdsWg7F9DHh3M0GYmW1UNqykP6OYU3xwGdLgX9Jk0qjEV5U2XQXsOUpkxwN/jYhzm721Kkn/EdbtfLUt6UXw0uLzPcAWzQZiZrbR2LCS/iLgiYbl2GGusBUwAVhcWr8YmDlsSNKrgcNI4501rUrrnR8AZ0s6BugFPglcHxGriu07Aw9sSDBmZhuFzWbC9KnVjpn45DxSs4BlDVtG66pbnldcw6yjeF/6DeCDEbGkWmDrqpL0P0Eq6V9GekLdBPxzw/Zg+CeamVlnWbYYBpqeGH1ZRCwdbVdgCalfU7lUvw1DS/8AzwWeDVymp6qRegAkrQF2jog/jyXMKp2zlgMHSJoKTCw+N24v102ZmXWmNrfTj4g+SfOAvVl3fvG9SSMelN0JvKi07mRgBvCvpBfAY1K5c1ZEeC5cM+tuPUpL1WOqOR24UNItpJqTDwGzgbMBJF0APBARxxZ59/bGgyX9DSAi1lm/Pi3rkWtm1jXGYbrEiLhY0pakFjnbkZL6/hFxb7HLbGBttSDWz0nfzKxMaqJ6p3qP3IiYC8wdYduc9Rx7cOUL4qRvZjaUJ0Y3M8tIFw+45qRvZlbmpG9mlhFX75iZZWTTbav3yO3pjNbsTvpmZmUr/gpRcWL0ld05MbqZWfdz9Y6ZWU6aGR/fL3LNzDqTS/pmZhlxk00zs4w46ZuZZWScxt6pg5O+mdkQKpaqx2z8nPTNzMpcvWNmlpFNt4bpm1Q7RqvWv89GwEnfzKxs5aOgisMwrPQwDGZmncnVO2ZmGXHnLDOzjLikb2aWEZf0zcwy4pK+mVlGnPTNzDLi6h0zs4xssgVMm1btmLUr2xNLiznpm5mVrX4cJlTsbLXaPXLNzDqTq3fMzDLiF7lmZjnxHLlmZvlw9Y6ZWUZcvWNmlhEnfTOzjDjpm5llpEdpqXpMB3DSNzMrm7o5TK3YI3fNlPbE0mJO+mZmZb1LYVJ/xWPcI9fMrDNJTdTpu3rHzKwzuZ2+mVlOVCxVj9n4OembmZW5yaaZWU5c0jczy4fr9M3MctJD9VEzXb1jZtaZXNI3M8vIlM1g6vRqx/R3RjrtjCjNzMZT3zLoHah4jCdGNzPrUN1bp98ZUZqZjafBOv2qS+XL6EhJCyStljRP0mtG2feDkm6Q9Hix/FzSK6pe00nfzKxsHJK+pAOAM4FTgJcCNwBXSJo9wiFzgG8BrwP2AO4DrpL0zCrXddI3Mxuip8mlkqOBcyPinIiYHxFHAfcDRwy3c0QcFBFzI+L3EXEn8MHiom+oemdmZtaozSV9SZOB3YCrSpuuAvYc42mmAZOAx8Z8Yfwi18xsqA0bWnmG1n0A9EZEb2nvrYAJwOLS+sXAzDFe8TTgAeDnVcJ0Sd/MbAg1uQCwCHiiYTl2lAvFMBcurxsanfQx4L3AOyNi9VjuaJBL+mZmQzTTGufJ/WcByxo2lEv5AEuAAYaW6rdhaOl/3atI/wH8J/DGiLi1YpBO+mZmQ0yZAVMq9sjte7LiZFlELB1t14jokzQP2Bv4fsOmvYEfjnScpGOATwD7RsQt1QJMnPRtg3TGaCM2moj11iZ0rKbvrW8F9DVxTDWnAxdKugW4CfgQMBs4G0DSBcADEXFs8fljwKeBfwIWShr8lrA8IpaP9aJO+rZBepz1O1+srTuC9mky6UtC1dvdV9o/Ii6WtCVwPLAdcDuwf0TcW+wyG2j8xzkSmAxcUjrVScCJY72uk34b/b+P7FN3CG33v75RbnFmnSYe+0vdIbRNLBtzAbhkfCZRiYi5wNwRts0pfX525QsMo9akHxFd/dXyxbtsUXcIZut3Zxc/uFdUatjyFA+t3B7NfIXqJNP+5dt1h2Ab6KsH7Vt3CG13+EVX1h1C2/QsXQqcXP1AJ30zs5x07yibTvpmZmUu6ZuZZcRJ38wsI5Omw6RNKx7TnlBazUnfzKysfxX0V6yj71/VnlhazEnfzKzM1TtmZhlRTxNDK7v1jplZB+uMkntVTvpmZmWu3jEzy4ird8zMMuKSvplZTsZnlM06OOmbmZW5esfMLCMTp6al0jFr2hNLiznpm5mVremFNRXHVVgz3PznGx8nfTOzIVynb2aWjy5uvdMZbx7MzKwlXNI3Myvr3todJ30zs6G6N+s76ZuZlXVxnb6TvpnZEC7pm5nlwyV9M7OMTJhSvUfuhP72xNJiTvpmZmUDfdV72A70tSeWFnPSNzMbwnX6Zmb5cJ2+mVluOiOJV+Wkb2Y2RBMl/Q55SDjpm5kN0b11+h5wzcwsIy7pm5mVSEIVq3eq7l8XJ30zsyG6t3rHSd/MrGzi5LRUPaYDOOmbmZUN9Kel6jEdwEnfzGwIV++YmeXDPXLNzHLikr6ZWT5c0jczy4lL+mZm+XBJ38wsJy7pm5nlY8LEtFQ9pgN4wDUzs4w46ZuZla0daG6pSNKRkhZIWi1pnqTXrGf/d0m6Q1Jv8fMdVa/ppG9mVjb4IrfqUukSOgA4EzgFeClwAyLJsQIAAArCSURBVHCFpNkj7L8HcDFwIbBr8fM7kl5Z5bpO+mZmQ6jJpZKjgXMj4pyImB8RRwH3A0eMsP9RwM8i4tSIuDMiTgWuLtaPmZO+mVlZm0v6kiYDuwFXlTZdBew5wmF7DLP/laPsP6zOeN1sZjaOli5fCT0Tqh+TzChNqNIbEb2l3bcCJgCLS+sXAzNHuMTMivsPy0nfzOwpfcDDz3re31VKpA2WA4tK604CThxh/yh91jDrNmT/IZz0zcwKEbFa0g5AK2dEKZfyAZYAAwwtpW/D0NL8oIcr7j8sJ30zswYRsRpY3eZr9EmaB+wNfL9h097AD0c47KZi+xkN6/YBbqxybSd9M7N6nA5cKOkWUkL/EDAbOBtA0gXAAxFxbLH/F4FfSPo/pAfD24A3An9f5aJO+mZmNYiIiyVtCRwPbAfcDuwfEfcWu8wG1jbsf6OkA4GTgU8DfwYOiIibq1zXSd/MrCYRMReYO8K2OcOsuwS4ZEOu6Xb6ZmYZcdI3M8uIk76ZWUac9M3MMuKkb2aWESd9M7OMOOmbmWXESd/MLCNO+mZmGXHSNzPLiJO+mVlGnPTNzDLipG9mlhEnfTOzjDjpm5llxEnfzCwjTvpmZhlx0jczy4iTvplZRpz0zcwy4qRvZpYRJ30zs4w46ZuZZcRJ38wsI076ZmYZcdI3M8uIk76ZWUac9M3MMuKkb2aWESd9M7OMOOmbmWXESd/MLCNO+mZmGXHSNzPLyMS6AzDbmE3KoFgUA2vqDqFtYmCg7hA2OrUm/ehbRfR173NHk6fVHYJtoEMvvLLuENruooP3rTuEtlnZ76RfVmvG1eRNujoxRqytO4S2kzIoCne5lf11R9A+q7r43prVvcXsjYATonWCbi6adPO9NctZycwsI7WW9JcuW1bn5c0MWNXF9d6ru/jemqWIGP+LSs8EFo37hc0sV7Mi4oG6g9gY1JX0BTwDGK+i/gzSQ2bWOF5zvHX7PXb7/YHvsZ3XfDDqSHYboVqqd4pf/rg9ddMzBoBlEbF0vK47nrr9Hrv9/sD32EZd+btsll/kmpllxEnfzCwjuST9XuCk4me36vZ77Pb7A9+jjYNaXuSamVk9cinpm5kZTvpmZllx0jczy4iTvpmNStJ1ks5s4fnOl/SD8bymPaWrk76k10q6TNKDkkLS2+uOqZUkHSvpN5KWSXpE0g8k7Vx3XK0k6QhJt0paWiw3SXpT3XG1S/FvGk541Ug6WNLf6o6jE3R10gemA38APlJ3IG2yF3AW8Cpgb1IP66skTa81qtZaBHwc2L1YrgF+KOmFtUbVBpJeDnwIuLXuWKx7dXXSj4grIuITEXFp3bG0Q0TsFxHnR8QfI+IPwCHAbGC3mkNrmYi4LCJ+EhF3F8txwHLSg65rSNoUuAj4IPB4jXFMl3SBpOWSHpL076XtTy+2Py5ppaQrJD2vYfuJkn5fOuYoSQuHudYJxTfUpZK+KmnyKHFNlvQ5SQ9IWiHpZklzim1zgPOApxXfkkLSiRvye+hmXZ30M/S04udjtUbRJpImSDqQ9A3uprrjabGzgB9HxM9rjuPzwOuAdwD7AHNYtxBxPukb11uBPQABP5E0qeJ13gDsUlzrvcX1Thhl//OAVwMHAi8Gvgv8tHjg3AgcRRpjZ7ti+ULFeLLhmbO6RDFy6enALyPi9rrjaSVJLyIl+amkUv47IuKOeqNqneJB9jLg5TXHsSlwGPD+iPhZse4DFMOgFwn2rcCrI+LGYt1BwP3A20mJeKz6gEMjYiXwR0nHA5+X9MkozTMq6bmkB8OsiHiwWP0FSfsBh0TEf0p6gjSW48NN3XxGnPS7x3+TSkB/X3cgbXAX8BJgc+BdwNcl7dUNiV/Ss4AvAvtExOqaw3kuMJmGb1ER8Ziku4qPuwBrgJsbtj9abN+l4rX+UCT8QTcBmwLPAu4t7fsy0jeKuxtG6QSYAjxa8brZc9LvApK+TCqBvTYium5ymojoA/5UfLyleOH5r8Dh9UXVMrsB2wDzGhLaBOC1kj4CTImI8Zr+SU1uFzA4nsvaYfarUvUz3LgwPcAA6XdV/l0sr3Buw0m/oxVVOl8m1YfOiYgFNYc0XkQq5XWDq4EXldadB9wJfHYcEz6kB2s/6SX5fZBe3AI7AdcDd5ByxitJ9ehI2rLYPr84x1+BmZLUMGnJS4a51q6SNomIVcXnV5ES+HCFlt+RHoTbRMQNI8TeV+xj69HVSb+oo9yxYdUOkl4CPBYR99UUViudBfwT8DZgmaSZxfonGv5n6miSPgNcQao3nkF6kTcH2K/GsFomIpYB67yDkbQCeHS8381ExHJJ55Lq1h8FFgOnkErvRMQ9kn4IfE3S4aSZr04jTYj0w+I01wFbAx+TdAnp3+lNDJ3IZDJwrqSTge1JI2/+d7k+v7ju3ZIuAi4oWhP9DtgKeD1wW0T8BFgIbCrpDaRm2itL1UdW6PbWO7uT/kB+V3w+vfjvT9UWUWsdQWqxcx3wUMNyQI0xtdq2wIWkev2rSaXM/QZfNFrLHQP8AvgR8HPgl8C8hu2HFJ8vJ9XDC9g/IvoBImI+cCTwL6Tk+wqGb0lzNXBPca3vAJcBJ44S1yHABcB/kf4WfkT6W7i/uO6NwNnAxaRvGx+rctM58dDKZmYZ6faSvpmZNXDSNzPLiJO+mVlGnPTNzDLipG9mlhEnfTOzjDjpm5llxEnfzCwjTvrWUsX8p4MTWfRLWizpZ5IOleS/N7Oa+X9Ca4efkiayeDZp3JVrScMHXy6ptvGelHT1eFNm6+Okb+3QGxEPR8QDEfHbiPgMaVC4NwEHA0iaLemHxbR8SyV9R9K2gycYnHZP0vskLZT0hKRvS5rRsM8USV8qptxbLemXxbDLg9vnFN849pV0C9ALvEbSrpKuVZpQfqmkeZJ2H69fjlmdnPRtXETENaQBuN5ZDAn9A2AL0uTue5Mm8Li4dNhzSTMyvblY9iJNkj7oc6RJVT5AmmjjT8CVkrYonedzwLGkiT5uJc1Fu4g0U9VupJEi+1txn2YbO3/VtfF0J2l2rzcWP3eIiPsBJL2PNG3eyyPiN8X+PcDBxfDDSLqQNLfqcZKmk0YZPTgirii2f5D0ADmMNNfroOMbR+WUNBv4fETcWay6py13a7YRcknfxtPgDEu7APcPJnyAYurDv7HutHsLBxN+4SHSLFOQvgVMAn7VcI5+4NcMnbrvltLn04FzJP1c0seLOVjNsuCkb+NpF2AB606v16i8vlzlEjz1N6uGdaOdA2DFOieJOBF4IfBj0kQcd0h6x/rDN+t8Tvo2LiS9njQt4PdI0+7NLiYFH9z+AtKEMPOHP8MQfyJNkffkRPCSJpEmzlnvOSLi7og4IyL2AS4lTdJh1vVcp2/tMKWYunECaear/UgvUi8nzX60luKFqqSjSH+Hc4HrI6JcFTOsiFgh6Sukqf0eI83p+jFgGnDuSMdJ2oRU338J6VvHLNIL3e81cZ9mHcdJ39phP1L9+xrgcVKrnY8CXx+cA1XS20mTuv+C9BD4KfC/K17n46RvqxeS5s+9Bdg3Ih4f5ZgBYEvSw2dbYAmppH9CxWubdSRPl2hmlhHX6ZuZZcRJ38wsI076ZmYZcdI3M8uIk76ZWUac9M3MMuKkb2aWESd9M7OMOOmbmWXESd/MLCNO+mZmGXHSNzPLyP8HL3wrbd5d/aIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 400x500 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from vireoSNP.plot.base_plot import heat_matrix\n",
    "\n",
    "fig = plt.figure(figsize=(4, 5), dpi=100)\n",
    "# assign_prob_comb = res['ID_prob']\n",
    "assign_prob_comb = np.append(res['ID_prob'], \n",
    "                             np.sum(res['doublet_prob'], axis=1, \n",
    "                                    keepdims=True), axis=1)\n",
    "im = heat_matrix(assign_prob_comb, \n",
    "                 cmap=\"Oranges\", alpha=0.8,\n",
    "                 display_value=False, row_sort=True)\n",
    "plt.colorbar(im, fraction=0.046, pad=0.04)\n",
    "plt.title(\"cell assignment prob\")\n",
    "plt.xlabel(\"Donors\")\n",
    "plt.ylabel(\"%d cells\" %(res['ID_prob'].shape[0]))\n",
    "plt.yticks([])\n",
    "plt.xticks([0, 1, 2, 3, 4], [1, 2, 3, 4, \"doublet\"])\n",
    "\n",
    "# plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "## donor id with highest probability\n",
    "donor_ids_best = np.argmax(res['ID_prob'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13 cells are called doublet\n",
      "4 cells are unassigned to singlet or doublets\n"
     ]
    }
   ],
   "source": [
    "## doublet rate\n",
    "\n",
    "doublet_threshold = 0.9\n",
    "is_doublet = np.sum(res['doublet_prob'], axis=1) > doublet_threshold\n",
    "print(\"%d cells are called doublet\" %(sum(is_doublet)))\n",
    "\n",
    "prob_threshold = 0.9\n",
    "is_unassigned = (np.max(res['ID_prob'], axis=1) < prob_threshold) & (~is_doublet)\n",
    "print(\"%d cells are unassigned to singlet or doublets\" %(sum(is_unassigned)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Allelic ratio per variant per donor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "## If ASE_mode is False\n",
    "AF_SNPs = np.tensordot(res['GT_prob'], res['theta_mean'][0, :], axes=[2, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAHqCAYAAAAAtunEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxbZdUH8N9JMpnJbAVKC13oXnZE2VQQqIUKKJZSikVLSwEVqCJYkRcQBAQRAXFhFUFKy1alWPZKEWQTZFF2ShdoS/dCl5lkMktmnvePpDIzzZObe9PMeTLz+/KZT2lyb+6Z25mcPNt5xBgDIiKiLULaARARkVuYGIiIqAMmBiIi6oCJgYiIOmBiICKiDpgYiIioAyYGIiLqgImBiIg6YGIgIqIOmBhKgIhMFRGT+RqV5XkRkcWZ5//Z9RFueyJymYiYTo/9M+j3l+3czP26LHCQ25iITBORqVkeH5KJdavniIohoh0A+VIP4HQA/+z0+OEAhmeep/x9GcAK7SDamQbgEwAzOj2+GulYl3R1QNQzscVQWmYDOEFEajs9fjqAlwAs7/qQSpcx5mVjTFESQ6YVF9sWr2WMacrEun5bvB6RFyaG0nJf5s9vb3lARHoBOAHAn7OdICJREblYRBaISJOIrBeRO0WkT6fjJorIkyKyWkSSIvK+iFwtIlWdjpshInERGSEij2f+/2MR+Y2IlHt9A/leJ1/5fn+Wc7fqShKRASJyW+Z7ahaRVSLygIjslMdr3SgiZ4rI+wCaAJySee5SEfm3iGwQkToR+Y+InC4i0u78pQD2AnB4u27DpZnnsnYlichXROQfIlIvIg0i8i8R+UY+940oF3YllZY6AA8AOA3AHzOPfRtAG9KtiXPbHywiIQAPATgUwDUA/gVgMIDLAfxTRA4wxiQzh48E8DiA3wFIANgdwP8BOAjA6E5xlAF4GMAdAH4D4DAAlwDYDOAXHt+Dn+vk5PP7y+f1BgB4Fenv7yoAbwHoDeAoANsDWOvxEuMysfwCwBoA6zKPD0H632tLi+5LAG4AMACf3a/jkf633Yx0lxKQTi62WA8HMD8T4+mZY6cBeEREvm2Mme31/RJZGWP45fgXgKkADIADAIzK/P9emedeAXBn5v/fAfDPduedlDl2fKfXOyDz+FmW6wnSHxoOyxz3uXbPzcg8dmKncx4DsMDn95XrOpelfzw7HP/PoN9f53MzjxkAl7X7+x0AmgHsEeDfyADYBGB7j+NCme/5EqTHE6Tdc+90jjHz+JDM609t99hLSCeq6naPhQG8DeDj9q/LL375/WJXUul5FulByNNEZB8AB8LSjQTgWKTfrB4RkciWLwBvIP2JdtSWA0VkmIjcKyJrALQCaMlcCwD26PS6BsAjnR57C+lP6zn5vI6XvL+/PB0D4BljzPs+z9viaWPMxs4PishoEXlKRDbjs+/5F0i3Rvr6vUim2+2LAB4wxsS3PG6MaQUwC8BAALsF+xaI2JVUcowxRkTuBPAjABUAFhpjnrccvhOA7ZD+FJzNjgAgItUAngfQCOBiAAsBNADYBcCDADoPojYYYxo7PdaUiccqwHW85PX9+dAHhc1SWt35ARE5CMCTSLdYvpd5/Waku51+Bv/fM5Du1pJs1wOwKvNn7wCvSwSAiaFUzUD6E+eZSL+52HwC4FMAR1ue3zK9dTSA/gBGGWO2fHqHiGxXcKQdbevr5Pv95Ws90p+2g8q2HeJJSLcQjm2fTEVkXAHX2Yj0uFK/LM/1z/z5SQGvTz0cE0MJMsasFJFrkR64vSvHoY8i/cYUNsb8O9dLZv7sPNh5RvAou+Q6+X5/+XoCwGQR2c0Y88E2eD0g/T2nkO5CAgBkprFOznJsE/JoQRhjEiLybwDjReQ8kxlgzwzGn4x0q2ThNoideigmhhJljLkgj8PuBzAJwOMi8nukB6pbkP5U/FUADxlj/ob0bJ6NAG4Vkcszx0wCsO82DntbXyff7y9fP0d6nOE5EbkK6YHc7ZBukVxvjFkQIMbHAEwHcK+I3IZ0F895yD7j6G0AJ4nIRAAfAmg0xrxted0LkZ6V9IyIXId099Q0AHsD+LYxhpu5U2AcfO7GMoORY5GeejkewN8AzAVwAdL9/G9njvsUwDeQ7u+/G+nB7DiAids4nm16nXy/Px+vtxLpabOPZl5jHtLTSnsB2BAwxqeRnl68D9ID9r9Eelrq1VkOvxTpgfg/IZ3kOg/wt3/dZ5Humksg3bV4fybOsYZTValAwg8WRETUHlsMRETUARMDERF1wMRAREQdMDEQEVEHTAxERNQBEwMREXWgssAtU4e+P7jjGBHpqwGwiosCP6O18rk/3NpSkYh6toEAVmoH4QqtxFAPAL9+5nFUVAfauKtHOGnYcO0QnNdiWrRDoBJWX1+PPYZ+DmDvRQeqtZIkVoNQrFozBKfV1NZoh+C8xtbO1b+ps0iIJdHIH9WfmAdu/Q8i0W2yX3q39L2bR2iH4LzycFQ7BCphYeH8m2xUE8PY7x2Iimq2GGxWJDgM46VXdFtvGdH9VEfYXWsTQlg7BCepJoaaaAtiUfYR2zy9ivfGy5SRheyrQz1deoIkdaaaGOY+3YBIBf9hbGafwW17vTSkGrRDcB7HGOya22y7wvZsqj8x40ZXIsZZSVYV4ZxbKBOAVQ2rvA/q4frG+mqHQCVGNTHcM3MJIuWVmiE4bcrVu2qH4LyqCMeoiLY13TZmPA40t3of10NtbNqkHYLzOCuJaNtTTQwnnfUFxDgryertDUu1Q3DeITvvrR2C85pa2Y9uk2rjB9NsOCrlsCX1TJpeDtlZOwL3cfDZLhLidNVsVH9iTh4xDLW1tZohUImbuWihdgjOmzKSY1U20RC7IrNRTQytJoVWk9IMwWn3LVmiHYLzRvfnrDYv/B2z473JTjUxhCWCsLCZa5NMsZnrhbOSvHGMwa6plYtIs1F9V775neUsiZFDLy5joG2A62HsmsNMmtnw47rD+lUmtUOgbiDF7hIr3pvsVBNDbUUKsQr+w9jc9BDvjZejf8AiehQcB5+zU00M99/yX5bdzmHO9Udrh+A8LgL0xkWAdg0ptsqzUU0Mk6btywVuOdRGeW+8rG5Yox1CCeDPkU0zB+azUk0M3xo2hOsYcuBsEm9Pr2J1VS+TRnAVoFW0TTsCJ6kmhseWL0ZlDT/N2Hy1X3/tEJx3YB/tCNzHDxh2nK6anWpiOHbQSLYYcljVsFo7BOe9ul47AvftsR0rGNukIpzgkY1qYqhrjsM0c6MemylXvKkdgvMevfKr2iE4r82wu8SG9yY71cQQDUc5YyKHmZfsqx2C8x74iPtie2GtJLuQhLRDcJLu4PP1/+FGPTlMnjhAOwTnjRvM3cm88FOxHe9Ndrorn9euAbiOwWrMgOHaITivqowfLLxsbq7TDsFZ9c312iE4STUx/PEXX0ZNLWcl2Uy5+BXtEJz3j+uP1w7BeUyedq1lHHzORjUxPPRRChXV/IexqRo5VDsE57UYTjf0EmY/uhXvTXaqiWHq7n1QW1ujGYLT9u7NLgAv65LrtENwHkuT29U3J7RDcJJqYtjYtBGpJn7is6ku46cZIup6qolh0eaNqGzjqkybuhZWRfcykBu4eeKUcLsm3pusVN95RvbanoPPOewQ3UE7BOfNW7lUOwTnHT2QP0fkj2pieG5NA2JxdpfYrIlzjrWXc/cdoR2C8xItLDRow+qq2akmhmMG9kENB5+tEqm4dgjOu2vhcu0QnHfabkO0Q3AXq6tmpZoY7l20GRXVrZohOG3lxjLtEJx37eFsMXjhZkZ2nJWUnWpi2D7WjFiMTTmbs/cZqB2C82YuWqgdgvNOGj5EOwRncYFbdpz24jABx1+8tBreI6JtTTUxfHv4CO7HkMN9Sz7QDsF5E4ftoh2C87gI0K4+yXG8bHQ36rnmRVZXzeGHk7k9mZdIiI1eL1z5bNcWMdohOEn3t2rjBiCaVA3BZUcPPFA7BOoGOLvNLpHi4HM2qonh5O/ugVg1P83YVEbYmvKyIsGNeoi2NdXE8LenGxCp4NaeNpNGcI61F3aTeGOLgfxSTQyzTt+dg885cIMVb/cs3qgdgvNO220n7RCcFSnnWqFsdPd8DkURDbGIlc3s5Uu1Q3BeLX98PMXC3CXRpiXM6s7ZqCaGTU2b0drE7hIKbsJQLgL0srJhpXYIzqpvYDdbNqqJ4alV6xCrY4Evm36VTdohOG9D0wbtEJzHcRg7TlfNTjUxTBg6kmMMObDcg7ePE43aIThv6sgB2iE4KxTl5JdsdBe4TX8AkSj7P20e/cNx2iFQN8CuJDt2JWWnmhgm/eDzXMeQw6eNn2iH4LynVrEIo5cpI4doh+CsqhRn/mWjmhgmDhvGrqQc7vjgQ+0QnDdh6I7aIThvdcMa7RCcxRZDdqqJYexFj7ErKYc5vx6tHYLzHlnGAnFepozcVTsEZ9W1ssWQjWpi+MsVR6KWO7hZzf5wlXYIzusVZT19L81t7G6z4b3JTjUx/OrlTYhW8Rfb5vW5b2mH4LxHrz1aOwSibkc1MQzvm0SsOqwZgtP2POtz2iE4j+sYvHEdg11DC6s7Z8Ni9g5LcVE4ESlQTQzl4TaUh/nuZzNpxDDtEJx3w9ssu+3l+3uyoJRNUzjYvRGRCgBBTm42xji/KlM1MRwzsA9qOPhs1WqYNKlwTa0cYLVpDnBvRKSiLFabbEkGmtG0RkSGup4cVBPDGTOXIVJRpRmC0+actZd2CM6bMJy9oV64H4NdwB3coi3JOmx3+P9BIuV5n2RSTdj07K93RrqlwcRgc9Y3wqis4eCzzdE/eVQ7BOfNv36sdgjOq2tmYrAppIielFVAyip8nFA6dZlUE8NBfYdxHUMOF52nHYH7Zi5aqh2C8yYO668dgrOCjjEAACSU/vJzfIlQTQwPL1vJWkk5TBqxt3YIzrtrEcuGeOEYg12QMYb/EfHXCmCLIT933/4+S2Lk8OmpnE3iZcJwH035Hqq8kE/F3VxhLQYmhqK4+7LDOCsphxPOe1I7BOdNuPpg7RCop2JiKI65yz5FrJq7lNnM+/047RCcd8/i5dohOG/cYHYl2RTUlRSS9Jef40uEamIYN7g3Www53L9kqXYIzhszoFI7BOdxuqpdwOmqAICRu/T21RWeak7ilcBX61qqiWHeivWIVbNWCQU3f6XT08GdcMpIrqC3qStgo55FKzZCyvJ//zItpfOzqpoYTOaLKKjjBvfWDsF5qxpYvt2moI16OMZQHMfswpIYucxdyq09vTy6fL12CM6bMHSgdgjOipSXBT+ZiaE4Ei0JhFpK52aRe8pCbHOSkhB8Dj4XLZJtTjUxPLO6EbF61rqxYd+wt/97brV2CM47brB2BN0UVz4Xx8kjhqG2tlYzBKed8dAS7RCcN3xQAV0BPURlhDO3bFKRAnaQZFdScdy9+EOWxMjhhm8O0g7BeT97nmMMXhpSDdohOKshVcCsSLYYimPS8KFsMeRwxA/u1w7BeTOvPVQ7BOexVpJdIQvcRg7cHpHy/FtjqaYGrmPIxyPLP0BlDVsMNnzT88Z1DN7GDebvWDEsWrUJUpZ/5QZTQvtLqyaGxRtrUNHCH1qbR14Nviqzp7ji69zoiZRwjKE4KiKtqIi0aobgtMP25SZGXp5exf5zLycN76sdgrNaywoZfOYYQ1E8+9ImRGItmiE4bebpQ7VDcN520V7aITjv3+vf1Q7BWYn6QlrlPlsMYIshL+NGVyJWza4Am83Nm7RDcN7cZVwd7mXcYK58tqmX+uAnsyupOCYO4zqGXO5evFg7BOdNGjFEOwTn/eeThdohOKugFgO7koojJCGESuhmdbWTR4zQDsF5R174hHYIznvqV8doh+CsuvLg1VXZYiiSY6950dc84J7m1JP6aYfgvDmXH6IdgvNYXdWukOqqIwf08rmOIcp1DPk4aXw/rnzOIcX6cJ4eWb5OOwTnsbqqXSHVVRetroNE8588Y5q5joG2gbC0aYfgvFbDrkhSwjGG4mARvdxG//hB7RCcN+fq0dohOG9D0wbtEJxV31TIRj3wOcYQ/FJdTTUxXPLCakSrCpgu1s3NvPIg7RCcVx6OaofgvA82c9qzTaKhkFlJHHwuij37JRCrLp2b1dUGVu2qHYLz7lzIKb1exgzYTjsEZ9W3FvAWyMRQHFWRFCoLWZLezbFcsreyEEuqkJKQ+NzBjYkhL7f9/h1EojHNEJz21d/21w7BeS1trCflpSrCmX82bZECpv6xxVAcD143loPPOURD7D/3MuvPL2qH4LxJV3HDJ5tCiuiN3LkXIhX5l/RJNZZxHUM+oqEo3/xymLmIpQy8XHT2LtohOG9dkms9bOqTwWclLVpbD4nm35Vpmkuna1g1MTSkkoikuJTC5sA+2hG47+V1XDnvZXeOPRcHu5KK49Z3NqCimtsO2kz//BDtEJz38rrl2iFQT8XB5+LY1BBBubDFYJNq44wtL++t5uQFL+MGs7S9TWGDz1z5XBQfvLsWkVgBKw+7udk7cWtPL/v0Z9kQL+Xh7bVDcFZTQQskuVFPUXzxwN6oYBE9q4nDemuH4LzZH7JyqBeWxLArrCSGQHwkBsMxhvyctkc1amprNENwWlMrx1+8jBu8o3YIzntmNZOnTUN98MSw687VPqerhvDvwFfrWqqJYVHdKlQZ9n/aLNhcrh2C8+LNHKPyMmkEF0ra1BdQq23R2jhC5fl3ZbY1BZuuKiLTAPwUQD8A7wI41xjzfI7jzwVwFoBBAD4B8ACAC40xjfleU/W36vJrFnDlcw7f+O4XtUNw3iH9WITRG8cYikFCAvEx08jPsf87R2QigN8BmAbgRQBnAHhCRPY0xmw1JU9EJgG4GsBpAP4FYFcAMzJP/zjf66omhtGTDuIYQw6n7cFuNi9zl+a/UUpPNbCKEzxsEqngEzy6aBnDdAB3GGNuz/z9XBE5CukWwYVZjv8ygBeNMfdm/r5URO4D4KtUs2pieHr2fxEpZ4vBZsEX99QOwXmj9iudKYBaBlZxBzebutbgez6Lz8FnP8dmjo8C2B/pFkB7TwI42HLaCwBOFpGDjDGviMgwAF8HcJefa6smhkd/dSxrJeUw+oLHtENw3vv999MOwX0jtQPongT+JqC2O7amU5JoMsY0ZTllRwBhAGs7Pb4WwM7ZrmGMuV9E+gB4QdIXiQC4xRjTObnkpJoY7l78Ifd8zuHU0/kb7eWEob20Q3Dey+ve0Q7BWYn6QrqSArcYVnR66nIAl+U4tfMqPMny2JZrjALwM6THJP4NYASA34vIamPMFfnGqpoYJg7bhS2GHFra2H/uhTu4eXt9PWf+2TTGg698LiAxDATQftZEttYCkJ5R1IqtWwd9sXUrYosrAMxqNybxtohUAbhNRH5pjMlrGpVqYnjgo2VsMeQwql+FdgjO+8cqrvXwctwQ/hzZ1NeV4WcBzx25U5XPdQyCNZnLGmM8BzeMMc0i8jqAMQD+1u6pMQAespxWCaDzm38rfPZ8cRK4wwZVs46+lykjWU/Ky8xFS7VDcFYynvfU/q0sXpdAqDz/FkfAdQzXA5glIq8BeAnA95Fen3ArAIjITAArjTFbZig9AmC6iPwXn3UlXQHgYWNM3jXCVRNDSAxCUkARq26u1fBNj8hVxZ6VBADGmNki0hvAz5Fe4PYOgK8bY5ZlDhmEji2EK5Eef7gSwAAA65FOFr4aRmJM178xi0gtgM0rPvmIYww5cIMVb1OuX6IdgvMevfBA7RCcVVdXjyF9RwBAr3y6d4DP3r8G/OSvCJXnvx9IW1MDVv7mRF/X0qLaYghLCOESKkXb1Z74OHgzt6eYPHGAdgjOYxE9u0KK6HVFi0GLcmKIIMz9GKy+t/sI7RCcx+42b9e9sVo7BGcVNiup227gppsY3vx0AaqaOSvJpm+M3Wxe+lVmXedD7fSrTmqH4Kwkgt+b9AZufjJD4Et1OdXEMLi6H2pqWA/IZu6yT7RDKAFLtQNw3pgB3Bfbpr4u+EZPbDEUyRkzlyES4+Ibm0lHs46Ul7+/XkK/bUpO3ZW1kmwKqZU0YsdqX+9fqaSUzMcY1cTwza9Wc4FbDkcN4CY0XsYP4adhL7MWLdQOwVnJePDB5yWfJhDysWVKW1PpbNWrmhjKQq0oC+W95qLHqY0yaXrhLnfeUoYz/2wKuTcCn11Jga/U9XRbDIMGopZbe1p9sIlz9L388MY12iE4j+sY7Orq6jE96Mk+p6uW0iCDamJ4Zd2HqExyjMHmsH67a4fgvJnT+fPjZd6KpdohOKuQPZ8llP7yc3ypUE0Mo/rvxZXPOfxr7VvaITjv87251sNLoiXYXsM9QTIV/C2QC9xIxcE7fU47BOetSHQubU+dNbWV0EfVLlbIveF01SI5dvoDiEQ5JdPm6Zu/rR2C87htpbdoiLOSbFpDhaxjYIuBFGgUOCw1pfTLRt3LsB2rUeZjHUNLUlAqKVo1MUw4cz+uY8jh/OdWaYfgvAHbc5c7L6fs1lc7BGfV18VwTsBzP/o0jnBF/h/eWhu5jiEv5eE2lIeDN+W6u2sPZ+VQLxubNmmH4DyWVrErZIEbu5KK5KRhwzkrKYc7Fy7WDsF5s2av1A7BefN+9mXtEJxVV1cXeB1Duohe/seb0skLuomhxbSgxbArwGbMAO7V62XWOtue6LRFqo2lyW1SbQVUXvCbGfwcq0w1Mcz5aDnHGHKYMnJX7RCcN/8P47VDcN4Nb3NKr01jQV1JnK5aFMcM7IMalsSwGv2D+7VDcN7kcw7QDsF5p+zGYow29XVhXBbwXMn85+f4UqGaGOYu+xSx6ibNEJw29dz9tEMoAZy84CWRCv6puLtLpILPFGKLoUieeDaOiI/pXj3NzO8O1Q6BuoFnVnPas00htZKG9a5CWSz/rvCWJPB24Kt1LdXEkFiwiCufczjhxwu0Q3DenN8eox2C8w7qw5l/NvXlwUtifLShAeGK/JsBrY2lU7NKNTFMmrYvB59zmHE/N3H3Mm8FZyV5OXrgTtohOKstErzHgusYiuSem95giyGH+TedqB2C88LCqi5e/rSA62FsClvgxjEGUjDmyhe1Q3DeDyf30Q7BeZNGDNEOwVl1dfU4P+C5IgLxsTaBLYY83f2rMZyumsOcpSxl4GXc4GHaITjvrkUfaofgLLYYstOdlbRiPWLVSc0QnHbvX9dph+C88Le0I3DfuMFcx2BTX1ceuCSGiCDk492+rYQyg2pimDhsGGsl5dA2sXR+kLTMmV+vHYLzxp3OxFAMHHwuktUNaxCPcPGNjSmlqltKEu9wSq+X7cu/oB2Cs8IFTFcdukMMZZWVeR/f0tCG1wJfrWupJobn1jQgFue2gza9oix+5oW73HmbuahUtofpeoWMMSzdmESkMf/3r1SydLrNVRNDv1gTKis5McpmdP+9tUNwHvdj8LYhGdUOwVmNBdwbCfmclcTqqvlZ0RBDLJR/U6ynuXsx55974ZuetwnD+eHLpr4uEryIHscYiuM7w4dw8DmHI698XjsE51X15s+Pl7P32Uc7BGdVpeoCn8vpqkVS11wP06wZgdv4puftxDH5b8beU61uWKMdgrPqG7i1ZzaqiaGyLIaqMnYl2cw8ndVVvbCktLe5S7lLok1jPPi94RhDkYz/2d9ZK4kKMvPKg7RDcN64IWXaITirvq4MFwc8l11JRXL3ZYexJEYOf3oveP9nT9E31lc7BOfds3i5dgjOSsYbA587eLsYoj7WMTSXt6FUqp+pJoa/r1yH2ObSqVHe1X64zy7aIThvXZJlQ0jH8s2NiDSH8z4+lQyehLqaamJImRBShgvcbGZ/yJ23vJwykkX0vISF055twhJ8a9iQwFetpBIaYtBNDA8/nUCEQwxWpo3bnnpJtbFyqJcJQ9ndZlNfF8M5Ac/lGEORjBtdiVg1pxvalIWYGLxMHMYWgxeW3bYrqOw2ZyUVR6sJoZVdSVatrdoRuC8k/PkhHQKf6xjAxJCXgZVJVFbxF9tm/ofcnczL7xqXaofgvAnDK7RDcFZ9HQtVZsNaSQ7bfaeEdgjO4+CzN3Yl2RUyXZVdSUUyfkg/1HIdgxVnJXk7cvpc7RCc9+i1R2uH4Ky6uvrAO7gN2i6GaGX+s2eao6XTN6yaGCKhMCIhVn60SbTkP0e6p5pz9WjtEJzHDxh2hQw+f7y5EZGW/N+/Ug1cx5CX8ec9zJIYOVTttbt2CM6rKiudT2FauOezXWF7PnO6alE89JtxLLudQ10zC8R5YRFGbz97fr12CM5qThRQecFnddVSygyqiSGEMEJgd4nNCdPnaYfgvMnnHKAdgvPOOYCzkmzq6wz+GPBclt0ukm9e9iQi5fzEZ9VvoHYE1A1URaq1Q3BWWyT4ItJ0SQx/x5cK1cRw0uSRiFXzh9amkDouPcWEof21Q3DeXz5aqR2Cs7jyOTvVxCBiIMKyDzZcFe6tMsIWp5ckZ7dZNRZwbwbVViBalX83XXOkdBbTqSaGReurUd7AdQw2e+20WTsE57UZtqq8REK8RzbhAu7Nx/WNKGvN/y20hdNV87Nbn3rEqtlisOlXWTo/SFq4wM0bF7jZ1dXV4/yA53bV4LOITAPwUwD9ALwL4FxjzPM5jt8OwC8BjAewPYCPAPzEGPN4vtdUTQyr62OoMFzHYMPqqt7m/eZY7RCcd9dC7uBmU9AYQxckBhGZCOB3AKYBeBHAGQCeEJE9jTFb/cOKSBTAfADrAEwAsALALgDq/VxXNTE8/eRSzkrKYc75+2qH4Dzu4OYtFuEiQKsC7o2E0l9+jg9gOoA7jDG3Z/5+rogcBeAsABdmOf40ADsAONgY05J5bJnfi6omhgnj+nFWUg6JFBe4eZlywb+0Q3DenOu+ph2Cs+rrgn8wLaDFUNPpvCZjTFOW46MA9gdwdaenngRwsOUyYwG8BOAmETkOwHoA9wL4tTEm7yyomhgeuHcJWww5nPLLI7RDcN7Mq22/H7TF3GWfaIfgrMK6kgKXxFjR6anLAVyW5ZQdAYQBrO30+FoAO1suMwzAaAD3APg6gJEAbkL6vf4X+caqmhjuvvAg1LC6qtXRv3xJOwTn/ePiw7RDcF5YFqFB6rYAACAASURBVGmH4KxC1goV0GIYiI59/lu1FjrpPNgoWR7bIoT0+ML3My2E10WkP9KD16WRGOatWI9YdVIzBLeVs5SBlxmLlmiH4LwxA9gqt6mvC54YBlZHUV5VnvfxTaHm/13WGFOXxymfAGjF1q2Dvti6FbHFagAtnbqN3gews4hEjTHNlvM68J0YRGQXAMYYsyLz94MAfAfAe8aY2/y81qE7V6Kmlj+0NuZ47QhKAefoe1mR2KQdgrMSDcE3w1qVaEbZ1kMDVi0Neb0n/48xpllEXgcwBsDf2j01BsBDltNeBPAdEQkZ879FPrsCWJ1vUgCCtRjuBXAbgFkisjPSU6PeBXCyiOxsjMm7uVIVqWIdlxwmDN1BOwTnNbX6+2Xria56mTsB2jQnygKfK/BXuUGsvT85XY/0e+1rSA8qfx/AIAC3AoCIzASw0hizZYbSLQDOBvB7EbkB6TGGiwD8wc9FgySGvQG8kvn/bwF4xxhziIh8LRNs3omhPFKBigi7S2y+9cf3tUNwXuKdBdohOG/+TSdqh+CsurpYAdVVi78fgzFmtoj0BvBzpBe4vQPg68aYLVNQB6Fds9kY83Hmvfi3AN4CsBLA7wH82s91gySGMnw2WHIkgIcz/78gE3jeTjz/UW7Uk8PYM76kHYL7vvJl7Qict7phjXYIzqpvCD4rKQSDkI9WgJ9j2zPG3AzgZstzo7I89hKAgt48giSGdwGcKSKPId3XdUnm8f4APvXzQt//0V6orGFXks1X+22nHYLzaqP8+fFywXO2cUpqTgSveMod3Dr6P6QHQn4K4C5jzJuZx8fisy6mvIwdvDt3cMth9BXPaofgvHk/Y4vBy579OPPPJhkPfm9EDEJ+xhhKqJK078RgjPmniOwIoNYYs7HdU7cB8LVP3iUvrEa0ylcJjx6Fb3reuNbD28zpw7VDcFZ9XQGlsCXz5ef4EuErMYjIF5FuGUQA/APppdkAAGPMUr8Xv/BL26GWC9ysbn5nlXYI7tvgq/eyR6qKsOaWTSE7uA2orkC5j/0YmqQb7scgIscD+CuARgApAOeJyE+MMb8LevGKcDkqwpyVZPOPFzZ6H9TDzbzyIO0QnPfMan7AsGmoDz74vDreiKjJ/7N1c6J0yuj7aTFcBGAGgDONMSkRuRjAxUiXhA1k7E8e5KykHIaPOVA7BOfNX1k6v2xaxgzgOJ5NfXnwXRL97kDZXccYdgMwyRizpT10LYDLRGRHY0ygKl0PXjeWg885hAPW6e1JOBXT2+9fK6HO7S5W0KwkdNshBl+JoRrA/9bWG2OaRCQJoBbpmh6+vbVhCapaqoKc2iO89Sm72bycNNxWZJK2GLA9q6vaNJa1eB9kEfI5K8nPsdr8zko6SkTab0QcAnCEiOy95QFjzMNbn5bdfr13Y4shh5UJrur1srmZdYBIB9cxfOauLI+1X1FukK4fnpdv3/4OIhVsMdjcMrWPdgjOm3L+C9ohOG/O9dzz2aa+Lpx1I4R8MDEAMMZs8w7vI77SCxXcwc3q7g9Kp+mpZf6NE7RDcN4Nb3feF4a2aCxgo56uKomhQXU/hmi4FeVh7kdrU1XNe+OFg8/eesdYgdYmWUB13p2rKlBRnf84YCO65zqGvLbKMsY8l+9rPnDrfzhdNZdw3r1yPdYLo76gHYLzbjl+iHYIzqqrq8OPAp67tqERUfGxjqGhdKZW+2kx/DPHc6bdn3m/5vfP2ZtF9HJYlWDS9PLgrPe0Q3Be4liOVdk0tASvlcRZSWnbWx6vBHAOgB8B+NDPxROpCNpaVHuznFYbDT6VrqeYeQnLPXiZu4zTVW2SBYwxcPAZgDGm/TRViEgIwGkALkV6o4gfIPusJatBVUlUVXMRl82XdtpTOwTnRUNR7RCc98Lbi7RDcFZLQ/BP8WwxdCIi4wFcBaAPgF8BuMEYH5ufZixPxBAT7vlsx24SL/vsMEQ7BOeNHMqFkjaN8cIGhEuoEeCL3+qqhyO9Rdw+yGwX17kl4ccJQwdxgVsO65LrtENwHrtJvA2oafM+qIdKSmH7MfT4Wkki8jiAIwDcCWCcMabgeYLjr3sVkXK2GGzmX3yIdgjO29zIOfpe9u5dpx2CsxpafXd0/E9I0l9+ji8VfloMRyNdbnsigG+JZSTFGLND3q+4dg3A6apW4255WzsE5805cy/tEJz36HJf+2f1KA0twaeE7xSrQEWlj3UMbd1wHQOAU7f1xbnnc243zlqvHYLzjr5ovnYIzpt31RjtEJxVVxe8NbWusRHl4fzfQpsau+E6BmOMrxlH+ahriaKlmbNKrBLBp9L1GE3BuwJ6ioiPRVg9TSH3hrOSLESkAumupSoA840xvubFLVxXhfIEi+jZzLmcYwxeXv9kuXYIzrtrka/lRT1KIesYujM/g8/XAogaY87J/D0K4CUAewFoAHCNiIwxxuS9O/seO9UjVl06WbSrbV/eTzsE532c4Kwk0sHB57RjkN7ec4tJAAYDGAlgOYA/I73V5zfyfUG/TbGeps1wmiGRqwQ+p6t20+qqg9BxxdXXADxgjFkGACLyewCP+7n4xGEjuY4hh9E/flA7BPft0Fs7AufNOZ9lQ2zq68oxPeC5LImR1oaOC/2+BOCKdn/fBHs9pazuW7IYMe7HYMc3PU9VvfnBgnSE4G+PhVIq/uMnMSwA8E0A14vIXki3IJ5p9/xgAGv9XLzVhNC67ff/6TYmTxygHYLzhtdw8NDLM6tXaYfgrIb64D8/fWMVqIjlvw6rsbV09lfxkxiuBXCfiHwD6QHnx40xH7V7/usAXvFz8ZNHDGNXUg5/WrBYOwTnLalni9NLryg36imG9Y2NKI/kv0Cuu65jmCMiX0d6cPlJADd0OqQBwM1+Lj7xT28hEuN0VQou8TF3cPMy75dHaIfgrEIWuLFWUoYx5ikAT1meu9zvxb/51WqOMeQwcVh/7RCcN/vDGu0QnHfPYq71sClkHQOnqxbJI8/EEYmVThbtag/iA+0QnDf3rH20Q3Denz9Yqh2Cs1IFzAhni6FIZn/vcxxjyGH0BY9ph+C81Q3cttJLWbh03pC6WqqAeyPwN9OohBoMuolhztKF7ErKgSUxvHE/hnyU0ltS1zIm+L1hi6FImltDCLdyuioRlZ4Q/LUYSumdruDEICJTAfwtyE5uAn6WyeWqlxPaIThvwbv12iE475ZT2N1mU18XfI+EHSsqEKvIfx1DMtU91zHY3Abg3wB8J4YThw3hGEMOd97xd+0QnDfzEpZ78DJ/ZenMn+9qyXjwe/NpYxIVkfzbAY2NwbcR7Wp+qqtuyPEaL4lIG+BvB7em1hY0tXLxjc3Y73B3Mi8LNtl+LOkz3D63GFgrKa0MwLMA/truMQFwO4BrAKz0e/FvXTgPEW7tafXI9XkXqu2xmtu20w7BeX9fzC5Jm+ZE8O6dEIzPWkndc/D5CwDuBTAawA+MMXEAEJE/AZhrjHkv18nZPHLNOHYl5XDED2drh+C8mdd8RTsE5+3Zj11JNsl48O4dthgAGGMWi8jBAH4J4A0ROcUY82IhF//5i2sQrWIRNJsZvz5YOwTnPbCkdDZY1zKolt21NmVlwe8Nt/bMMMakAPyfiPwdwL0icg8QvH10xVf6scWQw50LWUTPS68K7QjcF2/hns82yVTwe+N3VmUJNRiCzUoyxjwtIvsB+BOABIBAHXWvrn8fVY0somczqLp0prdp2dxcph2C8zY3R7VDcFYhZf/ZlZSFMeZTAOMLufhe2w9BbS2LoNlsaOKMGy/zVrD/3Eu8iS0Gm8bG4Pemd7nPdQwtpfNBz9ddEZEw0hv0LDPGtIlIOYDjkF7U94wxxtdGPd+64jlEyjmVzubpqzkryUtY2N3mpVcFx2Fsoqng92ZDUxIVZT7WMTR1z3UM+wKYB6AvgHcyG/Y8AWAo0uMMLSJylDHm1Xxfc9Dnh6GskrWSbDjG4G3KyCHaITjvohfWa4fgrKZE8AFhlsRIuwbACwAuB/BdAH8H8C6A/ZBODHcCuArAmHxf8NKje6GGXUlWtWW8N17GTH9YOwTnXXXhrtohOCtRn8CtAc9lEb20AwEcYox5X0QuBPBDAKcaY1oAQESuRnoBXN76VOyI2grOSrJZ1cC9er3sP+6L2iE47/1NpdOF0dWS8cK62UpoPNkXP4lBAGy5i53/BNIzk3y1lv760RKW3c5h/97d9cdu29mzXyk10Kk7YYsh7XWk1zBcCuB0AB8h3Wo4LfP82QDe8XPxowb0ZVdSDidMn6cdgvOq9tpNOwTnsbqqXSHVVTnGkHYh0oPPpwL4BMBXAfxZRFYDaAOwPYBv+rn4GTOXIRLjOgYblnvw9tQqdpN4eWZVg3YIzkrGg9+bHcpjqKzIf1ZlQ0sB+4h2MT8lMV4VkcEAdgPwgTEmLiKjAEwCEAMw3xjja5Pi+767F1c+53DkL1/SDsF5V32fRfS8LNzMBW42hXTucOVzRqZw3uvt/t4I4I6gF1/TsA6JCD/N2IwdO1A7BOct3Mw6QKRjU3MDmpryf7tPNgd7rxORaQB+CqAf0jNBzzXGPJ/HeScBuA/AQ8aYcX6u6WcdwwkAnjDGbLN38spIJSoj7Eqy6VXB/Yy91LEkhqfaaIt2CN1SVww+i8hEAL8DMA3AiwDOAPCEiOxpjFme47zBAK4D4JlAsvHTYvgrgLiI3A/gDmPMv4NcsL2LHl+Hskq2GGw2fMr+cy83nryjdgjOe2pV6ZRi6GqF1ErqosHn6Ui/396e+fu5InIUgLOQHvfdSqZCxT0ALgVwKADf/a1+C4VcC+B4AN8VkfeQ3qRnVqZukm83jRvBMYYc/rHqbe0QnPfUKtZK8rK5gHpA3V0htZKK3WIQkSiA/QFc3empJwHkqsn/cwDrjTF3iMihvi6a4feu/NEYc4WI7I/0lNVLAVwtIg8D+JMxZr6fF/vLR4u5jiGHP9+3UTsE51X15qdhL5yualdfF8FlAc8tYPC5RjqWWm0yxjRlOWVHAGEAnWvQrQWwc9ZriByC9Hvz532EtpWgZbdfB/C6iEwHcCLSaxnmicjHxpgh+b7OhKFD2WLIoWViWDsE54WldKYAapm7lGMMNo3x4PemgI16VnR66nIgZ37qfBHJ8hhEpAbA3QC+Z4wpaIDST2LYKpDMrKRZAGaJyAik1zjkrdW0odXwF9tm4rD+2iE4r6mVs5K8TLn9I+0QnJVqDL4f9nbRalSW5z95prz5f62EgQDq2z2VrbUApNeLtWLr1kFfbN2KAIDhAIYAeKRdiyQEACKSArCbMWZJPrH6LYlhZYxZDOBnPl4PDS1JhLm7lNXcZZyV5GVtvFw7BOed+DXO/LNJxg2eCXju5pY4mpvzbzEkW/6XhOqNMXVexxtjmkXkdaQLk/6t3VNjADyU5ZQFAPbp9NiVAGoAnAPg43xj9fOuPBTANq3fe/KF8xGJ5r/RRU9z01X7a4fgvPej9d4H9XDcwc2usFlJBiEfS+T8HNvO9Uj3yLwG4CUA30d6T5xbAUBEZgJYaYy5MNOD06EskYhsAgBjjK9yRX5WPi/z88L5+P45e6OyhoPPNntsx3LJXgZXc7qzl4m3+SpI0KOkksG7kgQ+t/YMcA1jzGwR6Y30TKN+SL/xf73d+/EgpEsSbVN+d3DbA8CXALxkjFkgIrsj3UQpB3C3MeZpP693262L2GLI4cZy9g17arJ1z9IWj157tHYIzqqrq8eQ/wt+fleUuTDG3AzgZstzozzOnRrkmn5WPh+NdL9WHECliBwPYCaAN5G+P3/P7OCWf3JIpYAQtx20+cGP9tAOwXmJFGdueZm3Yql2CM5qqI8HPpdlt9N+DuBaY8zFmRoc9wK4xRjzMwAQkV8CuABA3onh5DP34TqGHJKtpfODpKWQPuKegmMMdsmW4PeGZbfT9gIwJfP/f0F6muqcds/fh/TCirwZFFbdsLurDHMqr5fNTXzT81IVZavcppBP8WwxdGKMaRORRgCb2j1cD6CXn9dpg6CtpIrRdq36FKfyeqks48pnL8bwd8ymkHvTK1qNqmj+U4HLoqXz7+DnnWcpgBEAFmf+/mUA7av77QJgtZ+Ll4XaUBbip2KbJPvPPYVL53eNupm65jhSPtYxNDQHnwHV1fwkhluQrtsBIOu82GPgY3wh/RrpL8qOb3reUvxc4SlSSp3bJYRdSQCMMbd6PO9r1TMAzLrxTU5XzWHqj7nAzcusP/pat9MjzbzyIO0QnFXYns+CkI+ucD/HalPtxK7aY1fu+ZwDZ9x4mzptb+0QnPf8mtLpwuhqyXgBC9zE5wK30skLuokhsWI1IuX5b6bd08z6s3YE1B3MufwQ7RCcVV8X/P1HMv/5Ob5U6LYYBvZjiyGHE8fw3lDhHlm+TjsEZyXjhSxwY4uhKCZ9LYrKGlbHtBk3eIR2CM67a9GH2iE4773VHMezaU4UsI6BLYbiSKQiaGPZbatVDau0Q3DepBGDtENw3gnz39UOwVmF7MdQU1aNqrL8W/XhMiaGvBwzsA9qams0Q3Baryh3t/Ny3C1vaofgvKNG+d4LvsdojEcC78cQb4mjrcXHOoaW0pkEoJoYKstiqCrj4LMNWwx5KJ2p4WqquDrcKlTAvRERiI+BAz/HalNNDOPPe5jrGHIYe+aXtUNw3qnf4BiVl0QLVwHaRArYM5xF9IrkL9d8A7XsSrKa/SFbDF6OG8zNjLyMu/lt7RCclWoMvtETZyUVSUW4HBXhCs0QnMaBVW/jbuGbnpcjDt1eOwRnNcbLAo8xcFZSkYy9bD4XuOUw85J9tUNw3u579dUOwXk7xEpn0LOrJVubg5/ss8VQQnlBNzE8fNkY1NZy5o0NB5+9/Wh/9p97mTaTicEmlSygK4kthuIwaIPZ9vtYdxtVEe5u52X7ck7F9JJY/pZ2CM5KNQVPDNVlVb7WMUhZ4Et1OdXE0NTagqZCmnLdXHmYu5N52di0yfugnq6KHzCsCqhJnkglAB/FWROp0mm5qSaGb53/GKer5jDvxhO0Q3DexNs+0A7BeT+c3Ec7BGc11MfxytXBzmVXUpEcMPFQlFdxuqrNPYuXex/Uw7HQoLfldWx52jTGg3dlC/yNJ5dOWlBODOtW1aGskmMMNv++57/aIbhv537aEThv5vTh2iE4q74ugssCnsuVz0Vy8/EjOSsph9YJ/IX2MuZHD2qH4Lz5Kwdoh+CsZLwx8LlsMRRJQyqJSIrVVW0iId4bL5PP3k87BOc99yZrJdm0NLBWUjaq7zznPbICZZWcMWGz23Au/vNSEWHy9HIY10laJeNhPBDwXLYYioRdSbnNXLRQOwTnldIvm5a3VnGA3qY5EbzFUFVWiWof6xhQVjqlgFUTw8zFHyJWzRaDzem7skCcl9FXPKsdgvMmT+QYg00yngx8bkMqCUnlvw6iIRX8Wl1NNTGUSRvKCih7290dd+sb2iG4LxF8z96e4kt9+TtmE49xumo2qonhnpve4AK3HKb+eH/tEJy3IbmXdgjOe3mdj+W5PUwyHnzlMxe4FcmJZ+3HrqQcpowcoR2C89oMPw17OfLCJ7RDcFYhtZK4UU+RnDJyGAefcxg97T7tEJzHXe68TT5tN+0QnJWMx/HKDcHO5XTVIjHGwJjSGanvao/+4TjtEJzHIozeTrj0Re0QnFVIi4FjDEXSYlrQYlo0Q6ASl0hx8NkTq6vaFVBdlS2GIrny5Y0or2JisHnv7TXaIThv/JEswuiF01XtkvHg1VVj4RgqI/kvQm0Ll854mGpi+MXBO3GMIYc7dyqd+u1aHnyqXjsE5zF5FkdjaxJhH+sYGlu5jiEvl7ywGtEq/mLb7DuANW68zD1rH+0QnDfulre1Q3BWKlnAhy+fXUn+NojWpZoY9uyXQKy6dG5WVzt5BGeTUOESiz7SDsFZqebgn+I5+FwkVZEUKsu4+MZm9AWPaYfgvKmnj9QOwXnfPuNz2iE4KxmP45U/BjuXC9yK5LYZKxApZwVRK84m8TTjd//RDsF5M6/5inYIzqqvC/7BNCTpLz/HlwrVxHD3BQeippYDY1Yl9IOkxrCmtBeuY7ArbB0DWwxF8diKTxCrDr6DUnc38w/c2tMLN+rxxpXPdoWsfO7OowyqiSEEgxC48tmqlbOSiFwVC1cgFsm/CGhruHR+n1UTQ79YEyoruQOXDUtieCsLR7VDcN5R0x/WDsFZhcxKamxtRCQV9nV8qVB9V77uuvdZdjuHOb8dpB2C8ybe9oF2CM6bfBZnJdkUMisp3ZPkZx1DwOsoUE0Mj/xmPFc+53DEr17SDsF9Kz/WjsB5e0zYWzsEZyXCTYHP5eBzkbSYFFoM1zFYrV+rHYHz9v/WodohOO/9TaVTiqGrJePB33+679CzcmIYf97D7ErK4Yfncyqml4P6cPKCl0seL51Bz67W0lDIvem+qUE1MZz8g89zB7ccxg3mDm5eSmlAT8uyf72iHYKzCiuJ4bPsNhNDfkzmi7JbkVihHYLzEikmBi9VI4dqh+CsgorodWOqieG4wQNQy5XPVk98vEw7BOeN7t9POwTnJT59UzsEZxWy8rkiHENFOP+u8BTXMeSnMuJvo4ueZkNTmXYIzusV5aw2TwnucmdVQGJoam1EWWv+b6FNJdTtqZoYWk0KrZyVZNXQwsV/Xo688nntEJzHkhh2pVASQ0SmAfgpgH4A3gVwrjEm6w++iHwPwBQAW+Yovw7gImOMr4Em1XceQQiC4HuudncsSe6tqjdbDF42N/IDhk1jAfdGxOf6tgB5QUQmAvgdgGkAXgRwBoAnRGRPY8zyLKeMAnAfgH8BaARwPoAnRWQvY8zKfK+r+hMz9rL5LLudw4yLuTuZFxZ78FYV5QcMm1BB96ZLWgzTAdxhjLk98/dzReQoAGcBuLDzwcaYSR2umG5BTABwBICZ+V5UNTGcd8ZAVNZUaYbgtKkXvqwdgvNYXdUbWwx2jc0FtBgQOC3UdJrm2mSM2WoJtohEAewP4OpOTz0J4OA8L1sJoAzABh+h6iaGz/ceyP0YcphzHWslUeFOmD5POwRnFbyOIVhJjM7z0C8HcFmWU3YEEAbQuQTCWgA753nZqwGsBPBUnscDUE4MiVQCoVTpLProavNXls4sBi2z/swiel5O+/EXtENwVjIexyt3Bjs3GqlAeaQi7+ObI//rthoIoL7dU14Fmzov95Isj21FRM4H8G0Ao4wxvt5MVBPD/JWNiG1mM9emqTX/kr49VdWg/tohOK+ljT9HNoXcm+bWJjS15j+lvLn1f+//9caYujxO+QRAK7ZuHfTF1q2IDkTkPAAXATjSGPNW3kFmqL4rjxlQgZra/DNuT7ND+Q7aITivNtqiHYLz/jq/3vugHqqQlc/Frq5qjGkWkdcBjAHwt3ZPjQHwkPU6Ij8FcDGAo4wxr/m6aIZqYuhf2Q+1lZxuaJNsZVVML6P7c1ablxnvv6EdgrMKG2PokhJ61wOYJSKvAXgJwPcBDAJwKwCIyEwAK40xF2b+fj6AKwB8B8BSEdnS2ogbY/Je6aiaGP60YDkqWETPavft+EnPy1V/3qgdgvM4c8suGY/jlTuCnl381GCMmS0ivQH8HOkFbu8A+LoxZku9nEEA2tqdMg1AFMADnV7KNsCdlWpiiIZbUV5C9UO62uj+3GDFy+7TV2mH4LwpV7BWkk0htZK6YoEbABhjbgZws+W5UZ3+PiTYVTpSTQxPPBtHpIL1VW3u+8Ns7RCcN/bML2uH4DyWxLArhZIYGlQTw6579EG0il1JNudM7aMdgvOmnP+CdgjOm3zOAdohdEvdNy0oJ4Yzv2BQU8sWg81Zd63XDsF9hj8/XnapCt5d0t01tAUffI6Gy1Ee9rGOIVw6M+hUE8MZ17zNWkm51OUz1blnG/+DL2mH4LyPExzHs0k2tHkfZNHc2tx+bUJex5cK1cQw+vi9OSsph52q8/+h66lm3fCqdgjO++FPOInBpqysgDfrrhp9VqCaGPrVJBGr5qpMKkA1a215iXNfD6tkSqWInvNUf2KOGdiHRfRyOOHSF7VDcB8nL3iq5r4eVqFI8HtT7JXPmlQTQ020CrVR/mLbcJqhNxbR87a8brB2CM5qjAcfY+jObQbVxBBCGCGwK8lm1i2+a1/1OFUjh2qH4LxeFWwx2ERTBbQYuu8Qg25iaEMr2sAZEzb/uP547RCcN2PREu0QnDdrdt47OvY4hax8ZouhSOqbE0Az93y2mf5wzsq6BOCT9Zyj72XscexKsmmMx/FK5/3R8hQNRRENlfs4ntNV83LyVS9zHUMOp313hHYIzhvVj6vDvXChpF0hZbdb2lrQ0pb/m31LGxe45WXCpJGIcR2D1ah+3KvCy5SfPq8dgvNYT8quMV6GZwKem+5I8jMrqXSoJoZYpBWxCMcYbLi1p7dH/3CcdgjOO/bioG993V9BYwzdd4hBNzG0tIUR4baDVIBjf2TdyIoypv54f+0QnFVIdVWuYyiSbw0biNpa7uBmc/TZD2qH4LyZ1x6qHYLzOMZgl0oWsB8DE0Nx3P/hCo4x5DB+2he1Q3Aexxi8scVgl4ybwGMM7EoqksXrqlCeqNIMwWmv/nOhdgjOm/PbY7RDcN6UOz7SDsFZbDFkp5oYhu2YQKy6dG5WV/vRdE5X9ZJI5b2/eY+VWMTEYJNqDr4fQ1koirJQ1NfxpUI1MUzddRjHGHIYfcFj2iE474dnsiSGl8EH76UdgrNaGuJ45Y/Bzk21tSDlY22Cn2O1sR6vw6oG9dcOwXk3/uYd7RCcN/ns/bRDcFYyHsYDAc/txkMMuolhZWIV6sLcpYwKEOFnGy8RKaSCaPdW0L3pxlX0VH+rekV7oSbK/RhsEote0Q7BfTv01o7Aec2G9chsCrk3HHwukqpIJaojnJVkM/ZUdgF4YUlpbw8+Va8dgrMKqZXErqQi2di0CakmlsSw4ZsekcPYlVQcnHEFeAAAELNJREFU81d9glgd6wHZPPgUS0p7SXzKMSovkycO0A7BWcm4FFBEj11JRdFmBG2mdG5WV0u8v0g7BOdV7TFSOwTqoSISQUTyfwv1c6w21Ui/PXwE1zHk8PYEjr9Q4Z57k4sAbVoagndlt5oUWk3+3b1+jtWmmhiOvfTv3Kgnh8mn7aYdgvOmjhyuHYLzjpj+N+0QnFXIyufuPPysmhgmnzKMRfRymP0kZ5N4efCpN7RDcB5nt9k1xoOvfO7GY88cY3DZ+CO5xsNLr2jp7KOr5c7HN2qH4KzCpqty8LkoWk0IrVx8Y7WpsXQGq7RMHMayIV5ufPdh7RCcxa6k7FTfeXbr1YSqGr752Vz0K5bd9vIIVz57mnzOAdohOCsZj+OVO7WjcI/qu/Ll1yxAJBrTDMFpLH7m7dRdWZrcy+iz/6IdgrMKaTGICMTHwIGfY7WpJoavn7I/Kjj4nEPpTG/TMnrafdohOG/smV/WDsFZjfE4Xrkj2Llhn+sYwlzHkJ/tY82IxTh4aPP1XXbQDsF95x6oHYHzZs1eph2Cs1JNwasLcB1DkRw3eABqaznzxubYHz2kHYL7wmHtCJzHLkm7ZDyOV64Odi5nJRXJUyuXobKOXUlW/TjjxtOGT7UjcF6Y+zFYFXRvuvFCBtXEsKYhhooQB59txo8fpB2C8/7+LEuqeHl7FUur2DQnTOBzu+9kVeXEUBNtQSxaOvugdrW/zg+++KanSLz7gXYIztuTCyWtkvECZiWxK6k47rn9fU5XzYGlDLydcvrR2iE474RLX9QOwVmFDD4zMRTJpO/uwVpJOcz47UvaITjvYVbn9TT1dJYmt0nG43jlhoAnd+O+JNXEcNKwYSy7nUPruSwXQoWbcf9K7RCcVUiLISxhX2sTwlI6M+hUE4NICCJ887OZMnKIdgjOu+C5tdohuC/B/RisCkgMbaYVbT7WJrSZ0tnGWHeMYckSdiXlMPP+VdohuG/Nau0InMdaSXaFdCVxjKFI7puzhhv15DDn/H21Q3De3GXcz9jLrL+s0Q7BWanGAmb+CXyuYwh+qa6mO8Ywvh9bDDmc8Jt3tUNw3tMXHaIdgvNmrXpVOwR3FVJED9127Fk3MYSljasyc9j/kCHaITjvzoWLtUNwXtXeu2uH4Cxu1JOdamJoag0h1MrBZ5ujRrDcg5errmNi8LL/uC9qh+Cs5kQczwQ8l4mhSMrDbSgPs8Vgc9Uv2ZXkqbxcOwLn7bZzIbuUdW+FrHzuzn1JqonhsJ0rUVPLwWebGVzj4S3BsiFe4k2lsw9AV2ss4N7E65II+XgLjdeVToJW/YnpV7kzaiv55mf3pnYAzmNJaW8PPrVROwRnBRxjaAawZsSQ3XcOcO6azPlOU00Msz/8kLOScphzOWfceJm77BPtEJyXWPSRdgjOCrK1pzGmUUSGAogGuGSzMaYxwHldikX0HDb6yoO0Q3Dee6v58+OlauRQ7RCcFXRWUubN3fk3+KBUE8Ojv/omayXlwP2MvU3m1p6eFnAOA/kkxgTfqCLwRUVqAWw+6Ow/c+VzLlXsZvO0mmVDvJw2neMwNsl4HGcfMAoAehlj6pTDcQanKxB1c22mhOZJdjHem+x0S2JMHsnB5xymjhyuHYLzjrjyOe0QnMedAO1SyeDVVbsz1cSwcF0VyhPcj9Zm0+DN2iE4r6o3x6iItjXVxGDa2tDWxpXPNg8tX68dAnUDO/bhOJ5NSwPff7JRTQxnfsGgprbrB79LxSVP8N54+dYYtji93Hk/96ywKWQHt+5MNTGccdNHnJWUw8zpHGPw8tSq0ikzoIY7uNkxMWSlmhguPXUHVNXwE5/NlAv+pR2C81gSw1vVoP7aITirkLLb3Znu4HNdFLE2Vse0qmLSJKKup5oYJg0fzpXPOYRO5xiDl0QL9/Pwknh/kXYIzgpSK6knUE0Mx05/gLWScvGzn2wP9fRNJ2mH4Ly/tv5bOwR3tbZqR+Ak1cRw96/GoKa2RjMEp7FyqLdZixZqh0DU7agmhkQqgVCKn4opOHa25YEtTzvem6xUE8P0ez9BpIJ9fDaJxayj72XO1aO1Q3DejBDHYax4b7LS3ajne5/j4HMOo6ct0A7Bea9/slw7BPdxdptdhC2GbFQTw31LFrOIXi79OP/cS3UZP/ERbWu6O7jd9AZnJeXw6B+O0w7BeZe+yP2MvS3TDoBKjG7Z7Wn7scWQw+wPuQmNl/13SWmH4LzX67j/jBXXMWSlmhjun7WItZJy4Q5unlhPyttt5awuYCWsrpqNamKYfMowthhy4O5S3p5exSJonpqatCNwVzPvTTaqiWGnimZUxvgPY/P7uVyV6eWIQ7fXDsF9nMRgx+qqWakmhlsea0Ukxjc/m8S7nK7qJTLqS9ohOI+73NmlkmHtEJzEdQwOa0jtph2C84796TztENzHMQYrbtSTne501SVLOMaQw5z53GDFy+SzPqcdgvNm/fkD7RCoxOgucLv1La5jyGHm1Qdrh+C8KVe8qR2C87hRjx036slONTGgT1+A01Wt+sb6aofgPk7p9ZT4lOsYbNiVlJ1uYlizGmCLwequhawD5IUDq1QIDj5np5sYQiFWN8yhPMzFN15OHMMCcV5m/PZ17RCcxR3cstNNDP0GABX8xbapLmO5By83zlqvHYLzqvbeXTsEZ3GMITvdxNCaAlItqiG47M7HufjPy3lTe2uH4LzrrucudzZsMWSnmhj++IOhqKnl4KENZ9x4W5vkWg9Pzc3aEbiL9yYr1cTw3JoGxOIcY7BiVUxPKQ7DeKvhvupWTRx8zkY1MWxMRpEMRzVDcBr7hr1VcIDeGz9g2LErKSvVxPDsS5sQiXGMgaioWBLDjmW3s9IdfKacJn6Ni/+8zH6Ss0o87cABeisucMtKNTF886vVrJWUQws/zHjasS+nO3tJvL9IOwR3sSspK90d3G75L2sl5fDtsw/QDsF5y17jm56Xqj1GaofgLK5jyE41MTxy3XiW3c7hjoVLtENwH2sleUos+kg7BGdxHUN2ui2GD1l2O5d77+SnYS9jv7OXdgjOe/hPa7VDcFeK1QWy0d2P4aY32JVEBdkhxgVKnqo4DmMV4b7q2agmhkevn8CupBxu/2CxdgjOm/HbV7VDcB9/x8gn1cRw9+IP2ZWUw+BqTqXzxDc9om1Od/D5mTgiMaMZgtM4aOiNW3t6m/U7tqqsOPiclWpiOH50JWLV7P+0MWP21g7BeTPuX6kdgvu48tmOK5+zUk0Mw2tbUFXDwUObiy57WzsE5+0/8TDtEJz3+uxV2iG4q5ml7bNRTQxL6soQa2MRPZujv3eIdgjOm/coB+g9scVgxxZDVqqJ4Qs7AtWsCGx1zRvs//SUiGtH4D4mhhxatQNwku6spLerEOUca6tlC/hp2BNXPlMhItwPJhvVxPDf+W8hUs4KohTcgd/YTzsE57163z+1Q3AXZyVlpVt2u74eaOKSdJvv/YRvel7+9Jt/aofgPq71sGvizgPZqN6VL510KMrZFWDV0lavHQJ1B9zBzY4thqxUE8PLf3uFtZJyGDmN6xhoGwhzX2Mr3pusVBPDCScNQgUXuFnd/usXtEOg7qCGU/9sWttYRC8bMabrS1KIyAAAK7r8wkRE2Q00xnAZfYZWYhAA/QG40oleg3SiGgh3YnIN75E33qPcXL0/NQBWGY03Q0epdCVl/gGcyc7pPAUAqDfGcKQuC94jb7xHuTl8f1yKxQlc3UFERB0wMRARUQdMDGlNAC7P/EnZ8R554z3KjfenRKgMPhMRkbvYYiAiog6YGIiIqAMmBiIi6qDHJwYRmSYiH4lIo4i8LiKHasfkEhE5TEQeEZFVImJEZJx2TC4RkQtF5FURqReRdSIyV0R2047LJSJyloi8JSJ1ma+XROQY7bjIrkcnBhGZCOB3AH4J4AsAngfwhIgMUg3MLVUA3gTwQ+1AHHU4gJsAfAnAGKQXjT4pIiwC9pkVAC4AcEDm62kAD4nIXqpRkVWPnpUkIv8G8B9jzFntHnsfwFxjzIV6kblJRAyA440xc7VjcZWI9AGwDsDhxpjntONxlYhsAPBTY8wd2rHQ1npsi0FEogD2B/Bkp6eeBHBw10dE3USvzJ8bVKNwlIiEReQkpFuiL2nHQ9n15O2LdgQQBrC20+NrAezc9eFQqcsUh7wewAvGmHe043GJiOyDdCKoABBHuuX5nm5UZNOTE8MWnfvSJMtjRPm4EcDnAHxFOxAHfQDg8wC2A3ACgLtE5HAmBzf15MTwCYBWbN066IutWxFEOYnIDQDGAjjMGMO9RjoxxjQDWJz562siciCAcwCcoRcV2fTYMYbMD+rrSM8kaW8MgH91fURUiiTtRgDjAYw2xnykHVOJEADl2kFQdj25xQCk+4NnichrSPd/fh/AIAC3qkblEBGpBjCi3UNDReTzADYYY5YrheWSmwB8B8BxAOpFZEsLdLMxhjvNAxCRqwA8AeBjpDfFOQnAKABHK4ZFOfTo6apAeoEbgPMB9APwDoAfc5rhZ0RkFIBnsjx1lzFmatdG457MFN5sTjXGzOjKWFwlIncAOALp37HNAN4C8GtjzHzVwMiqxycGIiLqqMeOMRARUXZMDERE1AETAxERdcDEQEREHTAxEBFRB0wMRETUARMDERF1wMRAREQdMDEQEVEHTAxUVCIyI7NXtBGRFhFZKyLzReQ0EeHPH5GD+ItJXWEe0nVyhgA4BunaS78H8KiIqBVyzFRG7emFJIm2wsRAXaHJGLPGGLPSGPMfY8xVSFcjPQbAVAAQkUEi8pCIxEWkTkT+IiI7bXkBEblMRN4QkckislRENovI/SJS0+6YchH5g4isE5FGEXkhU/d/y/OjMi2XozIVdZsAHCoi+4rIMyJSn7n26yJyQFfdHCLXMDGQCmPM0wDeBDA+syXmXAA7ADgc6T0xhgOY3em04QDGATg283U4gAvaPX8N0ruDnQJgP6Q3hvm7iOzQ6XWuAXAhgD2QrvR5D4AVAA5Eeh/wqwG0bIvvk6gUsRlNmhYgvRXmkZk/hxpjPgYAEZkM4F0ROdAY82rm+BCAqcaY+swxs5Au5/wzEakCcFbm+Scyz38P6SRzOoBr21335+1LPovIIADXGmMWZB5aVJTvlqhEsMVAmrbsr70HgI+3JAUAyOwFvCnz3BZLtySFjNVIb8UKpFsTZQBebPcaLQBe6fQaAPBap79fD+B2EXlKRC4QkeHBvyWi0sfEQJr2APARPksQnXV+vHP3jsFnP8PS7rFcrwEAiQ4vYsxlAPYC8BiA0QDeE5HjvcMn6p6YGEiFiIwGsA+AOQDeAzBIRHZp9/yeAHoBeD/Pl1wMoBnAV9q9RhmAA/J5DWPMQmPMb40xXwPwIIBT87wuUbfDMQbqCuWZvZDDAHZCeq/fCwE8CmAmgDZkBoFF5Fykfy5vBvCsMaZzt09WxpiEiNwC4FoR2QBgOdJbtlYCuMN2nojEkB5/eADp1stApAeh5wT4Pom6BSYG6gpHIz0ekAKwEenZSD9Cet/oNgAQkXEAbgDwHNKJYh6As31e5wKkW8GzkN50/jUARxljNuY4pxVAb6QT1E4APkG6xXCpz2sTdRvc85mIiDrgGAMREXXAxEBERB0wMRARUQdMDERE1AETAxERdcDEQEREHTAxEBFRB0wMRETUARMDERF1wMRAREQdMDEQEVEHTAxERNTB/wN0wMMJ0B3APgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 400x500 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from vireoSNP.plot.base_plot import heat_matrix\n",
    "\n",
    "fig = plt.figure(figsize=(4, 5), dpi=100)\n",
    "im = heat_matrix(AF_SNPs, cmap=\"GnBu\", 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(\"Donors\")\n",
    "plt.ylabel(\"%d SNPs\" %(AF_SNPs.shape[0]))\n",
    "plt.yticks([])\n",
    "plt.xticks([0, 1, 2, 3])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choose the number of donors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unsually, the number of donors are known. In case you are not sure on this, you could try a set of numbers, e.g., from 2 to 6 here, and pick the one when the evidence lower bound (ELBO) stops increasing, e.g., \n",
    "[Figure 2A in Vireo paper](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1865-2/figures/2)\n",
    "\n",
    "In this data set, `n_donor=4` is a sensible choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[vireo] lower bound ranges [-42583.3, -34775.1, -33917.3]\n",
      "[vireo] allelic rate mean and concentrations:\n",
      "[[0.064 0.437 0.849]]\n",
      "[[24679.1 40819.5 22355.4]]\n",
      "[vireo] donor size before removing doublets:\n",
      "donor0\tdonor1\n",
      "483\t469\n",
      "[vireo] lower bound ranges [-33921.1, -27814.6, -25535.5]\n",
      "[vireo] allelic rate mean and concentrations:\n",
      "[[0.026 0.445 0.91 ]]\n",
      "[[25407.2 40997.2 21449.6]]\n",
      "[vireo] donor size before removing doublets:\n",
      "donor0\tdonor1\tdonor2\n",
      "241\t242\t469\n",
      "[vireo] lower bound ranges [-29456.1, -25592.6, -19155.3]\n",
      "[vireo] allelic rate mean and concentrations:\n",
      "[[0.017 0.462 0.941]]\n",
      "[[28199.7 38020.2 21634.1]]\n",
      "[vireo] donor size before removing doublets:\n",
      "donor0\tdonor1\tdonor2\tdonor3\n",
      "241\t246\t234\t231\n",
      "[vireo] lower bound ranges [-29997.3, -24381.9, -18709.2]\n",
      "[vireo] allelic rate mean and concentrations:\n",
      "[[0.019 0.457 0.929]]\n",
      "[[28304.3 37007.5 22542.2]]\n",
      "[vireo] donor size before removing doublets:\n",
      "donor0\tdonor1\tdonor2\tdonor3\tdonor4\n",
      "51\t233\t223\t236\t209\n",
      "[vireo] lower bound ranges [-28605.3, -23995.5, -19379.3]\n",
      "[vireo] allelic rate mean and concentrations:\n",
      "[[0.027 0.452 0.911]]\n",
      "[[28677.7 35500.7 23675.6]]\n",
      "[vireo] donor size before removing doublets:\n",
      "donor0\tdonor1\tdonor2\tdonor3\tdonor4\tdonor5\n",
      "240\t165\t76\t131\t105\t235\n"
     ]
    }
   ],
   "source": [
    "n_donor_list = np.arange(2, 7)\n",
    "ELBO_list_all = []\n",
    "for _n_don in n_donor_list:\n",
    "    res = vireoSNP.vireo_wrap(AD, DP, n_donor=_n_don, learn_GT=True,\n",
    "                              n_extra_donor=0, ASE_mode=False, fix_beta_sum=False,\n",
    "                              n_init=50, check_doublet=True, random_seed=1)\n",
    "    ELBO_list_all.append(res['LB_list'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFtCAYAAADFxlkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXgV5dnH8e+dlSULKKsiggImChUIi6Kt+tbdtuLSuhC3CqKttmpbrbav2k3batW69HWvFXCrWrW0dal1qyJIABUMEpEIyL5lIXvyvH/MHDwcsnPCnEl+n+s618mZuc+cO2nld56ZZ2bMOYeIiIgktqSgGxAREZGWKbBFRERCQIEtIiISAgpsERGREFBgi4iIhIACW0REJAQU2CIiIiGgwBYREQmBlKAb6KrMzIB9gLKgexERkUBlAmtcC1cyU2AHZx9gddBNiIhIQhgEfNFcgQI7OGUAq1atIisrK+heREQkAKWlpey3337Qir2tCuyAZWVlKbBFRKRFmnQmIiISAgpsERGREFBgi4iIhIACW0REJAQU2CIiIiGgwBYREQkBBbaIiEgIKLBFRERCQIEtIiISArrSmYhIwJxzbC4po+CDxZRX11FeVec9V9exvbref/ZebyvbzvZN6xg8ZDD9emWxV89U9uqZ5j16pJHdI43kJNux7ZycHHr06BHgbyfxosAWEYkD5xzba+opqaylpKKWkspaSqv850rvubFHZF35F0Ws+8uVce/rwluf5CuHjqF/Vjf6ZaZ7z1np7N0zjZRk7WQNE2vhbl7SQcwsCygpKSnRtcRFEoRzjrLqOkoqWgrbul0Ct7SylrqG9v972lBbRcPWL8hITyEjPYWe3VJ2/JyRnrzj55K1xdx745Wce+2tpPXZj83lNWytqGbL9lq2VdQQ20Lq3oNISu22y+clGeydkf5liGem0y861DPT6ZeVTp+MdFIV7B2mtLSU7OxsgGznXGlztQrsgCiwRTpGQ4MXus2NamPDNvr1bmQuAKnJRnb3VLK6p5Id9cjqFvM6en0P77lnWjJm1uz2FyxYQF5eHgUFBYwdO3andXX1DWzZXsP60mo2lFXt9LyxrIoNZdWsL61iU3kN9a38Rc1g755p9Mv0Rub9/efYcO+bqWBvj7YEtnaJi0jCqW9wlFW1LWwju6LLquvY3XFIWkrSTuHqBW5K42EbFbjZ3VPpntpy6HaUlOQkL0izugHZTdbVNzg2b69mgx/oG0qrd4T7hrJqNpT6z2XV1Dc4NpXXsKm8ho/XNv/5e/dMo29UiEd2v0fCPhLs6SnJ8f3FuwgFtojEXUVFBUuXLm22Zn1pJc/MXcGK4mKSsvpRWZ/E9povJ1rFamrXblO6pSY1OsLdJWxjAje7eyrdUjt3oCQnmReimc0He0ODY0tFDetLo4K8tHrHSD2ybGN5NbX1js3ba9i8vYal65q/tXPvHqn0z+q2a7hH7Zbvl6Vgj6XAFpG4W7p0KXl5eXHd5rHX/Zn9D9o/JmxTdgnbSCDrH/vdl5Rk9MnwjmMf0kxdQ4Nja0XNTkG+MfJzaTXr/VH8xrJqauob2FpRy9aK2haDvVeP1B1h3li4R5Z39i9YEQpsEYm7nJwcCgoKGl336YYybnhhCSWVtfSp28jiWb/hqt/czUE5Od7Eqshkq7QUUlOSdtqmTk9KTElJxt4Z6eydkU7uwKbn5Djn2FZRuyPAIwG/MSroI881dQ1sq6hlW0Uty9aXN/v52d1Td4zK+2d2o2/UsfYdI/fMbnRPC3ewK7BFJO569Oixy4QogPnFW/jVC+9Tlb0/4w/J5uoxqRwz6zfknzip0XrpXMyM3j3T6N0zjZwBTdc55yitrNsR7NFhHhvulZUVbFj3KRta+GxXV8OEPvVc++2v0r1792ZrE/XLoQJbRPaIt4s2csljBVTW1jNhyF48fOE4ij7+KOi2JAGZmXeoo0cqI/pnNlnnnOOtOe9z9BFntmq7fwf+/seW6xqbgZ8IFNgi0uFeWbKOyx9fSE19A0eN6Mt9+Xmh3z0pwTMzxo8e2eThlwjnHAs/WsK0iy5g5syZ5ObmNlufk5MTzzbjRoEtIh3qhUVfcPXTH1Df4Dhp5ADuPHu0JoRJ3DR1+CVW5FS73NzchBw9t0YoznI3syFm9rCZrTCzSjNbbma/MLO0mLpRZvamX/OFmd1gMSdEmtkZZvaxmVX7z6fFrDczu8nM1vjbecPMDomp6W1mM8ysxH/MMLNeHfcXEAmnx+eu5MqnFlHf4Dh97L7cfc4YhbVIO4UisIEcvF6nA4cAVwGXAjdHCvwrh70KrAHGA1cAPwaujqo5HHgKmAEc6j8/bWYToz7rGv89l/vbWQe8ambRB1IeB0YDJ/qP0f62RMT3wFvLuf5vH+EcnH/4/tx25qG6drXIbgjFLnHn3EvAS1GLPjOzg4DL8EIZYArQDbjQOVcNLDazEcDVZna7867BeiXwqnPuFv89t5jZUf7yc/zR+JXAb5xzzwGY2QXAeuBc4H4zy8UL6cOcc3P9mmnAHDM7yDn3SUf9HUTCwDnHHf8u4q7XigC47OgDueaEgwK7+pdIZxHmr7vZwJao14cDb/phHfEysA8wJKrmlZjtvAxM8n8eCgyIrvG392ZUzeFASSSs/Zr3gJKoml2YWbqZZUUeQNNTH0VCyjnHr2YX7gjrn5xwENeemKOwFomDUAa2mR2It8v7vqjFA/BGwtHWR61rrmZATF1LNY2d8rchqqYx1+GFeuSxuplakdCpb3Bc99xHPPLOCgB+8a1D+P4xwwLuSqTzCDSw/cldroXHuJj37IO3e/yvzrmHYjYZe8l/a2R5YzWxy1qqaezWAo1tJ9oteHsFIo9BzdSKhEptfQM/fHIhT76/iiSD2759KBdMGhJ0WyKdStDHsO8Bnmyhpjjygx/WrwNzgEti6tax6wi3n/+8voWa6PX4NWubqenfSJ992XVkvoO/a33H7nrtIpTOoqq2nu/PWsBrSzeQmmz88ewxnDxqYNBtiXQ6gQa2c24TsKk1tWa2L15YFwAXOecaYkrmADebWZpzrsZfdjzerPHiqJrjgDui3nc88K7/8wq8QD4OWOh/bhpwFHBt1DayzWyCc26eXzMRb9T8LiJdyPbqOqY9Np93l28mPSWJ+87L45iD+rX8RhFps1Acw/ZH1m8Aq/Bmhfc1swFmFj1afhxvBPuomY30z6++HojMEAf4I3C8mV1rZjlmdi1wLHAngF93J3C9mZ1mZiOBR4EKf/s45wrxdsk/aGaHmdlhwIPAbM0Ql66kpKKW/Ifn8u7yzWSkp/CX705QWIt0oKB3ibfW8cAw/xE7WcsAnHMlZnYccC8wH9gK3O4/8GveNbOzgV8DvwKWA2dFz/gGfg90B/4E9AbmAsc756LvAzcFuIsvZ5O/iHfetkiXsKm8mvMenkfh2lKyu6fy2HcncOh+unaQSEcKRWA75x7FG+m2VPcR8LUWap4BnmlmvQNu8h9N1WwB8lvqR6QzWrOtkvyH5vLZpu30yUhn5tQJ5Axo+paKIhIfoQhsEUkMn2/ezrkPzuWLbZXs26s7M6dOZGifnkG3JdIlKLBFpFWWrS8j/6G5bCirZmifnsycOpF9ezV/X2ERiR8Ftoi06MPV2zj/kXlsq6glZ0AmMy6eSN/M9KDbEulSFNgi0qx5K7bw3Uffp7y6jkP368VfLhpPrx5pLb9RROJKgS0iTXpz2Uamz5hPVW0Dhx2wFw9dMJ6MdP2zIRIE/ZcnIo16afFarnhiIbX1jmMO6sv/5efRLVX3shYJigJbRHbx3ILV/OSZD6lvcJwyaiB3nDWatJRQXGdJpNNSYIvITmbMKeZ/X1gCwLfzBvHbM75CcpKufS8SNAW2iOzwf28s53cvLQXgwklDuOEbB5OksBZJCApsEcE5x22vfMK9ry8H4Ir/GcbVx43QXeVEEogCW6SLa2hw/HL2xzz6bjEAPz0ph0uPOjDYpkRkFwpskS6svsHx02c/5K8FqzGDX546kvMO2z/otkSkEQpskS6qpq6Bq55axD8+WktyknHbt7/CaWMGBd2WiDRBgS3SBVXV1nPpzALe+GQjqcnG3eeM5cSRA1p+o4gERoEt0sWUV9dx8aPvM3fFFrqlJnH/eeM4akTfoNsSkRYosEW6kG0VNVzw5/f5YNU2MtNTeOSi8YwfslfQbYk0qqioiLKysrhsq7CwcKfn3ZWZmcnw4cPjsq3WUmCLdBEbyqo4/+F5LF1XRu8eqTz23YmMGpQddFsijSoqKmLEiBFx325+fn7ctrVs2bI9GtoKbJEu4IttleQ/NJcVm7bTLzOdmVMnMqJ/ZtBtiTQpMrKeOXMmubm5u729yspKiouLGTJkCN2779593AsLC8nPz4/b6L+1FNgindyKTduZ8uB7rCmpYt9e3Xl82kT237tn0G2JtEpubi5jx46Ny7aOOOKIuGwnKApskU5s6bpS8h+ax6byag7o25NZUycyMHv3RhciEgwFtkgntWjVNi54ZB4llbXkDsxixsUT6JORHnRbItJOCmyRTui9zzZz8aPvs72mnrGDe/HnCyeQ3SM16La6tHjNeO4Ms52lfRTYIp3M659s4NIZBVTXNTDpwL158Pxx9EzXf+pB6ogZz2Ge7Szto/+KRTqRf3y4liufWkhtvePY3H7cc+5YuqUmB91WlxfPGc+dYbaztI8CW6ST+Ov8VVz77Ic0OPjmoftw+3cOJTU5Kei2JEq8ZjyHfbaztI8CW6QTePSdFdz0948BOHv8fvzmtFEkJ3XMvax19SmRYCiwRULu3tc/5daXPwHg4iOH8vNTcjHruLDW1adEgqHAFgkp5xy/f/kT/u+N5QD88OvDufLY4R0W1qCrT4kESYEtEkINDY4bX1zCjPc+B+BnJ+cy7WsH7LHP19WnRPY8BbZIyNTVN3DNsx/y3IIvMIPfTB7FuRMHB92WiHQwBbZIiFTX1fPDJxbx0pJ1JCcZt3/nUE4dvW/QbYnIHqDAFgmJypp6ps8s4K1lG0lLTuKec8dw/CEDgm5LRPYQBbZICJRV1XLxo/OZV7yF7qnJPHj+OI4c3ifotkRkD1JgiyS4rdtruODP8/hwdQmZ3VJ49KLx5O2/V9BticgepsAWSWAbSqvIf3guy9aXs1fPNB777gRG7psddFsiEgAFtkiCWrWlgvyH5/L55gr6Z6Uza+pEhvXLDLotEQmIAlskAS3fWE7+Q3NZW1LFfnt15/Gph7HfXj2CbktEAqTAFkkwH68p5fxH5rKpvIZh/TKYefFEBmR3C7otEQmYAlskgSxYuZULH5lHaVUdh+yTxWPfncDeGelBtyWyx1ldFWMGJNF92zJYk1h3neu+bRljBiRhdVV79HMV2CIJ4t1PNzH1sflU1NQzbv/ePHLReLK6pQbdlkggupWvZMH0DHhrOrwVdDc7ywUWTM+gsHwlMGmPfa4CWyQBvFa4nstmLaCmroEjh/XhgfPz6JGm/zyl66rKGMzY+8uZNWsWuTk5Qbezk8KlS5kyZQoPn7xnLwmsfxFEAvb3D9Zw1VOLqGtwHHdwf+4+ZwzdUpODbkskUC6lGwvXNVDZawTsMzrodnZSua6BhesacCl7dm6JAlskQE+9v5KfPvcRzsHk0ftw67cPJTU5sY7XiUhiUGCLBOTh/67gV7M/BuDciYP59akjSUrquHtZS3ASdQJVUJOnpH0U2CJ7mHOOu//zKbe/ugyAS752ANedlIOZwrqzStQJVEFNnpL2UWCL7EHOOW7511IeeOszAH503Agu/59hCutOLlEnUAU1eUraR4Etsoc0NDh+/sJiHp+7EoD//cbBXHzk0IC7kj0hUSdQBTV5StoncQ6mNMPMhpjZw2a2wswqzWy5mf3CzNJialwjjxNjtnWUmRWYWZWZfWZmlzbyed/zP6vKr/1qzPp0M7vbzDaZ2XYze9HMBnXcX0DCrq6+gaufXsTjc1diBr87Y5TCWkTaJBSBDeTg9TodOAS4CrgUuLmR2mOBgVGP/0RWmNlQ4J/A28AY//13mdkZUTVnAXcCv/Fr3gb+ZWbR+4zuBE4DzgaOBDKA2Wamc3FkF9V19Xxv1gKeX7SGlCTjrrPHcNZ47YIUkbYJxS5x59xLwEtRiz4zs4OAy4Afx5Rvds6ta2JTlwIrnXNX+q8LzWycv41n/WVXAw875x7yX19pZif4n3WdmWUDFwPnOef+DWBm+cAqvC8LL7f395TOp6KmjukzCni7aBNpKUn835SxfD23f9BtiUgIhWWE3ZhsYEsjy180sw1m9o6ZnRmz7nDglZhlLwPjzCzV38We10jNK3w5hTIPSI2ucc6tARbTzDRLfzd6VuQB6D6JnVxJZS3nPzyPt4s20SMtmUcvGq+wFpF2C2Vgm9mBwBXAfVGLy/FGx2cCJwOvAU/5o9+IAcD6mM2tx9vT0Md/JDdRMyBqGzXOua3N1DTmOqAk6rG6mVoJuc3l1Zz74HvM/3wrWd1SmDl1IpMO7BN0WyISYoEGtpnd1MREsejHuJj37IO3e/yvUbutcc5tcs7d4Zyb55yb75y7AfgTcE3Mx7rYNhpZ3lhN7LJdfp0Wam7B2ysQeWiSWie1vrSKsx54jyVrStm7ZxpPXnI4Ywf3DrotEQm5oI9h3wM82UJNceQHP6xfB+YAl7Ri++8BU6Ner2PXUXA/oA7YjBe69U3UREbd64A0M+sdM8ruB7zbVCPOuWqgOup3aUX7EjartlQw5aG5rNxSwcDsbsycOpED+2YE3ZaIdAKBBrZzbhOwqTW1ZrYvXlgXABc55xpa8bYxwNqo13OAb8bUHA/Md87V+p9TABwH/C2q5jjgBf/nAqDWX/a0/56BwEh2Hc1LF/LphjKmPDSX9aXV7L93D2ZePJH99uoRdFsi0kkEPcJuFX9k/QawEm9Gd9/ICDUyI9zMLsAL0oVAA14w/wC4NmpT9wGXm9ntwIN4k9AuBs6JqrkdmGFm8/lyJD/Yfy/OuRIzexj4g5ltxpv4dhvwEfDvOP/qEhKLvyjh/EfmsWV7DSP6ZzDz4on0y9LFKEQkfkIR2Hij4GH+I3ayVvS+5Z8D++Pt1l4GfNc5NzOy0jm3wsxOBu4Avg+sAX7gnHs2quYpM9sbuAHvPO7FwMnOuc+jPucqvN3oTwPd8Sa4Xeicq4/D7yoJpKKigqVLlzZbs3ZbJT+cNY9tG77goOEHcv2Esaz+9ONmZxXm5OTQo0f4Rt+JehML0I0spPMLRWA75x4FHm2h5i/AX1qxrTeBsS3U/AlvwlpT66vwZqlf0dLnSbgtXbqUvLy8Vte/CxxzR8t1BQUFjB3b7P8NE1Ki3sQCdCML6fxCEdgiQcnJyaGgoKDRdZEbebzz6Sa6Vayj+K+/ZebMmeTm5rZqu2GUqDexAN3IQjo/BbZIM3r06NHkSPje1z9lfnkvMvbtzU2TvkL+X39Lbm5uKEfOrZWoN7EA3chCOj8Ftkg7/Gfpem575RMAfnnqSA5KadXJDiIi7ZZYs0ZEQuCzjeX88IlFOAdTJg7mnAnaBSsiHU8jbJE2KKuq5ZIZBZRV1zFu/97c+M1Dgm5JpFOqqKgAYMGCBXHZXmVlJcXFxQwZMoTu3bvv1rYKCwvj0lNbKbBFWqmhwXH10x/w6YZy+mel86f8saSlaCeVSEeInE45bdq0gDtpWmbmnr2HkwJbpJXu+k8Rr368nrTkJO4/bxz9MjW5SaSjTJ48GYjfNQsKCwvJz89v9ZkcLcnMzGT48OG7vZ22UGCLtMIrS9Zx57+LAPj1aSMZvV+vgDsS6dz69OnD1KlTWy5sozCfyaH9eSItKFpfxlVPLQLgwklD+M64/QLuSES6IgW2SDNKKr1JZttr6pk4dC9+dsru70oTEWkPBbZIE+obHFc+uZAVm7azT3Y37p0yltRk/ScjIsHQvz4iTbjj1WW8/slG0lO8SWZ9MtKDbklEujAFtkgj/vnRWu55/VMAfnvGKEYNyg64IxHp6hTYIjGWrivlx3/9AICpRw7ltDGDAu5IRESndUkXV1RURFlZ2Y7XZZW1XPX0IraWVDF6v14c379/q660FLnyUbyugBTEOZ4iktgU2NJlFRUVMWLEiCbXvwS89Ou2bTM/P3/3moqybNkyhXYnEc/LbHaGS2xK+yiwpcuKjKwjVz565J0VPFuwmvSUJG779qEc0Dej1duK9z+i+fn5O438JdwS/TKbe/oSm9I+Cmzp8nJzc1mdPIDZX6wlfcAw7j5nDN88dJ82b+eII47ogO6kM4jnZTY7wyU2pX0U2NLlLd9Yxv/+dz0Alx51YLvCWqQ5HXGZzTBfYlPaR7PEpcv79exCqmob+NqIvvzkhIOCbkdEpFEaYUuXVVffAMCGsmpGDO7B3WePITnJAu4qsekexSLBUWBLl/XIO8UAdEtN5sHzx5HdIzXYhkIg0SdPgSZQSeelwJYu6bkFq3lh0RcA/Oi44Yzor3/kW0P3KBYJjgJbupwPV2/jp899tOP1pGF9A+wmXHSPYpHgKLClS9lYVs30GQXU1DVw5JAMlg9Iovu2ZbAmceZfdt+2jDEDkrC6qqBbEZEEosCWLqO2voHvz1rA2pIqDujbk+tGdmMsGfDWdHgr6O6+lAssmJ5BYflKYFLQ7YhIglBgS5fxq9kfM694CxnpKTxw3jjKitMYe385s2bNIjcnJ+j2dihcupQpU6bw8MmDg25FRBKIAlu6hKffX8Vjcz4H4M6zRjOsXwYLVndj4boGKnuNgH1GB9zhlyrXNbBwXQMupVvQrYhIAkmcA3ciHWThyq38/PnFAFx17AiOPbh/wB2JiLSdAls6tQ2lVVw6s4Ca+gaOP7g/V/zPsKBbEhFpFwW2dFrVdfVcOrOA9aXVDO+Xwe1njSZJVzITkZBSYEunddOLH7Ng5TYyu6XwwPnjyEjXlA0RCS8FtnRKs+Z+zhPzVmIGd50zhqF9egbdkojIbtGQQzqd+cVbuOnFJQD85ISDOOagfgF3JCIdpaKiYsc17psTuTlMa24SE69L78abAls6lXUlVVw6cwG19Y5TRg3ksqMODLolEelAS5cuJS8vr9X1+fn5LdYUFBQk5KVyFdjSaVTV1jN9ZgGbyqvJGZDJrd/+CmaaZCbSmeXk5FBQUNBiXVtu5ZqTQBdSiqbAlk7BOcfPn1/MB6u2kd09lQfOG0ePNP3fW6Sz69GjR6tHw0cccUQHd9OxNOlMOoXH5nzOMwWrSTK459wxDN478Y4/iYjsDgW2hN6c5Zv55eyPAbjupFy+Oly3yxSRzkeBLaH2xbZKvv/4AuobHKeO3oepXx0adEsiIh1CgS2hVVVbz/QZ89myvYZD9snit6drkpmIdF4KbAkl5xzXPfcRi78oZa+eadx/Xh7d05KDbktEpMMosCWUHv7vCv628AuSk4x7zh3DoN6aZCYinVu7z3sxs17AMMABy51z2+LWlUgz/lu0iZv/6V2t6Oen5DLpwD4BdySxutLVp0T2lDYHtpkNAe4FTgAiBwydmb0EXO6cK45XcyKxVm2p4PInFtDg4Iyxg7hw0pCgW5JGdKWrT4nsKW0KbDPbD3gPqAX+FyjEC+1c4DJgjpmNd86tjnejIhU1dUx7bD7bKmo5dFA2vzltpCaZJaiudPUpkT2lrSPsXwCfACc456qilv/NzO4AXvJrLo5TfzuY2YvAaKAfsBX4N3Ctc25NVM0o4B5gArAFuB/4lXPORdWcAfwKOBBYDvzMOfe3qPUG3AhcAvQG5gLfd84tiarpDdwFfMtf9CJwhQ4LdBznHNc88yFL15XRJyON+87Lo1uqJpklqq509SmRPaWtk85OxAu4qtgVzrlKvFH3SfForBGvA98BDgLOwAvcZyIrzSwLeBVYA4wHrgB+DFwdVXM48BQwAzjUf37azCZGfc41/nsu97ezDnjVzDKjah7H+/Jwov8Y7W9LOsj9b33G7A/XkpJk/GlKHgOzmx+NiYh0Nm0dYe8NFDez/jO/Ju6cc3dEvfzczH4LPG9mqc65WmAK0A240DlXDSw2sxHA1WZ2uz/KvhJ41Tl3i7+dW8zsKH/5Of7o+krgN8655wDM7AJgPXAucL+Z5eKF9GHOubl+zTS8wwEHOec+6Yjfvyt745MN/O4lbwLTjd86hAlD9wq4IxGRPa+tgb0GOARo6hj1SGDtbnXUCma2F15Av+uHNcDhwJt+WEe8DNwCDAFW+DXRwR+pudL/eSgwAHglstI5V21mbwKT8HaxHw6URMLar3nPzEr8mkYD28zSgfSoRZmN1cnOijdt5wdPLMQ5OHv8fuRPHBx0SyIdoiNm1oNm13cmbQ3sF4BbzWyBc25j9Aoz6wf8Dng+Xs3FMrPf4e2q7oE3+e0bUasHsOvof33UuhX+8/pGagZE1dFEzf5RNRsaaW9D1Psbcx3esXFppfLqOi6ZMZ/SqjrGDO7FL049RJPMpNPqiJn1oNn1nUl7Jp2dDCw3s5lA5OvgwXi7jNcBv2ztxszsJloOsfHOufn+z7cCD+OF543AY2b2jahJZS7mvdbI8sZqYpe1VBO7vqntRLsFuD3qdSZN76no8pxz/PjpD1i2vpx+mencl59Hekp8J5lVVFQAsGDBgt3eVltmO7ektSMn6Vw6YmZ9ZLvSObQpsJ1zW/0JWjcDZwO9/FXb8CZi/cw5t6UNm7wHeLKFmuKoz98EbAKWmVkhsAo4DJiD92UhdoTbz3+OjJibqolej1+ztpma/o302ZddR+Y7+Lvqd+yu10ixefe+/ikvLVlHWnIS952XR/+sbnH/jMjux2nTpsV92/GQmamjJl2JZtZLS9p84RTn3FbgMjP7Hl5IAWyMPnWqDduKBHB7RBIvclx4DnCzmaU552r8ZcfjHXcvjqo5jp2PYx8PvOv/vAIvkI8DFgKYWRpwFHBt1DayzWyCc26eXzMRyI7ajuyG1wrX84dXlwHwy1MPYezg3h3yOZMnTwbic4yvsLCQ/Px8Zs6cSW5u7m73lpmZyfDhw3d7OyLSebT70qR+QG8AL9T8oCyPW2dRzGwC3rnV/8U7B/sAvF3vy/ECFLwR/o3Ao2Z2MzAcuB74ZdSXiT8Cb5nZtXjH465xP24AABjZSURBVE8FjgWOjPxOZnYncL2ZFQFF/jYq/O3jnCv0r+r2oJlN97f7ADBbM8R33/KN5Vz55CKcg/zDBnP2hI6bZNanTx+mTp0a123m5ubqeKGIdIg23/zDzC4ys7vNbIr/+hagDCgxs1fNrCNO66oETgdew5uF/QiwGDgqMivcOVeCNzIeBMwH/oR3zHjHcWPn3Lt4u/IvAj4ELgTOip7xDfweuNN//3xgX+B451xZVM0U4CO82eSv+Ns6L56/cFdUWlXLtMfmU1Zdx/ghvbnhG4cE3ZKISMJo66VJfwb8DG/X77lmdiQwGbgBaAB+APwa7zKlceOc+wj4n1bWfa2FmmeIuuBKI+sdcJP/aKpmC9C6KZrSKg0NjqufWsRnG7czIKsbf5qSR1qKbiYnIhLR1l3iFwIXO+eeMLNxeJftPMsPQcxsMXBffFuUruDO14r4d+EG0lKSuP+8PPpmprf8JhGRLqStQ5jBeMeR8U+1qsPbNRzxITAwPq1JV/HyknXc9VoRALecNopD9+vVwjtERLqetgZ2KlGnJgE1eHfuiqgDdEcGabWi9WVc/dQiAC6cNIQz8gYF3JGISGJqzyzxg80sci6zATlmluG/7hOftqQrKKn0Jpltr6nnsAP24men7P7pUCIinVV7Avs1vjwHGmC2/+xo+WpfIgDUNzh++ORCijdXsG+v7tx77lhSkzXJTESkKW0N7KEd0oV0OX945RPe+GQj3VK9SWZ7Z2iSmYhIc9p6adLPm1tvZr2BbwKP7U5T0rn948O1/OmN5QD87oyvMHLf7IA7EhFJfPHeBzkY+HOctymdSOHaUn781w8AuORrB3Dq6H0D7khEJBx00FD2mG0VNVwyYz6VtfUcOawP15xwUNAtiYiEhgJb9oi6+gaueGIhq7ZUst9e3bn7nDGkaJKZiEir6V9M2SN+//InvF20ie6pyTxw3jh690wLuiURkVBp67XEf9BCiQ5Iyi5eWPQFD7z1GQC3fftQcgdmBdyRiEj4tPW0rqtaUbOyPY1I57T4ixKueeZDAL539IGc8hVduVZEpD3aelqXzsOWVttcXs30GQVU1zVw9EF9+dHxmmQmItJebTqGbWb/NLPsqNc/M7NeUa/3NrOP49mghFNtfQOXP76QL7ZVMrRPT/549hiSk6zlN4qISKPaOunsRCD6klTXAntFvU4BNIwSbv5nIXM+20zPtGQeOC+P7O6pQbckIhJquztLXEMm2cUzBav58zvFAPzhO6MZ3j8z2IZERDoBndYlcfXBqm1c/zfvFuk/+PpwThw5oIV3iIhIa7Q1sB273o1Ld+cSADaWeZPMauoaODa3H1d+fXjQLYmIdBptPa3LgEfNrNp/3Q24z8y2+691y6Uuqqauge/NKmBdaRUH9u3JHWeNJkmTzERE4qatgf2XmNczG6nRnbq6oF/N/pj3i7eSmZ7CA+ePI7ObJpmJiMRTW8/DvqijGpHwenLeSma89zlmcOfZozmwb0bQLYmIdDqadCa7peDzrdzwwhIArj52BF/P7R9wRyIinZMCW9ptfWkVl80soKa+gRMPGcD3jxkWdEsiIp2WAlvapbqunktnFrChrJoR/TO47TuHapKZiEgHUmBLmznnuPGFJSxcuY2sbik8cN44MtLbOn9RRETaQoEtbTZr7kqefH8VSQZ3nzuWIX16Bt2SiEinp8CWNpm3Ygs3vehNMvvJCTkcNaJvwB2JiHQNCmxptbUllXxvVgF1DY5TvjKQS486IOiWRES6DAW2tEpVbT3TZxSwqbyGnAGZ3HrmVzDTJDMRkT1FgS0tcs7xs78t5sPVJfTqkcqD54+jR5ommYmI7EkKbGnRX94t5tkFq0kyuPfcsey3V4+gWxIR6XIU2NKsOcs386t/FAJw/cm5HDGsT8AdiYh0TQpsadLqrRV8//EF1Dc4Jo/eh4uPHBp0SyIiXZYCWxpVWeNNMtuyvYaR+2bx2zM0yUxEJEgKbNmFc46fPvchS9aUsnfPNO4/bxzdUpODbktEpEtTYMsuHnp7BS8sWkNyknHvlLHs26t70C2JiHR5OjenC6moqGDp0qXN1ixcuZUbnl1A7bb1XHrK4aRt+5wFCz5vsj4nJ4cePTrvrPHW/M0ACgsLd3puSWf/u4lI/JlzLugeuiQzywJKSkpKyMrK2u3tFRUVUVZW1mxNYWEh+fn5u/1Z0WbOnElubm6zNZmZmQwfPjyun7unLFiwgLy8vLhvt6CggLFjx8Z9uyISLqWlpWRnZwNkO+dKm6tVYAcknoFdVFTEiBEj4tNYB1m2bFkoQ7u1I+zKykqKi4sZMmQI3bu3fAhBI2wRgbYFtnaJdwKRkXVrRrut0dbwaU5kVN/S6D9R9ejRo9Uj4SOOOKKDuxGRrkyB3Ynk5ubGbTerwkdEJLFolriIiEgIKLBFRERCQIEtIiISAgpsERGREFBgi4iIhEBoZomb2YvAaKAfsBX4N3Ctc26Nv34IsKKRt57knHspajtHAbcDhwBrgN875+6L+azvAT8BBgJLgCudc29HrU8HbgPOAboDrwHfc86tjsfv2lZWV8WYAUl037YM1iTWd7Du25YxZkASVlcVdCsiIqEWmsAGXgduBtYC++IF5jPApJi6Y/FCNmJL5AczGwr8E3gQyAeOAP5kZhudc8/6NWcBdwLfA94BpgP/MrODnXMr/U3dCXwTOBvYDPwBmG1mec65+rj9xq3UrXwlC6ZnwFvT4a09/enNywUWTM+gsHwlu/5PJSIirRWawHbO3RH18nMz+y3wvJmlOudqo9Ztds6ta2IzlwIrnXNX+q8LzWwc8GPgWX/Z1cDDzrmH/NdXmtkJwGXAdWaWDVwMnOec+zeAmeUDq/C+LLy8e79p21VlDGbs/eXMmjWL3JycPf3xzSpcupQpU6bw8MmDg25FRCTUQhPY0cxsL2AK8G5MWAO8aGbdgCLgDufcM1HrDgdeial/GbjYzFIBA/KA38bUvMKXw8M8IDV6O865NWa22K9pNLD93ejpUYsym/0l22B7TQML1zXwzmflVPZq2O3txfVKZ2vrWbiuAZfSbbf7EhHpykIV2Gb2O+ByoAfwHvCNqNXleKPjd4AG4FvAU2Z2gXNupl8zAFgfs9n1eH+HPniBndxEzYCobdQ457Y2U9OY64Abm/v92ityretp06Z1xObjIjMzbt9PRES6pEAD28xuouUQG++cm+//fCvwMLC//77HzOwbzrMJiN5tPt/MegPXADOjlsfe7cSillszNS3dJaWlmlvwJrtFZAJxmaQ2efJkIH43lIhc/zte1yYP8926REQSRdAj7HuAJ1uoKY784IfyJmCZmRXiHTc+DJjTxHvfA6ZGvV7HrqPgfkAd3uQxA+qbqImMutcBaWbWO2aU3Q94t6lfwjlXDVRHXptZU6Vt1qdPH6ZOndpyYRvF89rkIiKyewIN7KgAbo9I4qU3UzMGb1Z5xBy82d3RjgfmR46Fm1kBcBzwt6ia44AX/J8LgFp/2dP+ewYCI/FG8yIiInEX9Ai7VcxsAjAB+C/eOdgHAL8EluOPrs3sArwgXYh3DPubwA+Aa6M2dR9wuZndjndq1+F4M77Piaq5HZhhZvP9bV8CDPbfi3OuxMweBv5gZpvxThu7DfgI79xwERGRuAtFYAOVwOnAL4CeeKPml4Cz/V3NET/HO75dDywDvhs14Qzn3AozOxnvWPf38S6c8oPIOdh+zVNmtjdwA96FUxYDJzvnPo/6nKvwdqM/zZcXTrkwiHOwRUSkazDnWppLJR3BzLKAkpKSErKysoJuZycLFiwgLy+PgoICHcMWEelApaWlZGdnA2Q750qbq02s61iKiIhIoxTYIiIiIRCWY9gSBxUVFTsustKcwsLCnZ6bE69zv0VEpHk6hh2QII5hR45Nx5OOc4uItF9bjmFrhN2F5OTkUFBQ0GJdW64lnpNgNxsREemsNMIOSCLPEhcRkT1Ds8RFREQ6GQW2iIhICCiwRUREQkCBLSIiEgIKbBERkRBQYIuIiISAAltERCQEFNgiIiIhoMAWEREJAQW2iIhICCiwRUREQkCBLSIiEgIKbBERkRBQYIuIiISAAltERCQEFNgiIiIhoMAWEREJAQW2iIhICCiwRUREQkCBLSIiEgIKbBERkRBQYIuIiISAAltERCQEFNgiIiIhoMAWEREJAQW2iIhICCiwRUREQkCBLSIiEgIKbBERkRBQYIuIiISAAltERCQEFNgiIiIhoMAWEREJAQW2iIhICCiwRUREQkCBLSIiEgIKbBERkRBQYIuIiISAAltERCQEFNgiIiIhELrANrN0M1tkZs7MRsesG2Vmb5pZpZl9YWY3mJnF1JxhZh+bWbX/fFrMejOzm8xsjb+dN8zskJia3mY2w8xK/McMM+vVcb+1iIh0daELbOD3wJrYhWaWBbzqrxsPXAH8GLg6quZw4ClgBnCo//y0mU2M2tQ1/nsu97ezDnjVzDKjah4HRgMn+o/R/rZEREQ6hDnngu6h1czsJOB24AxgCTDGObfIX3cZcAvQ3zlX7S/7KV5wD3LOOTN7Cshyzp0Utc2XgK3OuXP80fga4E7n3O/89enAeuBa59z9ZpYLfAwc5pyb69ccBswBcpxzn7Tyd8kCSkpKSsjKytrNv4yIiIRRaWkp2dnZANnOudLmakMzwjaz/sCDwHlARSMlhwNvRsLa9zKwDzAkquaVmPe9DEzyfx4KDIiu8bf3ZlTN4UBJJKz9mveAkqiaxvpPN7OsyAPIbKpWREQkVigC2x/5Pgrc55yb30TZALyRcLT1UeuaqxkQU9dSzYZGPn9DVE1jrsML9chjdTO1IiIiOwk0sP3JXa6Fxzi83dpZeLu8mxO7f98aWd5YTeyylmoaO47Q2Hai3QJkRz0GNVMrIiKyk5SAP/8e4MkWaoqBnwOHAdUxk77nm9ks59wFeJPDYke4/fznyIi5qZro9fg1a5up6d9In33ZdWS+g79rfcfu+pjfQ0REpFmBBrZzbhOwqaU6M/sBXmhH7IN37PksIHIseQ5ws5mlOedq/GXH400iK46qOQ64I2pbxwPv+j+vwAvk44CF/menAUcB10ZtI9vMJjjn5vk1E/FGze8iIiLSAYIeYbeKc25l9GszK/d/XO6cixwLfhy4EXjUzG4GhgPXA790X06F/yPwlpldC7wAnAocCxzpf44zszuB682sCCjyt1Hhbx/nXKE/s/xBM5vub/cBYHZrZ4iLiIi0VSgCuzWccyVmdhxwLzAf2Ip3CtjtUTXvmtnZwK+BXwHLgbOiZ3zjnefdHfgT0BtvBH+8c64sqmYKcBdfziZ/Ee+8bRERkQ4RqvOwOxOdhy0iIp3yPGwREZGuTIEtIiISAgpsERGREFBgi4iIhIACW0REJAQU2CIiIiGgwBYREQkBBbaIiEgIKLBFRERCQIEtIiISAgpsERGREFBgi4iIhIACW0REJAQU2CIiIiGgwBYREQkBBbaIiEgIKLBFRERCQIEtIiISAilBNyCJpb6+nrfffpu1a9cycOBAvvrVr5KcnBx0WyIiXZ5G2LLDc889x7BhwzjmmGM499xzOeaYYxg2bBjPPfdc0K2JiHR5CmwBvLA+88wzGTVqFHPmzKGsrIw5c+YwatQozjzzTIW2iEjAzDkXdA9dkpllASUlJSVkZWUF2kt9fT3Dhg1j1KhRPP/88yQlffk9rqGhgcmTJ7N48WKKioq0e1xEJI5KS0vJzs4GyHbOlTZXqxG28Pbbb1NcXMz111+/U1gDJCUlcd1117FixQrefvvtgDoUEREFtrB27VoARo4c2ej6yPJInYiI7HkKbGHgwIEALF68uNH1keWROhER2fN0DDsgOoYtIiI6hi1tkpyczB/+8Admz57N5MmTd5olPnnyZGbPns1tt92msBYRCZAunCIAnH766TzzzDP86Ec/YtKkSTuWDx06lGeeeYbTTz89wO5ERES7xAOSSLvEo+lKZyIie05bdolrhC07SU5O5uijjw66DRERiaFj2CIiIiGgwBYREQkBBbaIiEgIKLBFRERCQIEtIiISAgpsERGRENBpXQErLW32tDsREenE2pIBunBKQMxsX2B10H2IiEhCGOSc+6K5AgV2QMzMgH2AsqB7aUQm3peJQSRmf4lIf7P20d+t7fQ3a59E/rtlAmtcC4GsXeIB8f+HafbbVFC87xIAlLV0qTzx6G/WPvq7tZ3+Zu2T4H+3VvWjSWciIiIhoMAWEREJAQW2NKYa+IX/LK2jv1n76O/WdvqbtU/o/26adCYiIhICGmGLiIiEgAJbREQkBBTYIiIiIaDAFhERCQEFtgBgZteZ2ftmVmZmG8zseTM7KOi+Ep2ZXWZmH5pZqf+YY2YnBd1XmPj/33NmdmfQvSQyM7vJ/ztFP9YF3VeiM7N9zWymmW02swozW2RmeUH31R4KbIk4CrgXOAw4Du8qeK+YWc9Au0p8q4GfAuP8x3+AF8zskEC7CgkzGw9cAnwYdC8hsQQYGPUYFWw7ic3MegPvALXAScDBwI+AbUH21V66NKkA4Jw7Mfq1mV0EbADygLcCaSoEnHN/j1n0MzO7DO+Lz5IAWgoNM8sAZgHTgJ8H3E5Y1DnnNKpuvWuBVc65i6KWFQfUy27TCFuaku0/bwm0ixAxs2QzOxvoCcwJup8QuBf4h3Pu30E3EiLDzWyNma0wsyfN7ICgG0pw3wLmm9lf/UN9C81sWtBNtZcCW3bh30nsduC/zrnFQfeT6MxslJmV411B6T7gNOfcxwG3ldD8LzZjgeuC7iVE5gLnAyfg7ZUYALxrZnsH2lViOwC4DCjC+7vdB9xlZucH2lU76Upnsgszuxc4BTjSOad7drfAzNKAwUAv4AxgKnCUQrtxZrYfMB843jn3gb/sDWCRc+7KIHsLE39+yXLg986524PuJxGZWQ0w3zk3KWrZXcB459zhwXXWPhphy07M7G683UjHKKxbxzlX45z71Dk33zl3HfAB8MOg+0pgeUA/oMDM6sysDm/S4w/818nBthcOzrntwEfA8KB7SWBrgdgvzoV4X7BDR5POBNixG/xu4DTgaOfcioBbCjMD0oNuIoG9xq6zm/8MLAV+55yr3/MthY+ZpQO5wNtB95LA3gFiT08dAXweQC+7TYEtEfcC5wKnAmVmNsBfXuKcqwyurcRmZjcD/wJWAZnA2cDRwInNvK1Lc86VATvNjTCz7cBmzZlompndBvwdWIm3h+LnQBbwlyD7SnB34B3nvx54GpiAdxrhJYF21U4KbIm4zH9+I2b5RcCje7STcOkPzMA7J7YE73ziE51zrwbalXRGg4AngD7ARuA94DDnXChHi3uCc+59MzsNuAW4AVgBXOmcmxVsZ+2jSWciIiIhoElnIiIiIaDAFhERCQEFtoiISAgosEVEREJAgS0iIhICCmwREZEQUGCLiIiEgAJbRDqcmV1oZtuC7kMkzBTYIiIiIaDAFhERCQEFtojsYGZvmNldZvZ7M9tiZuvM7KZWvreXmT1gZuvNrMrMFpvZN5qpv8zMlptZjZl9Ymbnxax3ZjbVzP5mZhVmVmRm34qpOdjM/mlm5f7nzjCzPlHrzzSzj8ys0sw2m9m//ftIi4SOAltEYl0AbAcmAtcAN5jZcc29wcyS8O5aNgnIBw4Gfgo0eqtM/4YMfwT+AIwE7gf+bGbHxJTeiHeXpa8A/wRmmdle/jYGAm8Ci4BxeHdI6+/XR9Y/ATyCdxvKo4Hn8G5/KhI6uvmHiOxgZm8Ayc65r0Ytmwf8xzn302bedzxeYOc655Y1sv5C4E7nXC//9TvAEufcJVE1TwM9nXOn+K8d8Gvn3P/6r3sCZcDJzrmXzOyXwETn3AlR2xiEd6vTg4AMoAAYojtaSWegEbaIxPow5vVavPsvN2c0sLqxsG5CLvBOzLJ3/OWN9uKc244X2JFe8oBj/N3h5WZWDiz11x0IfAC8BnxkZn81s2lm1ruV/YkkHAW2iMSqjXntaPnfisp2fE7s7j1rZFlzvSQBf8f7shD9GA685ZyrB44DTgI+Bq4APjGzoe3oVSRwCmwRiYcPgUFmNqKV9YXAkTHLJvnLW2sBcAhQ7Jz7NOaxHcB53nHO3QiMAWqA09rwGSIJIyXoBkQk/Jxzb5rZW8CzZnY18CmQ461yLzXylluBp81sAd5u628CpwPHtuFj7wWmAU+Y2a3AJmAYcLa/fBzwdeAVYAPeJLq+tO1LgUjC0AhbROLlDOB9vJnZHwO/B5IbK3TOPQ/8EPgJsASYDlzknHujtR/mnFsDHOF/xsvAYryZ5yVAA1AKfA1vdvky4NfAj5xz/2r7ryYSPM0SFxERCQGNsEVEREJAgS0iLTKzKdGnT8U8lgTdn0hXoF3iItIiM8vEu4pYY2p1YRKRjqfAFhERCQHtEhcREQkBBbaIiEgIKLBFRERCQIEtIiISAgpsERGREFBgi4iIhIACW0REJAQU2CIiIiHw/xLXh3OcVymDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 500x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from vireoSNP.plot.base_plot import heat_matrix\n",
    "\n",
    "fig = plt.figure(figsize=(5, 4), dpi=100)\n",
    "plt.plot(n_donor_list - 1, np.max(ELBO_list_all, axis=1))\n",
    "plt.boxplot(ELBO_list_all)\n",
    "plt.xticks(n_donor_list - 1, n_donor_list)\n",
    "plt.ylabel(\"ELBO\")\n",
    "plt.xlabel(\"n_clones\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}