{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll take a dataset of documents in several different categories, and find topics (consisting of groups of words) for them. Knowing the actual categories helps us evaluate if the topics we find make sense.\n",
    "\n",
    "We will try this with two different matrix factorizations: **Singular Value Decomposition (SVD) and Non-negative Matrix Factorization (NMF)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "from sklearn import decomposition\n",
    "from scipy import linalg\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scikit Learn comes with a number of built-in datasets, as well as loading utilities to load several standard external datasets. This is a great resource, and the datasets include Boston housing prices, face images, patches of forest, diabetes, breast cancer, and more. We will be using the newsgroups dataset.\n",
    "\n",
    "Newsgroups are discussion groups on Usenet, which was popular in the 80s and 90s before the web really took off. This dataset includes 18,000 newsgroups posts with 20 topics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "categories = ['alt.atheism', 'talk.religion.misc', 'comp.graphics', 'sci.space']\n",
    "remove = ('headers', 'footers', 'quotes')\n",
    "newsgroups_train = fetch_20newsgroups(subset='train', categories=categories, remove=remove)\n",
    "newsgroups_test = fetch_20newsgroups(subset='test', categories=categories, remove=remove)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['/Users/princegrover/scikit_learn_data/20news_home/20news-bydate-train/comp.graphics/38816',\n",
       "       '/Users/princegrover/scikit_learn_data/20news_home/20news-bydate-train/talk.religion.misc/83741',\n",
       "       '/Users/princegrover/scikit_learn_data/20news_home/20news-bydate-train/sci.space/61092'],\n",
       "      dtype='<U100')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.filenames[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 3, 2])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.target[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2034,), (2034,))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.filenames.shape, newsgroups_train.target.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hi,\n",
      "\n",
      "I've noticed that if you only save a model (with all your mapping planes\n",
      "positioned carefully) to a .3DS file that when you reload it after restarting\n",
      "3DS, they are given a default position and orientation.  But if you save\n",
      "to a .PRJ file their positions/orientation are preserved.  Does anyone\n",
      "know why this information is not stored in the .3DS file?  Nothing is\n",
      "explicitly said in the manual about saving texture rules in the .PRJ file. \n",
      "I'd like to be able to read the texture rule information, does anyone have \n",
      "the format for the .PRJ file?\n",
      "\n",
      "Is the .CEL file format available from somewhere?\n",
      "\n",
      "Rych\n"
     ]
    }
   ],
   "source": [
    "# topic from computer graphics\n",
    "print('\\n'.join(newsgroups_train.data[:1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Seems to be, barring evidence to the contrary, that Koresh was simply\n",
      "another deranged fanatic who thought it neccessary to take a whole bunch of\n",
      "folks with him, children and all, to satisfy his delusional mania. Jim\n",
      "Jones, circa 1993.\n",
      "\n",
      "\n",
      "Nope - fruitcakes like Koresh have been demonstrating such evil corruption\n",
      "for centuries.\n"
     ]
    }
   ],
   "source": [
    "# topic about religion talk or miscellanious\n",
    "\n",
    "print('\\n'.join(newsgroups_train.data[1:2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " >In article <1993Apr19.020359.26996@sq.sq.com>, msb@sq.sq.com (Mark Brader) \n",
      "\n",
      "MB>                                                             So the\n",
      "MB> 1970 figure seems unlikely to actually be anything but a perijove.\n",
      "\n",
      "JG>Sorry, _perijoves_...I'm not used to talking this language.\n",
      "\n",
      "Couldn't we just say periapsis or apoapsis?\n",
      "\n",
      " \n"
     ]
    }
   ],
   "source": [
    "# topic about sci space\n",
    "# definition of *perijove* is the point in the orbit of a satellite of Jupiter nearest the planet's center \n",
    "\n",
    "print('\\n'.join(newsgroups_train.data[2:3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['alt.atheism', 'comp.graphics', 'sci.space', 'talk.religion.misc']"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.target_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['comp.graphics', 'talk.religion.misc', 'sci.space'], dtype='<U18')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# confirming from target names\n",
    "\n",
    "np.array(newsgroups_train.target_names)[newsgroups_train.target[:3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, scikit learn has a method that will extract all the word counts for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectorizer = CountVectorizer(stop_words='english')\n",
    "vectors = vectorizer.fit_transform(newsgroups_train.data).todense() # (documents, vocab)\n",
    "vectors.shape #, vectors.nnz / vectors.shape[0], row_means.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`columns` represent unique words  \n",
    "`rows` represent number of documents   \n",
    "`cells` represent number of occurence of word in document i  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all words unique words from corpus of docs\n",
    "vocab = np.array(vectorizer.get_feature_names())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['00', '000', '0000', ..., 'zware', 'zwarte', 'zyxel'], dtype='<U80')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26576,)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['cosmonauts', 'cosmos', 'cosponsored', 'cost', 'costa', 'costar',\n",
       "       'costing', 'costly', 'costruction', 'costs', 'cosy', 'cote',\n",
       "       'couched', 'couldn', 'council', 'councils', 'counsel',\n",
       "       'counselees', 'counselor', 'count'], dtype='<U80')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab[7000:7020]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Singular Value Decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would clearly expect that the words that appear most frequently in one topic would appear less frequently in the other - otherwise that word wouldn't make a good choice to separate out the two topics. Therefore, we expect the topics to be orthogonal.\n",
    "\n",
    "The SVD algorithm factorizes a matrix into one matrix with orthogonal columns and one with orthogonal rows (orthogonal topics) (along with a diagonal matrix, which contains the relative importance of each factor)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://research.fb.com/wp-content/uploads/2016/11/post00049_image0001.png\" alt=\"\" style=\"width: 60%\"/>\n",
    "\n",
    "(source: [Facebook Research: Fast Randomized SVD](https://research.fb.com/fast-randomized-svd/))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 8s, sys: 2.99 s, total: 1min 11s\n",
      "Wall time: 47.3 s\n"
     ]
    }
   ],
   "source": [
    "%time U, s, Vh = linalg.svd(vectors, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2034, 2034) (2034,) (2034, 26576)\n"
     ]
    }
   ],
   "source": [
    "print(U.shape, s.shape, Vh.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`s` are singular values in decreasing order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGUtJREFUeJzt3X2MHPd93/H3d2bvibzj84miSdqUbPpBCGJJph0msY3GahRZTkzlwYbStGJc1UIQu7XhFK5So2kKFKjVB7sREihVI8GUocR27BgiDCW2KstVClgPRz1LlEVKFk1SFHnio8h73N1v/5jf3i3vdnb2jru3O8vPC1jM7G9md7+cPX7mt7+dnTF3R0REulfU7gJERKS1FPQiIl1OQS8i0uUU9CIiXU5BLyLS5RT0IiJdTkEvItLlFPQiIl1OQS8i0uUK7S4AYN26db5ly5Z2lyEikit79ux5w92Hs9briKDfsmULIyMj7S5DRCRXzOxAI+tp6EZEpMsp6EVEupyCXkSkyynoRUS6nIJeRKTLKehFRLqcgl5EpMvlOugff/UEX/nBT5gqlttdiohIx8p10O85cJLbf7ifYllBLyKSJtdBb2Gq65uLiKTLd9CHpFfOi4iky3fQz/TpRUQkTa6DvsI1diMikirXQa+hGxGRbLkO+gp16EVE0uU66E1dehGRTPkO+jB1Jb2ISKp8B32lQ6+cFxFJle+gD1PlvIhIunwHvek4ehGRLA0HvZnFZvakmX0v3L/MzB41s/1m9k0z6w3tfeH+/rB8S2tKn6Xj6EVE0i2kR/85YG/V/duAr7r7O4CTwM2h/WbgZGj/alivJXTQjYhItoaC3sw2AR8D/ircN+AjwLfDKruAG8L8jnCfsPwaa9EYi05qJiKSrdEe/f8EvghUzge8Fjjl7sVw/xCwMcxvBA4ChOWnw/rNF/YfOrxSRCRdZtCb2a8Dx9x9TzNf2MxuMbMRMxsZHR1d3HNUZpTzIiKpGunR/zLwcTN7FfgGyZDNnwGrzKwQ1tkEHA7zh4HNAGH5SuD43Cd19zvdfZu7bxseHl5U8RqjFxHJlhn07v7H7r7J3bcANwI/dPffAx4CfiesthO4L8zvDvcJy3/oLTospnKaYo3Ri4iku5Dj6P8d8AUz208yBn9XaL8LWBvavwDcemElppvt0SvpRUTSFLJXmeXuPwJ+FOZfAT5QY50J4BNNqC2Tfi4lIpIt17+MrdDQjYhIulwHvb6MFRHJlu+gn/kyVlEvIpIm10GPTlMsIpIp10GvL2NFRLLlO+hNx9GLiGTJd9CHqY6jFxFJl++g19iNiEimXAd9hYZuRETS5TrodRy9iEi2fAe9jqMXEcmU76BXj15EJFOug75CHXoRkXS5DvrZS9Eq6UVE0uQ76MNUPXoRkXT5DnodRy8ikinXQV+hDr2ISLpcB72uGSsiki3fQa9rxoqIZMp30IepevQiIunyHfS68IiISKZcB32lT6+hGxGRdLkOevXoRUSy5Tvo212AiEgO5DroRUQkW66DXteMFRHJlu+gD1N9GSsiki7fQa8vY0VEMnVH0Le3DBGRjpbvoNelBEVEMuU66FGPXkQkU66DXue6ERHJlu+g15VHREQy5TroZ6lLLyKSJtdBr6EbEZFs+Q56fRkrIpIp30GvSwmKiGTKDHoz6zezx8zsaTN73sz+U2i/zMweNbP9ZvZNM+sN7X3h/v6wfEurip/9ZaySXkQkTSM9+kngI+7+XuBK4Doz2w7cBnzV3d8BnARuDuvfDJwM7V8N67XE7LluREQkTWbQe+JsuNsTbg58BPh2aN8F3BDmd4T7hOXXWKuOg9S5bkREMjU0Rm9msZk9BRwDHgBeBk65ezGscgjYGOY3AgcBwvLTwNoaz3mLmY2Y2cjo6OiiijddekREJFNDQe/uJXe/EtgEfAB494W+sLvf6e7b3H3b8PDwhT2XBm9ERFIt6Kgbdz8FPAT8IrDKzAph0SbgcJg/DGwGCMtXAsebUu0cpkF6EZFMjRx1M2xmq8L8APCrwF6SwP+dsNpO4L4wvzvcJyz/obfosBjlvIhItkL2KmwAdplZTLJj+Ja7f8/MXgC+YWb/GXgSuCusfxfwdTPbD5wAbmxB3YAuJSgi0ojMoHf3Z4CrarS/QjJeP7d9AvhEU6rLMPvLWCW9iEianP8yNqEevYhIunwHvc51IyKSKddBj46jFxHJlPOgT+hcNyIi6XId9Bq6ERHJlu+gr8wo6UVEUuU76CvH0SvpRURS5Tvow1RD9CIi6fId9DpNsYhIpnwHfeVSgm2uQ0Skk+U76HUYvYhIplwHfYWOoxcRSdcdQd/uAkREOliug15fxoqIZMt30OvSIyIimfId9OrRi4hk6o6gb28ZIiIdLd9Bjy4lKCKSJd9Br0sJiohkynfQt7sAEZEcyHXQV2joRkQkXa6DvnKa4rKSXkQkVa6DvhAlQV8qK+hFRNLkOujjEPRFBb2ISKpcB30hVo9eRCRLvoM+SspXj15EJF3Ogz706EvlNlciItK5ch30cawxehGRLLkOeh11IyKSLddBr6NuRESy5Tvoww+miiUFvYhImnwHfWXoRr+MFRFJleugNzMig7KGbkREUuU66CHp1atHLyKSLvdBH5mpRy8iUkfug74QmY66ERGpIzPozWyzmT1kZi+Y2fNm9rnQvsbMHjCzfWG6OrSbmd1uZvvN7Bkzu7ql/4DIdBy9iEgdjfToi8AfufsVwHbgM2Z2BXAr8KC7bwUeDPcBPgpsDbdbgDuaXnWVODKdj15EpI7MoHf3I+7+RJh/E9gLbAR2ALvCaruAG8L8DuAeTzwCrDKzDU2vPCioRy8iUteCxujNbAtwFfAosN7dj4RFrwPrw/xG4GDVww6FtpaITEEvIlJPw0FvZoPAd4DPu/uZ6mXu7sCC0tbMbjGzETMbGR0dXchDzxOrRy8iUldDQW9mPSQhf6+7/11oPloZkgnTY6H9MLC56uGbQtt53P1Od9/m7tuGh4cXW3/So9cYvYhIqkaOujHgLmCvu3+latFuYGeY3wncV9V+Uzj6ZjtwumqIp+kKsY6jFxGpp9DAOr8M/AvgWTN7KrT9e+DLwLfM7GbgAPDJsOx+4HpgPzAGfKqpFc8Rm46jFxGpJzPo3f3/AZay+Joa6zvwmQusq2GRDq8UEakr97+MjXXUjYhIXfkP+sjQJWNFRNJ1SdAr6UVE0uQ+6KPI0AWmRETS5T7oY114RESkrtwHfSGK9GWsiEgduQ/6KEJBLyJSR+6DXpcSFBGpL/dBr7NXiojUl/ug1/noRUTqy33Qx1Gkc92IiNSR+6DvifWDKRGRenIf9HGks1eKiNST+6DXGL2ISH25D/o4iijqHAgiIqlyH/Q9sVHUGL2ISKrcB70uDi4iUl/ug76gL2NFROrKfdBrjF5EpL7cB31BY/QiInXlP+g1Ri8iUldXBL3G6EVE0uU+6OMowl3npBcRSZP7oC/EBqBxehGRFPkP+igJevXoRURqy33Qx1GlR6+gFxGpJfdBX+nR61h6EZHach/0y/oKAJwZn25zJSIinSn3Qb9x1QAAR05PtLkSEZHOlPugH+pPevTnJottrkREpDPlPuiXh6Gbc1MKehGRWnIf9IMh6M+qRy8iUlPug36mR6+gFxGpKfdBv6wnBuDsZKnNlYiIdKbcB30UGetX9PHTN861uxQRkY6U+6AH2LR6GSfPTbW7DBGRjtQVQb+sN2ZMR92IiNSUGfRmdreZHTOz56ra1pjZA2a2L0xXh3Yzs9vNbL+ZPWNmV7ey+Iok6DVGLyJSSyM9+q8B181puxV40N23Ag+G+wAfBbaG2y3AHc0ps77lvQUdRy8ikiIz6N39YeDEnOYdwK4wvwu4oar9Hk88Aqwysw3NKjbNQG/MuHr0IiI1LXaMfr27HwnzrwPrw/xG4GDVeodC2zxmdouZjZjZyOjo6CLLSCzvK3BOh1eKiNR0wV/GursDCz5HsLvf6e7b3H3b8PDwBdWwrDdmfLqki4+IiNSw2KA/WhmSCdNjof0wsLlqvU2hraWW9SY/mhqfVq9eRGSuxQb9bmBnmN8J3FfVflM4+mY7cLpqiKdllvUmp0HQIZYiIvMVslYws78B/gmwzswOAf8R+DLwLTO7GTgAfDKsfj9wPbAfGAM+1YKa51nel/ToxyZLMLQUrygikh+ZQe/uv5uy6Joa6zrwmQstaqEGenSqYhGRNF3xy9hKj16HWIqIzNcVQV8Zo39TpyoWEZmnK4L+Lav6AXjt1HibKxER6TxdEfTrh/rpLUT87PhYu0sREek4XRH0UWRsXj3AAQW9iMg8XRH0AJcM9XP83GS7yxAR6ThdE/Q6VbGISG3dE/R9BQW9iEgNXRP06wZ7OXJ6nAmd70ZE5DxdE/S/8q5LmJgu8+OXj7e7FBGRjtI1Qf/zm1YC8PLo2TZXIiLSWbom6FcO9DDUV+BnJ3SIpYhIta4JejNj85plCnoRkTm6JugB3rNhBY/99ARlXWlKRGRGVwX99svXMDZV4vFX517LXETk4tVVQf+Rd18CwFMHT7W5EhGRztFVQb92sI8rNqzg23sOkVwDRUREuiroAT794cvYd+wsu59+rd2liIh0hK4L+o+/dyPv2bCC2/7+RV1xSkSELgz6ODL+w8few2unJ/iT+55rdzkiIm3XdUEP8EvvWMenP3QZf7vnEA+9eKzd5YiItFVXBj3AH137LrZeMsgf3vsE33/+9XaXIyLSNl0b9P09Mff+q1/gnZcO8a//+kkeeOFou0sSEWmLrg16gEtW9PO1338/lw8v59P3jHDHj15ud0kiIkuuq4MeYPXyXnZ/9oNce8V6bvuHF/kv9+9lbKrY7rJERJZM1wc9QG8h4o5//j5ufP9m/tfDr/Ar//1HfOvxg5R0ThwRuQhcFEEPyWGXX/7tn+fbf/CLbFg5wBe/8wwfu/0fefil0XaXJiLSUhdN0Fds27KG7/7hL/Hn/+wqzk0Vuenux7jp7sf48cvHKZbK7S5PRKTprBPOCbNt2zYfGRlZ8tedLJb4+o8PcPuD+zgzUWTt8l52XLmRT2zbxLsvHcLMlrwmEZFGmdked9+Wud7FHPQVZyeLPPzSKN975jUeeOEo0yVneKiPa69Yz/vetprtl69lw8p+Bb+IdBQF/SIdPzvJD144yj/uG+WhF0cZn07Ol7NusJf3bFgRbkO8+9IVvH14kN7CRTf6JSIdQkHfBKWys/fIGfYcOMmzh0/z4utneOnoWaaKyVh+T2y8fXiQreuH2LR6INyWsWn1ABtXDdDfE7f5XyAi3azRoC8sRTF5FUfGz21cyc9tXDnTViyVeeWNc+w9coa9R95k75EzPH3wFP/w3BGmS+fvNNcu7+XSlf1cuqKf9Sv7WT/UzyUr+li/oo9Lhvq5ZKiP1ct76Yn1qUBEWkdBv0CFOOKd64d45/ohdlw5214qO0fPTHD41DiHTo5x8MQ4R05P8PrpcV47PcGTB09x4txUzedc0V9g7WAfa5f3smZ5L2sHk+nqZb2sHOhhVZhWbisGCgz0xPrOQEQaoqBvkjgy3rJqgLesGuD9W9bUXGeqWGb07CTHzkxw9Mwko29OcPzcFCfOTSXTs1McOD7GEz87ycmx6bo/6OqJjcG+AoP9BQb7ehgK88v7Cgz2FRjqT3YGy3qT20BvIUxjlvXELOstMNAbJe09SXtfIdLOQ6QLKeiXUG8hYuOqZPw+S7nsvDlZ5PTYNKfHk9up8SnOjBc5M5HcPzdZ5M2J5HZ2cprRNyf56RvnODtZ5OxEceaL5EbFkdFfiBjojenvSYK/J47oK0T0hvneQkRvHNFTiOiL57SHZdXT85dZmMZhWXK/r8Zz98bJLYq04xG5UC0JejO7DvgzIAb+yt2/3IrX6WZRZDNDNYtVLjvj0yXGpkqMT5UYmy7Ozk+VGJsqzswn6xWZmC4zPl1iYqrEZLHMVKnMVHH2dnaymMxXtU9X5kvled9TXKie2M7fCaTseHrn7ZRmdyo9BaMvPLYQRxQioxAbhciIo2hmvrIsjoyeOFnWE+4n60TElfuREUU2s/5se0RsRhyeM7LZdUXapelBb2Yx8BfArwKHgMfNbLe7v9Ds15L6oshY3pcM5yyVctmZLlfvBDzsBEpMFX3eDqKyM5mu2nnMtNdYNm8HE+6PjRWZKjlTxVJ4jM97TDuZkewA5u0kIuKI83YiM+tYsoOJI5t5bPXy6p3LeTudGjuaOLxmbLPTOKJqPlm3er04YqYtrnru8x4z08a8turnrH6u6uecma88x7w27SCboRUJ8AFgv7u/AmBm3wB2AAr6i0AUGX1RTF+hsw4tdXeKZadUdqZL5TBN7hfLZYql+csr7aWyM112iqG9XPVcyXqz85X75Zn2MqUylMrl5L47pfBaM89TCu3nPbZc8/mniuXZx4baZh9bplxm3mMrtZTDenk7l19cY2cyd0cR1diJVnaU1Z+q4jlt1Z/WCrHRU5nGyafCQpx8qivE5y+vtPfEUfI9WF+Bwb7ku6/B0LlavaynY77zakXQbwQOVt0/BPxCC15HpGFmFoaBuOh/3+CehH2pKvxLnuwQZueZ2SmVZnYQszuc8nlt6c81s3ze8899TWbaSuU5y8NzVL9+7ZqT5ylV7ehmdnZhx1gsl5kozj5vZQdfDMOOxXIynS6VZ9Zf7HDkhpX9DDbwafrfXLOV33jvWxb1Go1q25exZnYLcAvAW9/61naVIXLRMTNiS3rKks2rPm1VdgDT4dPedKk8833X2ckSY5NFzk2VOHluiqcPnaLcwA9SL+R7uEa1IugPA5ur7m8Kbedx9zuBOyH5ZWwL6hARuWAWhnsKOf402IqfZD4ObDWzy8ysF7gR2N2C1xERkQY0vUfv7kUz+yzwfZLDK+929+eb/ToiItKYlozRu/v9wP2teG4REVkYnU1LRKTLKehFRLqcgl5EpMsp6EVEupyCXkSky3XEpQTNbBQ4sMiHrwPeaGI5zdKJdXViTdCZdammxnViXZ1YEzS/rre5+3DWSh0R9BfCzEYauWbiUuvEujqxJujMulRT4zqxrk6sCdpXl4ZuRES6nIJeRKTLdUPQ39nuAlJ0Yl2dWBN0Zl2qqXGdWFcn1gRtqiv3Y/QiIlJfN/ToRUSkjlwHvZldZ2Y/MbP9ZnbrEr7uZjN7yMxeMLPnzexzof1PzeywmT0VbtdXPeaPQ50/MbNfa2Ftr5rZs+H1R0LbGjN7wMz2henq0G5mdnuo6xkzu7oF9byrans8ZWZnzOzz7dhWZna3mR0zs+eq2ha8bcxsZ1h/n5ntbEFN/83MXgyv+10zWxXat5jZeNU2+8uqx7wvvO/7Q92LvqpISk0Lfr+a/f8zpa5vVtX0qpk9FdqXalulZUFb/67mcfdc3khOgfwycDnQCzwNXLFEr70BuDrMDwEvAVcAfwr82xrrXxHq6wMuC3XHLartVWDdnLb/Ctwa5m8Fbgvz1wN/DxiwHXh0Cd6z14G3tWNbAR8GrgaeW+y2AdYAr4Tp6jC/usk1XQsUwvxtVTVtqV5vzvM8Fuq0UPdHm1zTgt6vVvz/rFXXnOX/A/iTJd5WaVnQ1r+rubc89+hnLkLu7lNA5SLkLefuR9z9iTD/JrCX5Fq5aXYA33D3SXf/KbCfpP6lsgPYFeZ3ATdUtd/jiUeAVWa2oYV1XAO87O71fhzXsm3l7g8DJ2q83kK2za8BD7j7CXc/CTwAXNfMmtz9B+5eDHcfIblKW6pQ1wp3f8ST1Lin6t/RlJrqSHu/mv7/s15doVf+SeBv6j1HC7ZVWha09e9qrjwHfa2LkNcL25Ywsy3AVcCjoemz4SPZ3ZWPayxtrQ78wMz2WHJdXoD17n4kzL8OrG9DXZBcbaz6P2K7txUsfNssdX3/kqQHWHGZmT1pZv/XzD5UVeuhJahpIe/XUm+nDwFH3X1fVduSbqs5WdBRf1d5Dvq2M7NB4DvA5939DHAH8HbgSuAIyUfJpfZBd78a+CjwGTP7cPXC0ItZ8kOtLLms5MeBvw1NnbCtztOubZPGzL4EFIF7Q9MR4K3ufhXwBeCvzWzFEpXTce/XHL/L+Z2IJd1WNbJgRif8XeU56Bu6CHmrmFkPyRt7r7v/HYC7H3X3kruXgf/N7JDDktXq7ofD9Bjw3VDD0cqQTJgeW+q6SHY8T7j70VBf27dVsNBtsyT1mdnvA78O/F4ICsLwyPEwv4dkDPyd4fWrh3eaXtMi3q8lex/NrAD8FvDNqnqXbFvVygI67O8qz0HftouQh/HAu4C97v6Vqvbq8e3fBCpHB+wGbjSzPjO7DNhK8oVQs+tabmZDlXmSL/WeC69f+RZ/J3BfVV03hSMBtgOnqz5uNtt5Pa52b6sqC9023weuNbPVYfji2tDWNGZ2HfBF4OPuPlbVPmxmcZi/nGTbvBLqOmNm28Pf5k1V/45m1bTQ92sp/3/+U+BFd58ZklmqbZWWBXTa31WzvtVtx43kG+yXSPbWX1rC1/0gyUexZ4Cnwu164OvAs6F9N7Ch6jFfCnX+hAv4lj+jrstJjm54Gni+sk2AtcCDwD7g/wBrQrsBfxHqehbY1qK6lgPHgZVVbUu+rUh2NEeAaZIx0JsXs21Ixs33h9unWlDTfpLx2srf1l+GdX87vK9PAU8Av1H1PNtIwvdl4M8JP4ZsYk0Lfr+a/f+zVl2h/WvAH8xZd6m2VVoWtPXvau5Nv4wVEelyeR66ERGRBijoRUS6nIJeRKTLKehFRLqcgl5EpMsp6EVEupyCXkSkyynoRUS63P8HR2FZS5VVyFYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Topics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_top_words=8\n",
    "\n",
    "def show_topics(a):\n",
    "    top_words = lambda t: [vocab[i] for i in np.argsort(t)[:-num_top_words-1:-1]]\n",
    "    topic_words = ([top_words(t) for t in a])\n",
    "    return [' '.join(t) for t in topic_words]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['critus ditto propagandist surname galacticentric kindergarten surreal imaginative',\n",
       " 'jpeg gif file color quality image jfif format',\n",
       " 'graphics edu pub mail 128 3d ray ftp',\n",
       " 'jesus god matthew people atheists atheism does graphics',\n",
       " 'image data processing analysis software available tools display',\n",
       " 'god atheists atheism religious believe religion argument true',\n",
       " 'space nasa lunar mars probe moon missions probes',\n",
       " 'image probe surface lunar mars probes moon orbit',\n",
       " 'argument fallacy conclusion example true ad argumentum premises',\n",
       " 'space larson image theory universe physical nasa material']"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(Vh[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that topics are telling something. For example,   \n",
    "`Topic 10` is related to `physics and space`.  \n",
    "`Topic 3` is related to `computer graphics`.  \n",
    "`topic 2` is related to `image quality` etc.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non Negative Matrix Factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://sanghosuh.github.io/lens_nmf-icdm/img/nmf.png\" alt=\"NMF on documents\" style=\"width: 80%\"/>\n",
    "(source: [NMF Tutorial](http://perso.telecom-paristech.fr/~essid/teach/NMF_tutorial_ICME-2014.pdf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://cdn-images-1.medium.com/max/720/1*MLJVWz4EdOFsqhvBxEi9iA.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": []
   },
   "source": [
    "#### nmf on sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "m,n=vectors.shape\n",
    "d=5  # num topics. to be chosen "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = decomposition.NMF(n_components=d, random_state=1)\n",
    "\n",
    "W1 = clf.fit_transform(vectors)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def display_topics(model, n_top):\n",
    "    for topic_idx, topic in enumerate(model.components_):\n",
    "        print(\"Topic %d:\" % (topic_idx))\n",
    "        print(' '.join(vocab[i] for i in topic.argsort()[:-n_top:-1]), '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# no_top_words = 10\n",
    "# display_topics(nmf, tfidf_feature_names, no_top_words)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['00', '000', '0000', ..., 'zware', 'zwarte', 'zyxel'], dtype='<U80')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for each topic, words with highest weight represent that topic. We can visualize top k words of each topic and see whether that topic makes sense or not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([13816, 12642, 11163, 10286])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.components_[0].argsort()[:-5:-1] # last 6 except the last one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.12126083, 0.        , 0.        , ..., 0.00003202, 0.00006403,\n",
       "        0.00028996],\n",
       "       [0.11928575, 0.12085267, 0.00016592, ..., 0.0001165 , 0.000233  ,\n",
       "        0.05186412],\n",
       "       [0.05766391, 0.48160644, 0.0007993 , ..., 0.00028762, 0.00057525,\n",
       "        0.        ],\n",
       "       [0.        , 0.14794335, 0.        , ..., 0.00006585, 0.0001317 ,\n",
       "        0.        ],\n",
       "       [0.12742192, 0.19356417, 0.00050472, ..., 0.        , 0.        ,\n",
       "        0.        ]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Topic 0:\n",
      "jpeg image gif file color images format quality version \n",
      "\n",
      "Topic 1:\n",
      "edu graphics pub mail 128 ray ftp send 3d \n",
      "\n",
      "Topic 2:\n",
      "space launch satellite nasa commercial satellites year market data \n",
      "\n",
      "Topic 3:\n",
      "jesus god people matthew atheists does atheism said just \n",
      "\n",
      "Topic 4:\n",
      "image data available software processing ftp edu analysis images \n",
      "\n"
     ]
    }
   ],
   "source": [
    "display_topics(clf, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For all 5 topics, top 5 words with highest weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "683.3993486212418"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.reconstruction_err_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NMF Using TF-IDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Topic Frequency-Inverse Document Frequency](http://www.tfidf.com/) (TF-IDF) is a way to normalize term counts by taking into account how often they appear in a document, how long the document is, and how commmon/rare the term is.\n",
    "\n",
    "TF = (# occurrences of term t in document) / (# of words in documents)\n",
    "\n",
    "IDF = log(# of documents / # documents with term t in it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorizer_tfidf = TfidfVectorizer(stop_words='english')\n",
    "vectors_tfidf = vectorizer_tfidf.fit_transform(newsgroups_train.data) # (documents, vocab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectors_tfidf.shape # 2034 files and 26576 unique words in corpus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We basically normalized the number of occurences b/w 0 and 1 and this can give better topics as compares to BoW approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 232)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectors.min(), vectors.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1.0)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectors_tfidf.min(), vectors_tfidf.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "W1 = clf.fit_transform(vectors_tfidf)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Topic 0:\n",
      "people don think just like objective say morality know \n",
      "\n",
      "Topic 1:\n",
      "graphics thanks files image file program windows know format \n",
      "\n",
      "Topic 2:\n",
      "space nasa launch shuttle orbit moon lunar earth station \n",
      "\n",
      "Topic 3:\n",
      "ico bobbe tek beauchaine bronx manhattan sank queens vice \n",
      "\n",
      "Topic 4:\n",
      "god jesus bible believe christian atheism does belief faith \n",
      "\n"
     ]
    }
   ],
   "source": [
    "display_topics(clf, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.712926057952814"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.reconstruction_err_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use both tf-idf and countvectorizer to create topics using above NMF. TF-IDF is supposes to have richer information as we are normalizing the numbers and taking into account frequently occuring numbersm/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### End"
   ]
  }
 ],
 "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.4"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}