08_Document_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Document clustering with LSI, LDA, and HDP\n",
    "\n",
    "Practical course material for the ASDM Class 09 (Text Mining) by Florian Leitner.\n",
    "\n",
    "© 2017 Florian Leitner. All rights reserved.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This notebook discusses a use-case of LSI, using [gensim](http://radimrehurek.com/gensim/index.html), closely following [its own tutorial for topic modeling](https://nbviewer.jupyter.org/github/RaRe-Technologies/gensim/blob/develop/docs/notebooks/gensim_news_classification.ipynb).\n",
    "\n",
    "## Setup\n",
    "\n",
    "The tutorial corpus gensim provides is the [Lee Background Corpus](http://www.socsci.uci.edu/~mdlee/lee_pincombe_welsh_document.PDF), but shortened to 300 documents, from the Australian Broadcasing Coporation's news mail serivce, and from around the years 2000-2001."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hundreds of people have been forced to vacate their homes in ...\n",
      "Indian security forces have shot dead eight suspected milita ...\n",
      "The national road toll for the Christmas-New Year holiday pe ...\n",
      "Argentina's political and economic crisis has deepened with  ...\n",
      "Six midwives have been suspended at Wollongong Hospital, sou ...\n"
     ]
    }
   ],
   "source": [
    "import gensim\n",
    "import os\n",
    "\n",
    "# show text as HTML\n",
    "from IPython.display import HTML\n",
    "\n",
    "test_data_dir = os.sep.join([\n",
    "    gensim.__path__[0], 'test', 'test_data'\n",
    "])\n",
    "lee_train_file = (\n",
    "    test_data_dir + os.sep + 'lee_background.cor'\n",
    ")\n",
    "\n",
    "with open(lee_train_file) as f:\n",
    "    for n, l in enumerate(f):\n",
    "        if n < 5:\n",
    "            print(l[:60], \"...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to load the documents and do a simple tokenization with gensim.\n",
    "The `gensim.util.tokenize` method (used by `simple_preprocess`) split text into contiguous tokens of alphabetic characters (e.g., removing digits!), and `simple_preprocess` by default lower-cases all tokens.\n",
    "In a sense, this is quite some harsh preprocessing, but is useful for tutorials and/or very noisy texts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from gensim.utils import simple_preprocess\n",
    "\n",
    "def read_and_tokenize(file_path):\n",
    "    with open(file_path) as f:\n",
    "        for doc in f:\n",
    "            # strip accents and tokens\n",
    "            # with less than 3 characters\n",
    "            # NB: implicitly removes tokens\n",
    "            # *longer* than 15 characters\n",
    "            yield simple_preprocess(\n",
    "                doc,\n",
    "                deacc=True,\n",
    "                min_len=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "300"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raw_documents = list(read_and_tokenize(\n",
    "    lee_train_file))\n",
    "\n",
    "len(raw_documents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be looking at the [LSI](https://radimrehurek.com/gensim/models/lsimodel.html#gensim.models.lsimodel.LsiModel) and [LDA models](https://radimrehurek.com/gensim/models/ldamodel.html#gensim.models.ldamodel.LdaModel); To speed up computations a bit, we will be using gensim's [multicore implementation for LDA](https://radimrehurek.com/gensim/models/ldamulticore.html#module-gensim.models.ldamulticore)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from gensim.models import LsiModel, LdaModel\n",
    "# disable this line if your CPU has only one core:\n",
    "from gensim.models.ldamulticore import LdaMulticore as LdaModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation\n",
    "\n",
    "First, we will do a little trick we will actually only learn about in the future, but that will *vastly* improve performace: detecting collocations (roughly, idioms), lemmatize the words (\"am\" -> \"be\", \"going\" -> \"go\", \"men\" -> \"man\", etc.), and remove **stopwords** (ubiquitious words found almost in any text that are not *very* siginificant for the text's semantics)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from gensim.models import Phrases\n",
    "from gensim.models.phrases import Phraser\n",
    "from nltk.corpus import stopwords\n",
    "from nltk.stem import WordNetLemmatizer\n",
    "\n",
    "collocations = Phraser(Phrases(raw_documents))\n",
    "lemmatizer = WordNetLemmatizer()\n",
    "lemmatize = lambda t: lemmatizer.lemmatize(t, pos='v')\n",
    "stopwords_en = (frozenset(stopwords.words('english'))\n",
    "                | frozenset([\"also\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "united_states team monica seles jan michael gambill score decisive victory unseeded france first hopman cup match burswood dome perth pair runners million dollar mix team event last_year single encounter give unbeatable lead year_old seles currently rank"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def preprocess(doc):\n",
    "    \"\"\"\n",
    "    Stopword filtering,\n",
    "    collocation detection (joining),\n",
    "    and token lemmatization.\n",
    "    \"\"\"\n",
    "    \n",
    "    doc_filtered = filter(\n",
    "        lambda w: w not in stopwords_en,\n",
    "        doc)\n",
    "    doc_colloc = collocations[doc_filtered]\n",
    "    doc_lemmas = map(lemmatize, doc_colloc)\n",
    "    # journalist love the word \"say\" in all its \n",
    "    # forms - sayed, saying, says, ...\n",
    "    # so we filter it *after* lemmatizing\n",
    "    return [l for l in doc_lemmas if l != \"say\"]\n",
    "\n",
    "HTML(\" \".join(preprocess(raw_documents[6][:50])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "texts = list(map(preprocess, raw_documents))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we convert the documents (which by now are token sequences/lists) into numeric vectors, using the Bag-of-Word approach (i.e., ignoring token order, only measuring token *frequency*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(4, 1), (18, 2), (19, 1), (58, 1), (64, 3), (92, 1), (94, 1), (114, 1), (122, 1), (127, 1)]\n"
     ]
    }
   ],
   "source": [
    "from gensim.corpora import Dictionary\n",
    "\n",
    "dictionary = Dictionary(texts)\n",
    "corpus = [dictionary.doc2bow(doc) for doc in texts]\n",
    "\n",
    "print(corpus[6][:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LSI Clustering\n",
    "\n",
    "Latent semantic indexing ranks topics and `gensim` outputs them in a ranked order. However, LSI requires a `num_topics` parameter (set to 200 by default) to determine the number of latent dimensions after the SVD step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Topic 0</h3>\n",
       "<p>0.174*\"arafat\" + 0.172*\"palestinian\" + 0.161*\"israeli\" + 0.158*\"attack\" + 0.152*\"kill\" + 0.144*\"two\" + 0.142*\"force\" + 0.134*\"people\" + 0.123*\"government\" + 0.120*\"australia\"</p>\n",
       "<h3>Topic 1</h3>\n",
       "<p>-0.324*\"arafat\" + -0.311*\"palestinian\" + -0.269*\"israeli\" + -0.195*\"israel\" + -0.169*\"sharon\" + 0.145*\"australia\" + -0.142*\"hamas\" + 0.141*\"australian\" + -0.137*\"kill\" + -0.137*\"west_bank\"</p>\n",
       "<h3>Topic 2</h3>\n",
       "<p>0.377*\"fire\" + -0.255*\"afghanistan\" + 0.189*\"sydney\" + -0.176*\"force\" + -0.148*\"bin_laden\" + -0.142*\"pakistan\" + 0.139*\"firefighters\" + -0.131*\"taliban\" + 0.125*\"new_south\" + 0.125*\"wales\"</p>\n",
       "<h3>Topic 3</h3>\n",
       "<p>0.328*\"fire\" + -0.228*\"australia\" + 0.157*\"force\" + -0.153*\"match\" + -0.146*\"test\" + -0.136*\"play\" + 0.131*\"afghanistan\" + 0.121*\"sydney\" + 0.118*\"firefighters\" + -0.115*\"south_africa\"</p>\n",
       "<h3>Topic 4</h3>\n",
       "<p>-0.244*\"company\" + -0.188*\"qantas\" + -0.159*\"australian\" + -0.147*\"government\" + 0.138*\"test\" + -0.132*\"workers\" + 0.126*\"match\" + 0.114*\"south_africa\" + -0.113*\"meet\" + -0.111*\"unions\"</p>\n",
       "<h3>Topic 5</h3>\n",
       "<p>-0.254*\"force\" + 0.188*\"kill\" + 0.162*\"tell\" + -0.149*\"government\" + -0.135*\"afghanistan\" + 0.128*\"taliban\" + 0.127*\"one\" + 0.123*\"company\" + -0.121*\"australia\" + 0.112*\"people\"</p>\n",
       "<h3>Topic 6</h3>\n",
       "<p>0.511*\"pakistan\" + 0.310*\"india\" + -0.214*\"australian\" + 0.193*\"president\" + 0.172*\"indian\" + 0.147*\"musharraf\" + -0.141*\"afghanistan\" + 0.113*\"want\" + 0.104*\"tension\" + 0.097*\"attack\"</p>\n",
       "<h3>Topic 7</h3>\n",
       "<p>-0.426*\"metres\" + -0.247*\"event\" + -0.246*\"race\" + -0.175*\"pakistan\" + -0.149*\"water\" + -0.140*\"swim\" + -0.135*\"world\" + -0.115*\"women\" + 0.100*\"people\" + -0.098*\"ray\"</p>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lsimodel = LsiModel(\n",
    "    corpus=corpus,\n",
    "    num_topics=8,\n",
    "    id2word=dictionary)\n",
    "\n",
    "results = []\n",
    "\n",
    "for t in lsimodel.show_topics():\n",
    "    results.append(\n",
    "        \"<h3>Topic {}</h3>\\n<p>{}</p>\".format(*t)\n",
    "    )\n",
    "\n",
    "HTML(\"\\n\".join(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating cluster coherence\n",
    "\n",
    "Gensim can measure [topic coherence](https://radimrehurek.com/gensim/models/coherencemodel.html) following a meta-model by [Röder et al. (2015)](http://svn.aksw.org/papers/2015/WSDM_Topic_Evaluation/public.pdf) that combines several different approaches of measuring [word] topic/cluster coherence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.41634396699988108"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from gensim.models import CoherenceModel\n",
    "\n",
    "def get_topic_words(lsi_model):\n",
    "    topic_list = lsi_model.show_topics(formatted=False)\n",
    "    return [[word for word, prob in wp_pairs]\n",
    "            for topic_id, wp_pairs in topic_list]\n",
    "\n",
    "# window_size defines the \"context\"\n",
    "# tokens for the probability model\n",
    "# (see section 3.2, \"Boolean sliding\n",
    "#  window\" in Röder et al.)\n",
    "CoherenceModel(\n",
    "    topics=get_topic_words(lsimodel),\n",
    "    texts=texts,\n",
    "    dictionary=dictionary,\n",
    "    window_size=10).get_coherence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this, we can now estimate a good value for `num_topics` - simply by clustering with a range of topic numbers that seem realistic and picking the best result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def evaluate_graph(\n",
    "    dictionary,\n",
    "    corpus,\n",
    "    texts,\n",
    "    min_num=3,\n",
    "    max_num=15):\n",
    "    \n",
    "    n_topics = list(range(min_num, max_num + 1))\n",
    "    scores = []\n",
    "    \n",
    "    for num_topics in n_topics:\n",
    "        lm = LsiModel(\n",
    "            corpus=corpus,\n",
    "            num_topics=num_topics,\n",
    "            id2word=dictionary)\n",
    "        \n",
    "        cm = CoherenceModel(\n",
    "            topics=get_topic_words(lm),\n",
    "            texts=texts,\n",
    "            dictionary=dictionary,\n",
    "            window_size=10)     \n",
    "        \n",
    "        scores.append(cm.get_coherence())\n",
    "        \n",
    "    plt.plot(n_topics, scores)\n",
    "    plt.xlabel(\"N. topics\")\n",
    "    plt.ylabel(\"Coherence score\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd81fX1+PHXyQ4QAiEhjAAJSRgBgkBElmxURFHrqKh1\nthQVtY5a/NmvddTWqq1dilpbV2vROirLgQqIgspeCSMJI0AIYUMCmef3x73YiCG5kNx87jjPxyMP\n+Hzu53M/5yrJyXseUVWMMcaYuoQ4HYAxxhjfZ8nCGGNMvSxZGGOMqZclC2OMMfWyZGGMMaZeliyM\nMcbUy5KFMcaYelmyMMYYUy+vJgsRuUBENopIrohMq+O6K0RERSSrxrlMEVkiIutFZK2IRHkzVmOM\nMacm3lrBLSKhwCZgHLADWApMUtXsk66LAeYAEcBUVV0mImHACuBHqrpaRNoAB1W16lTPi4+P1+Tk\nZK98FmOMCVTLly/fq6oJ9V0X5sUYBgK5qpoPICIzgEuA7JOuewx4ErivxrnzgDWquhpAVffV97Dk\n5GSWLVvWGHEbY0zQEJFtnlznzW6ojkBBjeMd7nPfEpF+QCdVnX3Svd0AFZGPRGSFiNzvxTiNMcbU\nw5stC6nl3Ld9XiISAjwD3FjLdWHAMOBsoBT4VESWq+qn33mAyGRgMkDnzp0bJ2pjjDHf482WxQ6g\nU43jJGBXjeMYoDewQES2AoOAme5B7h3AQlXdq6qlwFyg/8kPUNUXVTVLVbMSEurtcjPGGHOGvJks\nlgLpIpIiIhHA1cDMEy+q6iFVjVfVZFVNBr4CJqrqMuAjIFNEmrkHu0fw/bEOY4wxTcRryUJVK4Gp\nuH7w5wBvqep6EXlURCbWc+8B4A+4Es4qYIWqzvFWrMYYY+rmtamzTS0rK0ttNpQxxpwe93hwVn3X\n2QpuY4wx9Qr6ZLHz4DF+9+EGig4fdzoUY4zxWUGfLErKKpm+II9PcoqcDsUYY3xW0CeL9LYt6NKm\nGZ9kW7IwxphTCfpkISKM7ZnIl3n7KCmrdDocY4zxSUGfLADG9kykvLKaRZuLnQ7FGGN8kiUL4Ozk\n1sRGhzMve4/ToRhjjE+yZAGEhYYwukdbPttQRGVVtdPhGGOMz7Fk4Ta2ZyIHSitYsf2g06EYY4zP\nsWThNqJ7AhGhITaF1hhjamHJwq1FZBiDUtswL7uIQNkCxRhjGoslixrG9WzLlr0l5BWXOB2KMcb4\nFEsWNYzNSASwrihjjDmJJYsa2sdG07tjS1vNbYwxJ7FkcZKxPRNZvv0Ae4+WOR2KMcb4DEsWJxnb\nMxFV+GyDLdAzxpgTLFmcpFeHlnSIjbKuKGOMqcGSxUlEhLEZiSzavJfjFVVOh2OMMT7BkkUtxvZM\n5FhFFV/m7nU6FGOM8QmWLGoxqGsbWkSG2RRaY4xxs2RRi4iwEEZ0T+CTnD1UV9tqbmOMsWRxCuN6\nJlJ8pIzVO2xjQWOMsWRxCqO6tyU0RKwryhhjsGRxSrHNwhmYHMc8m0JrjDGWLOoyNiORTUVH2bbP\nNhY0xgQ3SxZ1GNfzxMaCtprbGBPcLFnUoXObZnRPjGFe9m6nQzHGGEdZsqjH2Iy2LN16gIOl5U6H\nYowxjrFkUY9xGe2oqlYWbCx2OhRjjHGMJYt6ZHaMJSEm0mZFGWOCmiWLeoSECGN7tmXhpmLKKm1j\nQWNMcLJk4YFxGYkcLavk6/z9TodijDGOsGThgSGp8USHh1pXlDEmaFmy8EBUeCjnpsfzSU4Rqrax\noDEm+Hg1WYjIBSKyUURyRWRaHdddISIqIlknne8sIkdF5D5vxumJsRmJFB46zvpdh50OxRhjmpzX\nkoWIhALPAuOBDGCSiGTUcl0McCfwdS1v8wzwgbdiPB1jerRFBNtY0BgTlLzZshgI5KpqvqqWAzOA\nS2q57jHgSeB4zZMicimQD6z3Yowea9MikgGdW9u4hTEmKHkzWXQECmoc73Cf+5aI9AM6qersk843\nB34BPFLXA0RksogsE5FlxcXeXzQ3NiOR9bsOs+vgMa8/yxhjfIk3k4XUcu7b0WERCcHVzXRvLdc9\nAjyjqkfreoCqvqiqWaqalZCQ0KBgPTEuw7Wx4KfWFWWMCTJhXnzvHUCnGsdJwK4axzFAb2CBiAC0\nA2aKyETgHOAKEXkSaAVUi8hxVf2rF+OtV2pCC7rGN+fj7CJ+NDjZyVCMMaZJeTNZLAXSRSQF2Alc\nDVxz4kVVPQTEnzgWkQXAfaq6DDi3xvmHgaNOJ4oTxmYk8vKXWzhyvIKYqHCnwzHGmCbhtW4oVa0E\npgIfATnAW6q6XkQedbce/NK4jEQqqpTPN+11OhRjjGky3mxZoKpzgbknnXvoFNeOPMX5hxs9sAbo\n37k1rZuFMy97NxMy2zsdjjHGNAlbwX2aQkOE0T0S+WzDHiqqqp0OxxhjmoQlizMwLiORw8crWbb1\ngNOhGGNMk7BkcQbOTY8nIizEFugZY4KGJYsz0DwyjKGpbZiXs9s2FjTGBAVLFmdoXEY7CvYfY/Oe\nOtcNGmNMQLBkcYbG9GwLYF1RxpigYMniDCW2jKJvUqwlC2NMULBk0QBjeyayquAge44cr/9iY4zx\nY5YsGmBcrxMbC+5xOBJjjPEuSxYN0D0xhqTW0XxiXVHGmABnyaIBRISxPRP5IncvpeWVTodjjDFe\nY8migc7LSKSssppFm21jQWNM4LJk0UBnp8QRExVmXVHGmIBmyaKBwkNDGNW9LZ9t2ENVta3mNsYE\nJksWjWBcRiL7SspZVWAbCxpjApMli0YwonsCYSHCx9YVZYwJUJYsGkHLqHAGdW1j4xbGmIDlUbIQ\nkWgR6e7tYPzZuIxE8opLyC+2jQWNMYGn3mQhIhcDq4AP3cdnichMbwfmb05sLPhJjrUujDGBx5OW\nxcPAQOAggKquApK9F5J/SmrdjJ7tW/JJtm39YYwJPJ4ki0pVPeT1SALAuIxElm3bz/6ScqdDMcaY\nRuVJslgnItcAoSKSLiJ/ARZ7OS6/NK5nItUKn22w1oUxJrB4kizuAHoBZcAbwCHgZ94Myl/17tiS\ndi2jbFaUMSbghNX1ooiEAo+o6s+BB5smJP8lIozNaMu7K3ZyvKKKqPBQp0MyxphGUWfLQlWrgAFN\nFEtAGNszkdLyKpbk7XM6FGOMaTR1tizcVrqnyv4HKDlxUlXf9VpUfmxwahuaR4QyL6eIUT3aOh2O\nMcY0Ck+SRRywDxhd45wClixqERkWyvBuCXyaU0T1Jb0JCRGnQzLGmAarN1mo6k1NEUggGZeRyAfr\ndrN25yH6dmrldDjGGNNgnqzgThKR90Rkj4gUicg7IpLUFMH5q1Hd2xIitprbGBM4PJk6+zIwE+gA\ndARmuc+ZU2jdPIKs5Djm2RRaY0yA8CRZJKjqy6pa6f56BUjwclx+77yMRDbsPkLB/lKnQzHGmAbz\nJFnsFZHrRCTU/XUdrgFvU4cxPRMB64oyxgQGT5LFzcBVwG6gELjCfc7UISW+OWltW1iyMMYEhHqT\nhapuV9WJqpqgqm1V9VJV3ebJm4vIBSKyUURyRWRaHdddISIqIlnu43EislxE1rr/HH2qe33ZuIxE\nvs7fz6FjFU6HYowxDeLJbKhXRaRVjePWIvIPD+4LBZ4FxgMZwCQRyajluhjgTuDrGqf3Aherah/g\nBuD1+p7ni8b2TKSyWlmw0TYWNMb4N0+6oTJV9eCJA1U9APTz4L6BQK6q5qtqOTADuKSW6x4DngSO\n13jGSlXd5T5cD0SJSKQHz/QpZ3VqRXyLCD7JsWRhjPFvniSLEBFpfeJAROLwbOV3R6CgxvEO97lv\niUg/oJOqzq7jfS4HVqpqmQfP9CmhIcKYHoks2LiH8spqp8Mxxpgz5kmy+D2wWEQeE5HHcNWyeNKD\n+2rb50K/fVEkBHgGuPeUbyDSC/gd8NNTvD5ZRJaJyLLi4mIPQmp6YzMSOXK8km+27Hc6FGOMOWOe\nDHC/huu3+yJgD/ADVfVkDGEH0KnGcRKwq8ZxDNAbWCAiW4FBwMwag9xJwHvA9aqad4rYXlTVLFXN\nSkjwzaUfw9LiiQoPsVlRxhi/5skAdyqQp6p/BdYCY2sOeNdhKZAuIikiEgFcjWslOACqekhV41U1\nWVWTga+Aiaq6zP3+c4AHVPXL0/9YviM6IpRhaQnMyy5CVeu/wRhjfJAn3VDvAFUikga8BKTgqphX\nJ1WtBKYCHwE5wFuqul5EHhWRifXcPhVIA/5PRFa5v/x2v+9xGW3ZefAYOYVHnA7FGGPOiCcD1dWq\nWikiPwD+pKp/EZGVnry5qs4F5p507qFTXDuyxt9/Dfzak2f4g9E9EhFZy7zsIjI6tHQ6HGOMOW2e\ntCwqRGQScD1wYtZSuPdCCjwJMZEMTI7jrWUFVFTZrChjjP/xJFncBAwGHlfVLSKSAvzTu2EFnsnD\nu7Lz4DFmrtpV/8XGGONjPJkNla2qd6rqv93HW1T1Ce+HFlhG92hL98QYpi/Mo7raBrqNMf7Fk5aF\naQQiwq0jU8ndc9Sm0Rpj/I4liyZ0UWZ7OsVF89yCPJtGa4zxKx4nCxFp7s1AgkFYaAiTh6eyquAg\nX+Xbim5jjP/wZFHeEBHJxrVWAhHpKyLPeT2yAHXlgCTiW0Tw3IJcp0MxxhiPedKyeAY4H3d1PFVd\nDQz3ZlCBLCo8lJuHpbBo817W7jjkdDjGGOMRj7qhVLXgpFNVXoglaFw3qAsxkWE8v7DWLa+MMcbn\neJIsCkRkCKAiEiEi9+HukjJnpmVUONcN7sLcdYXkFx91OhxjjKmXJ8liCnA7rloUO4Cz3MemAW4e\nmkJ4aAgvLMx3OhRjjKmXJ4vy9qrqtaqa6K7BfZ2q7muK4AJZQkwkV2Ul8e7KHew+dLz+G4wxxkFe\nq8Ft6vfT4alUK7y0yFoXTjlaVul0CMb4BW/W4Db16BTXjIsz2/PGN9s5WFrudDhB54WFeQx4bB7r\ndtqsNGPq480a3MYDU0amUlpexauLtzkdSlDJ3XOE33+8ibLKah56f53t12VMPbxZg9t4oEe7lozp\n0ZZXFm+htNy6RJpCVbVy/9traBYZygPje7Bi+0HeWbHD6bCM8Wme1uC+gtOvwW08dNuoVA6UVjDj\nm5OXsxhveG3JVlZsP8hDF2Xwk3O7MqBLa574YAOHjlU4HZoxPsvTvaE2AO8C7wNHRaSz90IKPgO6\nxDEwOY6/LcqnvNKKI3lTwf5SnvxwIyO7J3BZv46EhAiPXtKLA6XlPDNvk9PhGeOzPJkNdQeuVsU8\nXJXy5vC/inmmkdw6KpXCQ8f576qdTocSsFSVB95dS4jAby7rg4gA0KtDLNcN6sJrS7aSveuws0Ea\n46M8aVncBXRX1V6qmqmqfVQ109uBBZuR3RLo2b4lz1txJK/5z7IdfJG7l2kX9qRDq+jvvHbvuO60\nbhbBr2aus+3jjamFR9t9ADa30MtOFEfKLy7h4+zdTocTcIoOH+exOdkMTInj2oHf70WNbRbOLy7o\nwdKtB3hvpbXujDmZJ8kiH1ggIg+IyD0nvrwdWDC6sHc7urRpxnQrjtSoVJVf/ncd5ZXV/O7yTEJC\npNbrrhiQxFmdWvGbuRs4fNwGu42pyZNksR3XeEUEEFPjyzQyV3GkrqzecYjFebajSmOZs7aQedlF\n3DOuGynxp67hFRIiPHZJb/aVlPHHeZubMEJjfF+9i+tU9RFwVcpT1RLvhxTcLu+fxB8/2cxzC3IZ\nmhbvdDh+b39JOb96fz2ZSbHcMiyl3uv7JMVyzcDOvLpkK1ednUSPdi29H6QxfsCT2VCDrVJe04kK\nD+XHw1L4MncfqwsO1n+DqdNjs7M5dKyCJ6/IJCzUs5niPz+/Oy2jwnjo/fXWHWiMmyffPX/EKuU1\nqWvO6UzLqDCmL7DiSA0xf8Me3lu5k9tGpZ1WC6FVswjuv6AH32zZz8zVu7wYoTH+wyrl+aCYqHCu\nH5zMR9m7yd1jxZHOxJHjFfy/99bSLbEFU0elnfb9P8zqRN+kWB6fk8MRG+w2xirl+aobhyYTERrC\nC1Z69Yw88cEGig4f58kr+hIR5ulGBf/jWtndm+KjZfz5UxvsNsYq5fmo+BaRXH12J95buZNdB485\nHY5fWZK3j399vZ2bh6ZwVqdW9d9wCn07teLqszvx8pdb2VR0pBEjNMb/1JksRCQU+JFVynPGT4Z3\nRYG/WXEkjx0rr+KBd9fQOa4Z957XvcHv9/Pze9AiKoyH3reV3Sa41ZksVLUKuKSJYjEnSWrdjEv6\ndmDGNwXsL7HiSJ545pNNbN1XyhOX9yE6IrTB7xfXPIL7zuvOV/n7mbWmsBEiNMY/edIN9aWI/FVE\nzhWR/ie+vB6ZAVzFkY5VVPHK4q1Oh+LzVhcc5KVF+Uwa2JkhqY23RmXSwM707tiSx+dkWxlWE7Q8\nSRZDgF7Ao7gKIf0eeNqbQZn/6ZYYw7iMRF5dvJUS+0F1SuWV1dz/9hraxkTxwIU9GvW9Q92D3UWH\ny/iLDXabIOVJ8aNRtXyNborgjMutI1M5dKyCf3+z3elQfNZzC3LZWHSExy/rTcuo8EZ///6dW3NV\nVhJ//2ILuXtssNsEH09WcCeKyN9F5AP3cYaI3OLJm4vIBSKyUURyRWRaHdddISIqIlk1zj3gvm+j\niJzvyfMCVf/OrRnU1VUcqazSlricbOPuIzw7P5dLzurAmJ6JXnvOLy7oQbOIUH4101Z2m+DjSTfU\nK8BHQAf38SbgZ/Xd5J5J9SwwHsgAJolIRi3XxQB3Al/XOJcBXI2r++sC4Dn3+wWtW0emUXS4jP/a\n9tnf4aqnvZqYqHB+dXEvrz6rTYtI7ju/O1/m7mPuWttG3gQXT5JFvKq+BVQDqGolnq3gHgjkqmq+\nqpYDM6h9ZtVjwJPA8RrnLgFmqGqZqm4Bct3vF7SGp8fTq0NLnl+YT5UVR/rWP77Ywuodh3h4Yi/i\nmkd4/XnXntOFjPYt+fWcbBtDMkHFk2RRIiJtAAUQkUF4VgypI67CSSfscJ/7loj0Azqp6sllWuu9\nN9iICLeNTGPL3hI+XGe/1QJs3VvC7+dtZGzPRC7ObN8kzwwNER67tBeFh47z1/m5TfJMY3yBJ8ni\nHmAmkCoiXwKvAXd4cF9tFWa+/ZVYREKAZ4B7T/feGu8xWUSWiciy4uJiD0Lybxf0bkdKfHOmL8wN\n+j7z6mpl2rtrCA8N4fHLen9bT7spDOgSx+X9k3hpUT55xbZ3lwkOnsyGWgGMwDWF9qdAL1Vd48F7\n7wA61ThOAmpu4RkD9MZVhW8rMAiY6R7kru/eE7G9qKpZqpqVkJDgQUj+LTRE+OnwrqzbeZhFm/c6\nHY6j/r10O1/l7+fBC3uS2DKqyZ8/bXwPosJDedgGu02Q8HSHtYFAX6A/roHq6z24ZymQLiIpIhKB\na8B65okXVfWQqsararKqJgNfARNVdZn7uqtFJFJEUoB04BuPP1UAu6x/RxJbRgb19uWFh47x27kb\nGJLahh+e3an+G7wgISaSe8Z1Y9HmvXy03roFTeDzZOrs67gW4Q0DznZ/ZdV5E98OhE/FNZMqB3hL\nVdeLyKMiMrGee9cDbwHZwIfA7e6tR4JeZFgoPx7WlSX5+1i5/YDT4TQ5VeXB99ZRVa088YPMJu1+\nOtmPBnWhR7sYHpudw7Fy++dpApvU14QWkRwgQ328rZ2VlaXLli1zOowmcbSskqFPfMbAlDj+dn29\neTug/HflTn725ir+76IMj8qkets3W/Zz1QtLmDoqjfvOb/jGhcY0NRFZrqr1/iDxpBtqHdCu4SGZ\nxtIiMowbBndhXnYRm4No6+y9R8t4ZNZ6+nVuxY1Dkp0OB4CBKXFc1q8jL36ez5a9VqLeBK5TJgsR\nmSUiM4F4IFtEPhKRmSe+mi5EU5sbh6YQHR7K9CAqjvTwzPWUlFXx5OWZhIY41/10sgfG9yAiLIRH\nZtlgtwlcYXW8ZpsF+rC45hFcPbATry/Zxj3jupHUupnTIXnVx+t3M3tNIfeO60Z6YozT4XxH25ZR\n/GxsOr+ek8O87CLO62UNcRN4TtmyUNWFJ76ADbimusYAOe5zxmE/PrcrAC8t2uJwJN516FgFv/zv\nOnq0i2HKyFSnw6nVDUOS6ZbYgkdnZ3O8wga7TeDxZDbUVbimrV4JXAV8LSJXeDswU7+OraK5tF9H\nZizdzr6jZU6H4zW/mZPDvpJynrqiL+Ghp19PuymEh4bw6CW92XHgGM8F8bRmE7g8+c57EDhbVW9Q\n1etxrbn4P++GZTw1ZURXyiqrA7Y40heb9/LmsgJ+cm5X+iTFOh1OnQZ1bcPEvh14fmEe2/bZYLcJ\nLJ4kixBV3VPjeJ+H95kmkNY2hvPcxZGOHK9wOpxGVVpeybR319A1vjk/G5vudDgeeXBCT8JDhEdn\nZTsdijGNypMf+h+6Z0LdKCI3AnOAD7wbljkdt45M4/DxSt74OrCKIz310UZ2HDjGE5dnEhXuHzvU\nJ7aM4q6x6Xy6YQ+f5hQ5HY4xjcaTvaF+DrwAZOLa8uNFVb3f24EZz53VqRVDUtvw0hdbAmZwdfm2\nA7yyeCvXD+7CwJQ4p8M5LTcNTSGtbQsemWWD3SZw1LXOIk1EhgKo6ruqeo+q3g3sExHfnJISxG4b\nmUbxkTLeXeH/xZFUlV/NXEeH2Gjuv6Bx62k3hfDQEB6d2Ivt+0t5YWG+0+EY0yjqaln8EahteXCp\n+zXjQ4amtSEzKZYXPs+jsqra6XAa5PPNe1m38zB3jUmnRWRdS4F815C0eCZktue5BbkU7C91Ohxj\nGqyuZJFc21bk7l1hk70WkTkjruJIqWzbV8qctYVOh9Mg0xfk0q5lFJf28+96V7+c0JPQEOHR2TbY\nbfxfXcmiriIB0Y0diGm48zLakd62Bc/Oz6XaT0uvrth+gK/y9/Pjc1OICPPvSXftY6O5Y3Q687KL\nrLqh8Xt1fTcuFZGfnHxSRG4BlnsvJHOmQkKE20elsanoKPP8dCbO8wvyiI0OZ9LAzk6H0ihuGZZC\n744t+dmbK1m+Lfi2lPeUqrKq4CCPzc5m5FPzmfzaMqtC6GNOuUW5iCQC7wHl/C85ZAERwGWq6lO/\nKgXTFuV1qayqZvTvFxIbHc7MqUMdrfdwunL3HGHsHz7nzjHp3DOum9PhNJriI2Vc+fxiDpRW8NZP\nB9O9nW/tbeUUVSWn8Aiz1uxi9ppdFOw/RnioMKhrG1ZsO0BZZTXXDerCnWPSiWse4XS4AcvTLco9\nqWcxClf5U4D1qvpZI8TX6CxZ/M+Mb7Yz7d21vHrzQEZ0859ys/f9ZzWz1+xi8bQxAffDoWB/KZdP\nX4wIvD1lCJ3iAnvjx7rk7jnK7DW7mLV6F3nFJYSGCEPT4rk4sz3n9WpHbHQ4xUfKeOaTTcz4ZjvN\nI8O4Y3QaNwxJJjLMP9bb+JNGSxb+wpLF/5RXVjPiqfkktY7mP1OGOB2OR3YdPMbwJ+dz3aAuPDyx\nl9PheMXG3Ue48vnFxDWP4O1bhxDfItLpkJpMwf5SZq3ZxazVheQUHkYEzkmJ4+K+HbigVzvanOK/\nxaaiI/xmbg4LNhaT1DqaaeN7MKFPe79qMfs6SxZB7uUvt/DIrGzenDyIc7q2cTqcej0yaz2vL9nG\nwvtH0bFV4M6fWL5tP9e+9DWpCS2YMXkQMVHhTofkNYWHjjFnTSGz1hSyuuAgAP07t+KizA5MyGxP\nYsu65tB816LNxTw+J4cNu4/Qv3MrHpyQwYAurb0VelCxZBHkjpVXce6Tn9GzfUtev+Ucp8Op0/6S\ncoY+8RkX9mnP76/q63Q4Xjd/4x5+8uoyBnRpzas3D/SbrUw8sfdoGR+sLWTW6kK+2bofgN4dW7oS\nRJ/2Dep+q6pW3l5ewNMfb6L4SBkTMtvzi/N70LlN8HbpNQZPk4V/rngy9YqOCOWWYV353YcbWF1w\nkL6dWjkd0im9ungrxyqqmDKiq9OhNIlR3dvy9JV9+dmbq7jz3yt57tr+hPno1uueOFhazkfrdzNr\ndSGL8/ZSrZDetgX3jOvGRZnt6ZrQolGeExoi/PDszlyU2YEXPs/nxc/zmLe+iBuGdGHq6HRiowO3\nleYLrGURwI4cr2DY7+YzMCWOv11f7y8Ojigpq2To7z7j7GTfjdFbTnQVXpWVxO8uz/Srfvgjxyv4\nJKeIWasLWbS5mIoqpUubZlyc2YGL+3Zokhlfuw8d5+mPN/LOih20ig7nrjHpXDuoi8/WPPFV1rIw\nxESFc+OQZP706WY27D5Mj3YtnQ7pe2YsLeBgaQW3+mgFPG+6aWgKB0rK+fNnubRuHsED43s6HVKd\nqqqVj9bvZuaqXczfuIeyymo6xEZx09AULs7sQO+OLZs04bWLjeLpK/ty09BkHp+Tw8OzsnltyTam\nje/BuIxEv0q+/sCSRYC7aWgyLy3K59n5efxlUj+nw/mO8spqXlqUzzkpcfTvHJyDlXeP68b+0nJe\nWJhPXLMIfjrCN5PmzoPHuPvNVXyzZT/xLSKZNLAzF2W2p3/n1oSEOPtDuVeHWP7143P4bMMefjM3\nh8mvL+eclDh+OSHD5wtm+RNLFgGuVbMIrhvchb99ns/dY9Mbrf+4Mby/aieFh47z2x/0cToUx4gI\nj0zszcHSCn77wQZaN4/gqqxOTof1HbPX7OL/vbuWqmrlySsyubx/EqEOJ4iTiQhjeiYyvFsCM77Z\nzjOfbObiv37BD/p15L7zu9MhgGfYNRXr3AsCPx7WlfDQEKb7UG3o6mrl+YV5ZLRv6VcLB70hNET4\nw1VncW56PNPeWcNH631jc4SjZZXc95/VTH1jJV0TWjD3rnO5KquTzyWKmsJDQ/jR4GQW/HwkU0ak\nMnttIaOeXsDTH23kaFml0+H5NUsWQSAhxtVt8N7Knew44BvbZX+cXURecQm3jky1vmUgIiyE568b\nQJ+kVtzx75UsydvnaDyrCg4y4c+LeGfFDqaOSuM/UwbTpU1zR2M6HS2jwpk2vgef3jOC83u146/z\ncxn51ALe+Hq732/h7xRLFkFi8vCuiOATxXhUlekL8+jSphnje7dzOhyf0TwyjFduPJvOcc34yWvL\nWLfzUJOvH3LYAAAUh0lEQVTHUFWtPDs/lyumL6aySpnxk0Hcd353v51h1CmuGX+e1I/3bhtCcptm\n/L/31nLhnxexYOMep0PzO/75L8Cctg6torm8fxJvLitgz+HjjsayJH8fqwsOMnl4V79eX+ANrZtH\n8PotA4mNDufGl79hy96SJnv2roPHmPS3r3jqo41c0Lsdc+861y9W/3uiX+fW/GfKYKZf25+yympu\nfHkpt/1rOYGydKAp2HdqEJkyIpXKqmpe+mKLo3FMX5BHfItILu+f5Ggcvqp9bDSv3TKQaoUf/f1r\nipoguc9ZU8gFf/yc9TsP8fSVffnLpH4Bt8hNRBjfpz0f3z2cqaPSmLt2N++v2uV0WH7DkkUQSY5v\nzsV9O/DPr7ZxoKTckRjW7TzEos17uWVYSkBtc9HYUhNa8MpNZ3OgpJzr//4NB0u98/+rpKySn/9n\nNbe/sYKUhBbMufNcrhiQFNDjSJFhodwzrht9OsbyxAcbKC23gW9PWLIIMrePSqO0vIqXv3SmdTF9\nYR4xkWFcOygwiht5U2ZSK168Poste0u4+ZWljf5DbbV7EPtt9yD221MGkxzvP4PYDRESIjw8MYPd\nh4/zvA/NEvRlliyCTLfEGM7vlcjLi7dy+HhFkz57y94SPlhbyHWDu9AygHdbbUxD0+L509Vnsarg\nILf9awUVjTCT58Qg9uXTF1NeWe33g9hnakCXOCb2de0zVbDfN2YJ+rLg+tdhAJg6Kp0jxyt5fcm2\nJn3ui5/nERYaws1DU5r0uf5ufJ/2PH5ZHxZsLOa+/6xuUH31XQePcY17EPv83u344K7hATOIfSam\nje+BCDzxwQanQ/F5liyCUJ+kWEZ0S+DvX2xpsv7aosPHeWf5Tq7KSiIhJniK/jSWSQM78/Pzu/P+\nql08Ojv7jGbxnBjEXrvzEE9dkclfJ/Ujtllwt/A6tIrm1hFpzFlbyNf5zq5t8XVeTRYicoGIbBSR\nXBGZVsvrU0RkrYisEpEvRCTDfT5cRF51v5YjIg94M85gdMfoNPaXlPPvbwqa5Hn/+GILldXVTD7X\nN/c+8ge3jUzllmEpvLJ4K3/5LNfj+04exJ5757lcmdUpoAexT8fk4V3pEBvFI7OyqWpAqy3QeS1Z\niEgo8CwwHsgAJp1IBjW8oap9VPUs4EngD+7zVwKRqtoHGAD8VESSvRVrMMpKjuOclDhe/DyPssoq\nrz7rUGkF//xqGxdldrBCNQ0gIjx4YU9+0L8jf5i3ide/qr8bMZgHsT0VHRHKAxf2JLvwMG8ta5pf\nnvyRN1sWA4FcVc1X1XJgBnBJzQtU9XCNw+bAibSuQHMRCQOigXKg5rWmEUwdnUbR4TLeXr7Dq8/5\n59fbKCmvYoqP7qjqT0JChN9dnsmYHm156P11zF5T+zqBqmrluQU2iO2pizLbc3Zya57+aGOTT/zw\nF978l9MRqJmmd7jPfYeI3C4iebhaFne6T78NlACFwHbgaVXdX8u9k0VkmYgsKy4ubuz4A96wtHj6\ndmrF9AV5jTLLpjbHK6r4xxdbGNk9gYwOvldPwx+Fh4bw7LX9ObtLHHe/uYrPN3333/6ug8e49qWv\nePJDG8T2lIjwq4t7sb+0nL98utnpcHySN5NFbR2i3+sQVNVnVTUV+AXwS/fpgUAV0AFIAe4Vke/V\n3FTVF1U1S1WzEhKCe+fSMyEiTB2Vxo4Dx5jppZWsby0rYF9JObeNTPPK+werqPBQ/nZDFmltY/jp\n68tZuf0AAHPXFjL+T4tYs8MGsU9X746xXDWgEy9/uZW84qNOh+NzvJksdgA1N+ZPAur6iTQDuNT9\n92uAD1W1QlX3AF8CwVVzs4mM6dGWHu1ieG5BboOmZNamoqqaFxbmM6BLa85ODs7iRt4UGx3Oqzef\nTduWkdz0ylLumrGS2/61guQ2zWwQ+wzdd353osJDeXxOjtOh+BxvJoulQLqIpIhIBHA1MLPmBSKS\nXuNwAnCi/bcdGC0uzYFBgE2E9oKQEOH2UWnkFZfwYSPXUZizppCdB49x6wjbhtxb2sZE8frN5xAe\nGsLM1bu4fVQqb986xAaxz1BCTCR3jknjsw17bGfak3gtWahqJTAV+AjIAd5S1fUi8qiITHRfNlVE\n1ovIKuAe4Ab3+WeBFsA6XEnnZVVd461Yg92FfdrTNb45f/kst9F24VRVpi/Io1tiC0b3aNso72lq\n17lNM/57+1Bm3zGMn5/fwwaxG+jGISkkt2nGY7OzvTaW54+8+q9KVeeqajdVTVXVx93nHlLVme6/\n36WqvVT1LFUdparr3eePquqV7tcyVPUpb8YZ7EJDhFtHppJTeJj5jfTb1PyNe9hYdIQpI1Idr9Ec\nDDq2iqZXB6s33RgiwkL45YQM8opLmnyXA19mv4IYAC7t15GOraIbrXUxfUEeHVtFc3HfDo0QnTFN\na0zPtpybHs8fP9nEfod2aPY1liwM4JqOOWVkKiu3H2xwSc+lW/ezdOsBfnJuinWJGL8kIjx0UQYl\n5VX8Yd5Gp8PxCfadbL515YAk2sZE8tf5nm8lUZvpC/KIax7BD8+2bciN/0pPjOFHg7rwxtfbySm0\nNcGWLMy3osJDmTy8K4vz9rF824Ezeo+cwsN8tmEPNw1JJjrCihsZ//azsem0jA7n0VlntnljILFk\nYb7jmnM607pZOM+eYevihYV5NI8I5frByY0bmDEOaNUsgnvHdWNJ/j4+Wl/kdDiOsmRhvqNZRBg3\nD03hsw17WLfz0GndW7C/lFlrCrnmnM62atgEjEkDO9M9MYbH52ZzvMK7m276MksW5nuuH5JMTGQY\nzy04vdbF3xblEyJwy7Dv7cxijN8KCw3hoYszKNh/jL9/4Uw5Yl9gycJ8T2x0ONcP6cIH63aTu+eI\nR/fsPVrGm0sL+EG/JNrFRnk5QmOa1tC0eM7LSOTZ+bkUHT7udDiOsGRhanXz0BSiwkJ5br5nxexf\n+XIr5VXVTB5hrQoTmB6c0JPKKuXJD4NzKq0lC1OrNi0iueaczry/ehfb99VdzP7I8QpeXbKV8b3b\nkZrQomkCNKaJdWnTnJuHpfDOih2sKjjodDhNzpKFOaXJw7sSKsL0hXW3Lt74ejtHjldacSMT8KaO\nTiMhJpJHZq0Puqm0lizMKSW2jOLKrCTeWb6D3Ydq76c9XlHFS19sYVhaPJlJrZo4QmOaVovIMO4/\nvzsrtx/kfS/VgPFVlixMnaaMSKVKlRc/z6/19fdW7qT4SBm3jrRWhQkOl/dPIjMplt9+kENJWaXT\n4aCqlFd6f3dcSxamTp3imnHpWR1545tt7D1a9p3XqqqVFxbmkZkUy5BUK9tpgkNIiPCrizMoOlzG\n8/V00Xrbht2HufrFr3jyQ++X+7FkYep126hUyiqr+cdJc8w/XLebrftKrbiRCToDusRxyVkdeOHz\nfAr21z0BxBsOHavg4ZnrmfDnL9hYdIT0RO9PLLFkYeqVmtCCC/u057Ul2zhUWgG4ixstzKVrfHPO\n69XO4QiNaXrTxvcgVIQnPmi6Ip7V1cpbSwsY/fQCXluylUkDOzH/3pFNsmmnJQvjkdtHpnG0rJJX\nl2wFYNHmvazbeZgpI1IJteJGJgi1j43m1pGpzFlbyFf5DdvW3xOrCw5y2fTF3P/OGpLjmzNz6jB+\nfWkfWjeP8PqzwZKF8VBGh5aM6dGWf3y5hZKySqYvyKNdyygu6WfFjUzwmjy8Kx1bRfPIrGyqqr0z\nlXbf0TKmvbOGS5/7kp0HjvGHq/ry9pTB9O7YtJURLVkYj90+Oo2DpRX84p01LMnfx4/PTSEyzLYh\nN8ErKjyUBy7sQU7hYd5cWtCo711ZVc1rS7Yy6ukFvL18B7cMTWH+fSP4Qf8kR8YIw5r8icZv9e/c\nmqFpbZi9ppDY6HCuHmjFjYyZ0Kc9ryVv4/cfb2RCZntioxu+4/I3W/bz0Pvr2LD7CEPT2vDwxb1I\nT4xphGjPnLUszGmZOiodgBsGd6FFpP2uYYyI8NDFGewvLecvn25u0HsVHT7Oz2as5KoXlnD4WAXP\nXduff95yjuOJAqxlYU7T4NQ2vPXTwfTt1LT9pcb4st4dY/lhVideWbyVSed0Pu090sorq3ll8Rb+\n9MlmKqqUO0ancevIVJpF+M6PaGtZmNM2MCXOxiqMOcm953UnOjyUx+fknNZ9izYXM/5Pn/ObuRsY\n1LUNH989nHvP6+5TiQKsZWGMMY0iISaSO8ek8/jcHOZv3MOo7m3rvH7HgVJ+PTuHD9fvpkubZvzj\nxixG90hsomhPnyULY4xpJDcMSeaNb7bz2OxshqXFEx76/c6b4xVVvLAwn+cW5CICPz+/O7cMSyEq\n3Ldb69YNZYwxjSQiLIRfTuhJfnEJry3Z9p3XVJV52UWMe2Yhz3yyibE9E/n03pHcPirN5xMFWMvC\nGGMa1egebRneLYE/frKJS8/qQJsWkeQXH+WRWdks3FRMetsWvPHjcxiSFu90qKfFkoUxxjQiEeH/\nJvTkgj8t4rcfbCAhJpKXFuUTGRbKLyf05IYhybV2T/k6SxbGGNPI0hNj+NGgLryyeCvgqoHxi/Hd\naRsT5WxgDWDJwhhjvODucd0AuLhvewZ0iXM4moazZGGMMV4QGx3OwxN7OR1Go/G/jjNjjDFNzpKF\nMcaYenk1WYjIBSKyUURyRWRaLa9PEZG1IrJKRL4QkYwar2WKyBIRWe++xn9Hhowxxs95LVmISCjw\nLDAeyAAm1UwGbm+oah9VPQt4EviD+94w4J/AFFXtBYwEKrwVqzHGmLp5s2UxEMhV1XxVLQdmAJfU\nvEBVD9c4bA6cKDV1HrBGVVe7r9unqlVejNUYY0wdvJksOgI1S0ftcJ/7DhG5XUTycLUs7nSf7gao\niHwkIitE5H4vxmmMMaYe3kwWtdX9+16RWlV9VlVTgV8Av3SfDgOGAde6/7xMRMZ87wEik0VkmYgs\nKy4ubrzIjTHGfIc3k8UOoFON4yRgVx3XzwAurXHvQlXdq6qlwFyg/8k3qOqLqpqlqlkJCQmNFLYx\nxpiTeXNR3lIgXURSgJ3A1cA1NS8QkXRVPVGHcAJw4u8fAfeLSDOgHBgBPFPXw5YvX75XRLbVdY0P\niAf2Oh1EIwmUzxIonwPss/gqX/8sXTy5yGvJQlUrRWQqrh/8ocA/VHW9iDwKLFPVmcBUERmLa6bT\nAeAG970HROQPuBKOAnNVdU49z/P5poWILFPVLKfjaAyB8lkC5XOAfRZfFSifxavbfajqXFxdSDXP\nPVTj73fVce8/cU2fNcYY4zBbwW2MMaZeliya1otOB9CIAuWzBMrnAPssviogPouofm82qzHGGPMd\n1rIwxhhTL0sWTUREQkVkpYjMdjqWhhCRViLytohsEJEcERnsdExnSkTudm9UuU5E/u1Pm1WKyD9E\nZI+IrKtxLk5E5onIZvefrZ2M0VOn+CxPuf+NrRGR90SklZMxeqK2z1HjtftEREXEvwpv12DJounc\nBeQ4HUQj+BPwoar2APrip59JRDri2l4mS1V745refbWzUZ2WV4ALTjo3DfhUVdOBT93H/uAVvv9Z\n5gG9VTUT2AQ80NRBnYFX+P7nQEQ6AeOA7U0dUGOyZNEERCQJ16LDl5yOpSFEpCUwHPg7gKqWq+pB\nZ6NqkDAg2r3LcTPq3mHAp6jq58D+k05fArzq/vur/G9HBJ9W22dR1Y9VtdJ9+BWuHSB82in+n4Br\nQfH91LLdkT+xZNE0/ojrH0u104E0UFegGHjZ3aX2kog0dzqoM6GqO4Gncf22VwgcUtWPnY2qwRJV\ntRDA/Wdbh+NpLDcDHzgdxJkQkYnAzhM7aPszSxZeJiIXAXtUdbnTsTSCMFx7dE1X1X5ACf7T1fEd\n7v78S4AUoAPQXESuczYqczIReRCoBP7ldCyny71d0YPAQ/Vd6w8sWXjfUGCiiGzFtVniaBHx15Xp\nO4Adqvq1+/htatng0U+MBbaoarGqVgDvAkMcjqmhikSkPYD7zz0Ox9MgInIDcBFwrfrnHP9UXL+M\nrHZ//ycBK0SknaNRnSFLFl6mqg+oapKqJuMaQP1MVf3yN1hV3Q0UiEh396kxQLaDITXEdmCQiDQT\nEcH1WfxysL6Gmbj3V3P/+b6DsTSIiFyAq2zBRPfO035HVdeqaltVTXZ//+8A+ru/j/yOJQtzuu4A\n/iUia4CzgN84HM8ZcbeO3gZWAGtxfS/4zUpbEfk3sAToLiI7ROQW4AlgnIhsxjX75gknY/TUKT7L\nX4EYYJ6IrBKR5x0N0gOn+BwBw1ZwG2OMqZe1LIwxxtTLkoUxxph6WbIwxhhTL0sWxhhj6mXJwhhj\nTL0sWRhzEvfuoL+vcXyfiDxczz3JInJNA5+7uCH3G+NNliyM+b4y4AenuZ10MtCgZKGq/r6C3AQw\nSxbGfF8lrgV6d5/GPU8A57oXkN0tIlEi8rKIrHVvujgKQERuFJH3ReRDEdkoIr868QYicrTG3+93\n37taRJ5wn7tTRLLdNR5mNM5HNcYzYU4HYIyPehZYIyJPenj9NOA+Vb0IQETuBVDVPiLSA/hYRLq5\nrx0I9AZKgaUiMkdVl514IxEZj2t78XNUtVRE4mo8I0VVy/yhGJAJLNayMKYWqnoYeA1XgaQzMQx4\n3f1eG4BtwIlkMU9V96nqMVwbGA476d6xwMsn9kRS1RM1Etbg2mrlOlytH2OajCULY07tj8AtwJnU\n7JA6Xjt5j52Tj6WWc+AqoPUsMABY7i7aZEyTsGRhzCm4f6N/C1fCqM8RXBvfnfA5cC2Au/upM7DR\n/do4d73saFzdTV+e9F4fAze76yGcqK0dAnRS1fm4Cmm1Alqc0Qcz5gxYsjCmbr8Hvp0VJSITReTR\nWq5bA1S6B6TvBp4DQkVkLfAmcKOqlrmv/QJXF9Uq4J2a4xUAqvohru3Gl4nIKuA+XDXC/+l+v5XA\nM35e0tb4Gdt11pgmJCI3AlmqOtXpWIw5HdayMMYYUy9rWRhjjKmXtSyMMcbUy5KFMcaYelmyMMYY\nUy9LFsYYY+plycIYY0y9LFkYY4yp1/8HeL3PawVzmVkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111f871d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 5.01 s, sys: 327 ms, total: 5.34 s\n",
      "Wall time: 3.95 s\n"
     ]
    }
   ],
   "source": [
    "%pylab inline --no-import-all\n",
    "%time evaluate_graph(dictionary, corpus, texts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If all went \"according to plan\", you should see that the best topic coherence is found at 3 topics, explaining our choice above (which was correct after finding this result).\n",
    "\n",
    "## LSI conclusions\n",
    "\n",
    "LSI clustering is relatively simple to implement.\n",
    "Its biggest (in fact, quite significant) drawback maybe is that it scales poorly with corpus size.\n",
    "Finally, note that a trained LSI model can be used to predict the topic/cluster distribution in new, unseen documents, just like this:\n",
    "\n",
    "`doc_topics = lsi_model[dictionary.doc2bow(preprocess(raw_doc))]`\n",
    "\n",
    "***\n",
    "***\n",
    "***\n",
    "\n",
    "## LDA clustering\n",
    "\n",
    "Similar to LSI, LDA can also be used to bootstrap clusters, and again a `num_topics` hyper-parameter is the \"hardest\" part about setting up this clustering algorithm.\n",
    "\n",
    "**WARNING**: LDA more expensive to run; expect long computations ahead if you are working with old hardware."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the \"best\" num_topics was estimated by \n",
    "# replacing `LsiModel` in `evaluate_graph` with `LdaModel`\n",
    "ldamodel = LdaModel(\n",
    "    corpus=corpus,\n",
    "    num_topics=11,\n",
    "    id2word=dictionary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now onto somthing cool - visually exploring the LDA model's results. Note that you need [`pyLDAvis`](https://nbviewer.jupyter.org/github/bmabey/pyLDAvis/blob/master/notebooks/pyLDAvis_overview.ipynb) installed for this step, e.g.:\n",
    "\n",
    "- `conda install pyLDAvis`\n",
    "- `pip3 install pyLDAvis`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pyLDAvis.gensim\n",
    "pyLDAvis.enable_notebook()\n",
    "# enable this line:\n",
    "# (disabled to save horizontal screen space)\n",
    "#pyLDAvis.gensim.prepare(ldamodel, corpus, dictionary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The circles' sizes are proportional to the proportions of the topics, using the \"default\" topic if no term is selected. The red bars represent the topic's weight for that term, the blue bars the overall term frequency. $\\lambda$ controls how term relevance is ranked. With a selected term, the circle sizes represent topic-term proportions. More information can be found in the original [R Vignette](https://cran.r-project.org/web/packages/LDAvis/vignettes/details.pdf).\n",
    "\n",
    "## LDA coherence score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.32927206726201924"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ldatopics = ldamodel.show_topics(formatted=False)\n",
    "topics = [[word for word, prob in topic]\n",
    "          for topicid, topic in ldatopics]\n",
    "\n",
    "CoherenceModel(topics=topics,\n",
    "               texts=texts,\n",
    "               dictionary=dictionary,\n",
    "               window_size=10).get_coherence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described before, the LDA model's `num_topics` was set up to use the value that produces the highest coherence score.\n",
    "Yet, the coherence score still is a distant second to what LSI produced.\n",
    "So how do we use LDA \"correctly\" to get good (better than LSI) results?\n",
    "\n",
    "\n",
    "## Cherry-picking the best LDA topics\n",
    "\n",
    "One of the problems with LDA is that if we train it on a large number of topics, the topics get \"lost\" among the numbers. With the below function, we sort each *topic* in a LDA model with 100 topics (the default number of topics for `LdaModel`) and then go cherry-picking the  model that has at least one very coherent topic, and then will select the most coherent topics among it.\n",
    "\n",
    "**WARNING:** LDA is very resource-hungry, computationally speaking; And this particular \"setup\" will take quite some time even on a modern laptop. If you have an old machine, this search can take dozens of minutes. On my machine MBP, this cell routinely takes around 2-5 minutes to complete.\n",
    "\n",
    "**HUGE WARNING**: This is LDA -- and as such, the loop below might not converege; If after around ten runs (print statement) you get no result (because no topic reached the required coherence), you are better of restarting this cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(43, 0.94240695001394914)\n",
      "(85, 0.94352607913652609)\n",
      "(74, 0.78718310305528516)\n",
      "(41, 0.87722206782657075)\n",
      "(74, 0.93087588176716329)\n",
      "(33, 0.86980675820426723)\n",
      "(16, 0.97137951353996921)\n",
      "CPU times: user 1min 23s, sys: 1.93 s, total: 1min 25s\n",
      "Wall time: 1min 25s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "import operator\n",
    "\n",
    "def find_top_topics():\n",
    "    top_topics = [(0, 0)]\n",
    "    \n",
    "    # create models until the top topic\n",
    "    # reaches the indicated coherence score:\n",
    "    while top_topics[0][1] < 0.95:\n",
    "        \n",
    "        lm = LdaModel( # num_topics=100 by default\n",
    "            workers=4, # from LdaMuliticore\n",
    "            corpus=corpus,\n",
    "            id2word=dictionary)\n",
    "        coherence_values = {}\n",
    "        topic_idx = lm.show_topics(\n",
    "            num_topics=-1,\n",
    "            formatted=False)\n",
    "        \n",
    "        for n, topic in topic_idx:\n",
    "            topic = [word for word, _ in topic]\n",
    "            \n",
    "            cm = CoherenceModel(\n",
    "                topics=[topic],\n",
    "                texts=texts,\n",
    "                dictionary=dictionary, \n",
    "                window_size=10)\n",
    "            \n",
    "            coherence_values[n] = cm.get_coherence()\n",
    "            \n",
    "        top_topics = sorted(\n",
    "            coherence_values.items(),\n",
    "            key=operator.itemgetter(1),\n",
    "            reverse=True)\n",
    "        print(top_topics[0])\n",
    "\n",
    "    return lm, top_topics\n",
    "\n",
    "lsildamodel, top_topics = find_top_topics()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('fire', 0.025204164466972451),\n",
       " ('line', 0.012416616962475485),\n",
       " ('north', 0.011203253029828619),\n",
       " ('burn', 0.0084761298283863428),\n",
       " ('south', 0.0081300712753220542),\n",
       " ('sydney', 0.0080226980343572445),\n",
       " ('containment', 0.0079644472525016736),\n",
       " ('rural_fire', 0.0076110078396065635),\n",
       " ('state', 0.0073798836120755728),\n",
       " ('west', 0.0072537481849166453)]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lsildamodel.show_topic(top_topics[0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Topic 16</h3><p>Coherence score: 0.971</p><p>\"fire\"=2.52e-02 + \"line\"=1.24e-02 + \"north\"=1.12e-02 + \"burn\"=8.48e-03 + \"south\"=8.13e-03 + \"sydney\"=8.02e-03 + \"containment\"=7.96e-03 + \"rural fire\"=7.61e-03 + \"state\"=7.38e-03 + \"west\"=7.25e-03</p>\n",
       "<h3>Topic 23</h3><p>Coherence score: 0.870</p><p>\"arafat\"=1.93e-02 + \"palestinian\"=1.89e-02 + \"israeli\"=1.42e-02 + \"west bank\"=9.75e-03 + \"attack\"=9.66e-03 + \"kill\"=9.16e-03 + \"israel\"=8.99e-03 + \"security\"=6.17e-03 + \"government\"=5.76e-03 + \"militants\"=5.52e-03</p>\n",
       "<h3>Topic 52</h3><p>Coherence score: 0.830</p><p>\"play\"=9.22e-03 + \"start\"=8.61e-03 + \"adelaide\"=8.34e-03 + \"macgill\"=8.26e-03 + \"innings\"=7.89e-03 + \"different\"=7.70e-03 + \"warne\"=7.18e-03 + \"australia\"=7.05e-03 + \"south africa\"=6.99e-03 + \"match\"=6.82e-03</p>\n",
       "<h3>Topic 45</h3><p>Coherence score: 0.822</p><p>\"afghanistan\"=1.65e-02 + \"force\"=1.35e-02 + \"united states\"=7.96e-03 + \"lade\"=7.91e-03 + \"still\"=7.47e-03 + \"believe\"=7.16e-03 + \"osama bin\"=7.12e-03 + \"bin laden\"=6.73e-03 + \"pentagon\"=6.58e-03 + \"mountains\"=5.68e-03</p>\n",
       "<h3>Topic 70</h3><p>Coherence score: 0.807</p><p>\"storm\"=1.49e-02 + \"sydney\"=1.31e-02 + \"fire\"=1.07e-02 + \"tree\"=9.09e-03 + \"around\"=7.67e-03 + \"service\"=7.31e-03 + \"damage\"=7.21e-03 + \"home\"=7.11e-03 + \"work\"=6.95e-03 + \"worst\"=6.88e-03</p>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HTML(\"\\n\".join([\n",
    "    \"<h3>Topic {}</h3><p>Coherence score: {:.3f}</p><p>{}</p>\".format(\n",
    "        topicid, c_v,\n",
    "        \" + \".join('\"{}\"={:.2e}'.format(w.replace('_', \" \"), p)\n",
    "                  for w, p in lsildamodel.show_topic(topicid))\n",
    "    )\n",
    "    for topicid, c_v in top_topics[:5]\n",
    "]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.86031333051745262"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "word_probs = lambda topicid: lsildamodel.show_topic(topicid)\n",
    "lda_lsi_topics = [[word for word, prob in word_probs(t_id)]\n",
    "                  for t_id, c_v in top_topics]\n",
    "\n",
    "CoherenceModel(\n",
    "    topics=lda_lsi_topics[:5], # top n using the results above\n",
    "    texts=texts,\n",
    "    dictionary=dictionary,\n",
    "    window_size=10).get_coherence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "With \"LDA as LSI\", we managed to beat our best LSI model, by selecting the four most coherent topics out of 100 - but at quite some \"computational cost\".\n",
    "\n",
    "## HDP (Hierarchical Dirichlet process) clustering\n",
    "\n",
    "Despite being beyond the content of the lecture, **nonparametric** [HDP clustering by Wang et al. (2011)](http://proceedings.mlr.press/v15/wang11a/wang11a.pdf) should be mentioned here for completeness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N. topics = 20\n",
      "Coherence = 0.559653747665\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<h3>Topic 0</h3>\n",
       "<p>0.005*arafat + 0.005*palestinian + 0.005*israeli + 0.004*sharon + 0.003*government + 0.003*kill + 0.003*group + 0.003*choose + 0.003*west_bank + 0.003*call + 0.003*make + 0.002*australia + 0.002*hamas + 0.002*attack + 0.002*team + 0.002*security + 0.002*target + 0.002*win + 0.002*suicide_attacks + 0.002*three</p>\n",
       "\n",
       "<h3>Topic 1</h3>\n",
       "<p>0.003*taliban + 0.003*airport + 0.002*bomb + 0.002*kill + 0.002*opposition + 0.002*fear + 0.002*kandahar + 0.002*target + 0.002*city + 0.002*end + 0.001*lali + 0.001*leave + 0.001*civilians + 0.001*night + 0.001*near + 0.001*gul + 0.001*report + 0.001*half + 0.001*area + 0.001*around</p>\n",
       "\n",
       "<h3>Topic 2</h3>\n",
       "<p>0.003*match + 0.002*israeli + 0.002*play + 0.002*australia + 0.002*team + 0.002*attack + 0.002*rafter + 0.002*ask + 0.002*double + 0.002*france + 0.002*territories + 0.002*tennis + 0.001*could + 0.001*world + 0.001*guarantee + 0.001*win + 0.001*day + 0.001*defeat + 0.001*still + 0.001*shanab</p>\n",
       "\n",
       "<h3>Topic 3</h3>\n",
       "<p>0.004*storm + 0.003*tree + 0.003*work + 0.002*ses + 0.002*areas + 0.002*get + 0.002*service + 0.002*damage + 0.002*sydney + 0.002*large + 0.002*home + 0.002*hit + 0.002*worst + 0.001*call + 0.001*go + 0.001*around + 0.001*hornsby + 0.001*people + 0.001*discover + 0.001*bring</p>\n",
       "\n",
       "<h3>Topic 4</h3>\n",
       "<p>0.002*rabbani + 0.002*president + 0.002*afghan + 0.002*green + 0.001*band + 0.001*northern_alliance + 0.001*draw + 0.001*know + 0.001*deteriorate + 0.001*later + 0.001*interim + 0.001*massive + 0.001*beatles + 0.001*survive + 0.001*beatle + 0.001*bonn + 0.001*factions + 0.001*international + 0.001*success + 0.001*object</p>\n",
       "\n",
       "<h3>Topic 5</h3>\n",
       "<p>0.003*arrest + 0.003*krishna + 0.003*ash + 0.002*hare + 0.002*benares + 0.002*gang + 0.002*harrison + 0.002*take + 0.002*ceremony + 0.002*bernadette + 0.002*river + 0.001*holy + 0.001*officials + 0.001*members + 0.001*boat_people + 0.001*fan + 0.001*smuggle + 0.001*presidency + 0.001*former + 0.001*rough</p>\n",
       "\n",
       "<h3>Topic 6</h3>\n",
       "<p>0.003*australian + 0.002*five + 0.002*afghanistan + 0.002*go + 0.002*head + 0.002*state + 0.001*troop + 0.001*world + 0.001*hiv + 0.001*today + 0.001*take + 0.001*united_nations + 0.001*observe + 0.001*renew + 0.001*force + 0.001*historian + 0.001*latest + 0.001*british + 0.001*hindu + 0.001*arrive</p>\n",
       "\n",
       "<h3>Topic 7</h3>\n",
       "<p>0.002*australian + 0.002*less + 0.002*day + 0.001*geoff + 0.001*australia + 0.001*lower + 0.001*future + 0.001*collapse + 0.001*new_zealand + 0.001*trend + 0.001*insurance + 0.001*bedside + 0.001*survey + 0.001*annan + 0.001*possible + 0.001*berlin + 0.001*party + 0.001*annual + 0.001*receive + 0.001*first</p>\n",
       "\n",
       "<h3>Topic 8</h3>\n",
       "<p>0.003*party + 0.002*president + 0.001*coalition + 0.001*weather_bureau + 0.001*two + 0.001*welsh + 0.001*congressmen + 0.001*australian + 0.001*third + 0.001*categorically + 0.001*national + 0.001*daryl + 0.001*reduce + 0.001*residents + 0.001*ford + 0.001*rule + 0.001*fight + 0.001*swan + 0.001*international + 0.001*stage</p>\n",
       "\n",
       "<h3>Topic 9</h3>\n",
       "<p>0.004*company + 0.002*austar + 0.002*receive + 0.002*entitlements + 0.002*staff + 0.002*million + 0.002*administrators + 0.001*meet + 0.001*leave + 0.001*cents + 0.001*cease + 0.001*brothers + 0.001*board + 0.001*celebrations + 0.001*save + 0.001*hundreds + 0.001*entitle + 0.001*morning + 0.001*operate + 0.001*today</p>\n",
       "\n",
       "<h3>Topic 10</h3>\n",
       "<p>0.002*guide + 0.002*canyoning + 0.002*adventure_world + 0.002*kill + 0.002*interlaken + 0.002*one + 0.002*three + 0.002*appear + 0.002*people + 0.002*swiss + 0.002*proceed + 0.002*court + 0.001*wall + 0.001*stuart + 0.001*first + 0.001*occur + 0.001*bulabarra + 0.001*massive + 0.001*deaths + 0.001*provision</p>\n",
       "\n",
       "<h3>Topic 11</h3>\n",
       "<p>0.003*harrison + 0.002*george + 0.002*liverpool + 0.002*beatle + 0.002*tonight + 0.002*memory + 0.002*die + 0.002*minute + 0.001*music + 0.001*cold + 0.001*man + 0.001*silence + 0.001*could + 0.001*shape + 0.001*portrayal + 0.001*world + 0.001*lord + 0.001*one + 0.001*lose + 0.001*sale</p>\n",
       "\n",
       "<h3>Topic 12</h3>\n",
       "<p>0.002*afghanistan + 0.002*call + 0.002*would + 0.002*agreement + 0.002*afghan + 0.001*offences + 0.001*two_years + 0.001*compose + 0.001*shot_dead + 0.001*bennett + 0.001*surrender + 0.001*strike + 0.001*stronger + 0.001*beset + 0.001*beretta + 0.001*structure + 0.001*assembly + 0.001*group + 0.001*unidentified + 0.001*govern</p>\n",
       "\n",
       "<h3>Topic 13</h3>\n",
       "<p>0.003*per_cent + 0.003*job + 0.002*olivier + 0.002*drop + 0.002*anz + 0.001*market + 0.001*advertise + 0.001*measure + 0.001*marines + 0.001*flank + 0.001*employment + 0.001*hole + 0.001*fall + 0.001*yassin + 0.001*motives + 0.001*internet + 0.001*sturesteps + 0.001*vulnerable + 0.001*eslake + 0.001*week</p>\n",
       "\n",
       "<h3>Topic 14</h3>\n",
       "<p>0.003*government + 0.002*help + 0.002*lew + 0.002*look + 0.002*truck + 0.002*attribute + 0.002*bid + 0.001*ask + 0.001*would + 0.001*could + 0.001*force + 0.001*dominance + 0.001*interim_government + 0.001*review + 0.001*alejandro + 0.001*assistance + 0.001*qantas + 0.001*spark + 0.001*dictator + 0.001*break</p>\n",
       "\n",
       "<h3>Topic 15</h3>\n",
       "<p>0.003*report + 0.003*company + 0.002*take + 0.002*ernst + 0.002*accident + 0.002*degree + 0.001*sarah + 0.001*hih + 0.001*apra + 0.001*downturn + 0.001*regime + 0.001*commission + 0.001*phenomena + 0.001*solvency + 0.001*tell + 0.001*untrue + 0.001*positive + 0.001*chancellor + 0.001*young + 0.001*cheaper</p>\n",
       "\n",
       "<h3>Topic 16</h3>\n",
       "<p>0.003*friedli + 0.002*director + 0.002*company + 0.002*reply + 0.002*know + 0.002*people + 0.002*make + 0.002*blake + 0.001*organise + 0.001*risk + 0.001*ask + 0.001*rage + 0.001*could + 0.001*murder + 0.001*question + 0.001*draconian + 0.001*appear + 0.001*think + 0.001*swiss + 0.001*mistake</p>\n",
       "\n",
       "<h3>Topic 17</h3>\n",
       "<p>0.002*two + 0.002*howard + 0.002*indonesia + 0.002*summit + 0.002*meet + 0.002*president + 0.002*megawati + 0.002*mogul + 0.001*afraid + 0.001*expect + 0.001*uneven + 0.001*talk + 0.001*palestinian + 0.001*beautiful + 0.001*australia + 0.001*bloomberg + 0.001*jakarta + 0.001*february + 0.001*taiba + 0.001*pad</p>\n",
       "\n",
       "<h3>Topic 18</h3>\n",
       "<p>0.002*virgin + 0.002*launceston + 0.002*morning + 0.002*second + 0.002*airline + 0.002*terminal + 0.001*flight + 0.001*indictment + 0.001*gould + 0.001*share + 0.001*ansett + 0.001*customers + 0.001*melbourne + 0.001*check + 0.001*pitch + 0.001*liabilities + 0.001*underway + 0.001*tell + 0.001*afternoon + 0.001*mean</p>\n",
       "\n",
       "<h3>Topic 19</h3>\n",
       "<p>0.002*sydney + 0.002*storm + 0.002*hit + 0.002*nitrous + 0.002*north + 0.002*damage + 0.002*come + 0.002*disaster + 0.001*carr + 0.001*tamworth + 0.001*volunteer + 0.001*natural + 0.001*areas + 0.001*apartment + 0.001*mental + 0.001*new_south + 0.001*sit + 0.001*could + 0.001*nicholas + 0.001*intimidation</p>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from gensim.models import HdpModel\n",
    "\n",
    "hdpmodel = HdpModel(\n",
    "    corpus=corpus,\n",
    "    id2word=dictionary)\n",
    "\n",
    "topic_words = []\n",
    "\n",
    "for topic in hdpmodel.show_topics():\n",
    "    topic_words.append(\n",
    "        \"<h3>Topic {}</h3>\\n<p>{}</p>\\n\".format(*topic)\n",
    "    )\n",
    "    \n",
    "hdptopics = hdpmodel.show_topics(formatted=False)\n",
    "topics = [[word for word, prob in topic]\n",
    "          for topicid, topic in hdptopics]\n",
    "\n",
    "print(\"N. topics =\", len(topics))\n",
    "print(\"Coherence =\", CoherenceModel(\n",
    "    topics=topics,\n",
    "    texts=texts,\n",
    "    dictionary=dictionary,\n",
    "    window_size=10).get_coherence())\n",
    "\n",
    "HTML(\"\\n\".join(topic_words))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HDP clustering was much faster that LDA, required no estimation of the number of topics, and has acceptable results in terms of coherence. At the very least a solid baseline to drive your LDA-based clustering work against!"
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}