2-svd-nmf-topic-modeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Topic Modeling with NMF and SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Topic modeling is a fun way to start our study of NLP. We will use two popular **matrix decomposition techniques**. \n",
    "\n",
    "We start with a **term-document matrix**:\n",
    "\n",
    "<img src=\"images/document_term.png\" alt=\"term-document matrix\" style=\"width: 80%\"/>\n",
    "\n",
    "source: [Introduction to Information Retrieval](http://player.slideplayer.com/15/4528582/#)\n",
    "\n",
    "We can decompose this into one tall thin matrix times one wide short matrix (possibly with a diagonal matrix in between).\n",
    "\n",
    "Notice that this representation does not take into account word order or sentence structure.  It's an example of a **bag of words** approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Latent Semantic Analysis (LSA) uses Singular Value Decomposition (SVD)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Motivation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the most extreme case - reconstructing the matrix using an outer product of two vectors. Clearly, in most cases we won't be able to reconstruct the matrix exactly. But if we had one vector with the relative frequency of each vocabulary word out of the total word count, and one with the average number of words per document, then that outer product would be as close as we can get.\n",
    "\n",
    "Now consider increasing that matrices to two columns and two rows. The optimal decomposition would now be to cluster the documents into two groups, each of which has as different a distribution of words as possible to each other, but as similar as possible amongst the documents in the cluster. We will call those two groups \"topics\". And we would cluster the words into two groups, based on those which most frequently appear in each of the topics. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started"
   ]
  },
  {
   "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": 138,
   "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": 139,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Additional Resources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Data source](http://scikit-learn.org/stable/datasets/twenty_newsgroups.html): 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.\n",
    "- [Chris Manning's book chapter](https://nlp.stanford.edu/IR-book/pdf/18lsi.pdf) on matrix factorization and LSI \n",
    "- Scikit learn [truncated SVD LSI details](http://scikit-learn.org/stable/modules/decomposition.html#lsa)\n",
    "\n",
    "### Other Tutorials\n",
    "- [Scikit-Learn: Out-of-core classification of text documents](http://scikit-learn.org/stable/auto_examples/applications/plot_out_of_core_classification.html): uses [Reuters-21578](https://archive.ics.uci.edu/ml/datasets/reuters-21578+text+categorization+collection) dataset (Reuters articles labeled with ~100 categories), HashingVectorizer\n",
    "- [Text Analysis with Topic Models for the Humanities and Social Sciences](https://de.dariah.eu/tatom/index.html): uses [British and French Literature dataset](https://de.dariah.eu/tatom/datasets.html) of Jane Austen, Charlotte Bronte, Victor Hugo, and more"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Look at our data"
   ]
  },
  {
   "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](http://scikit-learn.org/stable/datasets/), 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": 140,
   "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": 141,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2034,), (2034,))"
      ]
     },
     "execution_count": 141,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.filenames.shape, newsgroups_train.target.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at some of the data.  Can you guess which category these messages are in?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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",
      "\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",
      "\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": [
    "print(\"\\n\".join(newsgroups_train.data[:3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "hint: definition of *perijove* is the point in the orbit of a satellite of Jupiter nearest the planet's center "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['comp.graphics', 'talk.religion.misc', 'sci.space'], dtype='<U18')"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(newsgroups_train.target_names)[newsgroups_train.target[:3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The target attribute is the integer index of the category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 3, 2, 0, 2, 0, 2, 1, 2, 1])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.target[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_topics, num_top_words = 6, 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stop words, stemming, lemmatization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stop words"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [Intro to Information Retrieval](https://nlp.stanford.edu/IR-book/html/htmledition/dropping-common-terms-stop-words-1.html):\n",
    "\n",
    "*Some extremely common words which would appear to be of little value in helping select documents matching a user need are excluded from the vocabulary entirely. These words are called stop words.*\n",
    "\n",
    "*The general trend in IR systems over time has been from standard use of quite large stop lists (200-300 terms) to very small stop lists (7-12 terms) to no stop list whatsoever. Web search engines generally do not use stop lists.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### NLTK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['a',\n",
       " 'about',\n",
       " 'above',\n",
       " 'across',\n",
       " 'after',\n",
       " 'afterwards',\n",
       " 'again',\n",
       " 'against',\n",
       " 'all',\n",
       " 'almost',\n",
       " 'alone',\n",
       " 'along',\n",
       " 'already',\n",
       " 'also',\n",
       " 'although',\n",
       " 'always',\n",
       " 'am',\n",
       " 'among',\n",
       " 'amongst',\n",
       " 'amoungst']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.feature_extraction import stop_words\n",
    "\n",
    "sorted(list(stop_words.ENGLISH_STOP_WORDS))[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no single universal list of stop words."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stemming and Lemmatization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "from [Information Retrieval](https://nlp.stanford.edu/IR-book/html/htmledition/stemming-and-lemmatization-1.html) textbook:\n",
    "\n",
    "Are the below words the same?\n",
    "\n",
    "*organize, organizes, and organizing*\n",
    "\n",
    "*democracy, democratic, and democratization*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stemming and Lemmatization both generate the root form of the words. \n",
    "\n",
    "Lemmatization uses the rules about a language.  The resulting tokens are all actual words\n",
    "\n",
    "\"Stemming is the poor-man’s lemmatization.\" (Noah Smith, 2011) Stemming is a crude heuristic that chops the ends off of words.  The resulting tokens may not be actual words. Stemming is faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[nltk_data] Downloading package wordnet to\n",
      "[nltk_data]     /home/racheltho/nltk_data...\n",
      "[nltk_data]   Package wordnet is already up-to-date!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import nltk\n",
    "nltk.download('wordnet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nltk import stem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "wnl = stem.WordNetLemmatizer()\n",
    "porter = stem.porter.PorterStemmer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [],
   "source": [
    "word_list = ['feet', 'foot', 'foots', 'footing']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['foot', 'foot', 'foot', 'footing']"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[wnl.lemmatize(word) for word in word_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['feet', 'foot', 'foot', 'foot']"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[porter.stem(word) for word in word_list]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your turn!  Now, try lemmatizing and stemming the following collections of words:\n",
    "\n",
    "- fly, flies, flying\n",
    "- organize, organizes, organizing\n",
    "- universe, university"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "fastai/course-nlp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stemming and lemmatization are language dependent.  Languages with more complex morphologies may show bigger benefits.  For example, Sanskrit has a very [large number of verb forms](https://en.wikipedia.org/wiki/Sanskrit_verbs). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spacy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stemming and lemmatization are implementation dependent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spacy is a very modern & fast nlp library. Spacy is opinionated, in that it typically offers one highly optimized way to do something (whereas nltk offers a huge variety of ways, although they are usually not as optimized).\n",
    "\n",
    "You will need to install it.\n",
    "\n",
    "if you use conda:\n",
    "```\n",
    "conda install -c conda-forge spacy\n",
    "```\n",
    "if you use pip:\n",
    "```\n",
    "pip install -U spacy\n",
    "```\n",
    "\n",
    "You will then need to download the English model:\n",
    "```\n",
    "spacy download en_core_web_sm\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import spacy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "nlp = spacy.load(\"en_core_web_sm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note in the video we used\n",
    "```\n",
    "from spacy.lemmatizer import Lemmatizer\n",
    "lemmatizer = Lemmatizer()\n",
    "```\n",
    "which creates an *empty* lemmatizer. \n",
    "\n",
    "This has been replaced with an English specific lemmatizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "lemmatizer = nlp.Defaults.create_lemmatizer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['foot', 'foot', 'foot', 'foot']"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[lemmatizer.lookup(word) for word in word_list]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spacy doesn't offer a stemmer (since lemmatization is considered better-- this is an example of being opinionated!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stop words vary from library to library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['a',\n",
       " 'about',\n",
       " 'above',\n",
       " 'across',\n",
       " 'after',\n",
       " 'afterwards',\n",
       " 'again',\n",
       " 'against',\n",
       " 'all',\n",
       " 'almost',\n",
       " 'alone',\n",
       " 'along',\n",
       " 'already',\n",
       " 'also',\n",
       " 'although',\n",
       " 'always',\n",
       " 'am',\n",
       " 'among',\n",
       " 'amongst',\n",
       " 'amount']"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted(list(nlp.Defaults.stop_words))[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise: What stop words appear in spacy but not in sklearn?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ca',\n",
       " 'did',\n",
       " 'does',\n",
       " 'doing',\n",
       " 'just',\n",
       " 'make',\n",
       " 'quite',\n",
       " 'really',\n",
       " 'regarding',\n",
       " 'say',\n",
       " 'unless',\n",
       " 'used',\n",
       " 'using',\n",
       " 'various'}"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Exercise:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Exercise: And what stop words are in sklearn but not spacy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "frozenset({'amoungst',\n",
       "           'bill',\n",
       "           'cant',\n",
       "           'co',\n",
       "           'con',\n",
       "           'couldnt',\n",
       "           'cry',\n",
       "           'de',\n",
       "           'describe',\n",
       "           'detail',\n",
       "           'eg',\n",
       "           'etc',\n",
       "           'fill',\n",
       "           'find',\n",
       "           'fire',\n",
       "           'found',\n",
       "           'hasnt',\n",
       "           'ie',\n",
       "           'inc',\n",
       "           'interest',\n",
       "           'ltd',\n",
       "           'mill',\n",
       "           'sincere',\n",
       "           'system',\n",
       "           'thick',\n",
       "           'thin',\n",
       "           'un'})"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Exercise:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### When to use these?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<img src=\"images/skomoroch.png\" alt=\"\" style=\"width: 65%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "These were long considered standard techniques, but they can often **hurt** your performance **if using deep learning**. Stemming, lemmatization, and removing stop words all involve throwing away information.\n",
    "\n",
    "However, they can still be useful when working with simpler models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Another approach: sub-word units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[SentencePiece](https://github.com/google/sentencepiece) library from Google"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, scikit learn has a method that will extract all the word counts for us.  In the next lesson, we'll learn how to write our own version of CountVectorizer, to see what's happening underneath the hood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nltk\n",
    "# nltk.download('punkt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from nltk import word_tokenize\n",
    "\n",
    "# class LemmaTokenizer(object):\n",
    "#     def __init__(self):\n",
    "#         self.wnl = stem.WordNetLemmatizer()\n",
    "#     def __call__(self, doc):\n",
    "#         return [self.wnl.lemmatize(t) for t in word_tokenize(doc)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorizer = CountVectorizer(stop_words='english') #, tokenizer=LemmaTokenizer())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectors = vectorizer.fit_transform(newsgroups_train.data).todense() # (documents, vocab)\n",
    "vectors.shape #, vectors.nnz / vectors.shape[0], row_means.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2034 (2034, 26576)\n"
     ]
    }
   ],
   "source": [
    "print(len(newsgroups_train.data), vectors.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "vocab = np.array(vectorizer.get_feature_names())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26576,)"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "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": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab[7000:7020]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singular Value Decomposition (SVD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"SVD is not nearly as famous as it should be.\" - Gilbert Strang"
   ]
  },
  {
   "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** (along with a diagonal matrix, which contains the **relative importance** of each factor).\n",
    "\n",
    "<img src=\"images/svd_fb.png\" alt=\"\" style=\"width: 80%\"/>\n",
    "(source: [Facebook Research: Fast Randomized SVD](https://research.fb.com/fast-randomized-svd/))\n",
    "\n",
    "SVD is an **exact decomposition**, since the matrices it creates are big enough to fully cover the original matrix. SVD is extremely widely used in linear algebra, and specifically in data science, including:\n",
    "\n",
    "- semantic analysis\n",
    "- collaborative filtering/recommendations ([winning entry for Netflix Prize](https://datajobs.com/data-science-repo/Recommender-Systems-%5BNetflix%5D.pdf))\n",
    "- calculate Moore-Penrose pseudoinverse\n",
    "- data compression\n",
    "- principal component analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Latent Semantic Analysis (LSA) uses SVD.  You will sometimes hear topic modelling referred to as LSA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 19s, sys: 10 s, total: 1min 29s\n",
      "Wall time: 3.63 s\n"
     ]
    }
   ],
   "source": [
    "%time U, s, Vh = linalg.svd(vectors, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "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": [
    "Confirm this is a decomposition of the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([433.92698542, 291.51012741, 240.71137677, 220.00048043])"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s[:4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([433.92698542, 291.51012741, 240.71137677, 220.00048043])"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diag(np.diag(s[:4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Exercise: confrim that U, s, Vh is a decomposition of `vectors`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confirm that U, V are orthonormal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Exercise: Confirm that U, Vh are orthonormal\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "What can we say about the singular values s?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGXBJREFUeJzt3X2MHPd93/H3d2bvibzj85GiSdqUbNqWEMSSTDtMHBuNlSqynJjKgwqlacW6qoUgdmvDKRKlRtMUKFCrD3YqJFCqRoIpQ4nt2DFEGHJsQZaqFLAlkXqWKIuULIln0uSJz+Q97u63f8xv75Z3Ozt7x9uHWX1ewGJmfzO7++Xs8TO//e3sjLk7IiLSvaJ2FyAiIs2loBcR6XIKehGRLqegFxHpcgp6EZEup6AXEelyCnoRkS6noBcR6XIKehGRLldodwEA69at861bt7a7DBGRXNm3b9+b7j6ctV5HBP3WrVvZu3dvu8sQEckVM3u9kfU0dCMi0uUU9CIiXU5BLyLS5RT0IiJdTkEvItLlFPQiIl1OQS8i0uVyHfRPvHaCL33/x0wVy+0uRUSkY+U66Pe9fpI7fnCQYllBLyKSJtdBb2Gq65uLiKTLd9CHpFfOi4iky3fQz/TpRUQkTa6DvsI1diMikirXQa+hGxGRbLkO+gp16EVE0uU66E1dehGRTPkO+jB1Jb2ISKp8B32lQ6+cFxFJle+gD1PlvIhIunwHvek4ehGRLA0HvZnFZvaUmX0n3L/UzB4zswNm9nUz6w3tfeH+wbB8a3NKn6Xj6EVE0i2kR/9ZYH/V/duBL7v7NuAkcEtovwU46e7vAr4c1msKHXQjIpKtoaA3s83Ax4G/DvcN+CjwzbDKbuCGML8z3Ccsv8aaNMaik5qJiGRrtEf/58AfAZXzAa8FTrl7MdwfATaF+U3AIYCw/HRYf+mF/YcOrxQRSZcZ9Gb268Axd99X3VxjVW9gWfXz3mpme81s7+joaEPFznuO1GcXEZGKRnr0HwI+YWavAV8jGbL5c2CVmRXCOpuBw2F+BNgCEJavBE7MfVJ3v8vdt7v79uHh4UUVrzF6EZFsmUHv7n/i7pvdfStwE/ADd/894GHgd8Jqu4D7w/yecJ+w/AfepMNiKqcp1hi9iEi6izmO/o+Bz5vZQZIx+LtD+93A2tD+eeC2iysx3WyPXkkvIpKmkL3KLHd/BHgkzL8KfLDGOhPAjUtQWyb9XEpEJFuufxlboaEbEZF0uQ56fRkrIpIt30E/82Wsol5EJE2ugx6dplhEJFOug15fxoqIZMt30JuOoxcRyZLvoA9THUcvIpIu30GvsRsRkUy5DvoKDd2IiKTLddDrOHoRkWz5DnodRy8ikinfQa8evYhIplwHfYU69CIi6XId9LOXolXSi4ikyXfQh6l69CIi6fId9DqOXkQkU66DvkIdehGRdLkOel0zVkQkW76DXteMFRHJlO+gD1P16EVE0uU76HXhERGRTLkO+kqfXkM3IiLpch306tGLiGTLd9C3uwARkRzIddCLiEi2XAe9rhkrIpIt30EfpvoyVkQkXb6DXl/Giohk6o6gb28ZIiIdLd9Br0sJiohkynXQox69iEimXAe9znUjIpIt30GvK4+IiGTKddDPUpdeRCRNroNeQzciItnyHfT6MlZEJFO+g16XEhQRyZQZ9GbWb2aPm9kzZvaCmf3n0H6pmT1mZgfM7Otm1hva+8L9g2H51mYVP/vLWCW9iEiaRnr0k8BH3f19wJXAdWa2A7gd+LK7bwNOAreE9W8BTrr7u4Avh/WaYvZcNyIikiYz6D1xLtztCTcHPgp8M7TvBm4I8zvDfcLya6xZx0HqXDciIpkaGqM3s9jMngaOAQ8CrwCn3L0YVhkBNoX5TcAhgLD8NLC2xnPeamZ7zWzv6Ojoooo3XXpERCRTQ0Hv7iV3vxLYDHwQuLzWamFaK33n9bnd/S533+7u24eHhxutt3Z9GrwREUm1oKNu3P0U8AiwA1hlZoWwaDNwOMyPAFsAwvKVwImlKHYu0yC9iEimRo66GTazVWF+APhVYD/wMPA7YbVdwP1hfk+4T1j+A2/SYTHKeRGRbIXsVdgI7DazmGTH8A13/46ZvQh8zcz+C/AUcHdY/27gq2Z2kKQnf1MT6gZ0KUERkUZkBr27PwtcVaP9VZLx+rntE8CNS1JdhtlfxirpRUTS5PyXsQn16EVE0uU76HWuGxGRTLkO+tpHcoqISLWcB31C57oREUmX66DX0I2ISLZ8B31lRkkvIpIq30FfOY5eSS8ikirfQR+mGqIXEUmX76DXaYpFRDLlO+grlxJscx0iIp0s30Gvw+hFRDLlOugrdBy9iEi67gj6dhcgItLBch30+jJWRCRbvoNelx4REcmU76BXj15EJFN3BH17yxAR6Wj5Dnp0KUERkSz5DnpdSlBEJFO+g77dBYiI5ECug75CQzciIulyHfSV0xSXlfQiIqlyHfSFKAn6UllBLyKSJtdBH4egLyroRURS5TroC7F69CIiWfId9FFSvnr0IiLpch70oUdfKre5EhGRzpXroI9jjdGLiGTJddDrqBsRkWy5DnoddSMiki3fQR9+MFUsKehFRNLkO+grQzf6ZayISKpcB72ZERmUNXQjIpIq10EPSa9ePXoRkXS5D/rITD16EZE6ch/0hch01I2ISB2ZQW9mW8zsYTPbb2YvmNlnQ/saM3vQzA6E6erQbmZ2h5kdNLNnzezqpv4DItNx9CIidTTSoy8Cf+julwM7gE+b2RXAbcBD7r4NeCjcB/gYsC3cbgXuXPKqq8SR6Xz0IiJ1ZAa9ux9x9yfD/FlgP7AJ2AnsDqvtBm4I8zuBez3xI2CVmW1c8sqDgnr0IiJ1LWiM3sy2AlcBjwEb3P0IJDsDYH1YbRNwqOphI6GtKSJT0IuI1NNw0JvZIPAt4HPufqbeqjXa5iWxmd1qZnvNbO/o6GijZcwTq0cvIlJXQ0FvZj0kIX+fu/99aD5aGZIJ02OhfQTYUvXwzcDhuc/p7ne5+3Z33z48PLzY+pMevcboRURSNXLUjQF3A/vd/UtVi/YAu8L8LuD+qvabw9E3O4DTlSGeZijEOo5eRKSeQgPrfAj4l8BzZvZ0aPsPwBeBb5jZLcAbwI1h2QPA9cBBYAz45JJWPEdsOo5eRKSezKB39/9H7XF3gGtqrO/Apy+yroZFOrxSRKSu3P8yNtZRNyIideU/6CNDl4wVEUnXJUGvpBcRSZP7oI8iQxeYEhFJl/ugj3XhERGRunIf9IUo0pexIiJ15D7oowgFvYhIHbkPel1KUESkvtwHvc5eKSJSX+6DXuejFxGpL/dBH0eRznUjIlJH7oO+J9YPpkRE6sl90MeRzl4pIlJP7oNeY/QiIvXlPujjKKKocyCIiKTKfdD3xEZRY/QiIqlyH/S6OLiISH25D/qCvowVEakr90GvMXoRkfpyH/QFjdGLiNSV/6DXGL2ISF1dEfQaoxcRSZf7oI+jCHedk15EJE3ug74QG4DG6UVEUuQ/6KMk6NWjFxGpLfdBH0eVHr2CXkSkltwHfaVHr2PpRURqy33QL+srAHBmfLrNlYiIdKbcB/2mVQMAHDk90eZKREQ6U+6Dfqg/6dGfnyy2uRIRkc6U+6BfHoZuzk8p6EVEasl90A+GoD+nHr2ISE25D/qZHr2CXkSkptwH/bKeGIBzk6U2VyIi0plyH/RRZGxY0cdP3jzf7lJERDpS7oMeYPPqZZw8P9XuMkREOlJXBP2y3pgxHXUjIlJTZtCb2T1mdszMnq9qW2NmD5rZgTBdHdrNzO4ws4Nm9qyZXd3M4iuSoNcYvYhILY306L8CXDen7TbgIXffBjwU7gN8DNgWbrcCdy5NmfUt7y3oOHoRkRSZQe/ujwIn5jTvBHaH+d3ADVXt93riR8AqM9u4VMWmGeiNGVePXkSkpsWO0W9w9yMAYbo+tG8CDlWtNxLa5jGzW81sr5ntHR0dXWQZieV9Bc7r8EoRkZqW+stYq9FW8/zB7n6Xu2939+3Dw8MX9aLLemPGp0u6+IiISA2LDfqjlSGZMD0W2keALVXrbQYOL768xizrTX40NT6tXr2IyFyLDfo9wK4wvwu4v6r95nD0zQ7gdGWIp5mW9SanQdAhliIi8xWyVjCzvwX+CbDOzEaA/wR8EfiGmd0CvAHcGFZ/ALgeOAiMAZ9sQs3zLO9LevRjkyUYasUriojkR2bQu/vvpiy6psa6Dnz6YotaqIEenapYRCRNV/wyttKj1yGWIiLzdUXQV8boz+pUxSIi83RF0L9tVT8Ah0+Nt7kSEZHO0xVBv2Gon95CxBvHx9pdiohIx+mKoI8iY8vqAV5X0IuIzNMVQQ+wfqif4+cn212GiEjH6Zqg16mKRURq656g7yso6EVEauiaoF832MuR0+NM6Hw3IiIX6Jqg/5X3rGdiuswPXzne7lJERDpK1wT9z29eCcAro+faXImISGfpmqBfOdDDUF+BN07oEEsRkWpdE/RmxpY1yxT0IiJzdE3QA1y+cQWP/+QEZV1pSkRkRlcF/Y7L1jA2VeKJ1+Zey1xE5K2rq4L+o+9NrlH+9KFTba5ERKRzdFXQrx3s44qNK/jmvhGSa6CIiEhXBT3Apz5yKQeOnWPPM02/JrmISC50XdB/4n2buHzjCm7/7ku64pSICF0Y9HFk/MePX87h0xP86f3Pt7scEZG267qgB/ild63jUx++lL/bN8LDLx1rdzkiIm3VlUEP8IfXvodt6wf5g/ue5Hsv/Kzd5YiItE3XBn1/T8x9/+YXePclQ/zbv3mKB1882u6SRETaomuDHmD9in6+8q8+wGXDy/nUvXu585FX2l2SiEjLdXXQA6xe3suez/wy116xgdv/4SX+6wP7GZsqtrssEZGW6fqgB+gtRNz5L97PTR/Ywv9+9FV+5X88wjeeOERJ58QRkbeAt0TQQ3LY5Rd/++f55u//IhtXDvBH33qWj9/xjzz68mi7SxMRaaq3TNBXbN+6hm//wS/xF//8Ks5PFbn5nse5+Z7H+eErxymWyu0uT0RkyVknnBNm+/btvnfv3pa/7mSxxFd/+Dp3PHSAMxNF1i7vZeeVm7hx+2bee8kQZtbymkREGmVm+9x9e+Z6b+Wgrzg3WeTRl0f5zrOHefDFo0yXnOGhPq69YgPvf8dqdly2lo0r+xX8ItJRFPSLdPzcJN9/8Sj/eGCUh18aZXw6OV/OusFeLt+4ItyGeO8lK3jn8CC9hbfc6JeIdAgF/RIolZ39R86w7/WTPPfT07z0szO8fPQcU8VkLL8nNt45PMi2DUNsXj0QbsvYvHqATasG6O+J2/wvEJFu1mjQF1pRTF7FkfFzm1byc5tWzrQVS2VeffM8+4+cYf+Rs+w/coZnDp3iH54/wnTpwp3m2uW9XLKyn0tW9LNhZT8bhvpZv6KPDSv6WD/Uz/qhPlYv76Un1qcCEWkeBf0CFeKId28Y4t0bhth55Wx7qewcPTPBT0+NM3JyjEMnxjlyeoKfnR7n8OkJnjp0ihPnp2o+54r+AmsH+1i7vJc1y3tZO5hMVy/rZeVAD6vCtHJbMVBgoCfWdwYi0hAF/RKJI+NtqwZ426oBPrB1Tc11poplRs9NcuzMBEfPTDJ6doLj56c4cX4qmZ6b4vXjYzz5xklOjk3X/UFXT2wM9hUY7C8w2NfDUJhf3ldgsK/AUH+yM1jWm9wGegthGrOsJ2ZZb4GB3ihp70na+wqRdh4iXUhB30K9hYhNq5Lx+yzlsnN2ssjpsWlOjye3U+NTnBkvcmYiuX9+ssjZieR2bnKa0bOT/OTN85ybLHJuojjzRXKj4sjoL0QM9Mb09yTB3xNH9BUiesN8byGiN47oKUT0xXPaw7Lq6YXLLEzjsCy531fjuXvj5BZF2vGIXKymBL2ZXQf8LyAG/trdv9iM1+lmUWQzQzWLVS4749MlxqZKjE+VGJsuzs5PlRibKs7MJ+sVmZguMz5dYmKqxGSxzFSpzFRx9nZuspjMV7VPV+ZL5XnfU1ysntgu3Amk7Hh65+2UZncqPQWjLzy2EEcUIqMQG4XIiKNoZr6yLI6MnjhZ1hPuJ+tExJX7kRFFNrP+bHtEbEYcnjOy2XVF2mXJg97MYuAvgX8KjABPmNked39xqV9L6osiY3lfMpzTKuWyM12u3gl42AmUmCr6vB1EZWcyXbXzmGmvsWzeDibcHxsrMlVypoql8Bif95h2MiPZAczbSUTEERfsRGbWsWQHE0c289jq5dU7lwt2OjV2NHF4zdhmp3FE1XyybvV6ccRMW1z13Bc8ZqaNeW3Vz1n9XNXPOTNfeY55bdpBLoVmJMAHgYPu/iqAmX0N2Ako6N8Cosjoi2L6Cp11aKm7Uyw7pbIzXSqHaXK/WC5TLM1fXmkvlZ3pslMM7eWq50rWm52v3C/PtJcplaFULif33SmF15p5nlJov+Cx5ZrPP1Uszz421Db72DLlMvMeW6mlHNbL27n84ho7k7k7iqjGTrSyo6z+VBXPaav+tFaIjZ7KNE4+FRbi5FNdIb5weaW9J46S78H6Cgz2Jd99DYbO1eplPR3znVczgn4TcKjq/gjwC014HZGGmVkYBuIt//sG9yTsS1XhX/JkhzA7z8xOqTSzg5jd4ZQvaEt/rpnl855/7msy01Yqz1kenqP69WvXnDxPqWpHN7OzCzvGYrnMRHH2eSs7+GIYdiyWk+l0qTyz/mKHIzeu7GewgU/T/+6abfzG+962qNdoVDOCvtYubN6WMrNbgVsB3v72tzehDBGpxcyILekpSzav+rRV2QFMh09706XyzPdd5yZLjE0WOT9V4uT5KZ4ZOUW5gR+kXsz3cI1qRtCPAFuq7m8GDs9dyd3vAu6C5JexTahDROSiWRjuKeT402AzfpL5BLDNzC41s17gJmBPE15HREQasOQ9encvmtlngO+RHF55j7u/sNSvIyIijWnKcXfu/gDwQDOeW0REFkZn0xIR6XIKehGRLqegFxHpcgp6EZEup6AXEelyHXEpQTMbBV5f5MPXAW8uYTlLpRPr6sSaoDPrUk2N68S6OrEmWPq63uHuw1krdUTQXwwz29vINRNbrRPr6sSaoDPrUk2N68S6OrEmaF9dGroREelyCnoRkS7XDUF/V7sLSNGJdXViTdCZdammxnViXZ1YE7SprtyP0YuISH3d0KMXEZE6ch30Znadmf3YzA6a2W0tfN0tZvawme03sxfM7LOh/c/M7Kdm9nS4XV/1mD8Jdf7YzH6tibW9ZmbPhdffG9rWmNmDZnYgTFeHdjOzO0Jdz5rZ1U2o5z1V2+NpMztjZp9rx7Yys3vM7JiZPV/VtuBtY2a7wvoHzGxXE2r672b2Unjdb5vZqtC+1czGq7bZX1U95v3hfT8Y6l70VUVSalrw+7XU/z9T6vp6VU2vmdnTob1V2yotC9r6dzWPu+fyRnIK5FeAy4Be4Bngiha99kbg6jA/BLwMXAH8GfDva6x/RaivD7g01B03qbbXgHVz2v4bcFuYvw24PcxfD3yX5KpgO4DHWvCe/Qx4Rzu2FfAR4Grg+cVuG2AN8GqYrg7zq5e4pmuBQpi/vaqmrdXrzXmex4FfDPV+F/jYEte0oPerGf8/a9U1Z/n/BP60xdsqLQva+nc195bnHv3MRcjdfQqoXIS86dz9iLs/GebPAvtJrpWbZifwNXefdPefAAdJ6m+VncDuML8buKGq/V5P/AhYZWYbm1jHNcAr7l7vx3FN21bu/ihwosbrLWTb/BrwoLufcPeTwIPAdUtZk7t/392L4e6PSK7SlirUtcLdf+hJatxb9e9YkprqSHu/lvz/Z726Qq/8nwF/W+85mrCt0rKgrX9Xc+U56GtdhLxe2DaFmW0FrgIeC02fCR/J7ql8XKO1tTrwfTPbZ8l1eQE2uPsRSP4wgfVtqAuSq41V/0ds97aChW+bVtf3r0l6gBWXmtlTZvZ/zezDVbWOtKCmhbxfrd5OHwaOuvuBqraWbqs5WdBRf1d5DvqGLkLe1ALMBoFvAZ9z9zPAncA7gSuBIyQfJaG1tX7I3a8GPgZ82sw+UmfdltVlyWUlPwH8XWjqhG1VT1odrdxmXwCKwH2h6Qjwdne/Cvg88DdmtqJFNS30/Wr1+/i7XNiJaOm2qpEFqaumvH5Tt1eeg76hi5A3i5n1kLyx97n73wO4+1F3L7l7Gfg/zA45tKxWdz8cpseAb4cajlaGZML0WKvrItnxPOnuR0N9bd9WwUK3TUvqC1/G/Trwe2GIgTA8cjzM7yMZA393qKl6eGfJa1rE+9Wy99HMCsBvAV+vqrdl26pWFtBhf1d5Dvq2XYQ8jAfeDex39y9VtVePb/8mUDk6YA9wk5n1mdmlwDaSL4SWuq7lZjZUmSf5Uu/58PqVb/F3AfdX1XVzOBJgB3C68nGzCS7ocbV7W1VZ6Lb5HnCtma0OwxfXhrYlY2bXAX8MfMLdx6rah80sDvOXkWybV0NdZ81sR/jbvLnq37FUNS30/Wrl/89fBV5y95khmVZtq7QsoNP+rpbqW9123Ei+wX6ZZG/9hRa+7i+TfKx6Fng63K4Hvgo8F9r3ABurHvOFUOePuYhv+TPquozk6IZngBcq2wRYCzwEHAjTNaHdgL8MdT0HbG9SXcuA48DKqraWbyuSHc0RYJqkB3XLYrYNybj5wXD7ZBNqOkgyXlv52/qrsO5vh/f1GeBJ4Deqnmc7Sfi+AvwF4ceQS1jTgt+vpf7/Wauu0P4V4PfnrNuqbZWWBW39u5p70y9jRUS6XJ6HbkREpAEKehGRLqegFxHpcgp6EZEup6AXEelyCnoRkS6noBcR6XIKehGRLvf/AV2RYlNTShJWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1781c7aa58>]"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHwdJREFUeJzt3X2UXHWd5/H3tx76Mf1QlXRC6CRdAcJDQEnoSojguo4ooM4aR8dZOI4yLiO6E1yc45kd9MwenXE5x52jMjoiKwiKjoIcxCXjcmRA8cwiAqkQCAkBE0In6SSQDt2dp04/f/ePut2pJJ10d9KdW1X38zqnT9361a+qvlUn6U/f3+/e3zV3R0REoicWdgEiIhIOBYCISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKIUACIiEaUAEBGJqETYBZzMrFmzPJPJhF2GiEhJWbt27V53bxqvX1EHQCaTIZfLhV2GiEhJMbNtE+mnISARkYhSAIiIRJQCQEQkohQAIiIRpQAQEYkoBYCISEQpAEREIqosA6C7p59v/3ozG3buC7sUEZGiVdQngp2qeMz4pyf+wLA7lzQ3hF2OiEhRGncPwMyqzOw5M3vRzDaa2d8H7QvN7Fkz22xmPzOziqC9Mri/JXg8U/BaXwzaXzWza6brQ9VVJbngrHpybV3T9RYiIiVvIkNAfcB73P1SYAlwrZmtAP4XcLu7LwK6gBuD/jcCXe5+HnB70A8zWwxcB1wMXAt818ziU/lhCmVbUqzb3sXg0PB0vYWISEkbNwA872BwNxn8OPAe4KGg/T7gw8H2yuA+weNXmZkF7Q+4e5+7vw5sAZZPyacYQzaT4lD/EK+8cWC63kJEpKRNaBLYzOJm9gKwB3gceA3odvfBoEs70BxsNwM7AILH9wEzC9vHeE7he91kZjkzy3V0dEz+EwWymTQAubbOU34NEZFyNqEAcPchd18CzCP/V/tFY3ULbu0Ej52o/dj3usvds+6ebWoadzXTE2purObshipy2zQPICIylkkdBuru3cBvgRVAo5mNHEU0D9gVbLcD8wGCxxuAzsL2MZ4zLVozaXJtXbgflzMiIpE3kaOAmsysMdiuBt4LbAKeBP406HYD8EiwvTq4T/D4bzz/G3g1cF1wlNBCYBHw3FR9kLFkW1K8sb+Xnd2Hp/NtRERK0kTOA5gL3BccsRMDHnT3X5rZy8ADZvY/gXXAPUH/e4Afm9kW8n/5Xwfg7hvN7EHgZWAQWOXuQ1P7cY6WzaQAWLuti3mpmul8KxGRkjNuALj7emDpGO1bGeMoHnfvBT52gte6Dbht8mWemgvPqmdGZYI1bZ2sXHLcfLOISKSV5VIQI+IxY+mCRp0QJiIyhrIOAIBsS5pX3zzAvsMDYZciIlJUyj8AMincYd127QWIiBQq+wBYMr+ReMxYq/MBRESOUvYBUFuZYPHcetbojGARkaOUfQAAtLakeGFHNwNaGE5EZFQkAmBZJk3vwDAbd+0PuxQRkaIRiQAYOSFMC8OJiBwRiQCYU1/FvFS1JoJFRApEIgAgPwy0RgvDiYiMikwAtLak2Huwj+2dPWGXIiJSFCITAMtGLxCjYSAREYhQACyaPYP6qgS5bZoIFhGBCAVALGZc1pLSHoCISCAyAQD5YaDNew7S3dMfdikiIqGLVAC0thy5QIyISNRFKgAunddIIma6ULyICBELgOqKOJc0N+iMYBERIhYAkL9Q/Ivt++gbnNbLEYuIFL3oBUAmTf/gMBt2amE4EYm2yAXAyESwhoFEJOoiFwBNdZVkZtZoIlhEIi9yAQD5YaC127QwnIhEWzQDoCVF56F+tu49FHYpIiKhiWYABBeIWatlIUQkwiIZAOc2zSBVk9SF4kUk0sYNADObb2ZPmtkmM9toZrcE7V8xs51m9kLw84GC53zRzLaY2atmdk1B+7VB2xYzu3V6PtL4zIzWlpSWhBCRSEtMoM8g8AV3f97M6oC1ZvZ48Njt7v71ws5mthi4DrgYOBt4wszODx6+A3gf0A6sMbPV7v7yVHyQycpm0jyxaQ9vHexj5ozKMEoQEQnVuHsA7r7b3Z8Ptg8Am4DmkzxlJfCAu/e5++vAFmB58LPF3be6ez/wQNA3FNmR8wG0FyAiETWpOQAzywBLgWeDppvNbL2Z3WtmqaCtGdhR8LT2oO1E7aG4pLmBinhMw0AiElkTDgAzmwH8HPi8u+8H7gTOBZYAu4FvjHQd4+l+kvZj3+cmM8uZWa6jo2Oi5U1aVTLO2+c1aCJYRCJrQgFgZknyv/x/4u4PA7j7m+4+5O7DwN3kh3gg/5f9/IKnzwN2naT9KO5+l7tn3T3b1NQ02c8zKa2ZFBt27qN3QAvDiUj0TOQoIAPuATa5+zcL2ucWdPsTYEOwvRq4zswqzWwhsAh4DlgDLDKzhWZWQX6iePXUfIxTk21JMzDkrG/fF2YZIiKhmMhRQFcCnwBeMrMXgrYvAdeb2RLywzhtwGcA3H2jmT0IvEz+CKJV7j4EYGY3A48BceBed984hZ9l0kYWhlvT1snyhekwSxEROePGDQB3f4qxx+8fPclzbgNuG6P90ZM970xL11ZwblOtJoJFJJIieSZwoWWZNLm2ToaHtTCciERL5AOgtSXF/t5BtnQcDLsUEZEzKvIBkM3kx/5zWhhORCIm8gGQmVnDrBkVukKYiERO5ANgZGE4LQkhIlET+QCA/PkA2zt72LO/N+xSRETOGAUARy4Qo70AEYkSBQBw8dkNVCZimggWkUhRAAAViRhL5jeS26aJYBGJDgVAIJtJsXHXfnr6B8MuRUTkjFAABLItaYaGnRd2dIddiojIGaEACFy2IIWZTggTkehQAAQaapKcP7tORwKJSGQoAAq0ZlKs29bFkBaGE5EIUAAUWJZJcaBvkFffOBB2KSIi004BUCDbkl8Ybq0OBxWRCFAAFJiXqmZOfSVrNBEsIhGgAChgZmRb0rpCmIhEggLgGK0tKXZ2H2ZX9+GwSxERmVYKgGMsG7lAjPYCRKTMKQCOcdHcOmoq4qzVBWJEpMwpAI6RiI8sDKc9ABEpbwqAMWQzaTbt3s/BPi0MJyLlSwEwhmxLimGHddu1FyAi5UsBMIalCxqJGTofQETKmgJgDHVVSS48q15nBItIWVMAnEA2k2Ld9m4Gh4bDLkVEZFqMGwBmNt/MnjSzTWa20cxuCdrTZva4mW0OblNBu5nZt81si5mtN7PLCl7rhqD/ZjO7Yfo+1unLZtL09A+xabcWhhOR8jSRPYBB4AvufhGwAlhlZouBW4Ffu/si4NfBfYD3A4uCn5uAOyEfGMCXgcuB5cCXR0KjGGVb8qXpOsEiUq7GDQB33+3uzwfbB4BNQDOwErgv6HYf8OFgeyXwI897Bmg0s7nANcDj7t7p7l3A48C1U/ppptDZjdWc3VCl8wFEpGxNag7AzDLAUuBZYI6774Z8SACzg27NwI6Cp7UHbSdqP/Y9bjKznJnlOjo6JlPelMtm0uTaOnHXBWJEpPxMOADMbAbwc+Dz7r7/ZF3HaPOTtB/d4H6Xu2fdPdvU1DTR8qZFNpPizf19tHdpYTgRKT8TCgAzS5L/5f8Td384aH4zGNohuN0TtLcD8wuePg/YdZL2ojVygRjNA4hIOZrIUUAG3ANscvdvFjy0Ghg5kucG4JGC9k8GRwOtAPYFQ0SPAVebWSqY/L06aCtaF5xVR11lgpxOCBORMpSYQJ8rgU8AL5nZC0Hbl4CvAQ+a2Y3AduBjwWOPAh8AtgA9wKcA3L3TzL4KrAn6/YO7F/Wf1vGYsWRBoy4QIyJladwAcPenGHv8HuCqMfo7sOoEr3UvcO9kCgzbskya25/4A/sOD9BQnQy7HBGRKaMzgceRbUnhDs9rYTgRKTMKgHEsWdBIPGbkdIEYESkzCoBx1FQkuPjsek0Ei0jZUQBMQGtLihfbu+kf1MJwIlI+FAATsCyTpndgmI279oVdiojIlFEATMDIwnA6HFREyokCYAJm11cxP12teQARKSsKgAla1pImt00Lw4lI+VAATFBrJsXeg/1se6sn7FJERKaEAmCCRhaGW6PzAUSkTCgAJmjR7BnUVyU0ESwiZUMBMEGxmNHaktIVwkSkbCgAJiGbSbNlz0G6DvWHXYqIyGlTAEyCzgcQkXKiAJiES+c3koybhoFEpCwoACahKhnnkuYGrQwqImVBATBJ2ZYU69v30TswFHYpIiKnRQEwSa0tafqHhtmwUwvDiUhpUwBMUjaTnwjWPICIlDoFwCTNmlHJwlm1WhhOREqeAuAUtLakWKuF4USkxCkATsGyTIqungFe6zgUdikiIqdMAXAKWoOF4dZu0+GgIlK6FACn4NymWlI1SdZoHkBESpgC4BSYGa0taS0JISIlTQFwirKZFK/vPUTHgb6wSxEROSXjBoCZ3Wtme8xsQ0HbV8xsp5m9EPx8oOCxL5rZFjN71cyuKWi/NmjbYma3Tv1HObOWZbQwnIiUtonsAfwQuHaM9tvdfUnw8yiAmS0GrgMuDp7zXTOLm1kcuAN4P7AYuD7oW7IuaW6gIhHTRLCIlKzEeB3c/d/NLDPB11sJPODufcDrZrYFWB48tsXdtwKY2QNB35cnXXGRqEzEeXtzgyaCRaRknc4cwM1mtj4YIkoFbc3AjoI+7UHbidpLWjaTZuOufRzu18JwIlJ6TjUA7gTOBZYAu4FvBO02Rl8/SftxzOwmM8uZWa6jo+MUyzszsi0pBoacF9u7wy5FRGTSTikA3P1Ndx9y92Hgbo4M87QD8wu6zgN2naR9rNe+y92z7p5tamo6lfLOmFZdIUxEStgpBYCZzS24+yfAyBFCq4HrzKzSzBYCi4DngDXAIjNbaGYV5CeKV5962cUhVVvBebNn6AIxIlKSxp0ENrP7gXcDs8ysHfgy8G4zW0J+GKcN+AyAu280swfJT+4OAqvcfSh4nZuBx4A4cK+7b5zyTxOCbEuKR1/azfCwE4uNNdIlIlKcJnIU0PVjNN9zkv63AbeN0f4o8OikqisB2UyaB9bsYPOeg1xwVl3Y5YiITJjOBD5N2ZaRC8RoGEhESosC4DS1zKxh1owKXSBGREqOAuA0mRnZlrT2AESk5CgApkA2k2JH52He3N8bdikiIhOmAJgC2Uz+AjEaBhKRUqIAmAIXn11PVTKmYSARKSkKgCmQjMe4dF6j9gBEpKQoAKbIskyal3fv51DfYNiliIhMiAJgirRmUgwNOy/u0MJwIlIaFABT5LIFKczQ9QFEpGQoAKZIQ3WSC+bUaSJYREqGAmAKtbakWLe9m6HhMS91ICJSVBQAU2hZJs3BvkFeeWN/2KWIiIxLATCFdIEYESklCoApNC9VzZz6Sk0Ei0hJUABMITMjm0mzVlcIE5ESoACYYtmWFLv29bKz+3DYpYiInJQCYIplW0YWhtNegIgUNwXAFLtobh01FXFNBItI0VMATLFEPMbSBY2aCBaRoqcAmAbZljSvvrGfA70DYZciInJCCoBpkM2kGHZYt10Lw4lI8VIATIOlC1LETBPBIlLcFADTYEZlgovm1pPTRLCIFDEFwDTJtqR4YUc3A0PDYZciIjImBcA0ac2k6ekfYtNuLQwnIsVJATBNlmXyC8PpOsEiUqzGDQAzu9fM9pjZhoK2tJk9bmabg9tU0G5m9m0z22Jm683ssoLn3BD032xmN0zPxykecxuqaW6s1gViRKRoTWQP4IfAtce03Qr82t0XAb8O7gO8H1gU/NwE3An5wAC+DFwOLAe+PBIa5SybSZFr68JdF4gRkeIzbgC4+78Dx/4ZuxK4L9i+D/hwQfuPPO8ZoNHM5gLXAI+7e6e7dwGPc3yolJ1sS4o9B/p4aee+sEsRETnOqc4BzHH33QDB7eygvRnYUdCvPWg7UftxzOwmM8uZWa6jo+MUyysO7108h1kzKvj43c/ym1feDLscEZGjTPUksI3R5idpP77R/S53z7p7tqmpaUqLO9PmNlTzyM3vZMHMGm68L8d3f7tFw0EiUjRONQDeDIZ2CG73BO3twPyCfvOAXSdpL3vNjdU89Nkr+ODb5vKPv3qVWx54gcP9Q2GXJSJyygGwGhg5kucG4JGC9k8GRwOtAPYFQ0SPAVebWSqY/L06aIuE6oo4/3z9Uv7mmgv41/W7+Nj3nmaXLhgjIiGbyGGg9wO/By4ws3YzuxH4GvA+M9sMvC+4D/AosBXYAtwN/BWAu3cCXwXWBD//ELRFhpmx6o/O4+5PZGnb28OHvvM7rRUkIqGyYh6Tzmaznsvlwi5jym1+8wCf/lGOnd2H+erKS7hu+YKwSxKRMmJma909O14/nQkcgkVz6nhk1TtZcc5Mbn34Jb6yeqPWDBKRM04BEJKGmiQ/+Itl/OU7F/LDp9u44d7n6DrUH3ZZIhIhCoAQJeIx/u6PF/ONj11Krq2LD93xFK+8ocXjROTMUAAUgY+2zuNnn1lB38AwH/nu0/xqwxthlyQiEaAAKBJLF6T418+9k0Vz6vjsv6zlW09sZni4eCfoRaT0KQCKyJz6Kn520wo+srSZ25/4A6t++jyH+gbDLktEypQCoMhUJeN8488u5e8+eBGPbXyDj975NDs6e8IuS0TKkAKgCJkZf/kfzuEHn1rOzu7DfOg7T/H7194KuywRKTMKgCL2H89v4pFVV5KureAT9zzLj5/ZFnZJIlJGFABF7pymGfxi1ZW86/wm/sf/2cCXfvES/YM6aUxETp8CoATUVyW5+5NZ/uu7z+Wnz27nz7//LHsP9oVdloiUOAVAiYjHjL+99kK+dd0SXmzvZuV3fscGXWlMRE6DAqDErFzSzEOfvYJhd/70fz/NL9dH4rIKIjINFAAl6G3zGlh98zu55OwGbv7pOr7+2Ks6aUxEJk0BUKKa6ir5yacv5z9n5/OdJ7dw049zHOgdCLssESkhCoASVpmI87WPvo2//9DFPPlqBx/57tO07T0UdlkiUiIUACXOzLjhigw//i/L6TjYx8o7fsdTm/eGXZaIlAAFQJm44rxZrF71Ts6qr+KT9z7LPU+9TjFf7U1EwqcAKCMLZtbw8F9dwfsWz+Grv3yZv3loPX2DQ2GXJSJFSgFQZmorE9z58VZuuWoRD61t57q7nmHP/t6wyxKRIqQAKEOxmPHX7zufOz9+Ga/sPsB/+s5TvLCjO+yyRKTIJMIuQKbP+982l8ysWj79oxwfvuN3nNtUy4pzZvKOc2dy+cKZNNVVhl2iiITIinmiMJvNei6XC7uMktd1qJ8Hczt4ZutbrGnr4mBwkZnzZs9gxTlp3nHOLC4/J82sGQoEkXJgZmvdPTtuPwVAtAwODbNh136e2fpWPhBe7+RQf36ieNHsGQV7CGlmKhBESpICQCZkcGiYl3bu45mtncEeQic9QSCcPycfCCvOUSCIlBIFgJySgdFAeItntnaSKwiEC+bUseKcdD4QzplJurYi5GpFZCxnJADMrA04AAwBg+6eNbM08DMgA7QBf+buXWZmwLeADwA9wF+4+/Mne30FQPhGAuH3r+WHjHJtXRweyAfChWfVBXsIaS5fOJOUAkGkKJzJAMi6+96Ctn8EOt39a2Z2K5By9781sw8AnyMfAJcD33L3y0/2+gqA4jMwNMz69n2jcwhjB0J+yEiBIBKOMAPgVeDd7r7bzOYCv3X3C8zse8H2/cf2O9HrKwCKX//gMC/t7OaZrZ38/rW3yG3rpHdgGDO48Kz6I0NGC9M01igQRM6EiQbA6Z4H4MC/mZkD33P3u4A5I7/UgxCYHfRtBnYUPLc9aDthAEjxq0jEaG1J09qSZtUfnUf/4DDr27t5Zutb/H7rW9z/3HZ+8Ls2zOCis+pZfHY9tRVxqisS1FTEqamIUz1ymzy2LXFkOxknEdd5iyJT6XQD4Ep33xX8kn/czF45SV8bo+243Q8zuwm4CWDBggWnWZ6caRWJGNlMmmwmzc3vWUTf4FB+yOi1fCA8tXkvPf2DHB4YYmBocnufFfHYkbAIbmuSiePbKhJUJ+MFYXJ0kNRUJJifrtYeiUTeaQWAu+8KbveY2S+A5cCbZja3YAhoT9C9HZhf8PR5wHHXMwz2Iu6C/BDQ6dQn4atMxFmWSbMsk+ZzVy066rGBoWEODwxxuH+Inv6hfDCMbg9xeGAwf1vY1p9v6xl93iDdPf3s6h55Tr6td2B43NrOmz2DZZkU2ZZ8ffPT1eSPVRCJhlMOADOrBWLufiDYvhr4B2A1cAPwteD2keApq4GbzewB8pPA+042/i/lLxmPkYzHqK9KTvlrDw97EAZBUBSEyaG+QTbvOUiurZP/u3439z+XH5lsqqscDYRsJsXiufUadpKydjp7AHOAXwR/MSWAn7r7r8xsDfCgmd0IbAc+FvR/lPwRQFvIHwb6qdN4b5GTisWM2soEtZVj/xO/+uL87fCws3nPQda05c95yG3r4tGX3gCgpiLO0gWNtLakWZZJsXRBihkneD2RUqQTwUSOsXvfYXJtXaOBsGn3foYdYgaLz64f3UPItqQ5q6Eq7HJFjqMzgUWmyIHeAdZt7ya3LR8K67Z3j577MC9VzbLMkUBYNHsGsZjmESRcZ+owUJGyV1eV5F3nN/Gu85uA/OT1pt37WRPsJfy/zXv5xbqdANRXJYKjoPKB8PZ5DVQl42GWL3JC2gMQOU3uzvbOntFAyG3rYsueg0D+0NW3zWsg25Iim0nT2pLSGkoy7TQEJBKizkP9rN3WRW5bJ7m2Lta3d4+e93BuUy3LgjBobqwenayuq8rf1iTjGkaS06IAECkivQP5E+JGAiHX1sn+3sET9q+tiFNbmWBGEA61lfGC7aC9It8+EhyF7TOC59RWJqhMxHR+Q8RoDkCkiFQl4yxfmGb5wjSQP/x0696DdBzo51DfIIf6BznYN8ihvkEO9g1xsDfY7s/fHuobZGd37+j2wb5B+gbHP9kNIBEcElsYCqNBUZXfrkzE8j/J+Oh2RSJGZSK4n8xvV4z0K2iviB95XiJmCpsSogAQCUEsZpw3u47zZo/f90QGhobp6RsaDYnRAOk9sn2of6ggWEbCZIgDvYO8sa/3qDCZaKCc9HMZRwdF8khYHB8e8SA88u3VyTiNNUlSNRWkaytI1Vbkt2sqqKtKaFhsGigAREpUMh6joSZGQ83UnEnt7vQPDdMfhEHf4DB9A0Oj2/n2IfoGRh4fKuh7pL1/6OjnHdVvYJjunv6C1zvy3J6BIYaGxx6SjseMVE2SxiAQUrVJ0rUVBfcrSNcmj7pfX5XQ3sg4FAAiAoCZBX+dx6kL4f3dnYN9g3QdGqCzp5+uQ/109fTTOXo7QHdw//W9h1i7rZvunn4GTxAaiZiN7lGkao8Ex+geRnDbWJMc3eOoq4xWaCgARKQomBl1VUnqqpIsmFkzoee4Owf6Buk6lA+G7p6BgsDop6tnIP9YTz+vdRyka9sAXT39J9zTSMSM+uokVcfMhxzZjlOZjFEV3I60VSWPnhc58ni+rSp59FzKUf0TsdDWnFIAiEjJMjPqq5LUVyVpmVk7oee4O/t7B0f3MEb2Lkbu7zs8cNwQWO9Aft5k72B/wXDX0cNepyMRs6OCpioZ55LmBv75+qWn9brjvu+0vrqISJExMxqqkzRUJ8kwsdAYz/CwH5nPCOY6eguDIgiRwjmP3mPmSXqPCZV5qeopqe1kFAAiIqcpFjOqg4sOlRItdi4iElEKABGRiFIAiIhElAJARCSiFAAiIhGlABARiSgFgIhIRCkAREQiqqgvCGNmHcC203iJWcDeKSqn1Om7OJq+j6Pp+ziiHL6LFndvGq9TUQfA6TKz3ESuihMF+i6Opu/jaPo+jojSd6EhIBGRiFIAiIhEVLkHwF1hF1BE9F0cTd/H0fR9HBGZ76Ks5wBEROTEyn0PQERETqAsA8DMrjWzV81si5ndGnY9YTKz+Wb2pJltMrONZnZL2DWFzcziZrbOzH4Zdi1hM7NGM3vIzF4J/o28I+yawmRmfx38P9lgZvebWVXYNU2nsgsAM4sDdwDvBxYD15vZ4nCrCtUg8AV3vwhYAayK+PcBcAuwKewiisS3gF+5+4XApUT4ezGzZuC/AVl3vwSIA9eFW9X0KrsAAJYDW9x9q7v3Aw8AK0OuKTTuvtvdnw+2D5D/D94cblXhMbN5wAeB74ddS9jMrB54F3APgLv3u3t3uFWFLgFUm1kCqAF2hVzPtCrHAGgGdhTcbyfCv/AKmVkGWAo8G24lofon4L8Dp3cV7/JwDtAB/CAYEvu+mU3NRXJLkLvvBL4ObAd2A/vc/d/CrWp6lWMA2BhtkT/UycxmAD8HPu/u+8OuJwxm9sfAHndfG3YtRSIBXAbc6e5LgUNAZOfMzCxFfrRgIXA2UGtmfx5uVdOrHAOgHZhfcH8eZb4bNx4zS5L/5f8Td3847HpCdCXwITNrIz80+B4z+5dwSwpVO9Du7iN7hA+RD4Soei/wurt3uPsA8DBwRcg1TatyDIA1wCIzW2hmFeQncVaHXFNozMzIj/Fucvdvhl1PmNz9i+4+z90z5P9d/Mbdy/ovvJNx9zeAHWZ2QdB0FfByiCWFbTuwwsxqgv83V1Hmk+KJsAuYau4+aGY3A4+Rn8W/1903hlxWmK4EPgG8ZGYvBG1fcvdHQ6xJisfngJ8EfyxtBT4Vcj2hcfdnzewh4HnyR8+to8zPCtaZwCIiEVWOQ0AiIjIBCgARkYhSAIiIRJQCQEQkohQAIiIRpQAQEYkoBYCISEQpAEREIur/A8bAc5ovxt0fAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "hidden": true
   },
   "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": 115,
   "metadata": {
    "hidden": true
   },
   "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": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(Vh[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We get topics that match the kinds of clusters we would expect! This is despite the fact that this is an **unsupervised algorithm** - which is to say, we never actually told the algorithm how our documents are grouped."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We will return to SVD in **much more detail** later.  For now, the important takeaway is that we have a tool that allows us to exactly factor a matrix into orthogonal columns and orthogonal rows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-negative Matrix Factorization (NMF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Motivation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/face_pca.png\" alt=\"PCA on faces\" style=\"width: 80%\"/>\n",
    "\n",
    "(source: [NMF Tutorial](http://perso.telecom-paristech.fr/~essid/teach/NMF_tutorial_ICME-2014.pdf))\n",
    "\n",
    "A more interpretable approach:\n",
    "\n",
    "<img src=\"images/face_outputs.png\" alt=\"NMF on Faces\" style=\"width: 80%\"/>\n",
    "\n",
    "(source: [NMF Tutorial](http://perso.telecom-paristech.fr/~essid/teach/NMF_tutorial_ICME-2014.pdf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Idea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than constraining our factors to be *orthogonal*, another idea would to constrain them to be *non-negative*. NMF is a factorization of a non-negative data set $V$: $$ V = W H$$ into non-negative matrices $W,\\; H$. Often positive factors will be **more easily interpretable** (and this is the reason behind NMF's popularity). \n",
    "\n",
    "<img src=\"images/face_nmf.png\" alt=\"NMF on faces\" style=\"width: 80%\"/>\n",
    "\n",
    "(source: [NMF Tutorial](http://perso.telecom-paristech.fr/~essid/teach/NMF_tutorial_ICME-2014.pdf))\n",
    "\n",
    "Nonnegative matrix factorization (NMF) is a non-exact factorization that factors into one skinny positive matrix and one short positive matrix.  NMF is NP-hard and non-unique.  There are a number of variations on it, created by adding different constraints. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Applications of NMF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Face Decompositions](http://scikit-learn.org/stable/auto_examples/decomposition/plot_faces_decomposition.html#sphx-glr-auto-examples-decomposition-plot-faces-decomposition-py)\n",
    "- [Collaborative Filtering, eg movie recommendations](http://www.quuxlabs.com/blog/2010/09/matrix-factorization-a-simple-tutorial-and-implementation-in-python/)\n",
    "- [Audio source separation](https://pdfs.semanticscholar.org/cc88/0b24791349df39c5d9b8c352911a0417df34.pdf)\n",
    "- [Chemistry](http://ieeexplore.ieee.org/document/1532909/)\n",
    "- [Bioinformatics](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0485-4) and [Gene Expression](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2623306/)\n",
    "- Topic Modeling (our problem!)\n",
    "\n",
    "<img src=\"images/nmf_doc.png\" alt=\"NMF on documents\" style=\"width: 80%\"/>\n",
    "\n",
    "(source: [NMF Tutorial](http://perso.telecom-paristech.fr/~essid/teach/NMF_tutorial_ICME-2014.pdf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**More Reading**:\n",
    "\n",
    "- [The Why and How of Nonnegative Matrix Factorization](https://arxiv.org/pdf/1401.5226.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NMF from sklearn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use [scikit-learn's implementation of NMF](http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.NMF.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "m,n=vectors.shape\n",
    "d=5  # num topics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "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": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['jpeg image gif file color images format quality',\n",
       " 'edu graphics pub mail 128 ray ftp send',\n",
       " 'space launch satellite nasa commercial satellites year market',\n",
       " 'jesus god people matthew atheists does atheism said',\n",
       " 'image data available software processing ftp edu analysis']"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 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": 119,
   "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": 160,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[\"a\\n\\nWhat about positional uncertainties in S-L 1993e?   I assume we know where\\nand what Galileo is doing within a few meters.   But without the\\nHGA,  don't we have to have some pretty good ideas, of where to look\\nbefore imaging?  If the HGA was working,  they could slew around\\nin near real time (Less speed of light delay).  But when they were\\nimaging toutatis????  didn't someone have to get lucky on a guess to\\nfind the first images?   \\n\\nAlso, I imagine S-L 1993e will be mostly a visual image.  so how will\\nthat affect the other imaging missions.  with the LGA,  there is a real\\ntight allocation of bandwidth.   It may be premature to hope for answers,\\nbut I thought i'd throw it on the floor.\",\n",
       " \"I would like to program Tseng ET4000 to nonstandard 1024x768 mode by\\nswitching to standard 1024x768 mode using BIOS and than changing some\\ntiming details (0x3D4 registers 0x00-0x1F) but I don't know how to\\nselect 36 MHz pixel clock I need. The BIOS function selects 40 MHz.\\n\\nIs there anybody who knows where to obtain technical info about this.\\nI am also interested in any other technical information about Tseng ET4000\\nand Trident 8900 and 9000 chipsets.\\n\\n\\t\\t\\tthanks very much\",\n",
       " 'In-Reply-To: <20APR199312262902@rigel.tamu.edu> lmp8913@rigel.tamu.edu (PRESTON, LISA M)',\n",
       " \"\\n\\n\\n\\nI'm not sure, but it almost sounds like they can't figure out where the \\n_nucleus_ is within the coma. If they're off by a couple hundred\\nmiles, well, you can imagine the rest...\\n\",\n",
       " \"Hello,\\n     I am looking to add voice input capability to a user interface I am\\ndeveloping on an HP730 (UNIX) workstation.  I would greatly appreciate \\ninformation anyone would care to offer about voice input systems that are \\neasily accessible from the UNIX environment. \\n\\n     The names or adresses of applicable vendors, as well as any \\nexperiences you have had with specific systems, would be very helpful.\\n\\n     Please respond via email; I will post a summary if there is \\nsufficient interest.\\n\\n\\nThanks,\\nKen\\n\\n\\nP.S.  I have found several impressive systems for IBM PC's, but I would \\nlike to avoid the hassle of purchasing and maintaining a separate PC if \\nat all possible.\\n\\n-------------------------------------------------------------------------------\\nKen Hinckley (kph2q@virginia.edu)\\nUniversity of Virginia \\nNeurosurgical Visualization Laboratory\",\n",
       " '\\nIt was a test of the first reusable tool.\\n\\n\\nPointy so they can find them or so they will stick into their pants better, and\\nbe closer to their brains?',\n",
       " '\\nSize of armies, duration, numbers of casualties both absolute and as a\\npercentage of those involved, geographical area and numbers of countries\\ntoo, are all measures of size.  In this case I\\'d say the relevant\\nstatistic would be the number of combatants (total troops) compared to\\ntotal casualties from among the total civilian population in the\\naffected geographical area.\\n\\n\\nVietnam and Korea might make good comparisons.\\n\\n\\nWestern news in general, but in particular the American \"mass media\":\\nCBS, NBC, ABC, etc.  The general tone of the news during the whole\\nwar was one of \"those poor, poor Iraqis\" along with \"look how precisely\\nthis cruise missile blew this building to bits\".\\n\\n\\nI agree.\\n\\n\\nPerhaps so.  And maybe the atomic bomb was a mistake too.  But that\\'s easy\\nto say from our \"enlightened\" viewpoint here in the 90\\'s, right?  Back\\nthen, it was *all-out* war, and Germany and Japan had to be squashed.\\nAfter all, a million or more British had already died, hundreds of \\nthousands of French, a couple hundread thousand or so Americans, and \\nmillions of Russians, not to mention a few million Jews, Poles, and \\nother people of slavic descent in German concentration camps.  All \\nthings considered, the fire-bombings and the atomic bomb were\\nessential (and therefore justified) in bringing the war to a quick\\nend to avoid even greater allied losses.\\n\\nI, for one, don\\'t regret it.\\n\\n\\nSure.  And it\\'s the people who suffer because of them.  All the more\\nreason to depose these \"entrenched political rulers operating in their\\nown selfish interests\"!  Or do you mean that this applies to the allies\\nas well??\\n\\n\\nI make no claim or effort to justify the misguided foreign policy of the\\nWest before the war.  It is evident that the West, especially America,\\nmisjudged Hussein drastically.  But once Hussein invaded Kuwait and \\nthreatened to militarily corner a significant portion of the world\\'s\\noil supply, he had to be stopped.  Sure the war could have been\\nprevented by judicious and concerted effort on the part of the West\\nbefore Hussein invaded Kuwait, but it is still *Hussein* who is\\nresponsible for his decision to invade.  And once he did so, a\\nstrong response from the West was required.\\n\\n\\nWell, it\\'s not very \"loving\" to allow a Hussein or a Hitler to gobble up\\nnearby countries and keep them.  Or to allow them to continue with mass\\nslaughter of certain peoples under their dominion.  So, I\\'d have to\\nsay yes, stopping Hussein was the most \"loving\" thing to do for the\\nmost people involved once he set his mind on military conquest.\\n\\nI mentioned it.\\n\\nIf we hadn\\'t intervened, allowing Hussein to keep Kuwait, then it would\\nhave been appeasement.  It is precisely the lessons the world learned\\nin WW2 that motivated the Western alliance to war.  Letting Hitler take\\nAustria and Czechoslavkia did not stop WW2 from happening, and letting\\nHussein keep Kuwait would not have stopped an eventual Gulf War to\\nprotect Saudi Arabia.\\n\\n\\nSure.  What was truly unfortunate was that they followed Hitler in\\nhis grandiose quest for a \"Thousand Year Reich\".  The consequences\\nstemmed from that.\\n\\nWhat should I say about them?  Anything in particular?\\n\\n\\n\\nSo?  It was the *policemen* on trial not Rodney King!!  And under American\\nlaw they deserved a jury of *their* peers!  If there had been black\\nofficers involved, I\\'m sure their would have been black jurors too.\\nThis point (of allegedly racial motivations) is really shallow.\\n\\n\\nSo?  It\\'s \"hard to imagine\"?  So when has Argument from Incredulity\\ngained acceptance from the revered author of \"Constructing a Logical\\nArgument\"?  Can we expect another revision soon??  :)  (Just kidding.)\\n\\n\\nI have to admit that I wonder this too.  But *neither* the prosecution\\nnor the defense is talking.  So one cannot conclude either way due to\\nthe silence of the principals.  \\n\\n\\nOK.  It certainly seemed to me that there was excessive force involved.\\nAnd frankly, the original \"not guilty\" verdict baffled me too.  But then\\nI learned that the prosecution in the first case did not try to convict\\non a charge of excessive force or simple assault which they probably\\nwould have won, they tried to get a conviction on a charge of aggravated\\nassault with intent to inflict serious bodily harm.  A charge, which\\nnews commentators said, was akin to attempted murder under California\\nlaw.  Based on what the prosecution was asking for, it\\'s evident that \\nthe first jury decided that the officers were \"not guilty\".  Note, \\nnot \"not guilty\" of doing wrong, but \"not guilty\" of aggravated assault \\nwith the *intent* of inflicting serious bodily harm.  The seeds of the \\nprosecutions defeat were in their own overconfidence in obtaining a \\nverdict such that they went for the most extreme charge they could.\\n\\nIf the facts as the news commentators presented them are true, then\\nI feel the \"not guilty\" verdict was a reasonable one.\\n\\n\\nThanks mathew, I like the quote.  Pretty funny actually.  (I\\'m a \\nMonty Python fan, you know.  Kind of seems in that vein.)\\n\\nOf course, oversimplifying any moral argument can make it seem\\ncontradictory.  But then, you know that already.  \\n\\nRegards,',\n",
       " \"<stuff deleted>\\n\\nYou mean like: seconds, minutes, hours, days, months, years. . .  :-)\\n\\nRemember,  the Fahrenheit temperature scale is also a centigrade scale.  Some\\nrevisionists tell the history something like this:  The coldest point in a\\nparticular Russian winter was marked on the thermometer as was the body\\ntemperature of a volunteer (turns out he was sick, but you can't win 'em all).\\nThen the space in between the marks on the thermometer was then divided into\\nhundredths.\\n\\t\\t\\t\\t\\t\\t\\t\\t:-)\\n\\nFWIW,\\n\\nDoug Page\\n\",\n",
       " \"\\nIt wasn't especially prominent, as I recall.  However, quite possibly it's\\nno longer on display; NASM, like most museums, has much more stuff than it\\ncan display at once, and does rotate the displays occasionally.\",\n",
       " \"DM> Fact or rumor....?  Madalyn Murray O'Hare an atheist who eliminated the\\nDM> use of the bible reading and prayer in public schools 15 years ago is now\\nDM> going to appear before the FCC with a petition to stop the reading of the\\nDM> Gospel on the airways of America.  And she is also campaigning to remove\\nDM> Christmas programs, songs, etc from the public schools.  If it is true\\nDM> then mail to Federal Communications Commission 1919 H Street Washington DC\\nDM> 20054 expressing your opposition to her request.  Reference Petition number\\n\\nDM> 2493.\\n\\nFalse.  This story has been going around for years.  There's not a drop of\\ntruth.  Note that I don't care for O'Hare (O'Hair?) myself, but this\\nis one thing she's not guilty of.\\n\"]"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.data[10:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "W1 = clf.fit_transform(vectors_tfidf)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['people don think just like objective say morality',\n",
       " 'graphics thanks files image file program windows know',\n",
       " 'space nasa launch shuttle orbit moon lunar earth',\n",
       " 'ico bobbe tek beauchaine bronx manhattan sank queens',\n",
       " 'god jesus bible believe christian atheism does belief']"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f15cb8ae0f0>]"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuUFOWdN/DvTwzmrDFvNJDdHG8Yl2yWbBJNJsYkm2zuYi4QcxNzWc2a9U1OOMm+STaLGolBE42JJlExKwrGaBQxRkBAERCVi8AMtwFmGGYYbsMAM8AwA8y1Z37vH1091PRUdT9VXdXVVf39nMOhu7qm+nm6qn711FPPRVQVRESULKdFnQAiIgoegzsRUQIxuBMRJRCDOxFRAjG4ExElEIM7EVECMbgTESUQgzsRUQIxuBMRJdDpUX3xqFGjdMyYMVF9PRFRLK1fv/6wqo7Ot15kwX3MmDGoqqqK6uuJiGJJRPaYrMdqGSKiBGJwJyJKIAZ3IqIEYnAnIkogBnciogRicCciSiAGdyKiBDIK7iIyXkTqRKRBRKa4rPM1EakRkW0i8kSwySSiQm3c24at+9ujTgYVSd5OTCIyAsB0AJ8G0ASgUkTmq2qNbZ2xAG4E8GFVbRORt4SVYCLy56oHVgMAdt/5uYhTQsVgUnK/DECDqjaqai+A2QAmZq3znwCmq2obAKhqS7DJJCIiL0yC+7kA9tneN1nL7N4O4O0iskpE1ojIeKcNicgNIlIlIlWtra3+UkxERHmZBHdxWKZZ708HMBbAxwBcA+BhEXnTsD9SnaGqFapaMXp03nFviIjIJ5Pg3gTgfNv78wA0O6wzT1X7VHUXgDqkgz0llKrit4vrsPvwyaiTQkQOTIJ7JYCxInKRiIwEMAnA/Kx15gL4OACIyCikq2kag0wolZamti7cv7wB3/5TZdRJISIHeYO7qqYATAawGEAtgDmquk1EponIBGu1xQCOiEgNgOUA/ltVj4SVaCodff0DUSeBiBwYjeeuqosALMpaNtX2WgH8yPpHREQRYw9VIqIEYnAnorIzMKDo6O6LOhmhYnAnorLzu6U78O5bX0Tbyd6okxIaBnciKjsLqw8AAI52MrgTEVGMMLgTESUQgzsRUQIxuBMRJRCDOxFRAjG4ExElEIM7EVECMbgTESUQgzsRlS3NnnYoQRjciaj8OM0vlzAM7kRECcTgTuTB1v3t0CTfy1NiMLgTGXp1Rys+f99KPL5mT9RJIcqLwZ3I0J6jnQCA7QePR5wSovwY3ImIEojBnYgogRjciYgSiMGdiCiBGNyJiBLIKLiLyHgRqRORBhGZ4vD5dSLSKiKbrH/fCT6pRERk6vR8K4jICADTAXwaQBOAShGZr6o1Was+paqTQ0gjERF5ZFJyvwxAg6o2qmovgNkAJoabLCIiKoRJcD8XwD7b+yZrWbYvi0i1iPxVRM4PJHVEROSLSXB3Gj8te3CN5wCMUdV3A1gK4FHHDYncICJVIlLV2trqLaVERGTMJLg3AbCXxM8D0GxfQVWPqGqP9fYhAO9z2pCqzlDVClWtGD16tJ/0EhGRAZPgXglgrIhcJCIjAUwCMN++goi81fZ2AoDa4JJIRERe5W0to6opEZkMYDGAEQBmqeo2EZkGoEpV5wP4gYhMAJACcBTAdSGmmShSHPCX4iBvcAcAVV0EYFHWsqm21zcCuDHYpBGVljKYvIcShD1UiYgSiMGdiMpYcivZGNyJqOyUQxUbgzsRUQIxuBMRJRCDOxFRAjG4ExElEIM7EVECMbgTESUQgzsRUQIxuBMRJRCDOxFRAjG4ExElEIM7kUea3OFIKEEY3IkMSTkMSEKJweBORJRADO5ERAnE4E5ElEAM7kRECcTgTkSUQAzuREQJxOBORJRADO5EVLaS3CGNwZ2Iyo6UQY80o+AuIuNFpE5EGkRkSo71viIiKiIVwSWRiIi8yhvcRWQEgOkArgQwDsA1IjLOYb2zAPwAwNqgE0lERN6YlNwvA9Cgqo2q2gtgNoCJDuvdBuAuAN0Bpo+IiHwwCe7nAthne99kLRskIpcCOF9VFwSYNqISleCncJQYJsHd6cnD4NEtIqcB+B2AH+fdkMgNIlIlIlWtra3mqSQqAeJ4KhCVJpPg3gTgfNv78wA0296fBeBfALwsIrsBXA5gvtNDVVWdoaoVqloxevRo/6kmIqKcTIJ7JYCxInKRiIwEMAnA/MyHqtquqqNUdYyqjgGwBsAEVa0KJcVERJRX3uCuqikAkwEsBlALYI6qbhORaSIyIewEEhGRd6ebrKSqiwAsylo21WXdjxWeLCIiKgR7qBIRJRCDOxFRAjG4ExElEIM7EVECMbgTESUQgzsRUQIxuBMRJRCDOxFRAjG4ExElEIM7EVECMbgTUdlK8sj8DO5EVHbKYWR+BnciogRicCciSiAGd6IA9Q8oPnTHMszf3Jx/ZaIQMbgTBehkbwrN7d24+W9bok4KlTkGdyKiBGJwJyJKIAZ3IqIEYnAnCkGSO8dQPDC4EwWo2J1jbpm7FWOmLCzyt1IcMLgTxdhja/ZEnQQqUQzuREQJxOBO5JGyQp1iwCi4i8h4EakTkQYRmeLw+XdFZIuIbBKRlSIyLvikEkVLPFSoK68AFLG8wV1ERgCYDuBKAOMAXOMQvJ9Q1Xep6iUA7gJwT+ApJYoB8XIFIAqRScn9MgANqtqoqr0AZgOYaF9BVTtsb88EW4IREUXqdIN1zgWwz/a+CcAHslcSke8D+BGAkQA+4bQhEbkBwA0AcMEFF3hNKxERGTIpuTvdZw4rmavqdFW9GMD/APiZ04ZUdYaqVqhqxejRo72llIiIjJkE9yYA59venwcg13imswF8sZBEEcUd6yUpaibBvRLAWBG5SERGApgEYL59BREZa3v7OQD1wSWRKD74OJVKRd7grqopAJMBLAZQC2COqm4TkWkiMsFabbKIbBORTUjXu18bWoqJKBH2He3EzJW7Ik1DklusmjxQhaouArAoa9lU2+sfBpwuIkq4b85ciz1HOnHVpefinDNHFvW7y6HFKnuoEoUgySXCoBzvTgFgh6+wMLgTBagcSoQUDwzuREQJxOBOVIKa2jpRd/B41MkoClbKhMPogSoRFde//no5AGD3nZ+LOCXhYQ1WuFhyJ/LI5PmfsjxKEWNwJzJkL2lu3NuGituXoL2zL2sdlkepNDC4E/lw77J6HD7Ri/V7j0adlNhjS8hwMLgTUSTYbDRcDO5ERAnE4E4UgnKtajjZk8LJnlTUySCwKSRRoMq9quGdP18MwFsTTrYsCgdL7kQUkTK/EoaMwZ3yWll/GJf/ahm6evujTgoRGWJwp7x+tagWBzu6sbP1RNRJISJDDO5EIWAtsgf8sULB4E5EkSj3h89hY3AnIkogBnciIo8WVDfj4799GQMDpVunxHbuFEt9/QMAgNeNKH75xKRdNmsczEUZHv22sf/vp6vR1deP7lQ//m5kaYZRltwplt7/y6V4960vFvU7vdQRl255rnREeQHMjN7Zl1KkrIJC0jC4Uywd6+xDV1/07e7LdZiBpPjC/Svxmd+9GnUyQsHgTuSDsKlHYjQePhl1EkJhFNxFZLyI1IlIg4hMcfj8RyJSIyLVIrJMRC4MPqlE8XWiJ4WHVzRCWdQfJkk/ydrGI5i5clfUyQBgENxFZASA6QCuBDAOwDUiMi5rtY0AKlT13QD+CuCuoBNKFCtZAeu252pw+8JaLK9riSY9Ifn0Pa/g+j9V+vrbJN78XD1jDW5bUBN1MgCYldwvA9Cgqo2q2gtgNoCJ9hVUdbmqdlpv1wA4L9hkEsVbe1d6Or6evmQ9vKtvOYFl25N1wUoKk+B+LoB9tvdN1jI31wN4vpBEESVVgmogjDW0nMDDKxqjTkbZMWmg6XTz5HiMisg3AVQA+DeXz28AcAMAXHDBBYZJJIq/JFZBmLrqgVU43p3CdR8ag9Md+iVwPPdwmJTcmwCcb3t/HoDm7JVE5FMAbgYwQVV7nDakqjNUtUJVK0aPHu0nvUQlLUkPB4OSmZkpu4WRsKtXqEyCeyWAsSJykYiMBDAJwHz7CiJyKYAHkQ7srICjsjW7ci8AoNelY0w5B3+2FCquvMFdVVMAJgNYDKAWwBxV3SYi00RkgrXabwC8AcDTIrJJROa7bI4o0ZbWHnJcXs7VMm59AlgdEy6jQRFUdRGARVnLptpefyrgdBElCgut8dFt9Xx+/etGRJySwrCHKhEVRfb1rVTr3N9xywv44B3Lok5GwRjciQwVUvou62qZqBPgQ1tnn9F6pXxHxuBOvjQf64o6CSWpVEujFKw4XKwZ3H3o6u3HMpcHZ+Xi6hlrAACp/hIuugQsiBO6FB8iHmjvwpETjq2XA+VWyi3l0m+cMbj7cMu8rbj+0SrUNHdEnZRYuW9ZPeZt2h91MiJRyiX6D97xEt53+9LQtu92UYxD6TfOSnMKkRK350h6iNATVucMMnP3kh0AgImX5Bq9It7sAeuVHa3457eehbec9froEkSR23e0E/0DijGjzizq97LkXgB2yqBcrp21Dl/739eGLCvnQ6YUq6TCtnx7Cz5y13J87LcvF/27Gdx9KOVbbApOV28/nlnfVNBFfPeRzvwrJZzT+fJvv1mOA+3dAJI9mFqUw/8yuBO5mLagBj9+ejPWNB6NOimJYL9G7on4ohdFfX9Ncwfmbx42LFdoGNwL4FbiOHyiB9+auRZtJ3uLmh4KVktHumR5MohnK+V8s1fOebf57L0r8IMnNxbt+xjc/chzsM5cuQsr6g/jiXV7i5MeoiJRVazeebjsnzfFIfsM7gWIww6mwgW5mzPbUlWsaTwSuyD55Lp9+PpDa/Fc9YHAthnWb3DkRA/2HQ23+sdevZMZk6ZUMLj7kO8us5jn6+qdh/H9v2yILEiUYwsIP7KPmbmb9mPSjDV4ZkO82v3vOZpuBry/zbyHclS1MpffsQwfuWt5KNvucgjkN/1tSyjf5ReDewHyBbZiPLS5blYlFm454Dp+OBXOy250G942W+aB4t6QS5blrC+k3tMr6lsdl2/Y2xbK9/nF4O4De9ZRYGJWLVOIpGS17uDxqJNghME9JmoPdGDMlIXYtO+Y4+dJOXHKRan3lVBV/N/HqvDKDudSqpfquHyFIR674WBwN3T0ZC8eXtE4tG7bbSAkwwN/eV0LxkxZiIaWE3nXfWl7evbCxdsODv2gtGNEIuUKRknZHQMKLN52CN9+ZN2Q5faLkulzniQHb3veTKvkioVjyxj68ZxNWF7Xioox5xiXuvKtt9BqcbBhbxv+8S1vMNpmkk+UUue0P73uj+yAGOfdWcwOOeQdS+6GOrrTHVlStgeXriem1zPWYH3XkfU8flWpUFXc+LctJfcQyonT7nH73d12ZXapLspC3n3L6gPZTutxs2GCS6xAWzYY3H0I6mD1sxm3Kp+4lei7+vrx5Lq9+MZDa6NOSqSi2G+Z0TkLMXvdPs9/k5Rms6VW/eKGwT0iqxoO4+n1TQDMDnq3Kp5iHGepgfTdSpOHts2mknLCZ7jtjuzlW/a3h52UUGSONy9NOIN6ePzC1oM41lm6Q3qUWshncC9AvlJXrsA7c+WuYNMSYpDccSj9wPfmZ4PrpFHqrUWC5BSQltSkZ/JK2sUtLC0d3fju4+vxvcc3RJ2U2By5DO4hMDld7QeIp1vzrHXjHiTjVp2US6p/AMe6hk+sHFZnmrgpZF/3pNJ3j6XQ6SsmtTIM7mHKdQzYDxCji0FMDihTmfwkKez9bO5WbHboh5C0fedVvvx/5K7l6OgeflEM0smeFDp7U9h/rAvtDhfgQJTYfjYK7iIyXkTqRKRBRKY4fP5REdkgIikR+UrwySwdCvOTVQHMqdyHD9/5ksOn/o4Et2AYpxJw87Hg6+5LwXM+mgbGab8B4cWvjXudO+cNfm+BX/zOny/GJb9Ygg/f+RI+dc8r6B/w/8OXWAx3lTe4i8gIANMBXAlgHIBrRGRc1mp7AVwH4ImgE1gKpj1Xg/V7cjfZ29bcji/ctxKdvUPH/v7pM9XY7xDMhpTcTZpCui33caR94u6X8fEIpv3K6LOPgxOz4OaHSRXclqZ2PLZmT1HS42bWyl3YcSi8rvW5dnWm/f/9L9VjzJSFedfzIzP+UuvxHgzE7arqg0nJ/TIADaraqKq9AGYDmGhfQVV3q2o1gESOXjVrlfPDT/vDsDsWbceW/e1DLgK54u5ptg9venbLsLk23fT1D2Bh9QHXzjBff2gNrsvqVZitsfUkdh0+afR9dkFVL6jaq2W8n2TtneHewnsxdd5Wo/XcfrtM7r9w/0rcMtdsW2GZtqAGn793Zc51/IREk8Mms93fvujcTDPT/LC5vXvIxPSp/gHMqdyHAaskftuCGlzxu1d9pNIlXap4ua7F6KJSaiV6k+B+LgB7o9Yma5lnInKDiFSJSFVrq/OYFXHg9BAzUxJQdS9d1B7owBenr0Jnb2rYNtbtzj2VWyY4PLJqN77/xAa8sPWglZahVu88gpfrwvltgyzs+H0QvLP1BN4z7cXgEpKDW0Du6O7DMms4iGZrHtCk6O0fQFdvPy6+aVHUSRnCviuu/1Pl4OtZq3bhp89U46mqdIiauXIX6gK8+3h2435c90glnrS1609SO3ennPg6zVV1hqpWqGrF6NGj/WyipGSCXXp2miPp17bPs4+BXy6sxaZ9x7B+T1vBpeAjIU7h19Xbj3W7zOYNPdRh1kvRbk3jEWy0eqZ6vWDsah16x9FyvBtjpizEbGvWq288vAbvvnWx5zTls7yuBX1We//F2w55+lt7MIhD08e2HG3JCzlsc5Z+Pfws9rvjzHlwLKS7uczzof3H0q10Zq7chVddB1MrLSZjyzQBON/2/jwAZT2oRHZg9vNsxmtwL6TJ4/aDHRj1hjMw6g1nGK0/5W/VmLepGSt++nGcf87f+f5e9+2fai+vSJ+sZ54xAu/4hzd63tbuw+mT7pkNTbjkgjdhVcMRx/UWVDejrbMPq+oP4yNvH4VvfOBC4+9YvfMwHlm1G6PPMvv9ssWjnHdK0EHKpKSrUGz107HLltiUhzkN/Ja+b1tQ42n9Rh/Vn0ExKblXAhgrIheJyEgAkwDMDzdZ8eB0EpjUzakW3j5dka7mOdnbn/d7x/9+BT559yvG295+IH1bezLr4XAYd6P9A4ov/3E1xv9+ha+/t+c71zYmP7ERt8zdihe2HcTNz3qr2z58Il06dBtLZUV9Ky7/1TJ0p/IHF6/7/XuPr/e0/jPrmwoer8f+m4YxxaCbz983vL7/4RWN2He00+jY6wxxmju302vInXpo3+5P3uCuqikAkwEsBlALYI6qbhORaSIyAQBE5P0i0gTgqwAeFJFtYSa61AxpDQH3A8F+gHouuWetf+BYF678g3lADK1tb5Fl/w6vWrPihNmZK9+Wf7mwFgc7uo2a12VXy+QrCzy/9WDuFbL8+OnN+NIDqz39TdiMHqhq9ntFS0c3bl9Yi2sfWRdZZ7241K87MWrnrqqLVPXtqnqxqv7SWjZVVedbrytV9TxVPVNV36yq7wwz0aVurVVfnesusdCD5oGXd3r+mx88udHT+nFoLbbbmq4uzHM/364y2ZeFJu8rf1yNGz3M0WlvUeKV235/ua4FLYYjQRb6naqnqjtP9qTydvoLOwbH4FQYhj1U88hV3ZH5bMgaempQqMMnnE8ERfCxyOTgMx1/u1QKKxv2tmF5XcuQZcPSZmXcrcpk+faWnO2mg5Dv58rVp8HkAetTlXtRtacNT67bi/auPqO65YXVwT8Wu+6RSvzVGuzOj6AKC/0DWvBUd27nZiFKrZTP4J7H/3tq07BluXaiaWuI1xqHP/jrSfUP6wRVzr70wGp8+5HKnOtkfm+3dvtzN+13XD4woMYdYrwEb+e/F08nfnap+3+eOVVif88vXhzyQDrjeLdZ0I+EQdaz90T2++xNXPH7dFv2zU3pnq3p51jmwmpdU0oY3POYu8m9BJQvNLiN/HjtrHWOJc0r/7AC46Y6N+PL1zQxiFLRPUt24J9veaHwDRVRrnyv33MU81z239tuWmRczZEvMHstsOWrMrn/pYacn8/dOPyC9a5bX8SPn97sLSEuSqE6zvTCu6YxfV6EXShyfY4W6rcWhsHdh1w7NOf8mnmiQGOre7OpF2u8ta32495l9eiytTgohZM8W/aDtVzdyK9+cE3Obc2uPNUx5URPCivrD7t8Z26n5dmv9rs5VeC+l+rtHw7jtwTudiHzKrS2+LmauecdPzvYpJgqsZoWTxjcE6LSsNORm0LG7IhSrrpXLzn60VOb8M2Za30NamZy/mfuvF7Z0Yq+VDx/a79Mfp9HX9s95L3XX6gUxooptesAg7tHQ46h6I+nQd/5c1VBf3+w41Q3+j1H3MbMLrXD19ZaxoGXC1Z9S3pCki6nttKFVroDg4PHvVzXMmR4Wz+HUGpA8YFfLcVL2wu7m8s3SBeQO+tHDXtJm+RxbePQwkn2rov7vAVRYHD34VDH0PFEhnT6cDmSBwbUtdtyKbCn2zHAIT2mysU3LcKK+uLn4/2/XIr33bbE09+YBBWTAJUvsOQLO/bftqM7VVCLk4xDHT24bUGt0bpVu4+ioWX4HY7bKJSm18RcI0juc5hUI1Pdc9BhPJ7skre9auhQRw+OnMzduqUECu4lx2T4Acqy3aoKaG4ffgvvtAwAXtre4ri8VJjULfZaPTDzPfALQ+YBdHWTWRf1oyd7jU74ts5enHPmyEKSVlC9bCHVYaZ/+xXDEUe9f7/7Z91WAaG9q2/IxW/drqP42oNm6ak92DH4esozuR9+hx3bjeY5LrGbC5bcCzAwoLhl7tbB5lgAMHWec+fcHoOu6XERdm/Xqx98DfNcmjC+vMPsImk6FG/mYWgmUPo5P7fkueCEFXgK3a7bHUkQ6VUAy2oP4T2/eBEd3emWLKrpITOcZHfuVQUWVh8YfL/JYYar7PUL8fWHnB++57trmzpvW+j9KPxicC9AZ28/HluzB9fMWBt1UgrmdBC7lVa2Zz3EDHpG+rW7juKHs4f3LwDyz9iTYXoxPS0r2xPuXwUAmLdpP5bWpi8k+R7WpQqY1eehFbuwO6uNvunWCg1ormPMB1DHMaA62FPbz7Z3tp4Ytm9yUWhBnYgyo7p69cyGwqvYwsLg7lHL8eH1hXEYxtXJAlsvRrceniYumbYEPanwBm3yw3QatRe2HsTkJzYMvs+0QbdfXJ51aFcepP/0+TC8mJNFmwTlXtsF1a1Tkel18Ncv1OVtYmr34CuN+JefBzvU85gpCweHpo7jKc7g7tHkJ06NzzI4m1AJ7viHVzQ6LrfPYmTPy5IabwNUZevNKilv3NvmbwjXgJg+47jj+e1YYLv9D0O+kn+mpU5Gsapu3b7HzyxdwNDhLVSBmqwqGIWXu5LCSuLZ2/Ir07/k2Y37Pc0A1j+grlWLxcLgXoBMVYZJG1u/x2n/gOIXz3kfZPMPS+sdl7vOYlTgiZRdIrvqgdWOQ7iWo0dX7446CZ78bumpqe4yu9XpEM9elN35akVWp7Duvn7jQLv/WBfqA5pR6RbD5y8tHe6zarUc78F/PeU+8F52vr7+0BrXqsViYXAPQAFVrnllJorwyjRJzce68MJW55KrpwJPjnVP9qRwu8dJDorNXi21JODewE7NAkuZ34Ku/Txwqqq8+VnzUS0bW0+iKs+k9AOGJ97ja/YarffVPK14cs06lp0Up+cNxcamkAUIu659/uZmz8P02m3c2zY416qbLz2wGgc7uvHDT471/T1A+u5l7sb9aGg5gZ9c8U9DPvvZ3K0F1VtPey78C0Nm0hPAfx24m0KG3y3E+j1trq1TAPfhMJxipsmRbr+DdbpA1B44jg9dPMpgS2Z+v8z57jT9/d7PzT1HOn31UAbSd9ilVj1bdsF93a6jqLjwbJzm5VF8RJ5Y69zJxMSJnhSuMpi0IdMztdAWEpfaOhhlB3fTwH6ooxt//8bXD1s+a5XzAGxxkWl1U2xf/qO/STuyO+X1pgZOPVh00ZsaMDqGghwm4N4cwd0vP1WgQDpfpTAEgl1ZVcusqG/F1x58DTNcHjZ65aVLtHuXfme3L6gZHPHO1PFu/+3P73XomFS1+yg+62G2JydeLhof+NWygr4rKfb7LD3a5Rt7bE7VPuPv+cnTm3H1jOHtwO379tpZ64a0UHLa7Qc7ujGnat/wD0Kw49CJ/Cs5yNVCJ3vaSTvV0mtYUVbB/YDV7bmhxd+Oz+blSv3rF7Ybr/vbxXV42GW44Fz+stasbtHUrc/VDGvx4NVFNy4KKDXlw+vUek5uylO//dO/Vrt+lt2PwW2Sl+V1p4aheK3xSN46dwDYmWPk0yBlxnv3Kldwz1VA69fSaxAdu+B+rLMXS30+8MrsuKBun+543jxge3H/cn/d++8MKT3FNOH+8m5hc6lba6YYOGl7tpDdUiYuCmnVZtq3olhiF9y/9/gGfOfPVb6myTqthNulU5rp2DFJ1RbjGYLuXnKqCeVvFtdFmBL/vHScslNV405lpq18ChW74L7nSPq2buv+ds+dLewl923N7RgzZSHWOkx3B6THsni6SPWDSVPdZDZEAFGp6fdZ8vNSan/bTYuKUsqPXWuZTCuX66y5NXff+Tmjv9u879hgr8UBBVY1pG8bl9Yewgfe9uZh639x+qogkluWMuOzEMXNQp+9lU3b0mf09Q9gxGkjfH2XqdiV3EfkaMKY6h/AY2v2OE5TNnH6qsEHQwMDOvjwx6mtr9PY10REbuy9ek0UYwhwo+AuIuNFpE5EGkRkisPnZ4jIU9bna0VkTNAJzciuEzve3YdFWw6g5Xg3vvv4etwyd+vgxNSqil2HTw67Bero7husd5/xauNgt+hMXdin7vH3pJ2IyEQhA/WZylstIyIjAEwH8GkATQAqRWS+qtq7DV4PoE1V/1FEJgH4NYCrw0hwdsH9XbcOb12wad8xzN/cjBe3HcSC6gO49II3Dfl8Rf1h1DSfauL3jlteGHx96xfGBZtgIqIshQwTbcqk5H4ZgAZVbVTVXgCzAUzMWmcigEet138F8EkJaki3LCZPs5/fehA/eHKu7BLmAAAFVklEQVTj4Gh/TmOAH3GZXu3WInR1J6LyNt1nc2cvTIL7uQDszUaarGWO66hqCkA7gOFPKQOQq86diCgOTCcXL4RJcHeKptn3FCbrQERuEJEqEalqbfU3yfI3Lr/Q198REZWKu7/6ntC/w6QpZBOA823vzwOQ3R85s06TiJwO4P8AGDYwiqrOADADACoqKnxVOn3r8gvxLQZ4IqKcTErulQDGishFIjISwCQA87PWmQ/gWuv1VwC8pEFMxEhERL7kLbmrakpEJgNYDGAEgFmquk1EpgGoUtX5AGYCeExEGpAusU8KM9FERJSbUQ9VVV0EYFHWsqm2190Avhps0oiIyK/Y9VAlIqL8GNyJiBKIwZ2IKIEY3ImIEojBnYgogSSq5ugi0gpgj88/HwUgnvN4eVMO+WQek6Ec8giURj4vVNXR+VaKLLgXQkSqVLUi6nSErRzyyTwmQznkEYhXPlktQ0SUQAzuREQJFNfgPiPqBBRJOeSTeUyGcsgjEKN8xrLOnYiIcotryZ2IiHKIXXDPN1l3qROR3SKyRUQ2iUiVtewcEVkiIvXW/2dby0VE7rXyWi0i77Vt51pr/XoRudbt+4pBRGaJSIuIbLUtCyxPIvI+6zdrsP626NNxueTxVhHZb+3LTSLyWdtnN1rprRORK2zLHY9fa0jttVben7KG1y4qETlfRJaLSK2IbBORH1rLk7Yv3fKZqP0JVY3NP6SHHN4J4G0ARgLYDGBc1OnymIfdAEZlLbsLwBTr9RQAv7ZefxbA80jPdHU5gLXW8nMANFr/n229PjvCPH0UwHsBbA0jTwDWAfig9TfPA7iyRPJ4K4CfOKw7zjo2zwBwkXXMjsh1/AKYA2CS9fp/AXwvgjy+FcB7rddnAdhh5SVp+9Itn4nan3EruZtM1h1H9gnGHwXwRdvyP2vaGgBvEpG3ArgCwBJVPaqqbQCWABhf7ERnqOqrGD7zViB5sj57o6q+pukz5c+2bRWNSx7dTAQwW1V7VHUXgAakj13H49cqvX4C6cnlgaG/V9Go6gFV3WC9Pg6gFun5kZO2L93y6SaW+zNuwd1ksu5SpwBeFJH1InKDtezvVfUAkD7wALzFWu6W3zj8DkHl6VzrdfbyUjHZqpKYlamugPc8vhnAMU1PLm9fHhkRGQPgUgBrkeB9mZVPIEH7M27B3Wgi7hL3YVV9L4ArAXxfRD6aY123/Mb5d/Cap1LO6x8BXAzgEgAHANxtLY91HkXkDQCeAfBfqtqRa1WHZXHOZ6L2Z9yCu8lk3SVNVZut/1sAPIv0rd0h65YV1v8t1upu+Y3D7xBUnpqs19nLI6eqh1S1X1UHADyE9L4EvOfxMNJVGqdnLS86EXkd0gHvL6r6N2tx4valUz6Ttj/jFtxNJusuWSJypoiclXkN4DMAtmLoBOPXAphnvZ4P4N+tVgmXA2i3bosXA/iMiJxt3Tp+xlpWSgLJk/XZcRG53KrL/HfbtiKVCXiWq5Del0A6j5NE5AwRuQjAWKQfJDoev1b983KkJ5cHhv5eRWP9vjMB1KrqPbaPErUv3fKZtP1Z1Ke3QfxD+gn9DqSfUt8cdXo8pv1tSD9R3wxgWyb9SNfRLQNQb/1/jrVcAEy38roFQIVtW/+B9IOdBgDfjjhfTyJ9G9uHdGnm+iDzBKAC6RNtJ4D7YXW+K4E8PmbloRrpAPBW2/o3W+mtg61FiNvxax0b66y8Pw3gjAjy+K9IVx9UA9hk/ftsAvelWz4TtT/ZQ5WIKIHiVi1DREQGGNyJiBKIwZ2IKIEY3ImIEojBnYgogRjciYgSiMGdiCiBGNyJiBLo/wOjKVH1+iz6DQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(clf.components_[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.71292605795277"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.reconstruction_err_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NMF in summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Benefits: Fast and easy to use!\n",
    "\n",
    "Downsides: took years of research and expertise to create"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notes:\n",
    "- For NMF, matrix needs to be at least as tall as it is wide, or we get an error with fit_transform\n",
    "- Can use df_min in CountVectorizer to only look at words that were in at least k of the split texts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Truncated SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We saved a lot of time when we calculated NMF by only calculating the subset of columns we were interested in. Is there a way to get this benefit with SVD? Yes there is! It's called truncated SVD.  We are just interested in the vectors corresponding to the **largest** singular values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/svd_fb.png\" alt=\"\" style=\"width: 80%\"/>\n",
    "\n",
    "(source: [Facebook Research: Fast Randomized SVD](https://research.fb.com/fast-randomized-svd/))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Shortcomings of classical algorithms for decomposition:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Matrices are \"stupendously big\"\n",
    "- Data are often **missing or inaccurate**.  Why spend extra computational resources when imprecision of input limits precision of the output?\n",
    "- **Data transfer** now plays a major role in time of algorithms.  Techniques the require fewer passes over the data may be substantially faster, even if they require more flops (flops = floating point operations).\n",
    "- Important to take advantage of **GPUs**.\n",
    "\n",
    "(source: [Halko](https://arxiv.org/abs/0909.4061))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Advantages of randomized algorithms:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- inherently stable\n",
    "- performance guarantees do not depend on subtle spectral properties\n",
    "- needed matrix-vector products can be done in parallel\n",
    "\n",
    "(source: [Halko](https://arxiv.org/abs/0909.4061))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Timing comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 28s, sys: 9.3 s, total: 1min 37s\n",
      "Wall time: 4.06 s\n"
     ]
    }
   ],
   "source": [
    "%time u, s, v = np.linalg.svd(vectors, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import decomposition\n",
    "import fbpca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 18s, sys: 1.4 s, total: 1min 19s\n",
      "Wall time: 3.03 s\n"
     ]
    }
   ],
   "source": [
    "%time u, s, v = decomposition.randomized_svd(vectors, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Randomized SVD from Facebook's library fbpca:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 26 s, sys: 644 ms, total: 26.7 s\n",
      "Wall time: 1.19 s\n"
     ]
    }
   ],
   "source": [
    "%time u, s, v = fbpca.pca(vectors, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more on randomized SVD, check out my [PyBay 2017 talk](https://www.youtube.com/watch?v=7i6kBz1kZ-A&list=PLtmWHNX-gukLQlMvtRJ19s7-8MrnRV6h6&index=7).\n",
    "\n",
    "For significantly more on randomized SVD, check out the [Computational Linear Algebra course](https://github.com/fastai/numerical-linear-algebra)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}