LDA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "  and should_run_async(code)\n"
     ]
    }
   ],
   "source": [
    "import scanpy as sc\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.neighbors import BallTree\n",
    "\n",
    "import re\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pprint import pprint\n",
    "\n",
    "# Gensim\n",
    "import gensim\n",
    "import gensim.corpora as corpora\n",
    "from gensim.utils import simple_preprocess\n",
    "from gensim.models import CoherenceModel\n",
    "\n",
    "import pyLDAvis\n",
    "import pyLDAvis.gensim  # don't skip this\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.cluster import KMeans\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "home = '/Users/aj/Dropbox (Partners HealthCare)/Data/PCA/F8 image analysis/'\n",
    "adata = sc.read(str(home) + 'Z147_1_750.h5ad')\n",
    "x_coordinate='X_centroid';y_coordinate='Y_centroid';phenotype='phenotype_2';method='radius';radius=30;knn=10;imageid='imageid'\n",
    "data = pd.DataFrame({'x': adata.obs[x_coordinate], 'y': adata.obs[y_coordinate], 'phenotype': adata.obs[phenotype]})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "kdt = BallTree(data[['x','y']], leaf_size= 2) \n",
    "ind = kdt.query_radius(data[['x','y']], r=radius, return_distance=False)\n",
    "phenomap = dict(zip(list(range(len(ind))), data['phenotype'])) # Used for mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(ind)):\n",
    "    ind[i] = [phenomap[letter] for letter in ind[i]]\n",
    "texts = ind.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Create Dictionary\n",
    "id2word = corpora.Dictionary(texts)\n",
    "\n",
    "# Term Document Frequency\n",
    "corpus = [id2word.doc2bow(text) for text in texts]\n",
    "\n",
    "# Build LDA model\n",
    "lda_model = gensim.models.ldamodel.LdaModel(corpus=corpus,\n",
    "                                           id2word=id2word,\n",
    "                                           num_topics=15, \n",
    "                                           random_state=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0,\n",
      " '0.053*\"Mast cells\" + 0.053*\"CD163+ CD11C-  macrophages\" + '\n",
      " '0.053*\"Myofibroblast\" + 0.053*\"Unknown\" + 0.053*\"Patially Exhausted T cells\" '\n",
      " '+ 0.053*\"Melanocytes\" + 0.053*\"Blood Vessels\" + 0.053*\"CD163+ CD11C+ '\n",
      " 'macrophages\" + 0.053*\"CD11C+ CD163- PDL1+ cells\" + 0.053*\"Regulatory T '\n",
      " 'cells\"')\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Coherence Score:  0.5509164526578277\n"
     ]
    }
   ],
   "source": [
    "# Print the Keyword in the 10 topics\n",
    "pprint(lda_model.print_topics()[0])\n",
    "doc_lda = lda_model[corpus]\n",
    "\n",
    "# Compute Coherence Score\n",
    "coherence_model_lda = CoherenceModel(model=lda_model, texts=texts, dictionary=id2word, coherence='c_v')\n",
    "coherence_lda = coherence_model_lda.get_coherence()\n",
    "print('\\nCoherence Score: ', coherence_lda)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "Object of type complex is not JSON serializable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    339\u001b[0m                 \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    340\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 341\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mprinter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    342\u001b[0m             \u001b[0;31m# Finally look for special method names\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/pyLDAvis/_display.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(data, kwds)\u001b[0m\n\u001b[1;32m    309\u001b[0m     \u001b[0mformatter\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mip\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdisplay_formatter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformatters\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'text/html'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    310\u001b[0m     formatter.for_type(PreparedData,\n\u001b[0;32m--> 311\u001b[0;31m                        lambda data, kwds=kwargs: prepared_data_to_html(data, **kwds))\n\u001b[0m\u001b[1;32m    312\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    313\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/pyLDAvis/_display.py\u001b[0m in \u001b[0;36mprepared_data_to_html\u001b[0;34m(data, d3_url, ldavis_url, ldavis_css_url, template_type, visid, use_http)\u001b[0m\n\u001b[1;32m    176\u001b[0m                            \u001b[0md3_url\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md3_url\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    177\u001b[0m                            \u001b[0mldavis_url\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mldavis_url\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 178\u001b[0;31m                            \u001b[0mvis_json\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_json\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    179\u001b[0m                            ldavis_css_url=ldavis_css_url)\n\u001b[1;32m    180\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/pyLDAvis/_prepare.py\u001b[0m in \u001b[0;36mto_json\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    470\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    471\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mto_json\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 472\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mjson\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdumps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_dict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNumPyEncoder\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/json/__init__.py\u001b[0m in \u001b[0;36mdumps\u001b[0;34m(obj, skipkeys, ensure_ascii, check_circular, allow_nan, cls, indent, separators, default, sort_keys, **kw)\u001b[0m\n\u001b[1;32m    236\u001b[0m         \u001b[0mcheck_circular\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcheck_circular\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mallow_nan\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mallow_nan\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindent\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mindent\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    237\u001b[0m         \u001b[0mseparators\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mseparators\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdefault\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdefault\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msort_keys\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msort_keys\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 238\u001b[0;31m         **kw).encode(obj)\n\u001b[0m\u001b[1;32m    239\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    240\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/json/encoder.py\u001b[0m in \u001b[0;36mencode\u001b[0;34m(self, o)\u001b[0m\n\u001b[1;32m    197\u001b[0m         \u001b[0;31m# exceptions aren't as detailed.  The list call should be roughly\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    198\u001b[0m         \u001b[0;31m# equivalent to the PySequence_Fast that ''.join() would do.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 199\u001b[0;31m         \u001b[0mchunks\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0miterencode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_one_shot\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    200\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mchunks\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    201\u001b[0m             \u001b[0mchunks\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mchunks\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/json/encoder.py\u001b[0m in \u001b[0;36miterencode\u001b[0;34m(self, o, _one_shot)\u001b[0m\n\u001b[1;32m    255\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkey_separator\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitem_separator\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msort_keys\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    256\u001b[0m                 self.skipkeys, _one_shot)\n\u001b[0;32m--> 257\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_iterencode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    258\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    259\u001b[0m def _make_iterencode(markers, _default, _encoder, _indent, _floatstr,\n",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/pyLDAvis/utils.py\u001b[0m in \u001b[0;36mdefault\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    148\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    149\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 150\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mjson\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mJSONEncoder\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdefault\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/anaconda3/envs/scimap/lib/python3.7/json/encoder.py\u001b[0m in \u001b[0;36mdefault\u001b[0;34m(self, o)\u001b[0m\n\u001b[1;32m    177\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    178\u001b[0m         \"\"\"\n\u001b[0;32m--> 179\u001b[0;31m         raise TypeError(f'Object of type {o.__class__.__name__} '\n\u001b[0m\u001b[1;32m    180\u001b[0m                         f'is not JSON serializable')\n\u001b[1;32m    181\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: Object of type complex is not JSON serializable"
     ]
    },
    {
     "data": {
      "text/plain": [
       "PreparedData(topic_coordinates=                        x                   y  topics  cluster       Freq\n",
       "topic                                                                    \n",
       "7      0.092821+0.000000j -0.276135+0.000000j       1        1  46.303000\n",
       "1      0.079493+0.000000j -0.213885+0.000000j       2        1  10.809078\n",
       "11     0.166249+0.000000j -0.113758+0.000000j       3        1  10.035292\n",
       "9     -0.294451+0.000000j  0.111979+0.000000j       4        1   8.201760\n",
       "6     -0.064494+0.000000j  0.312759+0.000000j       5        1   7.728194\n",
       "2      0.253907+0.000000j  0.302327+0.000000j       6        1   5.139819\n",
       "13     0.253907+0.000000j  0.302326+0.000000j       7        1   5.103914\n",
       "5      0.196559+0.000000j  0.096789+0.000000j       8        1   1.908339\n",
       "14    -0.288400+0.000000j -0.044144+0.000000j       9        1   1.275266\n",
       "10    -0.391684+0.000000j  0.055260+0.000000j      10        1   0.732715\n",
       "12    -0.391069+0.000000j  0.052689+0.000000j      11        1   0.710952\n",
       "3      0.189587+0.000000j  0.082949+0.000000j      12        1   0.534665\n",
       "8      0.065858+0.000000j -0.223052+0.000000j      13        1   0.505669\n",
       "4      0.065858+0.000000j -0.223052+0.000000j      14        1   0.505669\n",
       "0      0.065858+0.000000j -0.223052+0.000000j      15        1   0.505669, topic_info=                   Term          Freq         Total Category  logprob  loglift\n",
       "10                Tumor  3.669468e+06  3.669468e+06  Default  19.0000  19.0000\n",
       "0               Unknown  1.729136e+06  1.729136e+06  Default  18.0000  18.0000\n",
       "2   CD11C+ CD163- cells  1.866604e+06  1.866604e+06  Default  17.0000  17.0000\n",
       "1         Blood Vessels  9.303250e+05  9.303250e+05  Default  16.0000  16.0000\n",
       "12          Melanocytes  1.062437e+06  1.062437e+06  Default  15.0000  15.0000\n",
       "..                  ...           ...           ...      ...      ...      ...\n",
       "1         Blood Vessels  3.502841e+03  9.303260e+05  Topic15  -2.9444  -0.2949\n",
       "12          Melanocytes  3.502841e+03  1.062437e+06  Topic15  -2.9444  -0.4277\n",
       "0               Unknown  3.502842e+03  1.729137e+06  Topic15  -2.9444  -0.9148\n",
       "2   CD11C+ CD163- cells  3.502841e+03  1.866605e+06  Topic15  -2.9444  -0.9913\n",
       "10                Tumor  3.502840e+03  3.669468e+06  Topic15  -2.9444  -1.6672\n",
       "\n",
       "[304 rows x 6 columns], token_table=      Topic          Freq           Term\n",
       "term                                    \n",
       "1         1  1.074892e-06  Blood Vessels\n",
       "1         2  9.887018e-01  Blood Vessels\n",
       "1        12  1.074892e-06  Blood Vessels\n",
       "1        13  3.765347e-03  Blood Vessels\n",
       "1        14  3.765347e-03  Blood Vessels\n",
       "...     ...           ...            ...\n",
       "0        11  2.461286e-02        Unknown\n",
       "0        12  5.783233e-07        Unknown\n",
       "0        13  2.025867e-03        Unknown\n",
       "0        14  2.025867e-03        Unknown\n",
       "0        15  2.025867e-03        Unknown\n",
       "\n",
       "[116 rows x 3 columns], R=19, lambda_step=0.01, plot_opts={'xlab': 'PC1', 'ylab': 'PC2'}, topic_order=[8, 2, 12, 10, 7, 3, 14, 6, 15, 11, 13, 4, 9, 5, 1])"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyLDAvis.enable_notebook()\n",
    "vis = pyLDAvis.gensim.prepare(lda_model, corpus, id2word)\n",
    "vis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "  and should_run_async(code)\n"
     ]
    }
   ],
   "source": [
    "topic_weights = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# n-1 rows each is a vector with i-1 posisitons, where n the number of documents\n",
    "# i the topic number and tmp[i] = probability of topic i\n",
    "topic_weights = []\n",
    "for row_list in lda_model[corpus]:\n",
    "    tmp = np.zeros(15)\n",
    "    for i, w in row_list:\n",
    "        tmp[i] = w\n",
    "    topic_weights.append(tmp)\n",
    "arr = pd.DataFrame(topic_weights).fillna(0).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEJCAYAAABGw1qNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABCIklEQVR4nO3de1xUdf7H8dcMM8DgoIAyysXMLPMalKObrQuZpSDgNS01qaxMs7I2SVQStcgyNtlS3O3XbTOtsE0sF9DU0FJTMS9peEnFQBSGi+Igl2GY3x/kJII6msOAfJ6PB4+c78yZ+Zwjzdtzvt/z/SosFosFIYQQwg6Uji5ACCHEjUtCRgghhN1IyAghhLAbCRkhhBB2IyEjhBDCblSOLqAxKS8vZ9++fXh7e+Pk5OTocoQQotEzm80YDAZ69OiBq6trneclZC6wb98+xo0b5+gyhBCiyVm2bBl6vb5Ou4TMBby9vYGag9WuXTsHVyOEEI3fqVOnGDdunPX782ISMhc4f4msXbt2+Pv7O7gaIYRoOi7VxSAd/0IIIexGQkYIIYTdSMgIIYSwGwkZIYQQdiMd/9dB+s5sPknNpKC4jDaeGiJDu3Jvr/aOLksIIRxOQuZPSt+ZzaIVe6gwmQEwFJexaMUeAAkaIUSzJ5fL/qRPUjOtAXNehcnMJ6mZDqpICCEaDwmZP6mguOyq2oUQojmxa8gYjUbCw8PJyckB4IcffmDIkCGEh4fz8ssvU1lZCUBubi7jxo0jJCSEyZMnU1paCkBJSQkTJ04kNDSUcePGYTAYAKisrCQqKorQ0FCGDx/OkSNHALBYLLz55puEhIQwePBgdu7cac/dA6CNp+aq2oUQojmxW8js2bOHMWPGkJWVZW2bNWsWCxcuZPXq1ZSXl7Nq1SoA5s6dy9ixY0lLS6NHjx4kJiYCkJCQgF6vJzU1lVGjRhEXFwfA0qVL0Wg0pKamMnPmTKKjowFYs2YNR44cISUlhcWLFxMdHU1VVZW9dhGAyNCuuKhr3+nqonYiMrSrXT9XCCGaAruFTFJSErGxseh0Omub2WzGaDRiNpupqKjAxcUFk8nEjh07GDRoEAAjRowgLS0NgPT0dCIiIgAIDw9n06ZNmEwm0tPTGTJkCAC9e/emuLiY3NxcNm7cyODBg1EqlXTs2BFfX1927dplr10Eajr3nx0VgPfvZy5KBUx5MEA6/YUQAjuOLjt/1nGhOXPmMH78eLRaLf7+/oSEhFBcXIxWq0WlqinF29ubvLw8APLz862TrqlUKrRaLUVFRbXaz29z6tQp8vPza4Xa+XZ7u7dXe+7t1Z70ndn8Y/lPtG3tZvfPFEKIpqDBOv4NBgPx8fGsXr2aH374gYCAAObPn4/FYqnzWoVCccn3USrrL1mpVNb7Xpd6vT306d4OtUrJD3tyG+wzhRCiMWuwb+CMjAw6d+7MTTfdhFKpZPTo0Wzfvh0vLy/rJTSoCaPzZyM6nY6CggIAqqqqMBqNeHh4oNPprIMALtymbdu29bY3FDdXNb266Ni85wTV1XUDTwghmpsGC5nOnTuzd+9ea2isX7+enj17olar0ev1pKSkAJCcnExQUBAAwcHBJCcnA5CSkoJer0etVhMcHGwdNJCRkYGLiwu+vr4EBQXxzTffYDabOX78OFlZWfTs2bOhdhGAvwX6UVRSQWZWUYN+rhBCNEYNdsd/p06dmDp1KpGRkTg5OdGhQwfmzZsHQGxsLNHR0SxZsgQfHx/efvttAKZOnUp0dDRhYWG4u7sTHx8PwPjx45k9ezZhYWE4OzuzYMECAEJCQti7d691UEBcXFy9y4HaU+9u7XBWKflh9wm639K6QT9bCCEaG4Wlvo6MZionJ4cBAwawfv36P7Vo2fz/bCfzWBEfzR6Ek/LS/UtCCNHUXel7U+74t4N+AX4Un63gl2OFji5FCCEcSkLGDnp3bYuz2onvd59wdClCCOFQEjJ24Oqione3tmzdexKzudrR5QghhMNIyNjJ3wL8OG2sYN9RuWQmhGi+JGTspFdXHS7OTnJjphCiWZOQsRNXZxV9urVjy95cuWQmhGi2JGTs6G+BvpSUVvLzkQJHlyKEEA4hIWNHd3Vpi8ZFLpkJIZovCRk7clE70aebD1v2nqRKLpkJIZohCRk76xfoy9lzlez9VS6ZCSGaHwkZO7vrdh0aFxU/yI2ZQohmSELGzpzVTvylRzu2/nwSU5VcMhNCNC8SMg3gbwF+GMtM7DlsuPKLhRDiBiIh0wDuvN0bN1cVP+yRS2ZCiOZFQqYBqFVO3N3Dhx/lkpkQopmRkGkg/QJ8KS2vYvehfEeXIoQQDUZCpoEEdtbRQqOWGzOFEM2KhEwDUauU9O3hw4/7TmKqMju6HCGEaBB2Dxmj0Uh4eDg5OTkA7Nq1i9GjRxMWFsbf//53KisrAcjMzGTkyJEMGjSIWbNmUVVVBUBubi7jxo0jJCSEyZMnU1paCkBJSQkTJ04kNDSUcePGYTDUjNyqrKwkKiqK0NBQhg8fzpEjR+y9izbrF+jLufIqdh2UUWZCiObBriGzZ88exowZQ1ZWFlATOM899xzz5s3jf//7HwBffvklAFFRUbzyyiusWbMGi8VCUlISAHPnzmXs2LGkpaXRo0cPEhMTAUhISECv15OamsqoUaOIi4sDYOnSpWg0GlJTU5k5cybR0dH23MWrEnCbN1qNWlbMFEI0G3YNmaSkJGJjY9HpdABs3ryZwMBAunTpAkBMTAwPPPAAJ06coLy8nMDAQABGjBhBWloaJpOJHTt2MGjQoFrtAOnp6URERAAQHh7Opk2bMJlMpKenM2TIEAB69+5NcXExubmNox9E5aSkb08ftu0/SYVJLpkJIW58dg2ZuLg49Hq99fHx48dxc3NjypQpRERE8O6779KyZUvy8/Px9va2vs7b25u8vDyKi4vRarWoVKpa7UCtbVQqFVqtlqKionrf69SpU/bczavSL9CPsgozPx2QUWZCiBtfg3b8m81mfvjhB6Kjo0lOTqasrIz33nsPi8VS57UKheKS7ZeiVNa/O5dqd4Q7bm2Du5uz3JgphGgWGvTbt02bNgQEBNC+fXucnJwIDQ1l7969tG3bloKCP2YpNhgM6HQ6vLy8MBqNmM3mWu0AOp3Ouk1VVRVGoxEPDw90Op11EMDF2zQGKicl99zhw/b9p+SSmRDihtegIdOvXz/279/PyZMnAfjuu+/o3r07fn5+uLi4sHPnTgCSk5MJCgpCrVaj1+tJSUmp1Q4QHBxMcnIyACkpKej1etRqNcHBwaxatQqAjIwMXFxc8PX1bcjdvKJ+Ab6UV5rZmZnn6FKEEMKuGjRkfHx8mDdvHpMmTSIkJIQzZ87w9NNPAxAfH8/8+fMJDQ2lrKyMyMhIAGJjY0lKSmLw4MFkZGTwwgsvADB16lR2795NWFgYy5cvZ/bs2QCMHz+eyspKwsLCiIuLY8GCBQ25izbp2akNrbTOcmOmEOKGp7DU1/HRTOXk5DBgwADWr1+Pv7+/XT8r8cs9bNiZzadzQ3B1Vtn1s4QQwl6u9L3ZeHrEm5l+gb5UVJrZmSmjzIQQNy4JGQfpfksbPLQucmOmEOKGJiHjIE5KBffc4cOOzDzKKqocXY4QQtiFhIwD9Qv0o9JkJuMXGWUmhLgxScg4ULeOrfF0d+F7uTFTCHGDkpBxICelgr/e4cvOzDzOlZscXY4QQlx3EjIO1i/Qj8qqanbIJTMhxA1IQsbBut7shVdLV5nLTAhxQ5KQcTClUkG/AF92HsiXS2ZCiBuOhEwj0C/AD1NVNdv3N54lCYQQ4nqQkGkEbu/gSZtWrny/W+YyE0LcWCRkGgGlUsFfA/z46WA+xjK5ZCaEuHFIyDQS/QJ9qTJXs33/SUeXIoQQ141M/9tI3H6TJ+5uahav2EPCZ7to46khMrQr9/Zq7+jShBDimknINBIbf8rhXHkV5uqalRcMxWUsWrEHQIJGCNFkyeWyRuKT1ExrwJxXYTLzSWqmgyoSQog/T0KmkSgoLruqdiGEaAokZBqJNp6aq2oXQoimQEKmkYgM7YqL2qlWm4vaicjQrg6qSAgh/jy7hozRaCQ8PJycnJxa7cuWLWP8+PHWx7m5uYwbN46QkBAmT55MaWkpACUlJUycOJHQ0FDGjRuHwWAAoLKykqioKEJDQxk+fDhHjhwBwGKx8OabbxISEsLgwYPZuXOnPXfvurq3V3ueHRWA9wVnLkF3+kmnvxCiSbNbyOzZs4cxY8aQlZVVq/3XX3/l3//+d622uXPnMnbsWNLS0ujRoweJiYkAJCQkoNfrSU1NZdSoUcTFxQGwdOlSNBoNqampzJw5k+joaADWrFnDkSNHSElJYfHixURHR1NV1XRWnby3V3s+jBnI1/FD6HyTB7sPGzBVmR1dlhBCXDO7hUxSUhKxsbHodDprW2VlJbNnz2bq1KnWNpPJxI4dOxg0aBAAI0aMIC0tDYD09HQiIiIACA8PZ9OmTZhMJtLT0xkyZAgAvXv3pri4mNzcXDZu3MjgwYNRKpV07NgRX19fdu3aZa9dtBuFQsG4QV0xFJexdttvji5HCCGumd3ukzl/1nGhf/zjH4wcORJ/f39rW3FxMVqtFpWqphRvb2/y8mrWVsnPz8fb27umUJUKrVZLUVFRrfbz25w6dYr8/PxaoXa+vSm683Zvut7sRdK6Q9zf56Y6/TVCCNEU2HQmU1BQwPr164Ga8IiMjOTAgQNX9UGbN2/m5MmTjBw5sla7xWKp81qFQnHJ91Eq6y9ZqVTW+16Xen1jp1AoeCS0C0Ul5aRtzXJ0OUIIcU1s+gaOjo4mOzubrVu3sm3bNoYNG8Zrr712VR+0evVqDh8+zNChQ4mJiWHfvn288MILeHl5YTQaMZtr+h4MBoP1bESn01FQUABAVVUVRqMRDw8PdDqddRDAhdu0bdu23vam6o5bvbnj1jZ8uf4w5RVNp29JCCHOsylkTp8+zWOPPcamTZsIDw9nxIgRlJVd3U2C8+fPJzU1lVWrVvHaa6/Ro0cPEhISUKvV6PV6UlJSAEhOTiYoKAiA4OBgkpOTAUhJSUGv16NWqwkODmbVqlUAZGRk4OLigq+vL0FBQXzzzTeYzWaOHz9OVlYWPXv2vKo6G5txIV04bazgf5uPOboUIYS4ajaFjMlkwmQy8f3333PPPfdQVlbGuXPnrlsRsbGxJCUlMXjwYDIyMnjhhRcAmDp1Krt37yYsLIzly5cze/ZsAMaPH09lZSVhYWHExcWxYMECAEJCQrjtttsYMmQIzzzzDHFxcbi6ul63Oh2hW8fW3HW7jv9+96usnCmEaHIUlvo6Mi7yzjvv8Mknn9C1a1eWLl1KeHg4oaGhTJkypSFqbDA5OTkMGDCA9evX1xqc4GiHfivmpX9u4pGQLjz0wO2OLkcIIayu9L1p0+iy559/ntGjR9OuXTsA4uPj6dKly/WtVFxS55s86dOtHSs3HiGs3y1oNWpHlySEEDax6XJZdXU1//vf/4iOjsZoNLJx40ZrR71oGONCulBaZiJ546+OLkUIIWxmU8gsWLCAgwcPsmdPzfom33//PfPnz7drYaK2W/xacc8dPny96SglpZWOLkcIIWxiU8hs3bqVN954AxcXF7RaLR9++CGbN2+2d23iImMHdaG8soqvvjvs6FKEEMImNoWMSqWqdVOjs7Oz9Q590XA6tGvJ3wL9WL35GMVnyx1djhBCXJFNIdO5c2eWLVuG2Wzm6NGjzJ49Wzr+HWTMwNsxmcz8d4P0zQghGj+bQmbWrFns37+fwsJCxowZQ2lpKTNnzrR3baIe/jp37u3VntQtxyg8I6tmCiEaN5uueWm1Wl5//XV71yJsNGbg7Wz8KYcV6w8zacQdji5HCCEu6bIhExcXx6xZs5g0aVK9z//rX/+yS1Hi8tq1bsH9fW5izY/HGdH/VnSebo4uSQgh6nXZkOnbty+Ada0X0XiMvr8z63dkk7TuEM+OCnR0OUIIUa/L9sncd999ABw/fpzhw4fX+tm/f3+DFCjqp/N0Y9DdHVi3/TdOFZY6uhwhhKjXZc9k3nnnHUpKSkhJScFoNFrbTSYTGzZsICYmxu4FiksbNeA2vt12nM/WHuTFMXc5uhwhhKjjsmcyAQEBeHh4oFQq8fDwsP60a9eOd999t6FqFJfQupWG0Hs6kr4zm5z8s44uRwgh6rjsmUxwcDDBwcEEBQVxxx0yiqkxevC+20j7MYvP1h4k6hG9o8sRQohabBpdlpiYWO/zMrrM8TzcXQj/a0e+Sv+V0QM608GnpaNLEkIIKxlddgMY0f82UrZksXztAWY82sfR5QghhNVlQ+b86LLhw4dz6tQpDh48SL9+/cjPz8fHx6dBChRX1rKFM0OCbuGLbw9x9MQZbvFr5eiShBACsHFamY0bN/Lwww8zd+5cCgsLGTx4MOvWrbN3beIqDAu+lRYaNcvXHHB0KUIIYWVTyLz77rskJSXRsmVLdDody5cv55133rHpA4xGI+Hh4eTk5ADwxRdfEB4eTkREBDNmzKCysmZtlMzMTEaOHMmgQYOYNWsWVVVVAOTm5jJu3DhCQkKYPHkypaU194SUlJQwceJEQkNDGTduHAaDAYDKykqioqIIDQ1l+PDhHDly5OqOSBOl1agZHtyJbftPcei3YkeXI4QQwFWsjKnT6ayPu3btikKhuOJ2e/bsYcyYMWRlZQFw7NgxPvjgAz7//HO+/vprqqurWb58OQBRUVG88sorrFmzBovFQlJSEgBz585l7NixpKWl0aNHD+sghISEBPR6PampqYwaNYq4uDgAli5dikajITU1lZkzZxIdHW370WjiIv52C+5uziyTsxkhRCNhU8hoNBpyc3OtwZKRkYGLi8sVt0tKSiI2NtYaUM7OzsyZMwetVotCoaBz587k5uZy4sQJysvLCQwMBGDEiBGkpaVhMpnYsWOHdeDB+XaA9PR0IiIiAAgPD2fTpk2YTCbS09MZMmQIAL1796a4uJjc3NyrOCRNl5urmpH9b+WnA/mMn5PGkJdWMeG1taTvzHZ0aUKIZsqmWZinTZvGhAkTMBgMPPTQQ2RlZdl0M+b5s4vz/Pz88PPzA6CoqIhly5Yxf/588vPz8fb2tr7O29ubvLw8iouL0Wq11gXSzrcDtbZRqVRotVqKiorqfa9Tp07h6+try642eS3dnAE4fbYCAENxGYtW1CybfW+v9g6rSwjRPNkUMnfeeSdJSUns2rWL6upqAgIC8PLyuuYPzcvL48knn2TkyJH85S9/4aeffqrzGoVCgcViqbf9Ui5cvdOW9hvRZ+sO1mmrMJn5JDVTQkYI0eBs/vY9cOAAlZWVmEwmMjIyWLt27TV94JEjRxgzZgzDhw9nypQpALRt25aCggLrawwGAzqdDi8vL4xGI2azuVY7gE6ns25TVVWF0WjEw8MDnU5nHQRw8TbNQUFx/QuZXapdCCHsyaYzmZdffpkff/yRDh06WNsUCgUDBw68qg8zGo088cQTvPjiiwwdOtTa7ufnh4uLCzt37qRXr14kJycTFBSEWq1Gr9eTkpJCRESEtR1qprxJTk5m0qRJpKSkoNfrUavVBAcHs2rVKvR6vbXvqLlcKgNo46nBUE+gtPHUOKAaIURzZ1PIZGRkkJqaSosWLf7Uh3355ZcUFBTw4Ycf8uGHHwI1N3xOnTqV+Ph4YmJiKC0tpVu3bkRGRgIQGxtLdHQ0S5YswcfHh7fffhuAqVOnEh0dTVhYGO7u7sTHxwMwfvx4Zs+eTVhYGM7OzixYsOBP1dzURIZ2ZdGKPVSYzNY2lZOCyNCuDqxKCNFcKSz1dXxc5JFHHuHTTz9tiHocKicnhwEDBrB+/Xr8/f0dXc41S9+ZzSepmRQUl6FSKTGbq/nnS/25WeY1E0JcZ1f63rTpTOauu+7ixRdfpH///ri6ulrbr/ZymWgY9/Zqb+3kLz5bzvP/SGfB0gzefiEIV2eb/sqFEOK6sOkbZ9euXQCsWLHC2nYtfTKi4Xm6u/LimLuIfW8rH369n2ceDHB0SUKIZsSmkFm6dKm96xB2dNftOkbceytfpf9KYGdv7rmj+QyEEEI4lk1DmA0GAxMnTmTQoEEUFhbyxBNP1BomLBq/R0K7cmt7D95J2k1+8TlHlyOEaCZsCpm5c+dy//334+LiQsuWLenSpQuzZs2yd23iOlKrlEQ90ovq6mreXv4TZnO1o0sSQjQDNoXMiRMnGD16NEqlErVaTVRUFCdPnrR3beI6822jZfLIAPYfLSRp3SFHlyOEaAZsChmFQkF19R//8jUajbUei6ajf6/29O/lz+ffHmT/0UJHlyOEuMHZFDIDBw5k2rRpnD17ls8//5xHH32U0NBQe9cm7GTSiDto27oF8ct2cvZcpaPLEULcwGwKmUmTJhEUFETPnj3ZsmULDz30kHXeMdH0uLmqiXqkF8Ul5bybtLveiUiFEOJ6sPnOvGHDhjFs2DA7liIa0m3tPYkc3I2PVu8n7cfjhPa92dElCSFuQDaFTPfu3ev0wbi6utK5c2def/11OnXqZJfihH0NC+7E7kP5vJ/8M906etGhnUw7I4S4vmy6XBYZGcmUKVP48ccf2b59O1FRUYSFhfHkk08yZ84cO5co7EWpVPDimLtwc1Xz1tKMWpNqCiHE9WBTyGzbto1nn32WVq1a4e7uzoQJEzhw4AAPPPAAJSUl9q5R2JFnS1deGHMnx0+d5cOv9zm6HCHEDcamkCkrKyM/P9/6OD8/n4qKmuV9zy8oJpquXl3aMiy4Eylbstj6s9z/JIS4fmzqk3nqqacYPnw4/fr1o7q6mq1btzJjxgwWLVrEXXfdZe8aRQOIHNyNfUcKeOeLXdzq74G3LHImhLgObAqZESNGEBAQwPfff49arWbKlCncfPPN5OTk0K5dO3vXKBpAzbQzeqa+nc4/lu8kbvJfcVIqHF2WEKKJs+lyGcDBgwc5cOAAu3btYvfu3QD4+/ujUsn6JDcKX28tk0fewf6jhaxYL9POCCH+PJtC5oMPPuDf//43Xbp0oXv37nz00UcsWbLE3rUJB+jfqz333uXPZ2sO8MsxmXZGCPHn2BQyycnJLFu2jMcee4zHH3+cTz/9lK+//tqmDzAajYSHh5OTkwPAli1biIiIYODAgSxcuND6uszMTEaOHMmgQYOYNWsWVVVVAOTm5jJu3DhCQkKYPHkypaWlAJSUlDBx4kRCQ0MZN26cdemByspKoqKiCA0NZfjw4Rw5csT2oyFQKBRMHnkHOi83XvtwG4/PW8OQl1Yx4bW1pO/MdnR5QogmxubLZVqt1vpnd3d3my6T7dmzhzFjxpCVlQVAeXk5M2fOJDExkZSUFPbt28fGjRsBiIqK4pVXXmHNmjVYLBaSkpKAmmUGxo4dS1paGj169CAxMRGAhIQE9Ho9qampjBo1iri4OKBmgTWNRkNqaiozZ84kOjra1l0Uv3NzVXOfvj1nz5koOFOOBTAUl7FoxR4JGiHEVbEpZPz8/PjPf/6DyWTCZDLx8ccf4+t75dUVk5KSiI2NRafTAbB37146dOhA+/btUalUREREkJaWxokTJygvLycwMBCoGWiQlpaGyWRix44dDBo0qFY7QHp6OhEREQCEh4ezadMmTCYT6enpDBkyBIDevXtTXFxMbm7u1R0Vwbfbf6vTVmEy80lqpgOqEUI0VTYvWrZu3ToCAwMJDAxk7dq1xMbGXnG7uLg49Hq99XF+fj7e3t7Wxzqdjry8vDrt3t7e5OXlUVxcjFartZ41nW+/+L1UKhVarZaioqJ63+vUqVO27Ka4QEFx2VW1CyFEfWwaGrZu3TqWLl1KWVkZ1dXVtGjR4po+rL7ZfhUKxVW3X4pSWX9mXqpdXFobTw2GegKljdw/I4S4CjZ9+3722WcAaDSaaw4YgLZt21JQUGB9nJ+fj06nq9NuMBjQ6XR4eXlhNBqtswqcb4eas6Dz21RVVWE0GvHw8ECn01kHAVy8jbBdZGhXXNROddr79vBxQDVCiKbKppDp2LEjMTExJCcns3btWuvP1QoICODYsWMcP34cs9nM6tWrCQoKws/PDxcXF3bu3AnUjGYLCgpCrVaj1+tJSUmp1Q4QHBxMcnIyACkpKej1etRqNcHBwaxatQqAjIwMXFxcbOo/ErXd26s9z44KwNtTgwJo4+FKOy83UrZkseeQ4YrbCyEE2Hi57PTp05w+fZrjx49b2xQKBQMHDryqD3NxceGNN97gueeeo6KiguDgYEJCQgCIj48nJiaG0tJSunXrRmRkJACxsbFER0ezZMkSfHx8ePvttwGYOnUq0dHRhIWF4e7uTnx8PADjx49n9uzZhIWF4ezszIIFC66qRvGHe3u1595e7a2Pz56rZGbiZl77aBuvPn0PXW72cmB1QoimQGGRZRGtcnJyGDBgAOvXr8ff39/R5TRKxSXlTF/0AyXnKpn/zF/p6NvK0SUJIRzoSt+bNl0uMxgMTJw4kUGDBlFYWMgTTzxRq99DNB+eLV15ddI9aJydmP3vreQajI4uSQjRiNk8hPn+++/HxcWFli1b0qVLF2bNmmXv2kQj1dbLjXlP34MFCzH/3kJ+8TlHlySEaKRsCpkTJ04wevRolEolarWaqKgoTp6UdUeas/Zt3Zn7VF9Ky0zM/vcWTp+tcHRJQohGyKaQUSgUVFdXWx8bjcZaj0Xz1Mnfg9lP3I3hdDmx723FWGZydElCiEbGppAZOHAg06ZN4+zZs3z++ec8+uij1lFhonnrfktrZj3Wh9/ySpj3/o+UV1Q5uiQhRCNiU8hMmjSJoKAgevbsyZYtW3jooYd47rnn7F2baCLu6qJj2iN6Dh4vIu7j7ZiqZEluIUQNm0Jm+fLlDBs2jISEBN555x1Gjx7Ne++9Z+/aRBPy1zt8eW70new+ZOCtT3diNsvlVCHEFW7G/OyzzygvL+fjjz+mouKPjl2TycTSpUuZOHGi3QsUTcf9fW7iXLmJ/1u1j3eSdjP1oTtRyhLOQjRrlw0ZlUrFoUOHKC8v59ChP5bjdXJy4pVXXrF7caLpGRLUidLyKpavOYCbq4qJw3pedlJTIcSN7bIhM2rUKEaNGsW6deu4//77G6om0cQ9/EBnzpWbSN54hBYaNY+EdHV0SUIIB7Fp7rK7776buXPncvToUf75z3/y9ttvM3369D81I7O4cSkUCiZEdKe0zMQX3x5i9Q/HOFdmoo2nhsjQrrXmQxNC3Nhs6viPi4ujZcuWFBYW4uLigtFoZPbs2fauTTRhCoWCHp3aoFQoKC0zyRLOQjRTNoVMZmYmL774IiqVCo1GQ3x8PJmZsgyvuLxP0zKpvmj+VVnCWYjmxaaQuXhlSbPZLKtNiiuSJZyFEDb1yfTu3Zu33nqL8vJyvv/+ez799FP69Olj79pEE3epJZzdNGoHVCOEcASbTkemTZuGm5sb7u7uJCQk0KVLF6Kjo+1dm2ji6lvC+XwfzRffHnRQVUKIhmRTyOzdu5cdO3ZQXl5OZWUlGzZsYOTIkfauTTRxFy/h7O2pYepDgfTv5c+naQf4JOUXZM08IW5sNl0ue+WVVxg9ejRdu3aVG+vEVbl4Cefzbc5qJ1asP0yFycyTQ3rI75UQNyibQsbZ2ZnHHnvsun3oqlWrrHOfBQUFMX36dDIzM4mJicFoNKLX65k7dy4qlYrc3FyioqIoLCykY8eOxMfH06JFC0pKSpg2bRrZ2dl4eXmRkJCAt7c3lZWVzJo1i3379uHq6kp8fDydOnW6brWLP0+pVDDlwQCc1U58vekolaZqJo+4Q6agEeIGZNPlsltuuYWff/75unxgWVkZcXFxLF26lFWrVpGRkcGWLVuIiorilVdeYc2aNVgsFpKSkoCaVTnHjh1LWloaPXr0IDExEYCEhAT0ej2pqamMGjWKuLg4AJYuXYpGoyE1NZWZM2dK31EjpVAoeGpoDx687zbStmbxzy92yaSaQtyALhsyERERREREsHv3bsaMGUNoaKi1LSIi4po+0Gw2U11dTVlZGVVVVVRVVaFSqSgvLycwMBCAESNGkJaWhslkYseOHQwaNKhWO0B6erq1hvDwcDZt2oTJZCI9PZ0hQ4YANaPiiouLyc3NvaZahX0pFAoiB3dlXEgXNmRkE79sJ1USNELcUC57ucwek2BqtVqmTp1KaGgorq6u9OnTB7Vajbe3t/U13t7e5OXlUVxcjFarRaVS1WoHyM/Pt26jUqnQarUUFRXVaj+/zalTp/D19b3u+yL+PIVCwcMP3I6zyomPVu/HVFXN9Eg9apXTlTcWQjR6lw0Ze9wLc+DAAf773//y3Xff4e7uzrRp09i8eXOd1ykUinpHHl2ug/hSN4jKjaON34j+t+KiVvKvlT/z6gfbmPl4H1ydbeoyFEI0Yg3+7fvDDz/Qt29fWrdujbOzMyNGjGDbtm0UFBRYX2MwGNDpdHh5eWE0GjGbzbXaAXQ6nXWbqqoqjEYjHh4e6HQ6DAZDnfcSjV9Yv1t4fnQguw8bmPf+Ns6VmxxdkhDiT2rwkOnSpQtbtmzh3LlzWCwWNmzYQJ8+fXBxcWHnzp0AJCcnExQUhFqtRq/Xk5KSUqsdIDg4mOTkZABSUlLQ6/Wo1WqCg4NZtWoVABkZGbi4uMilsibkgb904KWxvdh/rJDY97ZiLJOgEaIpU1gccDfce++9x1dffYVaraZnz57ExsZy7NgxYmJiKC0tpVu3bsyfPx9nZ2dOnDhBdHQ0hYWF+Pj48Pbbb9OqVStOnz5NdHQ02dnZuLu7Ex8fj7+/PxUVFcyePZt9+/bh7OzMa6+9Rvfu3W2qKycnhwEDBrB+/Xr8/f3tfBTE5WzZm8tbn2bQwacl8ybeQ8sWzo4uSQhRjyt9bzokZBorCZnGJSMzj9c/3o5vmxa8OukePN1dHV2SEOIiV/relJ5V0Wjpu7Yl9om7efWjbTz/j+9wUiopOlMui58J0YTIsCvRqAV09mZYcCdOn62k8Ey5LH4mRBMjISMavQ0ZdcNEFj8TommQkBGNnix+JkTTJSEjGr02npr6n1DA97tPyHIBQjRiEjKi0atv8TO1Som3h4YFSzOI+2g7BaflrEaIxkhCRjR69S1+9vzoQN6bcT8TIrqz65CBZxZsIHXLMaqr5axGiMZEhjCLJqG+xc8Aht97K3f38GHxl7tJ/O9eNu46wbOjAvDXuTugSiHExeRMRjR5Pm1a8OrT9zD1oUCOnyzh+X+k88W6g7JsgBCNgISMuCEoFAru79OBxJfvo0/3dnyaeoAXF27k0G/Fji5NiGZNQkbcUDxbuhId2ZtZj/ehpLSSqHc28cHX+yivqHJ0aUI0S9InI25Id/fwoWenNvznf7+QvPEIW38+yZQHAzhjrOCT1EwKistkehohGoCEjLhhtdCoeebBAILu9GPRit3Mfm8rSqXCOgLt/PQ0gASNEHYil8vEDa9Hpza881J/3FxVdYY4y/Q0QtiXhIxoFpzVTpSV198vI9PTCGE/EjKi2bjU9DQaVxXllTIwQAh7kJARzUZ909MoFQrOlVcxZcEGftx3UuZBE+I6k45/0Wyc79y/eHRZGw8N//pqL3EfbUfftS1PD+9Ju9YtHFytEDcGh4TMhg0bWLRoEefOnaNfv37ExMSwZcsW5s+fT0VFBaGhobz44osAZGZmEhMTg9FoRK/XM3fuXFQqFbm5uURFRVFYWEjHjh2Jj4+nRYsWlJSUMG3aNLKzs/Hy8iIhIQFvb29H7KZohC41PU3C3+9l9Q9HWb7mAM8s2MCo+25j5H234XzRmY8Q4uo0+OWy7OxsYmNjSUxM5JtvvuGXX35h48aNzJw5k8TERFJSUti3bx8bN24EICoqildeeYU1a9ZgsVhISkoCYO7cuYwdO5a0tDR69OhBYmIiAAkJCej1elJTUxk1ahRxcXENvYuiCVI5KRkWfCtLpg/g7h4+LF97kGff+o6MzDxHlyZEk9bgIfPtt98yePBg2rVrh1qtZuHChWg0Gjp06ED79u1RqVRERESQlpbGiRMnKC8vJzAwEIARI0aQlpaGyWRix44dDBo0qFY7QHp6OhEREQCEh4ezadMmTCZTQ++maKJat9Lw8ng9rz7dF6VSwdz3f+T1j7eTX3zO0aUJ0SQ1eMgcP34cs9nME088wZAhQ1i+fDn5+fm1LmnpdDry8vLqtHt7e5OXl0dxcTFarRaVSlWrHai1jUqlQqvVUlRU1IB7KG4EgZ11vDutP5GDu/LTwXwmv7mBFesPYaqSSTeFuBoN3idjNpvJyMhg6dKluLm58cwzz6DR1B1aqlAo6h3pc7n2S1EqZRCduHpqlZJRAzoTfJc/76/axycpmazfkc3kEXdQfLZcpqcRwgYNHjJt2rShb9++eHl5ATBgwADS0tJwcvqjgzU/Px+dTkfbtm0pKCiwthsMBnQ6HV5eXhiNRsxmM05OTtZ2qDkLKigooF27dlRVVWE0GvHw8GjQfRQ3Fp2nGzMf60NGZh7vrfyZmH9vQalQUG2R6WmEuJIG/yd+//79+eGHHygpKcFsNvP9998TEhLCsWPHrJfSVq9eTVBQEH5+fri4uLBz504AkpOTCQoKQq1Wo9frSUlJqdUOEBwcTHJyMgApKSno9XrUanVD76a4Aem7tmVR1O/T01hkehohbNHgZzIBAQE8+eSTjB07FpPJxF//+lfGjBnDLbfcwnPPPUdFRQXBwcGEhIQAEB8fT0xMDKWlpXTr1o3IyEgAYmNjiY6OZsmSJfj4+PD2228DMHXqVKKjowkLC8Pd3Z34+PiG3kVxA7vc9DSG4jL+t/kYAbe1wc9be9lLuEI0FwqL3OJslZOTw4ABA1i/fj3+/v6OLkc0UhNeW4uhnvnOLpzhuU0rV+64zZuA27wJuK0NrVvVP6WNEE3dlb435Y5/Ia5SZGhXFq3YQ4XJbG1zUTsx5cEAbr/Zkz2HC9hz2MCOX/LYkJENgL9Oaw2cnp3aoHVzBiB9Z7YMIBA3NAkZIa7SpaanOd/u20ZLaN+bqa62kHWyhN2HDOz51cC6Hb/xv83HUCqgk78Hnu4u7DpksA6LlgEE4kYkISPENbjU9DQXUioV3OLXilv8WjGi/62Yqqo59Fsxew4b2HPYwPZf6s4mUGEy85///SIhI24YEjJCNBC1Skn3W1rT/ZbWjB3UhSEvraK+DtGCM+U8/4/v6NaxNd06etGtY2vaeEifjmiaJGSEcJA2npp6BxC4uapo1cKF9b9fXgPQeWpqhU77tu4olX+MXpO+HdFYScgI4SCXGkAwecQd3NurPWZzNcdyS/glq5BfjhWx57CB9J9yAGihUdP1Zi+6dfSiwmQmOf2I9X2kb0c0JhIyQjjIlQYQODkpubW9B7e292DI3zphsVjIKzrHL8dqQueXY4WXnCW6wmTmk5RMCRnhcBIyQjiQLQMIzlMoFLRr3YJ2rVtwn/4mAM4YK3gkNq3e1xtOl/FiwkZu8a0ZfHCLbytu9m2JxkX+txcNR37bhGjCWmld8L5E347GRYXWVc3Wn0+ydttxABQK8G3Tgo7ng+f38PFs6QpI3464/iRkhGjiLtW388zImr4di8VC4Zlyjp44w9HcMxw9cYbD2af5YU+u9fWe7i601DqTk2fEXC0Tf4rrR0JGiCbuSn07CoWCNh4a2nho6NO9nXU7Y5mJY7lnOHbiDEdOnGHjrhxrwJxXYTKz6Ms9FJWUc1O7ltzs05LWrVxlXjZhMwkZIW4AV9O3c55Wo6Znp5ppbgC++30KnItVVJr5aPUv1sctNGo6tHOng09N6HRo15IOPi3Rav6Y7Vwuu4nzJGSEEMCl79vx9tTwz7/fy/GTJRw/WULWqbMcP1nCxp9ySL1gRuo2rVzp4NMShQJ2HzJQZZbLbkJCRgjxu0v17USGdsXdzZkendrQ4/ezHgCLxYLhdBm/nTpL1u8BdPxUCcdyS+q8d4XJzKIVeziWW4LOy422v/94e2pwdb7015CcETV9EjJCCODKfTsXUygU6Dzd0Hm6oe/a1tp+qelyKkxmvv7+KFXm6lrtHloX2nq5ofNyQ+ep+T2AWpB18gzL1xyUm0ybOAkZIYTVtfTtXOxyl93en/kAxWfLyS8qI6/4HPlF58grqvnvrzmn2fpzrvUyW30qTGb+vfJn3FzVeHtq0Hm60UJz5ZVv5YzIcSRkhBDX1eUuuymVClq30tC6lYauHb3qbGuutlBcUk5e0TmiF/9Q7/sby0y8+uE262M3VxU6z5pLb94eGuvZlbenBm9PDXsPF7D4yz1yRuQgEjJCiOvqai+7XchJ+cdw60vdZNq6lSszHu2N4XQZ+UVlGIrP1fy5+By/HCuitMx0xc+pMJl5L/ln3DRqPLQuNT/uLjirnS65jZwNXRuHhsybb75JcXExb7zxBpmZmcTExGA0GtHr9cydOxeVSkVubi5RUVEUFhbSsWNH4uPjadGiBSUlJUybNo3s7Gy8vLxISEjA29ubyspKZs2axb59+3B1dSU+Pp5OnTo5cjeFaHaux2W3S50RPRbWjds7eHF7h/q3O1duwlBcEzqG02Us+e/eel939pyJVz/YVqtN46LCw/2P0Gn1ewAZikvZuCvX2p/0Z86GmltYOSxktm7dysqVK7n33nsBiIqK4rXXXiMwMJCZM2eSlJTE2LFjmTt3LmPHjiUsLIzFixeTmJhIVFQUCQkJ6PV63nvvPZKTk4mLiyMhIYGlS5ei0WhITU1lx44dREdHs2LFCkftphDiGl3rGZGbq5oOPmo6+LQE4MsNh+s9I/Jq6cqsx/tw+mwFp40VnD5bwZnf/3vaWMEJg5FfjhVSUlqJpZ5uogqTmYTPd7F22294tXTFq5UrXi1daf37n1u3csWzpSsuF5wdpe/MrhWczeHSnUNC5vTp0yxcuJBJkyZx4MABTpw4QXl5OYGBgQCMGDGCd955h1GjRrFjxw4WL15sbX/kkUeIiooiPT2dZcuWARAeHs68efMwmUykp6czdepUAHr37k1xcTG5ubn4+vo6YleFEH+CPc+IHg/vRuebPK+4vdlczbCXv6n/uWoL5upqDv5WRNGZciqrquu8RqtRWwMoM6uoVh1w7TNmN5UzIoeEzOzZs3nxxRc5efIkAPn5+Xh7e1uf9/b2Ji8vj+LiYrRaLSqVqlb7xduoVCq0Wi1FRUX1vtepU6ckZIRopv5MHxHULLlwqf4hb08Nbz77N6DmvqHSMhOFJeUUnSmn8Ew5RSUX/Jwpp6LSXOc9oGbG7Mg5abRu5UrrVpqaM6GWNWdDXq001vYWrioUCsV1PSOyd1g1eMisWLECHx8f+vbty1dffQXU/OVcTKFQXLL9UpRK5VW1CyGahz97RnS5EXPnKRQKtG7OaN2c6dCuZb3vM+G1tZdcDVXftS2Fv4+s++VYEWfPVdZ5nYuzE14tXSk4XYbporOmCpOZD77eT/u27rTQqNFq1Ghc1TgpL/2d2RCX7xo8ZFJSUjAYDAwdOpQzZ85w7tw5FAoFBQUF1tcYDAZ0Oh1eXl4YjUbMZjNOTk7WdgCdTkdBQQHt2rWjqqoKo9GIh4cHOp0Og8FAhw4dar2XEEJcqz97NnTelVZDvVClyUxRye9nRGfKKSwpo/D3M6STBaX1vv9pYwUvLNxYq83NVUULjZoWrmpr+LT4/Wf9jt/qv3yXev0WvGvwkPnoo4+sf/7qq6/Yvn078+fPJzw8nJ07d9KrVy+Sk5MJCgpCrVaj1+tJSUkhIiLC2g4QHBxMcnIykyZNIiUlBb1ej1qtJjg4mFWrVqHX68nIyMDFxUUulQkh/rTr0T90NWHlrHayLlJ3sQPHi+o9I/LQuvDMgwGUlpk4V26itMyE8ff/1vxUYSgu49jJkt9fU1XnPQAK6nnva9Vo7pOJj48nJiaG0tJSunXrRmRkJACxsbFER0ezZMkSfHx8ePvttwGYOnUq0dHRhIWF4e7uTnx8PADjx49n9uzZhIWF4ezszIIFCxy2T0IIcTF7DmZ4Ykh3+vb0sfl9Jry6FsPpuoHSxlPzp+q7kMJSX8dHM5WTk8OAAQNYv349/v7+ji5HCCEu6Xp02F/cJwM1YfXsqACb3+tK35uN5kxGCCGE7Rr68t21kpARQohm7HqE1eXI2F4hhBB2IyEjhBDCbiRkhBBC2I2EjBBCCLuRjv8LmM01w/hOnTrl4EqEEKJpOP99ef7782ISMhcwGAwAjBs3zsGVCCFE03LhdF4XkpsxL1BeXs6+ffvw9vbGyenSK+QJIYSoYTabMRgM9OjRA1dX1zrPS8gIIYSwG+n4F0IIYTcSMkIIIexGQkYIIYTdSMgIIYSwGwkZIYQQdiMhI4QQwm4kZIQQQtiNhIwdLFq0iLCwMMLCwupd/nnRokX079+foUOHMnToUJYtW+aAKmuLjIwkLCzMWtOePXtqPb9lyxYiIiIYOHAgCxcudFCVf1ixYoW11qFDh9KrVy/mzZtX6zWN6TgbjUbCw8PJyckBbDueubm5jBs3jpCQECZPnkxpaanD6v3iiy8IDw8nIiKCGTNmUFlZWWeb5ORk+vXrZz3eDf17cnHNM2bMYODAgdZ6vv322zrbZGZmMnLkSAYNGsSsWbOoqqp/zfuGqHnjxo21fqfvvvtunn766TrbOPo4XzWLuK42b95seeihhywVFRWWyspKS2RkpGXt2rW1XvP0009bfvrpJwdVWFd1dbXlr3/9q8VkMtX7fFlZmSU4ONjy22+/WUwmk2XChAmW9PT0Bq7y0g4dOmR54IEHLIWFhbXaG8tx3r17tyU8PNzSvXt3S3Z2ts3Hc+LEiZbVq1dbLBaLZdGiRZYFCxY4pN6jR49aHnjgAcvZs2ct1dXVlpdfftny0Ucf1dlu3rx5lm+++aZBarzYxTVbLBZLeHi4JS8v77LbhYWFWXbt2mWxWCyWGTNmWJYtW2bvUq3qq/m8/Px8y4ABAyzHjh2rs50jj/O1kDOZ68zb25vo6GicnZ1Rq9V06tSJ3NzcWq/Zt28f//d//0dERATz5s2joqLCQdXWOHr0KAqFgqeeeoohQ4bw6aef1np+7969dOjQgfbt26NSqYiIiCAtLc1B1dY1Z84cXnzxRby8vGq1N5bjnJSURGxsLDqdDrDteJpMJnbs2MGgQYMAGDFiRIMd84vrdXZ2Zs6cOWi1WhQKBZ07d67zOw3w888/k5yczJAhQ5g2bRpnzpxpkHrrq/ncuXPk5ubyyiuvEBERwTvvvEN1dXWtbU6cOEF5eTmBgYFAwx7j+mq+0IIFC3j44Ye5+eab6zznyON8LSRkrrPbbrvN+kublZVFSkoKwcHB1udLS0vp2rUr06dPZ+XKlZSUlJCYmOigamuUlJTQt29fFi9ezMcff8znn3/O5s2brc/n5+fj7e1tfazT6cjLy3NEqXVs2bKF8vJyQkNDa7U3puMcFxeHXq+3PrbleBYXF6PValGpauaw9fb2brBjfnG9fn5+3HPPPQAUFRWxbNkyBgwYUGc7b29vnnvuOVatWoWPj0+dy5cNWXNhYSF33303r7/+OklJSWRkZPDll1/W2ubiv4eGPMb11XxeVlYW27dvJzIyst7tHHmcr4WEjJ0cPnyYCRMmMH369Fr/GmnRogX/93//R4cOHVCpVEyYMIGNGzc6rlDgzjvvZMGCBbi5ueHl5cWDDz5YqyZLPdPbKRSKhizxkj7//HMef/zxOu2N8TifZ8vxbIzHPC8vj0cffZSRI0fyl7/8pc7zixcvJiAgAIVCwZNPPsmmTZscUGWN9u3bs3jxYlq3bo1Go2H8+PF1/v4b4zGGmv6vsWPH4uzsXO/zjek420JCxg527tzJY489xksvvcTw4cNrPZebm1vrX1QWi8X6r1VHycjIYOvWrdbHF9fUtm1bCgoKrI/z8/PrPcVvaJWVlezYsYP77ruvznON8TifZ8vx9PLywmg0WtfoMBgMDj3mR44cYcyYMQwfPpwpU6bUef7s2bN8/PHH1seOPt4HDx5kzZo1l63n4r8HRx/j89avX8/gwYPrfa6xHWdbSMhcZydPnmTKlCnEx8cTFhZW53lXV1feeustsrOzsVgsLFu2jAceeMABlf7h7NmzLFiwgIqKCoxGIytXrqxVU0BAAMeOHeP48eOYzWZWr15NUFCQAyuucfDgQW6++Wbc3NzqPNcYj/N5thxPtVqNXq8nJSUFqBlR5KhjbjQaeeKJJ5g6dSoTJkyo9zVubm68//771lGJn376qUOPt8Vi4fXXX+fMmTOYTCa++OKLOvX4+fnh4uLCzp07Acce4/OKioooLy+nffv29T7f2I6zLSRkrrMPPviAiooK3njjDesQw88++4ynnnqKn3/+GS8vL+bNm8fkyZMJCQnBYrHUe7mnIfXv35/g4GCGDRvGyJEjGTlyJHfeeSdDhw4lLy8PFxcX3njjDZ577jkGDx7MLbfcQkhIiENrBsjOzqZdu3a12hrzcT7vcsdz1qxZrF+/HoDY2FiSkpIYPHgwGRkZvPDCCw6p98svv6SgoIAPP/zQ+jv9z3/+s1a9Tk5OJCQkMGfOHEJDQ9m/fz9RUVEOqRegS5cuTJw4kTFjxhAWFkbXrl0JDw8H/vgdAYiPj2f+/PmEhoZSVlZ2yX6QhpKTk1Pndxoa73G2hawnI4QQwm7kTEYIIYTdSMgIIYSwGwkZIYQQdiMhI4QQwm4kZIQQQtiNhIy44eXk5HD77bezYsWKWu0ffPAB0dHR1+1z7rvvPuvQWHszGo08/PDDhIWF1brp8HLGjx9/zXNznT171uHDe0XT1LhvFRXiOlEqlbz55pvo9Xo6duzo6HL+tMzMTAoLC+udvt4ezpw502ABKm4sEjKiWXB1deXxxx/npZde4vPPP68zL1R0dDS33XYbTzzxRJ3H9913H+Hh4aSnp3P69Gmee+45fvrpJ/bv349KpWLJkiW0bdsWgOXLl3PgwAEqKyt5/PHHefDBBwHYsGEDS5YswWQy4erqyvTp07nzzjt599132b17N/n5+dx+++3Ex8fXqmvdunUsWrQIs9mMVqtlxowZaLVaZs6cSV5eHkOHDuWLL77A1dXVuo3BYCA2NpajR4+iVCp5+OGHa52F5OTkEBERwa5du+o8NhgMTJ8+neLiYgCCg4N54YUXmDFjBuXl5QwdOpSvvvqKrKws4uLiOH36NGazmfHjx/Pggw+ybds24uLicHNz49y5cyxbtoxZs2Zx/PhxlEol3bt3Z968eSiVchGluZCQEc3G5MmT2bJlCwsXLmT69OlXtW1FRQVff/01KSkpvPTSS6xcuZIuXbowZcoUVq5cyaRJk4Cau/lXrlxJXl4ew4YNIyAgALVazcKFC/nkk0/w9PTk8OHDPP7446xduxaomXJ+9erVdeagOnLkCLGxsXz++ee0b9+erVu38swzz5CWlsZrr73Gq6++yqpVq+rUOnfuXG6++WYSExM5e/YsY8aMqTUT+OUkJSXh7+/Phx9+yLlz55g1axZnz55l/vz5REREsGrVKqqqqnj++edZsGAB3bt35+zZszz00EPceuutQM3ksOvWrcPPz4/k5GRKS0tZtWoVZrOZ2NhYsrOz6dChw1Udf9F0SciIZkOpVPLWW28xfPhw+vXrd1XbDhw4EKiZ3bdNmzZ06dIFgJtuuqnWeh4PP/wwUDP5Yr9+/di6dStOTk7k5+fz2GOPWV+nUCj47bffAAgMDKx3ksMff/yRu+++2zqPVd++ffHy8mLfvn2XnS14y5Yt1qlG3N3dWb16tc37+be//Y2JEydy8uRJ7rnnHl566SXc3d1r7WNWVha//fYbM2fOtLaVl5fzyy+/0KlTJ3x8fPDz8wOgV69eLFy4kPHjx3PPPffw6KOPSsA0MxIyolnx9fVlzpw5TJ8+nWHDhlnbFQpFranfTSZTre0uvLymVqsv+f4XXgY6P0Ou2Wymb9++JCQkWJ87efIkOp2Ob7/9tt4JPs9vX19bVVXVZWtQqVS1Qig7OxtPT0/r48vt6x133MH69evZunUrP/74I6NGjWLx4sW1Zic2m820bNmy1llUQUEB7u7u7N69u9b+tG/fnm+//ZZt27bx448/8vjjjxMTE9Mo5r4TDUMujIpmJzQ0lKCgIP7zn/9Y2zw9Pdm3bx9QMxNuRkbGNb33ypUrgZqlBrZs2ULfvn25++672bx5M0eOHAFg48aNDBky5IordZ7fLjs7G4CtW7dy8uRJAgICLrtd3759+e9//wvUjAp79NFHycrKsj7fsmVLTCYTv/76K0CtwQPx8fEkJiZy//33M2vWLG699VaysrKsYWmxWOjYsSMuLi7WkDl58iTh4eHW43eh5cuXM2PGDPr160dUVBT9+vXj8OHDl61f3FjkTEY0SzExMdYp3qFmeO+0adMYNGgQ/v7+9OnT55ret6KiguHDh2MymYiJibGOZJs3bx5///vfrWc3S5YsueQZzHm33norsbGxPPvss5jNZlxdXfnXv/6Fu7v7ZbebPXs2c+bMISIiAovFwtNPP02PHj2sz7u7uxMVFcVTTz2Fl5dXrbOKRx99lOjoaMLDw3F2dub2228nPDwcJycnunXrRmhoKJ999hmJiYnExcXx/vvvU1VVxdSpU+nVqxfbtm2rVcuwYcPYvn07gwcPRqPR4OvrK0OhmxmZhVkIIYTdyOUyIYQQdiMhI4QQwm4kZIQQQtiNhIwQQgi7kZARQghhNxIyQggh7EZCRgghhN1IyAghhLCb/wcYWQydXPiquwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "heterogeneities = []\n",
    "for i in range(2, 20):\n",
    "    km = KMeans(n_clusters = i, random_state = 0)\n",
    "    #km.fit(adata.uns['spatial_count_radius'].fillna(0))\n",
    "    km.fit(arr)\n",
    "    heterogeneities.append(km.inertia_)\n",
    "\n",
    "plt.plot(range(2, 20), heterogeneities, marker = 'o')\n",
    "plt.xlabel('Number of clusters')\n",
    "plt.ylabel('heterogeneities')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "kmeans = KMeans(n_clusters=15, random_state=100).fit(arr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "cluster_labels = list(map(str,kmeans.labels_))\n",
    "cluster_labels = pd.DataFrame(cluster_labels, index = data.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sub_data = data.merge(cluster_labels, left_index=True, right_index=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Vizualising using plotly\n",
    "import pandas as pd\n",
    "import plotly.express as px\n",
    "import plotly.io as pio\n",
    "pio.renderers.default = 'browser'\n",
    "\n",
    "\n",
    "def plotly (adata,phenotype,size=2, **kwargs):\n",
    "    data = pd.DataFrame({'x':adata['x'], 'y':adata['y'],'col': adata[phenotype]})\n",
    "    data = data.sort_values(by=['col'])\n",
    "    fig = px.scatter(data, x=\"x\", y=\"y\", color=\"col\", **kwargs)\n",
    "    fig.update_traces(marker=dict(size=size),selector=dict(mode='markers'))\n",
    "    fig.update_yaxes(autorange=\"reversed\")\n",
    "    fig.update_xaxes(showgrid=False,zeroline=False)\n",
    "    fig.update_yaxes(showgrid=False,zeroline=False)\n",
    "    fig.update_layout({'plot_bgcolor': 'rgba(0, 0, 0, 0)','paper_bgcolor': 'rgba(0, 0, 0, 0)'})\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "plotly (sub_data,phenotype=0,size=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "adata.obs['lda'] = sub_data[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import scimap as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/scimap/lib/python3.7/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sm.pl.stacked_barplot (adata,x_axis='lda',y_axis='phenotype_2', method='percent', plot_tool='plotly')"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}