NLP Lesson 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"width: 100%; overflow: hidden;\">\n",
    "    <div style=\"width: 150px; float: left;\"> <img src=\"data/D4Sci_logo_ball.png\" alt=\"Data For Science, Inc\" align=\"left\" border=\"0\"> </div>\n",
    "    <div style=\"float: left; margin-left: 10px;\"> \n",
    "        <h1>Natural Language Processing From Scratch</h1>\n",
    "        <h2>Topic Modeling</h2>\n",
    "        <p>Bruno Gonçalves<br/>\n",
    "        <a href=\"http://www.data4sci.com/\">www.data4sci.com</a><br/>\n",
    "        @bgoncalves, @data4sci</p></div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import string\n",
    "from collections import Counter\n",
    "from pprint import pprint\n",
    "import gzip\n",
    "import matplotlib.pyplot as plt \n",
    "import numpy as np\n",
    "from numpy.linalg import svd\n",
    "from numpy import diag\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load a standard list of stopwords from a file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "stopwords = set([word.lower().strip() for word in open(\"data/nltk_stopwords.txt\", \"rt\").readlines()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to extract the words from our input file. Here we have to be (slightly) more careful than in the previous lesson since we want to keep #tags and @mentions and remove stopwords."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_words(text, stopwords):\n",
    "    temp = text.split() # Split the text on whitespace\n",
    "    text_words = []\n",
    "\n",
    "    punctuation = set(string.punctuation)\n",
    "    \n",
    "    #Keep #tags and @mentions\n",
    "    punctuation.remove(\"#\")\n",
    "    punctuation.remove(\"@\")\n",
    "    \n",
    "    for word in temp:\n",
    "        # Remove any punctuation characters present in the beginning of the word\n",
    "        while len(word) > 0 and word[0] in punctuation:\n",
    "            word = word[1:]\n",
    "\n",
    "        # Remove any punctuation characters present in the end of the word\n",
    "        while len(word) > 0 and word[-1] in punctuation:\n",
    "            word = word[:-1]\n",
    "\n",
    "        # Simple rule to eliminate (most) URLs\n",
    "        if len(word) > 0 and \"/\" not in word:\n",
    "            # If it's not a stopword\n",
    "            if word.lower() not in stopwords:\n",
    "                # Append this word into our list of words.\n",
    "                text_words.append(word.lower())\n",
    "        \n",
    "    return text_words"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Process the data from the file. File was obtained from crowdflower (https://www.crowdflower.com/data-for-everyone/) and contains a list of tweets about Apple from 2014 and some information about their sentiment. For our purposes, we only care about the text of each tweet that is located in the last column of the csv file.\n",
    "\n",
    "We treat each tweet as an individual document, and limit ourselves to only the first 200 tweets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "tweets = []\n",
    "line_count = 0\n",
    "\n",
    "for line in open(\"data/Apple-Twitter-Sentiment-DFE.csv\", \"rt\"):\n",
    "    fields = line.strip().split(',')\n",
    "    \n",
    "    line_count += 1\n",
    "    \n",
    "    # Skip the first line of the file which contains the header\n",
    "    if line_count == 1:\n",
    "        continue\n",
    "    \n",
    "    text = \",\".join(fields[11:])\n",
    "    \n",
    "    if len(text) == 0:\n",
    "        continue\n",
    "    \n",
    "    words = extract_words(text, stopwords)\n",
    "    \n",
    "    if len(words) > 0:\n",
    "        tweets.append(words)\n",
    "        \n",
    "    if len(tweets) == 200:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the function to calculate the Inverse Document Frequency for each word and the TFIDF matrix. These functions were introduced in the previous lesson."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inv_doc_freq(corpus_words):\n",
    "    number_docs = len(corpus_words)\n",
    "    \n",
    "    document_count = {}\n",
    "\n",
    "    for document in corpus_words:\n",
    "        word_set = set(document)\n",
    "\n",
    "        for word in word_set:\n",
    "            document_count[word] = document_count.get(word, 0) + 1\n",
    "    \n",
    "    IDF = {}\n",
    "    \n",
    "    for word in document_count:\n",
    "        IDF[word] = np.log(number_docs/document_count[word])\n",
    "        \n",
    "    \n",
    "    return IDF\n",
    "\n",
    "def tf_idf(corpus_words):\n",
    "    IDF = inv_doc_freq(corpus_words)\n",
    "    \n",
    "    TFIDF = []\n",
    "    \n",
    "    for document in corpus_words:\n",
    "        TFIDF.append(Counter(document))\n",
    "    \n",
    "    for document in TFIDF:\n",
    "        for word in document:\n",
    "            document[word] = document[word]*IDF[word]\n",
    "            \n",
    "    return TFIDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the TFIDF matrix for our corpus of tweets. We should note that while we call it a matrix, this is effectively a list of dictionaries, which we can consider to be a sparse representation of a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "TFIDF = tf_idf(tweets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Counter({'#applesnob': 5.298317366548036,\n",
       "          '@apple': 0.2169130015635736,\n",
       "          'cat': 5.298317366548036,\n",
       "          'chews': 5.298317366548036,\n",
       "          'cords': 5.298317366548036}),\n",
       " Counter({'#aapl': 1.9310215365615626,\n",
       "          'expects': 5.298317366548036,\n",
       "          'inquisition': 5.298317366548036,\n",
       "          'nobody': 5.298317366548036,\n",
       "          'spanish': 5.298317366548036}),\n",
       " Counter({'#aapl:5': 5.298317366548036,\n",
       "          'apple': 2.5902671654458267,\n",
       "          'buy': 3.506557897319982,\n",
       "          'december': 4.605170185988092,\n",
       "          'gains': 5.298317366548036,\n",
       "          'rocket': 5.298317366548036,\n",
       "          'stocks': 5.298317366548036})]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TFIDF[:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, extract the word_dict and word_list from the TFIDF matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_vocabulary(TFIDF):\n",
    "    words = set()\n",
    "    \n",
    "    for document in TFIDF:\n",
    "        words |= set(document.keys())\n",
    "    \n",
    "    word_list = list(words)\n",
    "    word_dict = dict(zip(word_list, range(len(word_list))))\n",
    "    \n",
    "    return word_dict, word_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "word_dict, word_list = build_vocabulary(TFIDF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We have 927 words in our vocabulary\n"
     ]
    }
   ],
   "source": [
    "vocabulary_size = len(word_dict)\n",
    "print(\"We have\", vocabulary_size, \"words in our vocabulary\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['amendment',\n",
       " 'ios8',\n",
       " '#hft',\n",
       " '1',\n",
       " 'never',\n",
       " 'protocol',\n",
       " 'dominating',\n",
       " '@remuuusss',\n",
       " '6.6',\n",
       " 'anymr']"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "word_list[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use the TFIDF matrix and our vocabulary to generate the Term Document matrix. This is just a matter of rearranging the values in our (sparse) TFIDF matrix into the full TD matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def term_document_matrix(TFIDF, word_list, word_dict):\n",
    "    vocabulary_size = len(word_dict)\n",
    "    number_documents = len(TFIDF)\n",
    "    \n",
    "    TDM = np.zeros((vocabulary_size, number_documents))\n",
    "    \n",
    "    for doc in range(number_documents):\n",
    "        document = TFIDF[doc]\n",
    "        \n",
    "        for word in document.keys():\n",
    "            pos = word_dict[word]\n",
    "            \n",
    "            TDM[pos, doc] = document[word]\n",
    "            \n",
    "    return TDM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our dataset has:\n",
      "927 unique words\n",
      "200 documents\n"
     ]
    }
   ],
   "source": [
    "TDM = term_document_matrix(TFIDF, word_list, word_dict)\n",
    "print(\"Our dataset has:\\n%u unique words\\n%u documents\"%(TDM.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explicit Semantic Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In ESA we use the TD matrix of our corpus as a knowledge base that we can use to look up related documents. To demonstrate how this works, let's create a new document:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_tweet = ['#macbook', 'mini', 'rocket']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not an actual sentence since, for clarity, we are not including stopwords.\n",
    "\n",
    "In the TD formulation, the values in the word (row) vector representa how much each document is \"activated\" by that specific word. Naturally, documents in which the word does not appear have weight zero. We can use these word vectors to find a vector representation of our new document where each value is our related it is to every other document in our knowledge base. You can think of this as the set of results one might obtain after a query in a search engine!\n",
    "\n",
    "The **find_related** simply calcuates the vector corresponding to the new \"document\" and returns a list of the corresponding weights sorted in decreasing order:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_related_docs(tweet, TDM):\n",
    "    new_vector = np.zeros(TDM.shape[1])\n",
    "    \n",
    "    for word in tweet:\n",
    "        pos = word_dict[word]\n",
    "        new_vector += TDM[pos, :]\n",
    "        \n",
    "    # Now the entries of new_vector tell us which documents are activated by this one.\n",
    "    # Let's extract the list of documents sorted by activation\n",
    "    doc_list = sorted(zip(range(TDM.shape[1]), new_vector), key=lambda x:x[1], reverse=True)\n",
    "    \n",
    "    return doc_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "related = find_related_docs(new_tweet, TDM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look a the top 5 results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "166 ipad mini unboxing via @youtube @apple #ipadmini #ipad #macbook #macbookpro #startup #hipster #unboxing\n",
      "167 ipad mini first time startup via @youtube @apple #ipadmini #ipad #macbook #macbookpro #startup #hipster #unbox\n",
      "190 ipad mini unboxing via @youtube @apple #ipadmini #ipad #macbook #macbookpro #startup #hipster #unboxing\n",
      "191 ipad mini first time startup via @youtube @apple #ipadmini #ipad #macbook #macbookpro #startup #hipster #unbox\n",
      "2 #aapl:5 rocket stocks buy december gains apple\n"
     ]
    }
   ],
   "source": [
    "for tweet, score in related[:5]:\n",
    "    print(tweet, \" \".join(tweets[tweet]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With such small documents, the results are easy to interpret. The first two result are due to mentions of **macbook** and **mini**, the third due to the presence of the word **rocket** but the remaining two results don't include any words in common with our original tweet, just words that appear in similar contexts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def similarity(vec1, vec2):\n",
    "    sim = np.dot(vec1, vec2)\n",
    "    norm1 = np.sqrt(np.dot(vec1, vec1))\n",
    "    norm2 = np.sqrt(np.dot(vec2, vec2))\n",
    "    \n",
    "    return sim/(norm1*norm2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_similar_words(tweet, TDM):\n",
    "    new_vector = np.zeros(TDM.shape[1])\n",
    "    \n",
    "    for word in tweet:\n",
    "        pos = word_dict[word]\n",
    "        new_vector += TDM[pos, :]\n",
    "    \n",
    "    sim = [similarity(new_vector, TDM[i, :]) for i in range(TDM.shape[0])]\n",
    "    \n",
    "    sim_words = sorted(zip(range(TDM.shape[0]), sim), key=lambda x:x[1], reverse=True)\n",
    "    \n",
    "    return sim_words"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "similar = find_similar_words(new_tweet, TDM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#macbook 0.842346298979266\n",
      "@youtube 0.8273747381258668\n",
      "#startup 0.8273747381258668\n",
      "#hipster 0.8273747381258668\n",
      "mini 0.8227410508913036\n"
     ]
    }
   ],
   "source": [
    "for word, score in similar[:5]:\n",
    "    print(word_list[word], score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This approach also returns words that are not included in the document but that are used in a similar context, such as #hipster. One might use this extra words (above some finetuned threshold of similarity) to expand our query and retrieve furhter related documents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Latent Semantic Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LSA analysis, sometimes refered to as Latent Semantic Indexing, relies on the well known Singular Value Decomposition. Singular Values (and Singular Vectors) can be thought of as being a generalization of the more common eigenvalues (and eigenvectors) for the case of non-square matrices. \n",
    "\n",
    "The first step in any LSA analysis is to apply SVD to the Term-Document matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, sigma, vt = svd(TDM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting matrices have the expected shapes, except the singular value matrix, **sigma** that is simply a vector. This is an optimization of the numpy function meant to save memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(927, 927) (200,) (200, 200)\n"
     ]
    }
   ],
   "source": [
    "m, n = TDM.shape\n",
    "print(u.shape, sigma.shape, vt.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "927 200\n"
     ]
    }
   ],
   "source": [
    "print(m, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now reduce the internal dimention of these matrices to, say, 10. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 10\n",
    "\n",
    "# Convert the vector of singular values into a diagonal matrix\n",
    "sigma_k = sigma[:k]\n",
    "Sk = diag(sigma_k)\n",
    "\n",
    "# Drop the extraneous dimensions in the other two matrices.\n",
    "uk = u[:,:k]\n",
    "vtk = vt[:k, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this way, we preserve the most significant (higher sigma values) latent dimensions of the dataset while discarding the rest. The idea is simple, SVD projects the original dataset into an internal space which is a linear combination of the original one. Each dimension in this space can be thought of as an underlying topic in the data. By preserving only the most dominant topics in effect we are simplifying the data and, potentially, removing noise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(927, 10) (10, 10) (10, 200)\n"
     ]
    }
   ],
   "source": [
    "print(uk.shape, Sk.shape, vtk.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are no easy ways to determine the number of topics to choose, but a good rule of thumb is to identify kinks in the values of sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAocAAAHSCAYAAACev2E7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df2zc933f8ddbJ1vkteWcLKpEW56lFlKTOICUgnMDZK3aSOkko6tDrWgUCbW2ZqBFN+uP1YNTFAJScQXSzG26AQpjNcksDVmiLqUdI0uGaklalUAplf4hx6qsyjopjS2KUtpmjSFRno7v/fG57/e+PN2RR/Luvnffez6Awx2/PIqfXI7ky5/P5/3+mLsLAAAAkKQVaQ8AAAAA7YNwCAAAgBjhEAAAADHCIQAAAGKEQwAAAMQIhwAAAIitbOU3e9vb3ubr169v5bcEAABAFc8999x33X115fWWhsP169drcnKyld8SAAAAVZjZt6tdZ1kZAAAAMcIhAAAAYoRDAAAAxAiHAAAAiC0YDs2sx8xOmdlpMztjZr9Tuv6UmV00sxdLty3NHy4AAACaqZ5q5ZuS3ufub5jZHZLGzexrpc/9R3f/UvOGBwAAgFZaMBy6u0t6o/ThHaWbN3NQAAAASEddew7NLGdmL0q6Kum4u58sfep3zewlM/ukma1q2igBAADQEnWFQ3cvuvsWSeskPWBm75L0W5LeLumfS3qrpMerfa2ZDZnZpJlNXrt2rUHDBgAAQDMsqlrZ3b8n6c8k7XD3KQ9uSvpvkh6o8TWH3X3A3QdWr77thBYAAAC0kXqqlVeb2V2lx72Stkt6xcz6S9dM0gckvdzMgQIAAKD56qlW7pd0xMxyCmHyj939K2b2DTNbLckkvShpfxPHCQAAgBaop1r5JUnvrnL9fU0ZEQAAAFLDCSkAAACIEQ4BAAAQy2w4nJqStm6VrlxJeyQAAACdI7PhcGREGh+XDh5MeyQAAACdI3PhsLdXMpNGR6XZ2XBvFq4DAABgfpkLh4WCtGePlM+Hj/N5ae9e6eLFdMcFAADQCTIXDvv7pb4+aWZG6ukJ93190tq1aY8MAACg/WUuHErS9LS0f780MRHuKUoBAACoTz0npHScsbHy40OH0hsHAABAp8nkzCEAAACWhnAIAACAGOEQAAAAMcIhAAAAYoRDAAAAxAiHAAAAiBEOAQAAECMcAgAAIEY4BAAAQIxwCAAAgBjhEAAAADHCIQAAAGKEQwAAAMQIhwAAAIgRDgEAABAjHAIAACBGOAQAAECMcAgAAIAY4RAAAAAxwiEAAABihEMAAADECIcAAACIEQ4BAAAQIxwCAAAgRjgEAABAjHAIAACAGOEQAAAAMcIhAAAAYoRDAAAAxAiHAAAAiBEOAQAAECMcAgAAIEY4BAAAQIxwCAAAgBjhEAAAADHCIQAAAGKEQwAAAMQIhwAAAIgRDgEAABAjHAIAACC2YDg0sx4zO2Vmp83sjJn9Tun6BjM7aWbnzeyYmd3Z/OECAACgmeqZObwp6X3uvlnSFkk7zOw9kn5P0ifdfaOkf5D04eYNEwAAAK2wYDj04I3Sh3eUbi7pfZK+VLp+RNIHmjJCAAAAtExdew7NLGdmL0q6Kum4pAuSvufut0pPeU3SPTW+dsjMJs1s8tq1a40YMwAAAJqkrnDo7kV33yJpnaQHJL2j2tNqfO1hdx9w94HVq1cvfaQAAABoukVVK7v79yT9maT3SLrLzFaWPrVO0uXGDg0AAACtVk+18mozu6v0uFfSdklnJX1T0i+UnrZP0pebNUgAAAC0xsqFn6J+SUfMLKcQJv/Y3b9iZn8t6Ytm9p8kvSDps00cJwAAAFpgwXDo7i9JeneV6wWF/YcAAADICE5IAQAAQIxwCAAAgBjhEAAAADHCIQAAAGKEQwAAAMQIhwAAAIgRDgEAABDrinA4NSVt3SpduZL2SAAAANpbV4TDkRFpfFw6eDDtkQAAALS3TIfD3l7JTBodlWZnw71ZuA4AAIDbZTocFgrSnj1SPh8+zuelvXulixfTHRcAAEC7ynQ47O+X+vqkmRmppyfc9/VJa9emPTIAAID2lOlwKEnT09L+/dLERLinKAUAAKC2lWkPoNnGxsqPDx1KbxwAAACdIPMzhwAAAKgf4RAAAAAxwiEAAABihEMAAADECIcAAACIEQ4BAAAQIxwCAAAgRjgEAABAjHAIAACAGOEQAAAAMcIhAAAAYoRDAAAAxAiHAAAAiBEOAQAAECMcAgAAIEY4BAAAQIxwCAAAgBjhEAAAADHCIQAAAGKEQwAAAMQIhwAAAIgRDgEAABAjHAIAACBGOAQAAECMcAgAAIAY4RAAAAAxwiEAAABihEMAAADECIcAAACIEQ4BAAAQIxwCAAAg1lXhcGpK2rpVunIl7ZEAAAC0p64KhyMj0vi4dPBg2iMBAABoT10RDnt7JTNpdFSanQ33ZuE6AAAAyroiHBYK0p49Uj4fPs7npb17pYsX0x0XAABAu1kwHJrZvWb2TTM7a2ZnzOzXStc/Zmavm9mLpduDzR/u0vT3S3190syM1NMT7vv6pLVr0x4ZAABAe1lZx3NuSfpNd3/ezH5I0nNmdrz0uU+6+xPNG17jTE9L+/dLQ0PS4cOhOAUAAABzLRgO3X1K0lTp8ffN7Kyke5o9sEYbGys/PnQovXEAAAC0s0XtOTSz9ZLeLelk6dJHzOwlM/ucmb2lxtcMmdmkmU1eu3ZtWYMFAABAc9UdDs3sByX9iaRfd/d/lDQq6UclbVGYWfz9al/n7ofdfcDdB1avXt2AIQMAAKBZ6gqHZnaHQjD8vLuPSZK7T7t70d1nJf2RpAeaN0wAAAC0Qj3Vyibps5LOuvsfJK73J542KOnlxg+veTgtBQAA4Hb1zBy+V9IvSXpfRduaT5jZt8zsJUk/I+k3mjnQRuO0FAAAgNuZu7fsmw0MDPjk5GTLvl81vb2hz2Glnh7pxo3WjwcAACANZvacuw9UXu+KE1KSOC0FAACgtq4Lh5yWAgAAUFvXhUOpfFrKxES4pygFAAAgqOf4vMzhtBQAAIDqunLmEAAAANURDgEAABAjHAIAACBGOAQAAECMcAgAAIAY4RAAAAAxwiEAAABihEMAAADECIcAAACIEQ4BAAAQIxwCAAAgRjgEAABAjHAIAACAGOEQAAAAMcKhpKkpaetW6cqVtEcCAACQLsKhpJERaXxcOngw7ZEAAACkq6vDYW+vZCaNjkqzs+HeLFwHAADoRl0dDgsFac8eKZ8PH+fz0t690sWL6Y4LAAAgLV0dDvv7pb4+aWZG6ukJ93190tq1aY8MAAAgHV0dDiVpelrav1+amAj3FKUAAIButjLtAaRtbKz8+NCh9MYBAADQDrp+5hAAAABlhEMAAADECIcVaIgNAAC6GeGwAg2xAQBANyMcltAQGwAAgHAYoyE2AAAA4TBGQ2wAAADC4Rw0xAYAAN2u65tgJ1U2xI4ql48dYwYRAAB0B2YO50HlMgAA6DaEwypqVS739KQ9MgAAgOYiHFZRWbmcy4X73bvTGxMAAEArEA6riCqXr18PHxeL4f7IEXofAgCAbCMc1jA9Le3bJ+3cKa0sle3k89KuXdKWLVQyAwCAbCIc1jA2Jj31lHTffWHfYdT78Nw56dSpUKTCOcwAACBrCIcLiHofzs6G25kz5SKVu++WTpygmhkAAGQH4XABY2Oh5+GlS3OLVJKiauYVK5hFBAAAnY1wWKfk8XqrVoVryb2IGzeGxyw3AwCATkY4XIRoifnkSen++6Vbt8L169el8+cld5abAQBAZzN3b9k3GxgY8MnJyZZ9v2batSvMJg4OSo8+GnojRi1vKq1aFWYcAQAA2oWZPefuA5XXOVt5iZLnMG/bJl24EELgzZthufnWrdA8u1ikeTYAAOgcLCs3QK3lZppnAwCATsPMYQMkZxE3bZIGBqSrV6Xjx0NQzOelHTuky5dDkcratemNFQAAYD7MHDZYPc2zAQAA2tWC4dDM7jWzb5rZWTM7Y2a/Vrr+VjM7bmbnS/dvaf5wO8d8zbPpiQgAANpVPTOHtyT9pru/Q9J7JP2Kmb1T0kclfd3dN0r6euljlNRqnl3ZExEAAKCdLBgO3X3K3Z8vPf6+pLOS7pH0kKQjpacdkfSBZg2ykyWbZ0u390RkFhEAALSTRe05NLP1kt4t6aSkNe4+JYUAKemHa3zNkJlNmtnktWvXljfaDhUtMR8/HmYNc7lwnZNVAABAu6m7CbaZ/aCkP5f0u+4+Zmbfc/e7Ep//B3efd99hlppgL9XwsHT4cNh/ON9zPvWp1o0JAAB0n1pNsOuaOTSzOyT9iaTPu3vUuGXazPpLn++XdLVRg82yWrOISSw3AwCAtNRTrWySPivprLv/QeJTz0raV3q8T9KXGz+87IkKVbZvDyeruIeTVaRwsopE0QoAAEhPPTOH75X0S5LeZ2Yvlm4PSvq4pPeb2XlJ7y99jEWodbJKtaKVnp50xwoAALrDgiekuPu4JKvx6W2NHU53qTxZZetWaXBQevRRqVAIx+9xPjMAAGgljs9rE8mguG1bmDmU5p7PHJ3RfPkyR/ABAIDm4Pi8NjQ9Le3bJ+3cyT5EAADQWoTDNlR5PrNE82wAANAahMM2Vm/zbAAAgEapuwl2I9AEe+kWap7NXkQAALAYy2qCjfRxBB8AAGgFZg47EEfwAQCA5WLmMEMWcwRfTw+ziQAAoH70OexAlT0RL1wIR/DdvBla39y6Nbd59siIND4elp2ZTQQAAPNh5rDD1TqCL9k8e3Q0LEHTAgcAACyEcNjhxsakQ4ekzZvDEXyVzbNzufqKV1h6BgAAEgUpmRQVrNx5pzQzU9/zJenJJ6VHHmHpGQCAbkBBSheJlponJqQNG8JtoeKVyqXnnp7WjxsAAKSPgpQMShasFArlx/MVr0hhn2KykAUAAHQfZg67SLXilZ6eEAajApZkIQvFKwAAdB/CYRepLF559NHy0vO9984tZOHkFQAAuhMFKYhx8goAAN2DghQsaDEnr7DcDABANhEOEYuWnbdvD8Ur7qF4RWK5GQCAbkE4RFW1Tl65fl06fz4Ex9FR6e67pRMnpMcfJyQCAJAFtLJBVcl2OJs2heA3OBiKWAqFclVz5OjRcL9unfSd74RWOMeOSWvXtm7MAABg+QiHWFAyKFb2SqxULIbZRCksO1O8AgBAZ2FZGYtSudwsheKUaqLild7e1o0PAAAsDzOHWJRqy81Xr0pf+lIIgu7lk1fyeWnHDuny5bAXkSVmAADaHzOHWLKourlYDHsRX3hh7skrMzPSuXPSqVNUNgMA0Clogo2G2rVL6u+XPvMZ6c03qz9neFg6cICiFQAA0lSrCTbhEE0xNSU99pj0zDOh/U0tDz8sXbpESAQAoNU4IQUt1d8v9fWFpeXKRtpJR4+GPonr1s1ddmYJGgCAdBAO0TTVGmlHQbFS1ALnxImwP3FkRBofZ68iAACtxrIyWiLaizg0JO3dK505E1rgzM7W/28MD9M3EQCARmHPIdpGFBRrtcDJ5cLzKk9hiZiF9jjsUQQAYOnYc4i2sVALnGIx3Cr3Kubz0saN4fHBg+mMHQCArCMcIjVRSNy8OTTUfvRRaWJC2rAh3JJ7FaVQ9Xz+fJhljE5f6emhkAUAgEbihBS0heTJK4VC+XF0CsvgYAiPhUKYVczlwv3u3XOLV6TyY/YnAgCweOw5RMcYHpY+/enFfQ37EwEAqI49h+h409PSvn3Szp3lfYi5XLmAJfm4cn8iy80AANSHcIiOMTYmPfWUdN99oQVOsngl+Vi6fX9i1EPx8ccJiQAAzIdwiI4TNddOFq8kHx8/HmYNo1nEpFonsgAAgICCFHScWsUrycfbtkkXLoR2ODdv3v5vRCeySGHZ+cCBUNzCGc8AgG7HzCEyqfLoPimcyFJNctmZ/YkAgG5HtTIyb6ETWWrhuD4AQJZRrYyuNd+JLJWnsCRFjbZXrJBOn2Y2EQDQHQiH6BrVTmRJnsIy33F9e/eWm2uz7AwAyDKWldH1omXnoaEQAs+cqe/rhocpZAEAdK5ay8qEQyAhCorRcX0XLoSeigthfyIAoNPUCoe0sgESkm1yonY40TnOK1aEoFitkGV0tLxHkeP6AACdjD2HQA1RO5ytW8O+xJ/+6fr2J7IvEQDQyVhWBhZhsfsTH35YunQp7El0Z38iAKB9sKwMNEBy2XnTpjA7GO1PLBTKZztHjh4N9+vWhUAZVTyzPxEA0K4WXFY2s8+Z2VUzezlx7WNm9rqZvVi6PdjcYQLtJ2qNs3172J/oXl5urlQshj2Js7Nz+ydeucISNACgvdSz5/ApSTuqXP+ku28p3b7a2GEBnWWh4/pyuXCTbt+fODJSnlEEACBtde05NLP1kr7i7u8qffwxSW+4+xOL+WbsOUQ3qDyuL6pyrhcVzwCAVmjG8XkfMbOXSsvOb5nnGw+Z2aSZTV67dm0Z3w7oDJXH9T3/vLRhQ7gdPx5mDaNZxPlmFFluBgCkYakzh2skfVeSSxqR1O/uv7zQv8PMIRAaZh8+LN15pzQzs/DzqXgGADRDQ2cO3X3a3YvuPivpjyQ9sNwBAt0i2p84MVF7RjHp6FHpxIlQ8Zzcn8jMIgCgGZYUDs2sP/HhoKSXaz0XwFzRsvPmzaH9TaGwtIrnu+8OoZGgCABopHpa2XxB0l9K+jEze83MPizpE2b2LTN7SdLPSPqNJo8TyLzFVDwnERQBAI3ECSlAG5qv4nnVKunmzepnPCcl9yqyPxEAUKkZ1coAmmS+iudoZrHaGc9Jyb2KydlEZhYBAPPh+DygjSWP6ysUyo+jo/uSZzxHM4qVisWw7CyVG21HRS0HDlD9DACYi2VloMNFS9DJoLjYxtvDwwRFAOg2tZaVCYdAhlTuVTQLFdDR/sSooKVYrP1vsFcRALoD4RDoItVmE3t6yk236ylqyeXmL3gBAHQ2ClKALpLspbhpUyhqSTbdrlbUUqlYDDOPK1ZQvAIA3YRwCGRctabbydBYq68iZz0DQHdiWRnAbXsV5zM8LH3qU60ZFwCgeVhWBlBTZV/F+c56Hh1luRkAsoxwCCAWhcRqZz1HjbZZbgaAbCMcAqiq8qznqHL5+nXp/PkQHJPnOj/+OCERALKAcAigqmoVz/MtN3NcHwBkA8fnAVhQ8hi/bdukCxeWdlwfhSwA0P6YOQSwKJXLzdLcFjhJo6PhNjt7eyELM4oA0J6YOQSwKMlZxE2bQsCr57i+fF665x7p1VeZUQSAdkafQwDLttBxffVYtUq6eFHavZtznQGgFehzCKBpFjqur7KQJZeb+1gKoXBkpDybSFELAKSDmUMALTE8LB0+LN155+JmFIeHw/2TT0qPPCIdOMDsIgA0Qq2ZQ8IhgJZILj0PDoY9ie96V5hVrNyfWI+HH5YuXSIkAsBS1QqHFKQAaIlkIUuhEO6Hh0Mlc3J/YtQip1pRS9LRo+F+3TrpO99hNhEAGoU9hwBSE7XFSe5PTJ7I0tMTQmGxWD7Gr1LUV/HEiXIVNABg6VhWBtB2KpegJenpp8uV0CtWhBnHWsyky5eZRQSA+bCsDKBjVFuClhbuq1jZS5HiFQBYPJaVAXSMqGVOsRja5bzwQnkJWpKuX5fOnw+BcXR07nJz1A7n9Gla5ADAfFhWBtDRoiXowcEQGAuF+Sue779fOns2tMWRqrfIcWfGEUD20coGQOZFvRTvuGNuxfNi/w2pHBo52g9AVnFCCoDMi6qfkxXPUZVz1BInkjylJWl0NNxmZ8O9WSiAYekZQLcgHALIjGrH+EVBsVice2RfZYuclSvLn4uel8+HY/8k6fHH2Z8IoDuwrAwg86J9ia+8EmYX16yRLl4Mn0u2yEk2464ll6PpNoBsoJUNgK6VbI1TTdQiJ9lX8fDhMPN4/vzc50ZNtyXa5QDIJmYOAaCG4eFQmGI2f9Pt6LkERQCdhIIUAFik6ekQ+nbtCh+bhfuVVdZcavVVZI8igE5DOASAGuZrul1ZyJKUDIoUsgDoNOw5BIAFJPcsJvcnRoUsq1bV7qt49Gi4X7du8T0XASANhEMAWIR6g2KlYjEsS/f0SDdutG68ALBYhEMAWKKFguKKFeVClnxe2rFDunw5LDFTsAKgXbHnEAAaoFoD7qiQZcWK0D/x3Dnp1CkKVgC0N1rZAECTRM23P/MZ6c03qz9neJjzmwGkg1Y2ANBi0WzipUvSnj1hablSdH5zb2/LhwcAVREOAaDJ+vulvr6wtFzZAiefDzOMW7aEJeZoufn0aZadAaSDcAgALTA9Le3fL508We6VGJ3lnNyLODIijY+Hopbx8XANAFqJPYcA0GL17EWsZBYqnd3LR/QlH1P9DGCx2HMIAG2i1l7EXC5UNif19kobN4bHyZnFysdUPwNoFGYOASBFw8PS4cPSnXeGJWYphMRicen/HtXPAOrBzCEAtKFoL+LEhLRhQ7ht3Rr2Jb73vWHWMJcLz83lqj9OiqqfV6yYW9TCzCKAenFCCgCkKHnKSqFw++eHh6ULF8rFK9Lcx5XnOufz0j33SK++Gopazp4tF7VES9AHDrBXEUBtLCsDQBuLileGhqTBwXDt6afnPo6O61sslqCB7lZrWZlwCAAdLgqQg4Ph2L4LF8pnOkvl5eda+xirVUIzowhkX61wuOCyspl9TtLPSbrq7u8qXXurpGOS1ku6JOkX3f0fGjlgAEB9kkvT27aFcBgVtSSLW+Zbgq5cemZGEeheC84cmtlPSXpD0tFEOPyEpL9394+b2UclvcXdH1/omzFzCADNFc0ivvJKKHZZs0a6eDF8brFL0KtWlfc2AsieJc8cuvsJM1tfcfkhST9denxE0p9JWjAcAgCaKzmLWM2mTaFqOVqCLhTKM4zS3NnG3btDlTNLzUB3WWormzXuPiVJpfsfbtyQAADNEjXg3r49LEG7h+rnYrG8/BzdHzki3X23dOIEjbaBbtL0PodmNmRmk2Y2ee3atWZ/OwBAnSp7LN57r7RzZ9iTWGl0lKAIdIulhsNpM+uXpNL91VpPdPfD7j7g7gOrV69e4rcDADRaNIu4eXNYXv7bv5Xuuy9UOq9aFZ5TT1AEkC1LDYfPStpXerxP0pcbMxwAQJqi2cSTJ8MpLbduLRwUzcIZ0ACyYcFwaGZfkPSXkn7MzF4zsw9L+rik95vZeUnvL30MAOhwydnETZtC0cp8QTGfDxXSW7ZwXB+QFfVUK3+oxqe2NXgsAIA2kqx8jqqch4bK7XCiY/zOnQvH9NU6ro+eiUBn4YQUAMCiRL0UP/MZ6c036/uanh7pxo3mjgvA4tTqc9j0amUAQLZES8+XLkl79ty+3zCXK/dNTC47s8QMdAbCIQBgSfr7pb6+cCRfFAajBtrF4txl51Onbm+BEz1O7lUEkD6WlQEASzbfcX2vv1572Xl4ONw/+aT0jneEvYqPPML+RKCVai0rEw4BAE0xNSU99pj0zDPS9ev1f52ZdPkyx/UBzcaeQwBAS0XLzjMzt7fAyeWkFRV/gXp7pY0bw+NaS9AsPQPNRzgEADRNtaba0VnOs7PlvYpSqGY+fz6c91x5CsvISLk1DoDmYlkZANAS0f7EoSFpcDBc27AhBMi77pKuXg3H+BWLC/9bZtILL0i/+qvSsWMsQQNLUWtZecEm2AAANEKyqXahcPvnh4elCxfCEvTNm2EJ+tat8uxisRha49xzj/Tqq3ObblPIAjQOy8oAgLYw3xJ0NJt4/Xp56fnMmbA0zfnOQGMxcwgAaAu1juuLlqAPHw5nPV+4EEJhpLdX2rkzVDhfucISM7BczBwCANpOdArL5s1hCbpQkLZvl7ZtC59PNt2+eXPhRttUOQP1IxwCADpGtPS8dWtYenYPs4jJJeZaVc4ERaA+VCsDADrWUhptDw9TwAJINMEGAGTQQo22k30UI1EBS09P68YJdBLCIQCgoy1U5VwtNErS7t3pjBdod4RDAEBHSxavbNoUKponJkKD7Q0b5oZGqdwW58iRMIO4YgX7EIEkWtkAADKjVqPtTZukgYFwCsvx4yEoJhtq00gbKGPmEACQeWNj0lNPSffdV+6RmGyoHe1DZBYRIBwCALpItD/x+HFp48by/sN8PnwshVlEoJvRygYA0JWGh8OpK8nTVpLMwqkrnLiCrKKVDQAACfXOIkbNs0+fpok2ugMzhwCArrfQLKIUKp7PnpUeeUQ6cCC0wjl2jJlFdC5mDgEAqKHWLGJSrSP6OJYPWUMrGwBA10u2wNm2TbpwITTPvnkzBMWoN2Kl0dFwk2iHg+xg5hAAgITKE1eKxfJM4orSX82VVaZWonY4vb2tGyvQDMwcAgCQkJxF3LQpLBm/8koIjWvWhPszZ8oziytXlptqDw5KTzyR3tiBRmDmEACAGqKj+b7+denll8N9dERf5VnOMzNhhvGDHwz7D6lyRqdi5hAAgEWoNrM4NBSqnb/2Nenb3y430h4fl/buDVXO7ElEp6CVDQAAy9TbG2YO69HTI9240dzxAPWglQ0AAE1SKEh79oR9h1JYXl5R8Re2tzfMIl682PrxAYtBOAQAYJn6+6W+vjB72NMTKpxnZ8tVzrlcKF7p66NpNtofew4BAGiAqAXO0FCoWpakDRvKVc5vf3soUgHaHeEQAIAGSBaqFArpjQNYLpaVAQBokeRRexy7h3bFzCEAAC0yMhLa2yRb3dDiBu2GVjYAADRZPa1uaHGDVqOVDQAAKanW6iaqZM7npV27pC1bWGJGeyAcAgDQZNVa3RSL5WP3zp2TTp0qLzcDaSIcAgDQAlGrm4mJ0OJmw4bQC3F2VjpzJtyPjkpmoYE2s4hIC+EQAIAWGBuTDh2SNm8Oy8yFgnTp0tzl5nxe2rgxPGYWEWkhHAIAkJLkcrMkXb8unT8vuZdnEXt60h0jug/hEACAFEXLzcePh1nD5JF7krR7N/0R0Vr0OQQAIEXJk1W2bQszh1IoWJGkI0fCTaI/IlqDPocAALSJXbvCMvPVq2Em8dathb+G/ohYqlp9Dpk5BACgTUSziMPDoXp51Srp5k1p5coQFKOl5mIxFK/s2CFdvhyWmPUgaOoAAAzoSURBVNeuTW/cyBb2HAIA0GaifYgnT0r33x+CIf0R0SosKwMA0MZ27QpVzUND0uBguPb669Kbb97+XJaYsRgsKwMA0IGSBSuFQrifmpIee0x65pnQ/iafD8HxiSfSGSOyZVnh0MwuSfq+pKKkW9XSJwAAaKzK4/hmZsLH7DtEIzRiz+HPuPsWgiEAAK2TPI5v/376HqJxWFYGAKADJZebDx1KbxzInuXOHLqkPzWz58xsqBEDAgAAQHqWGw7f6+4/LmmnpF8xs5+qfIKZDZnZpJlNXrt2bZnfDgAAVMOxemiUZYVDd79cur8q6WlJD1R5zmF3H3D3gdWrVy/n2wEAgBpGRsrH6gHLseRwaGY/YGY/FD2W9LOSXm7UwAAAwMJ6eyUzaXQ0nKoyOho+7u1Ne2ToVMuZOVwjadzMTks6Jel/ufv/bsywAABAPQoFac+e0OtQCvd790oXL6Y7LnSuJVcru3tB0uYGjgUAACxSZc/DGzekb3wj7VGhk3G2MgAAHS7Z8/Cd7wzFKew9xFJxtjIAABnQ2xtmDyuZSZcvc3oKblfrbGVmDgEAyIBqew83bgyPmUXEYhAOAQDIgOTeQ0m6fl06f15yp4IZi0M4BAAgI6K9h8ePh1nDXC5cp4IZi8HZygAAZETyvOVt26QLF0IF88xMmFVk3yHqwcwhAAAZlKxg3r+fY/VQP2YOAQDIoOQs4qFD5bOXjx1jBhHzY+YQAIAukDx7OQqKV65Uf3z6dPkaug99DgEAyLBa/Q8laXg43D/5pPTII+XH73iHdPZsuHbggLR7NzOOWVSrzyHhEACADJuakh57THrmmdDeZqmGhwmKWUMTbAAAulCy/+GqVeHaylLFQS5XbneTy0kr5kkFo6PS3XdLJ07QVDvrCIcAAGRcVLl88qR0//3SrVuhxU2xGG7R49nZcliMguLKKqWrNNXONqqVAQDIuGTl8qZNodhkaEgaHAzXnn66/HjDhhAm16wJ92fOhBnHmzdDULx1KzTVHhyUnnii9f9b0HyEQwAAukgyKBYK1R9Hdu0qB8m9e0NQpKl29hEOAQBAVbVmHA8fDkfx0Tcxm9hzCAAAFjQ2Fpppb94c7tevL/dNRLYQDgEAQN16e0MxyuhoKGChOCV7CIcAAKBuhYK0Z08oSpHC/d69YZkZ2UA4BAAAdUv2TaQ4JZsIhwAAYFGivokTE+GeM5izhWplAACwKMkq5kOH0hsHmoOZQwAAsCxTU6GtDTOI2UA4BAAAyzIyQlubLCEcAgCAJanV1qanJ+2RYTkIhwAAYEkq29rkcuF+9+70xoTlIxwCAIAlidraXL8ePi4Ww/2RIzTG7mSEQwAAsGTT09K+fdLOndLKUg8UGmN3NsIhAABYsrEx6amnpPvuC/sOo8bYuZz0wQ9SwdyJCIcAAGDZKhtj/8VfUMHcqczdW/bNBgYGfHJysmXfDwAAtFZvb5g5rNTTI9240frxoDYze87dByqvM3MIAAAaprKCOZ+Xdu2StmxhiblTEA4BAEDDRBXMMzPl/YfnzkmnTrHE3CkIhwAAoKGi/Yezs+F25szcJtm0uGlvhEMAANBQY2PSoUPSpUu3LzHT4qb9EQ4BAEBTVFti7uuT1q5Ne2SYD+EQAAA0TWWLG4pS2t/KtAcAAACya2ys/PjQofTGgfoxcwgAAIAY4RAAAAAxwiEAAABihEMAAADECIcAAACIEQ4BAAAQIxwCAAAgRjgEAABAjHAIAACAGOEQAAAAMcIhAAAAYssKh2a2w8zOmdmrZvbRRg0KAAAA6VhyODSznKRDknZKeqekD5nZOxs1MAAAkD1TU9LWrdKVK+XHp0/ffq2ex1n9urStXMbXPiDpVXcvSJKZfVHSQ5L+uhEDAwAA2TMyIo2PSwcPho/Hx6W9e6WzZ+deq+dxVr/uU59a/OvaSObuS/tCs1+QtMPd/13p41+S9BPu/pFaXzMwMOCTk5NL+n4AAKBz9fZKMzNpj6Jz9PRIN24093uY2XPuPlB5fTl7Dq3KtduSppkNmdmkmU1eu3ZtGd8OAAB0qkJB2rNHyufDx7mctKIiheRy4bbQ4yx/XT4fZhovXlRqlhMOX5N0b+LjdZIuVz7J3Q+7+4C7D6xevXoZ3w4AAHSq/n6pry/MHvb0SMWiNDs7NyAVi+EWfb7W4yx/3cxMeJ3Wrm3t/z9JywmHfyVpo5ltMLM7Je2W9GxjhgUAALJmelrav1+amJA2bAi3rVul++8P99G15OdrPc7q1+3fn35RypL3HEqSmT0o6Q8l5SR9zt1/d77ns+cQAACgPdTac7icamW5+1clfXU5/wYAAADaByekAAAAIEY4BAAAQIxwCAAAgBjhEAAAADHCIQAAAGKEQwAAAMQIhwAAAIgRDgEAABAjHAIAACBGOAQAAECMcAgAAIAY4RAAAAAxwiEAAABihEMAAADECIcAAACImbu37puZXZP07ZZ9Q+ltkr7bwu/XKXhdquN1qY7XpTpel+p4XarjdamN16a6Vrwu97n76sqLLQ2HrWZmk+4+kPY42g2vS3W8LtXxulTH61Idr0t1vC618dpUl+brwrIyAAAAYoRDAAAAxLIeDg+nPYA2xetSHa9Ldbwu1fG6VMfrUh2vS228NtWl9rpkes8hAAAAFifrM4cAAABYhEyGQzPbYWbnzOxVM/to2uNJi5nda2bfNLOzZnbGzH6tdP1jZva6mb1Yuj2Y9ljTYGaXzOxbpddgsnTtrWZ23MzOl+7fkvY4W8nMfizxvnjRzP7RzH69G98zZvY5M7tqZi8nrlV9f1jwX0u/c14ysx9Pb+TNVeN1+c9m9krpf/vTZnZX6fp6M7uReN98Or2RN1eN16Xmz42Z/Vbp/XLOzP5lOqNuvhqvy7HEa3LJzF4sXe+m90utv89t8Tsmc8vKZpaT9DeS3i/pNUl/JelD7v7XqQ4sBWbWL6nf3Z83sx+S9JykD0j6RUlvuPsTqQ4wZWZ2SdKAu383ce0Tkv7e3T9e+g+Lt7j742mNMU2ln6XXJf2EpH+rLnvPmNlPSXpD0lF3f1fpWtX3R+mP/r+X9KDC6/Vf3P0n0hp7M9V4XX5W0jfc/ZaZ/Z4klV6X9ZK+Ej0vy2q8Lh9TlZ8bM3unpC9IekDS3ZL+j6RN7l5s6aBboNrrUvH535f0f939YJe9X2r9ff43aoPfMVmcOXxA0qvuXnD3NyV9UdJDKY8pFe4+5e7Plx5/X9JZSfekO6q295CkI6XHRxR+WLvVNkkX3L2VjevbhrufkPT3FZdrvT8eUvjj5+4+Iemu0i//zKn2urj7n7r7rdKHE5LWtXxgKavxfqnlIUlfdPeb7n5R0qsKf7syZ77XxcxMYbLiCy0dVBuY5+9zW/yOyWI4vEfSdxIfvyYCkUr/RfZuSSdLlz5Smpr+XLctnSa4pD81s+fMbKh0bY27T0nhh1fSD6c2uvTt1txf2rxnar8/+L1T9suSvpb4eIOZvWBmf25mP5nWoFJU7eeG90vwk5Km3f184lrXvV8q/j63xe+YLIZDq3ItW2vni2RmPyjpTyT9urv/o6RRST8qaYukKUm/n+Lw0vRed/9xSTsl/Upp+QOSzOxOST8v6X+WLvGemR+/dySZ2W9LuiXp86VLU5L+mbu/W9J/kPQ/zKwvrfGloNbPDe+X4EOa+x+gXfd+qfL3ueZTq1xr2nsmi+HwNUn3Jj5eJ+lySmNJnZndofDG+7y7j0mSu0+7e9HdZyX9kTK6nLEQd79cur8q6WmF12E6mqov3V9Nb4Sp2inpeXeflnjPJNR6f3T97x0z2yfp5yTt9dJm9tKy6d+VHj8n6YKkTemNsrXm+bnh/WK2UtIuSceia932fqn291lt8jsmi+HwryRtNLMNpdmP3ZKeTXlMqSjt5/ispLPu/geJ68l9CoOSXq782qwzsx8obQKWmf2ApJ9VeB2elbSv9LR9kr6czghTN+e/6HnPxGq9P56V9HCpovA9Chvsp9IYYBrMbIekxyX9vLtfT1xfXSpskpn9iKSNkgrpjLL15vm5eVbSbjNbZWYbFF6XU60eX8q2S3rF3V+LLnTT+6XW32e1y+8Yd8/cTaGa528U/qvjt9MeT4qvw79QmHZ+SdKLpduDkv67pG+Vrj+rUDGV+nhb/Nr8iKTTpduZ6H0i6Z9K+rqk86X7t6Y91hRem7ykv5P0TxLXuu49oxCOpyT9P4X/av9wrfeHwpLPodLvnG8pVMGn/r+hha/Lqwr7oaLfM58uPfdfl36+Tkt6XtK/Snv8LX5dav7cSPrt0vvlnKSdaY+/la9L6fpTkvZXPLeb3i+1/j63xe+YzLWyAQAAwNJlcVkZAAAAS0Q4BAAAQIxwCAAAgBjhEAAAADHCIQAAAGKEQwAAAMQIhwAAAIgRDgEAABD7/xEm6rIS5+U7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 792x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(200), sigma, 'b*')\n",
    "plt.gcf().set_size_inches(11, 8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also explore the way in which the topics are defined by considering the words for which they have stronger weights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort the list of words by the weight they have in a specific topic\n",
    "def top_words(vector, word_list):\n",
    "    doc_list = sorted(zip(word_list, vector), key=lambda x:x[1], reverse=True)\n",
    "    \n",
    "    return doc_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, for topic 5, we have strong mentions of AIDS day and the Produc (RED) iPhone launch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "together 0.2562979288804484\n",
      "ur 0.23461856127285463\n",
      "like 0.2139227443779649\n",
      "shit 0.17985072359458534\n",
      "ios 0.1636338889646671\n",
      "8 0.16013230849583537\n",
      "let 0.1400017310800703\n",
      "version 0.13320119712005074\n",
      "words 0.131666470491568\n",
      "give 0.1316664704915677\n"
     ]
    }
   ],
   "source": [
    "topic_words = top_words(uk[:,4], word_list)\n",
    "for word, weight in topic_words[:10]:\n",
    "    print(word, weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, if we were trying to map a specific sentence or document into the most relevant topic we would simply have to project it into the singular space. First we define the vector in the \"regular\" word space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "document = ['ipad', 'mini', 'price']\n",
    "\n",
    "doc_vector = np.zeros(vocabulary_size)\n",
    "for word in document:\n",
    "    doc_vector[word_dict[word]] += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and now we project it into the singular space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([ 1.62243400e-04,  1.44900491e-02, -1.09973213e-03,  4.30257980e-05,\n",
      "        9.69106965e-04, -4.54646796e-03,  7.58613138e-04, -2.40772326e-03,\n",
      "        3.82599561e-03, -1.95924079e-03])\n"
     ]
    }
   ],
   "source": [
    "doc_singular = 1/sigma_k*np.dot(uk.T, doc_vector)\n",
    "pprint(doc_singular)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From where we see that hte most relevant topic is topic number 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#macbook 0.3072355938213901\n",
      "@youtube 0.30540145699726673\n",
      "#startup 0.30540145699726673\n",
      "#hipster 0.30540145699726673\n",
      "#macbookpro 0.29809606289080237\n",
      "via 0.2918795721098773\n",
      "mini 0.2094966623380363\n",
      "ipad 0.20381766001212515\n",
      "#ipadmini 0.19891025664321696\n",
      "startup 0.19327488007614926\n"
     ]
    }
   ],
   "source": [
    "topic = np.argmax(doc_singular)\n",
    "\n",
    "topic_words = top_words(uk[:,topic], word_list)\n",
    "for word, weight in topic_words[:10]:\n",
    "    print(word, weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "topic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which indeed seems relevant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Non-Negative Matrix Factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea behind applying NMF to a Term-Document Matrix is similar to the goals behind using SVD to perform Latent Semantic Analysis. NMF is a general matrix factorization algorithm and, as such, the factor matrices it finds can also be thought of as internal or latent dimensions. \n",
    "\n",
    "Similarly to the case of SVD, we can perform a partial reconstruction of the original matrix by reducing the internal dimensions. Naturally, as the number of latent dimensions increases, the reconstruction error decreases since we are removing less and less information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NMF(V, k):\n",
    "    n, m = V.shape\n",
    "\n",
    "    W = np.random.rand(n, k)\n",
    "    H = np.random.rand(k, m)\n",
    "\n",
    "    error = 1\n",
    "    err = 1\n",
    "\n",
    "    while err > 1e-3:\n",
    "        old_error = error\n",
    "\n",
    "        hn = np.dot(W.T, V)\n",
    "        hd = np.dot(W.T, np.dot(W, H))\n",
    "\n",
    "        H = H*hn/hd\n",
    "\n",
    "        wn = np.dot(V, H.T)\n",
    "        wd = np.dot(W, np.dot(H, H.T))\n",
    "\n",
    "        W = W*wn/wd\n",
    "\n",
    "        error = np.sum(np.power(V-np.dot(W, H), 2.0))\n",
    "        err = np.abs(old_error-error)\n",
    "\n",
    "    return W, H, error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how the error changes with the number of components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 32018.60718378711\n",
      "3 31373.949867759144\n",
      "4 30774.367750722042\n",
      "5 30270.93203672034\n",
      "6 29820.781934073224\n",
      "7 29342.239005796724\n",
      "8 28901.919381172876\n",
      "9 28480.27493196867\n",
      "10 28042.691136101912\n",
      "11 27603.713891196814\n",
      "12 27234.76434115258\n",
      "13 26836.858225522406\n",
      "14 26407.808172868652\n",
      "15 26113.311406986068\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/bgoncalves/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:21: RuntimeWarning: divide by zero encountered in true_divide\n",
      "/Users/bgoncalves/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:21: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16 nan\n",
      "17 25271.928994232498\n",
      "18 24949.893217070265\n",
      "19 24530.746318234076\n",
      "20 nan\n",
      "21 23868.032504894025\n",
      "22 23482.011785337152\n",
      "23 23068.251083740837\n",
      "24 22785.393946004777\n",
      "25 22472.593647756155\n",
      "26 22088.90583026834\n",
      "27 21790.763445551387\n",
      "28 21522.41082872108\n",
      "29 21161.84906238318\n",
      "30 20770.869930907265\n",
      "31 20528.728198808945\n",
      "32 20195.26835168722\n",
      "33 19991.502745198934\n",
      "34 19611.52410313329\n",
      "35 19298.706206915\n"
     ]
    }
   ],
   "source": [
    "error = []\n",
    "\n",
    "for comps in range(2, 35+1):\n",
    "    W, H, err = NMF(TDM, comps)\n",
    "    print(comps, err)\n",
    "    error.append(err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nan values found above are due to numerical instabilities. In practice you should use a well tested numerical library such as sklearn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Non-Negative Matrix Factorization')"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xUVfrH8c83IaFjKAHpIKIIKiARUBGxIYsF/VnWjm1ZsZfdtWyxbVF3bdh7RdEVCyqKrAsCShGQ3qVIE0Gq9PL8/rgn6xhDMoRMZpI879drXrlzbnvuTTLP3HPuPUdmhnPOOVcUackOwDnnXOnlScQ551yReRJxzjlXZJ5EnHPOFZknEeecc0XmScQ551yReRJx5Y6kpyT9Odlx7A1JF0j6NNlxFBdJ3SRNT8B2e0v6uLi3634if06k/JG0EKgM7GdmG0PZFcCFZtYtAfsbDnQGWprZ4lB2AvCcmTUr7v3l2fclwBVm1iWR+wn7Wgg0ABqY2aqY8klAW6C5mS0sZBvNgAVAhpntSFCcBmwCcv/5d5hZ1l5sr0R+l3HEsT8w18yUzDjKG78SKb8qANeX4P42AqX623+cFgDn5b6RdAhRwi42kioUw2bamlm18CpyAikOxXQ8Lkk8iZRf/wR+JynfDxBJR0r6StK68PPImHnDJd0j6QtJGyR9KqlOIfvrB5wXvi3mt78GkgZKWilpgaTrYuZVlvSypDWSZkr6g6QlMfNvlfRNiGWGpDNC+UHAU8ARkn6UtDaUvyTpr2F6pqRTYrZVQdIqSYeF950lfSlpraTJkroVcpyvAhfHvO8NvJLnWE+W9LWk9ZIWS7ozZvaI8HNtiPkISZeEc/2QpNXAnaFsVNjekSHmxuF92xBvq0Ji/RlJtSUNDr+DNZI+kNQwz/yXJC0P8wdK2gf4AGgS4v1RUl1JlST1C8sulfSgpMywnRMkLZR0u6TvgGdzy8L8C2K29aOkrZL+E+adJmlS+F1/m6dackRYJne9wyVdEa6Ec4+hi6Tx4e96nKROMfNGSbor/L43SPpEUq09OYflkSeR8ms8MBz4Xd4Z4R/nI6IP/trAg8BHkmrHLHY+cClQF8jMbzt5LAWeBe7MZ39pRB9Ek4GGwPHADZJOCovcATQD9gNOBC7Ms4lvgKOBfYC7gNck1TezmcCVwOgCvnG/QcyVA3ASsMrMJoYP0I+AvwK1wjEOlJRdwHGOAWpIOkhSOvBr4LU8y2wkSjRZwMlAX0mnh3ldw8+sEPPo8L4TMJ/ofP8tdmNm9iXwNPCypMpEiexPZjargDjzk0b0O2oCNAW2A4/EzH+d6HfdGqgHPGJm64BTgW9jrmy+B/4C5ACHAu2Bo4DbYrbVCKgW9nVVnuPpn7utsNxCot8TwI9Ev/99wn6vj/kS0DWsnxvHV7HbDV90PgIeIPq77gcMllQzZrHziRJ/PaAqcFPhp6188yRSvv0FuDafD8WTieqWXzWzHWb2BjCL6J8214tmNsfMNgNvAe3i2N8/gFMltclTfjiQbWZ3m9k2M5tP9GF2bph/DvB3M1tjZkuI/vn/x8z+bWbLzGyXmb0JzAU6xhEPRB+Mp0mqEt6fH8og+rAabGaDw7aHEiXfnoVsM/dq5ESi87Y0T7zDzWxq2OYUog/IYwrZ5jIzezT8PjbnM/9Oog/WccAy4PFCtjcxXK2sldQvxLXSzN41s81mth74e25c4SrneKBv+D1sM7MRu988FwB3hm1+D9wNXBQzf0eYv203x5P75eIN4FMzez7E+F8zmxbO3WRgAIWfu1ynAtPN7I1wHl8jSswnxyzzvJnNNbNNwL+J7++6XPMkUo6Z2TTgQ+DWPLMaAIvylC0iukrI9V3M9Caib5W5dz7lVifcnmd/K4HHiD5QYjUFGsR8qK0Fbif6Npgbz+KY5WOnkXRxqOLIXfdgoLDqtdyY5gEziZJbFeA0fkoiTYGz88TVBahfyGZfJUpGl5CnKivE20nSsFBttI7oaqmweBcXNNPMtgMvER37A1b4HTOHmVlWeF0X4qoq6blQTbQe+G9MXI2JrtDWFbLdXPX5+d9Q3r+fFWa2rZBt3AdUBG7MLQjVe8Njzt0VxPm7Zi/+rt3ueRJxdwC/4ef/SMuIPkBjNSHPN+r8mNmVMdUJf89nkX8CxwIdYsoWAwtiPtSyzKy6meV+419OVK2Rq3HuhKSmRFct1wC1Q5XVNCD3Dp14bj/MrdLqBcwIiSU3rlfzxFXVzO4t5BwsImpg7wm8k88irwODgMZmtg9Ru01h8RZ4HKHq7Q7gReABSRULWn43/gA0BzqaWQ3guJh5i4E6kmrEGdtyfv43lPfvp7DjuQA4Ezg7z11qA4CB/HTuniP+33WR/67d7nkSKefCB+abwHUxxYOBAySdr6ih+ddE9eAfFsP+1hLVSf8hpngcsF7SLYoa0dMlHSzp8DD/LeA2STXDh+U1MetWJfrwWAkg6VKib+O5VgCNcht1d2MA0B3oy09XIRC1ZZwq6aQQUyVFzzM0yncrP3c5cFzuLdR5VAdWm9kWSR2JrlpyrQR2EbX/xEWSiK5Cng/7XQ7cE+/6eeLaBKwJ7V9/yZ0Rbs3+D/C4pCxJGZJy229WECWY6jHbegP4i6Q6obr0z/yybWh3x5MDPAz0MrMf8okx99x15qcqT4DvAZO0u3P3IdBG0q/D3/X5wP5Ef++uiDyJOIiql6rmvgn/uKcANwM/EH3gnxL77MNeegTYGbO/nUT11e2IvsGvIvqGuU9MfEvCvP8AbwNbw7oziJLSaKIPs0OAL2L29V9gOvCdpHzjN7PlYf0jiRJqbvlioquT24k+3BcDvyeO/xsz+8bMxu9m9lXA3ZI2EH1QvxWz3iaihvMvQhVa58L2RfQFoB7w51CNdSlwqaSj41g31oNE5/wH4Esg70N6uTc0zCE619eGmKcRXR0sDDHXJbrBYTIwFZgCjCVqE4vH6UBNYHRM1egHYV5f4B/h3N3Oz8/dhrCPsSGOnNiNhurU04BbwjHeSPR3vTrOuFw+/GFDV+pI6guca2bxNqg65xLEr0RcypNUX9JRktIkHUh0hfRusuNyzkVPLTuX6jKJnoNoDqwlasN4IqkROecAr85yzjm3F7w6yznnXJGVu+qsOnXqWLNmzZIdhnPOlSoTJkxYZWa/6PKn3CWRZs2aMX787u68dM45lx9JeZ/2B7w6yznn3F7wJOKcc67IPIk455wrMk8izjnnisyTiHPOuSJLWBIJPZ6OUzSk6HRJd4Xy/pJmS5om6QVJGaFciobTnCdpisLwpGFeb0lzw6t3THkHSVPDOv1Cb6bOOedKSCKvRLYSdYXdlqh31h6hR9L+QCui3lYrEw0qA/AroGV49QGehP8N1XoH0fCgHYE7YoazfDIsm7tejwQej3POuTwSlkQs8mN4mxFeFoYatdBl9Th+GmyoF/BKmDUGyJJUn2jM66FmttrM1gBDiRJSfaCGmY0O23qFqAvphHhj3Ld8PmdlojbvnHOlUkLbRMJAPpOIBosZamZjY+ZlEI25/EkoasjPhwBdEsoKKl+ST3mx27ZjF6+OXsSVr05g4rdrErEL55wrlRKaRMxsp5m1I7ra6CgpdsS5J4ARZjYyvM+vPcOKUP4LkvpIGi9p/MqVe341kVkhjZcv60i9GhW59MWvmLNiwx5vwznnyqISuTsrDIk6nNBmIekOIBu4KWaxJcSMnU2UeJYVUt4on/L89v+MmeWYWU529i+6folLdvWKvHp5JypWSOOi58eyePWmIm3HOefKkkTenZUtKStMVwZOAGZJuoKoneM8M9sVs8og4OJwl1ZnYF0YtnQI0D2Mr12TaCzsIWHeBkmdw11ZFwPvJ+p4ABrXqsIrl3dk87adXPzCOFb9uDWRu3POuZSXyCuR+sAwSVOAr4jaRD4EniIaD3q0pEmS/hKWHwzMB+YBzxKNQ00Y//iesI2vgLtjxkTuSzQW9zzgG345JnSxa7VvDV689HCWr9tM7xfGsWHL9kTv0jnnUla5G5QqJyfHiqMX32Gzv+c3L4+nQ9OavHxZRyplpBdDdM45l5okTTCznLzl/sR6ER17YF0eOKctYxes5ro3vmbHzl2Fr+Scc2WMJ5G90KtdQ+48tTWfzljBbe9Mpbxd1TnnXLkblKq4XXJUc1Zv2k6/z+ZSq2omt/U8KNkhOedcifEkUgxuPKElazZu4+kR88mqkknfbi2SHZJzzpUITyLFQBJ3ntaGtZu3c98ns8iskMblXZonOyznnEs4TyLFJD1NPHhOW7bv2MU9H84gXVFVl3POlWXesF6MMtLTePT89nRvXY87P5jBK6MXJjsk55xLKE8ixSwjPY3Hzj+MEw6qx1/en85rYxYlOyTnnEsYTyIJkFkhjScuOIzjW9XlT+9N4/Wx3yY7JOecSwhPIgmSWSGNJy48jGMPzOb2d6cyYJwnEudc2eNJJIEqVkjnyQs7cMwB2dz27lTe+mpx4Ss551wp4kkkwSplpPP0RR3osn8dbnlnCm9PWFL4Ss45V0p4EikBlTLSefbiHLrsX4ffvz2ZdyZ6InHOlQ2eREpIpYx0nrkohyNb1Obmf0/muZHzva8t51yp50mkBFXOTOe5iw+nR5t9+etHM7n93Wls995/nXOlmCeRElY5M53Hzz+Mvt1a8Ma4b7n0xa9Yt9kHtnLOlU6eRJIgLU3c0qMV9591KGMX/MD/PfEFi37YmOywnHNuj3kSSaJzchrz6uWdWPXjNs544kvGL1xd+ErOOZdCPIkkWef9avPuVUeyT+UMzn92LO99vTTZITnnXNw8iaSA/bKr8e5VR9K+SRY3vDmJBz+d7XduOedKBU8iKSKrSiavXt6Jszs0ot9/53HtG1+zZfvOZIflnHMF8vFEUkhmhTTuP+tQ9suuxn2fzGLe9z/yxAWHsV92tWSH5pxz+UrYlYikSpLGSZosabqku0J5c0ljJc2V9KakzFBeMbyfF+Y3i9nWbaF8tqSTYsp7hLJ5km5N1LGUJEn07daCFy85nO/Wb+G0x77goynLkx2Wc87lK5HVWVuB48ysLdAO6CGpM3Af8JCZtQTWAJeH5S8H1pjZ/sBDYTkktQbOBdoAPYAnJKVLSgceB34FtAbOC8uWCce2qstH1x1Ny3rVuPr1idzx/jS27vDqLedcaklYErHIj+FtRngZcBzwdih/GTg9TPcK7wnzj5ekUD7AzLaa2QJgHtAxvOaZ2Xwz2wYMCMuWGQ2zKvNmnyO4vEtzXh69iHOeGs3i1ZuSHZZzzv1PQhvWwxXDJOB7YCjwDbDWzHaERZYADcN0Q2AxQJi/DqgdW55nnd2VlymZFdL48ymteerCDsxftZGT+43kPzNWJDss55wDEpxEzGynmbUDGhFdORyU32Lhp3Yzb0/Lf0FSH0njJY1fuXJl4YGnoB4H78uH13ahca0qXPHKeP4xeKb3u+WcS7oSucXXzNYCw4HOQJak3LvCGgHLwvQSoDFAmL8PsDq2PM86uyvPb//PmFmOmeVkZ2cXxyElRdPaVRnY90gu7NyEp0fM57xnxvDdui3JDss5V44l8u6sbElZYboycAIwExgGnBUW6w28H6YHhfeE+f+16Im7QcC54e6t5kBLYBzwFdAy3O2VSdT4PihRx5MqKmWk89fTD+GRc9sxY/l6Tnl0lHeX4pxLmkReidQHhkmaQvSBP9TMPgRuAW6SNI+ozeP5sPzzQO1QfhNwK4CZTQfeAmYAnwBXh2qyHcA1wBCi5PRWWLZc6NWuIe9dfRTVKqZz3rNjeH2sj+HunCt5Km/da+Tk5Nj48eOTHUaxWbdpO9cN+JrP56zkgk5NuOPUNmRW8I4InHPFS9IEM8vJW+6fNqXcPlUyeOGSw7nymBb0H/stFzw3hpUbtiY7LOdcOeFJpAxITxO3/qoV/c5rz9Sl6zj10VFMXrw22WE558oBTyJlyGltGzCw75Gkp4mznx7NwAlLkh2Sc66M8yRSxrRpsA8fXNuFDk1qcvO/J3P3BzPY4c+TOOcSxJNIGVSraiavXN6RS49qxgtfLODC58eyYr0/T+KcK36eRMqojPQ07ji1DQ+c3ZbJi9fR85GRfD6ndD6t75xLXZ5EyrgzOzRi0DVHUadaRXq/MI77Ppnl3aU454qNJ5FyoGW96rx39VGc17ExTw7/hnOfGcPStZuTHZZzrgzwJFJOVM5M5x//dyj9zmvP7O820PORkXw6/btkh+WcK+U8iZQzp7VtEHoDrkyfVydw1wfTfbAr51yReRIph5rViXoDvuTIZrz4xULOenI0i37YmOywnHOlkCeRcqpihXTuPK0NT1/UgUU/bOTkfqN4Z+ISyltfas65veNJpJw7qc2+DL7+aA6qX52b3prMdQMmsW7z9mSH5ZwrJTyJOBrVrMKAPkfwu+4HMHjqcno+MpKx839IdljOuVKgwCQiKU3StJIKxiVPepq45riWvH3lEVRIF+c+O4Z/DvFnSpxzBSswiZjZLmCypCYlFI9LsvZNavLRdUdzdodGPD7sG8568ksWrPJGd+dc/uKpzqoPTJf0maRBua9EB+aSp1rFCtx/VlueuOAwFv6wiZP7jeTNr771Rnfn3C9UiGOZuxIehUtJPQ+pT/smWdz05mRuGTiVYbNWcu+Zh5BVJTPZoTnnUkShVyJm9jkwC6geXjNDmSsH6u9Tmf5XdOK2X7Xis1kr6PnISMYvXJ3ssJxzKaLQJCLpHGAccDZwDjBW0lmJDsyljrQ08dtjWvD2lUdSIT2NXz8zhseHzWPXLq/ecq68U2H13JImAyea2ffhfTbwHzNrWwLxFbucnBwbP358ssMotdZv2c7t70zlwynL6bJ/HR78dVvqVq+U7LCccwkmaYKZ5eQtj6dhPS03gQQ/xLmeK4NqVMrg0fPac+//HcL4Ravp+chIRs71cUqcK6/iSQafSBoi6RJJlwAfAYMTG5ZLZZI4t2MTBl3ThVpVM7n4hXHc7+OUOFcuxdOw/nvgaeBQoC3wjJndUth6khpLGiZppqTpkq4P5e0kjZE0SdJ4SR1DuST1kzRP0hRJh8Vsq7ekueHVO6a8g6SpYZ1+krTnp8AV1QH1qvP+1V049/DGPBHGKVmyZlOyw3LOlaAC20QkpQNDzOyEPd6wVB+ob2YTJVUHJgCnAw8DD5nZx5J6An8ws25h+lqgJ9AJeMTMOkmqBYwHcgAL2+lgZmskjQOuB8YQXR31M7OPC4rL20QS44PJy7j9nalI8M+z23JSm32THZJzrhgVqU3EzHYCmyTts6c7NLPlZjYxTG8AZgINiRJBjbDYPsCyMN0LeMUiY4CskIhOAoaa2WozWwMMBXqEeTXMbLRFmfAVoiTlkuDUtg346LqjaVanKr99dQJ3vD+NLdt9nBLnyrp4HjbcAkyVNBT4X/8XZnZdvDuR1AxoD4wFbgCGSPoXURI7MizWEFgcs9qSUFZQ+ZJ8yvPbfx+gD0CTJt6DS6I0qV2Ft688kvs/mcVzoxYwftEaHj2vPftlV0t2aM65BImnYf0j4M/ACKKqpNxXXCRVAwYCN5jZeqAvcKOZNQZuBJ7PXTSf1a0I5b8sNHvGzHLMLCc7Ozve0F0RZFZI40+ntOb53jksXbuZUx8dxXtfL012WM65BCmsF990omdEXs77imfjkjKIEkh/M3snFPcGcqf/DXQM00uAxjGrNyKq6iqovFE+5S4FHH9QPT6+/mjaNNiHG96cxO//PZlN23YkOyznXDGLp00kW9Ied5YU7pR6nqiblAdjZi0DjgnTxwFzw/Qg4OJwl1ZnYJ2ZLQeGAN0l1ZRUE+hO1Ni/HNggqXPY18XA+3sap0uc+vtU5vXfdOLa4/bn7YlLOO2xL5j93YZkh+WcK0bxtIksBL4IPffGtok8uNs1IkcBFxG1p0wKZbcDvwEekVSBqL2lT5g3mOjOrHnAJuDSsJ/Vku4BvgrL3W1muZ039QVeAioDH4eXSyEV0tO4ufuBdN6vNtcPmMRpj43ijlPbcF7Hxvgd2c6VfvF0e3JHfuVmVip79/VbfJNn5Yat3PTWJEbOXUX31vW4q1cb6u9TOdlhOefisLtbfAtNIjEbqGpmpX50Ik8iybVrl/HsyPk8OHQOFdLEzd0PpPeRzUhP86sS51JZkfvOknSEpBlEz3kgqa2kJxIQoysHcnsEHnrjMeQ0q8XdH86g1+OjmLpkXbJDc84VQTy3+D5M9MDfDwBmNhnomsigXNnXpHYVXrr0cB47vz0r1m+l1+OjuHPQdDZs2Z7s0JxzeyCu3njNbHGeIn8U2e01SZxyaAM+u/kYLuzclJdHL+SEBz/n46nLfShe50qJeJLIYklHAiYpU9LvCFVbzhWHGpUyuLvXwbx71VHUrlqRvv0ncvnL41m82jtzdC7VxZNErgSu5qduRtqF984Vq3aNsxh0zVH86eSDGDP/B056eASvjlnkIyg6l8LivjurrPC7s0qHpWs3c+vAKYycu4ojW9TmvjMPpXGtKskOy7lya29GNnSuxDXMqswrl3XkH/93CFOWrKPHwyN4fey33lbiXIrxJOJSliTO69iET244mnZNsrj93alc/MI4lq7dnOzQnHOBJxGX8hrVrMJrl3fir6cfzIRFazjpoRG8+ZVflTiXCgrtO0tSReBMoFns8mZ2d+LCcu7nJHFh56Ycc0A2f3h7CrcMnMrgqd9x75mHeNcpziVRPFci7xONOriDqAPG3JdzJa5xrSr0v6ITd/dqw7gFq+n+0Ag+mOwjADiXLPH04tvIzHokPBLn4pSWJi4+ohnHHJDNjW9O4to3vmbEnJXceVobqlaM50/aOVdc4rkS+VLSIQmPxLk91LR2Vd767RH/G6/klEdHMW2p98HlXEmKJ4l0ASZImi1piqSpkqYkOjDn4pE7Xskbv+nMlu07OeOJL3h2xHx/QNG5EhLPeCJN8ys3s0UJiSjB/GHDsmvtpm3cMnAKQ6av4OiWdXjgnLbUrV4p2WE5VyYU+WHDkCyygFPDK6u0JhBXtmVVyeSpCzvwtzMOZtyC1fR8ZCTDZn+f7LCcK9PiGU/keqA/UDe8XpN0baIDc64oJHFBp6Z8cG0X6lSryKUvfsVdH0xn07YdyQ7NuTIpnuqsKcARuaMaSqoKjDazQ0sgvmLn1Vnlx5btO/nH4Jm8PHoRtatm0qfrflx0RFOqZPodXM7tqb3pO0v8fPyQnaHMuZRWKSOdu3odzMC+R9C6QQ3+8fEsutw3jCeHf8PGrX5l4lxxiOdK5CagN/BuKDodeMnMHk5wbAnhVyLl14RFa3jks7mMmLOSmlUyuOLo/eh9ZDOq+bMlzhVqd1cicXUFL+kwolt9BYwws6+LP8SS4UnEff1tlEyGz15JVpUMrujSnN5HNqN6pYxkh+Zcytrj6ixJNcLPWsBC4DXgVWBRKCtsh40lDZM0U9L00ECfO+/a8NzJdEn3x5TfJmlemHdSTHmPUDZP0q0x5c0ljZU0V9KbkjILPROu3GvfpCYvXdqR964+isOa1ORfn86hy33DeHWM33To3J4q6Dr+deAUYAIQe7mi8H6/Qra9A7jZzCZKqk70wOJQoB5RX1yHmtlWSXUBJLUGzgXaAA2A/0g6IGzrceBEopEVv5I0yMxmAPcBD5nZAElPAZcDT8Z57K6ca9c4ixcuOZwpS9bS77O5yQ7HuVJpt0nEzE4JP5sXZcNmthxYHqY3SJpJNMTub4B7zWxrmJd7I38vYEAoXyBpHtAxzJtnZvMBJA0AeoXtHQecH5Z5GbgTTyJuDx3aKIvneh/uXcs7VwTxPCfyWTxlhWyjGdAeGAscABwdqqE+l3R4WKwhsDhmtSWhbHfltYG1ZrYjT3l+++8jabyk8StXrtyT0F05IvlNh87tqd1eiUiqBFQB6kiqyU+39dYgqm6Ki6RqwEDgBjNbL6kCUBPoDBwOvCVpP/K/bdjIP9FZAcv/stDsGeAZiBrW443dOedcwQpqE/ktcANRwpjATx/a64naKAolKYMogfQ3s3dC8RLgHYvqDsZJ2gXUCeWNY1ZvBOQOFJFf+SogS1KFcDUSu7xzzrkSsNvqLDN7JLSH/M7M9jOz5uHV1sweK2zDiuoGngdmmtmDMbPeI2rLIDScZxIlhEHAuZIqSmoOtATGAV8BLcOdWJlEje+DQhIaBpwVttubaAAt55xzJSSeJ9Z3ScrKfSOppqSr4ljvKOAi4DhJk8KrJ/ACsJ+kacAAoLdFpgNvATOAT4CrzWxnuMq4BhgCzATeCssC3ALcFBrhaxMlLeeccyUknifWJ5lZuzxlX5tZ+4RGliD+sKFzzu25vek7K00xt61ISieqgnLOOVfOxdNp0BCiO6ieIrr76Uqi6ibnnHPlXDxJ5BaiO7X6Et2h9SnwXCKDcs45VzoUmkTMbBfRU+D+JLhzzrmfKTSJSFpAPg/xmVlhfWc5V+6MmruKtycs5ubuB9K4VpVkh+NcwsVTnRXbGl8JOBsotBdf58qjhT9s5ONp3zF46ndcdERTrjl2f2pW9ftQXNkV13giv1hJGmVmXRIQT8L5Lb4u0Zav28xDQ+fw9oQlVK1Ygb7dWnDZUc2plJGe7NCcK7IiD0oVBqTKlUZ0ZdLXzNoWb4glw5OIKymzv9vA/Z/M4rNZ37NvjUrceGJLzurQmPQ07+jRlT57k0SGxbzdASwAHjCz2cUbYsnwJOJK2tj5P/CPj2cxafFaWtatxi09WnH8QXW912BXqhQpiUhKA842szcTGVxJ8iTiksHM+Hjad/xzyGwWrNpIx+a1+PsZh7B/3WrJDs25uBTpifVwe+/VCYvKuXJCEj0Pqc+nN3blntMPZs6KDZzy6EheG7PIB8NypVo83Z4MlfS7MGZ6rdxXwiNzrgzKSE/jos5NGXJDVw5vVos/vTeN37wynlU/bk12aM4VSTxtIgvyKbbS+pyIV2e5VLFrl/HSlwu595NZ1KhUgX+e1ZZjW9VNdljO5WtvOmA8KGYskeZhjJHWxR+ic+VLWpq4rEtzBl1zFHWqVeTSl77iz+9NY/O2nckOzbm4xZNEvoyzzDlXBK32rcF7Vx/F5V2a8+qYRZz62CimLV2X7LCci0RiJvUAABozSURBVMtuk4ikfSV1ACpLai/psPDqRjT2unOumFTKSOfPp7Tmtcs7sWHLds544gue+vwbdu7yRneX2grq9uQk4BKiscsf4Kcx1jcAtyc2LOfKpy4t6/DJ9V25/d2p3PvxLDZu3cHN3Q9MdljO7VY8DetnmtnAEoon4bxh3ZUGZsY7E5fS9YBssqtXTHY4zu1Vw3ojSTUUeU7SREndExCjcy6QxJkdGnkCcSkvniRymZmtB7oDdYFLgXsTGpVzzrlSIZ4kktsW0hN40cwmx5Q555wrx+JJIhMkfUqURIZIqg7sSmxYzjnnSoN4ksjlwK3A4Wa2CcgkqtIqUOgmZZikmZKmS7o+z/zfSTJJdcJ7SeonaZ6kKbFd0EvqLWluePWOKe8gaWpYp5+8W1TnnCtRcY2xLmkF0FpSPCMh5toB3GxmE8PVywRJQ81shqTGwInAtzHL/wpoGV6diMZ07xT66bqDaBwTC9sZZGZrwjJ9gDHAYKAH8PEexOicc24vxDPG+n3Ar4EZQG5/DAaMKGg9M1sOLA/TGyTNBBqG7TwE/AF4P2aVXsArFt1zPEZSlqT6QDdgqJmtDvEMBXpIGg7UMLPRofwV4HQ8iTjnXImJ58ridOBAMytyN6OSmgHtgbGSTgOWmtnkPLVPDYHFMe+XhLKCypfkU57f/vsQXbHQpEmToh6Gc6XSjp27SJNI8xEVXQLE0yYyH8go6g4kVQMGAjcQVXH9EfhLfovmU2ZFKP9lodkzZpZjZjnZ2dlxxe1cWfH0iPmc+8wY5q/8MdmhuDIoniSyCZgk6enQeN1PUr94Ni4pgyiB9Dezd4AWQHNgsqSFRF2qTJS0L9GVROOY1RsBywopb5RPuXMuRoOsSsxesYEej4zkyeHfsGOn31zpik88SWQQcA9Rz70TYl4FCndKPQ/MNLMHAcxsqpnVNbNmZtaMKBEcZmbfhf1cHO7S6gysC+0qQ4DukmpKqkn00OOQMG+DpM5hXxfz8zYW5xxwRvtGDL2pK8cdWJf7PpnF6U98wYxl65MdlisjCu07C0BSJnBAeDvbzLbHsU4XYCQwlZ+eK7ndzAbHLLMQyDGzVSERPEZ0h9Um4FIzGx+Wu4yfOn38m5m9GMpzgJeAykQN6tdaIQfkfWe58uzjqcv58/vTWbtpG327teCa4/anYoX0ZIflSoHd9Z0VTweM3YCXgYVE7RCNgd5mVuDdWanKk4gr79Zu2sY9H85k4MQl7F+3GvedeSgdmtZMdlguxe1NB4wPAN3N7Bgz60rURfxDxR2gc65kZFXJ5IFz2vLyZR3ZvG0nZz31JXd9MJ1N23YkOzRXCsWTRDLMbHbuGzObw17creWcSw3HHJDNkBu7clHnprz4xUJOengE4xasTnZYrpSJJ4mMl/S8pG7h9SxxNKw751JftYoVuLvXwbz12yNIk/j1M6O59+NZbN3h47y7+MSTRPoC04HrgOuJnji/MpFBOedKVsfmtRh83dGce3hjnvr8G05//Etmf7ch2WG5UiCehvWqwBYz2xnepwMVQ2eMpY43rDtXsKEzVnDrwCls2LqDP5x0IJcd1dyfdnd71bD+GdEttLkqA/8prsCcc6nlxNb1GHJjV7q2rMNfP5rJhc+PZdnazckOy6WoeJJIJTP7X38JYbpK4kJyziVbnWoVefbiHO79v0OYtHgtJz08gvcnLU12WC4FxZNENuYZ26MD4F9LnCvjJHFuxyYMvu5oWtatxvUDJnHN6xNZsX5LskNzKSSeNpHDgQH81C9VfeDXZlYq79DyNhHn9tyOnbt46vNvePg/c0mTOCunEVd2bUGT2l4pUV4U+Yn1sHIGcCDRE+uz4un2JFV5EnGu6L79YRNPjfiGt8cvYceuXZzWtgF9u+3PgftWT3ZoLsH2ptuTKsBNQFMz+42klkTji3yYmFATy5OIc3tvxfotPDdyPv3HfsumbTs5sXU9rj52f9o1zkp2aC5B9iaJvEn0cOHFZnawpMrAaDNrl5hQE8uTiHPFZ83Gbbz05UJe+nIh6zZvp8v+dbjq2BYcsV9t8gw650q5vbnFt4WZ3Q9sBzCzzeQ/IJRzrpypWTWTG088gC9uPY7be7Zi9ooNnP/sWM55erR3oVJOxJNEtoWrDwOQ1AIo8lC5zrmyp1rFCvTp2oKRfziWe3q1YdEPmzjn6dFc8uI4pi9bl+zwXALFU511IvAnoDXwKXAUcImZDU94dAng1VnOJd7mbTt5efRCnhz+Des2b+eUQ+tzc/cDaV6narJDc0VUpDaRMFBUI6JBojoTVWONMbNViQo00TyJOFdy1m3ezrMj5vP8qAVs27mLc3Iacd3xLam/T+XCV3YpZW8a1ieYWYeERVbCPIk4V/JWbtjK48Pm0X/sIiTR+4imXNVtf2pWzUx2aC5Oe9OwPiY8cOicc0WSXb0id57Whv/e3I1TD23A86MWcOwDw3ln4hLieVbNpa54rkRmEI2vvgjYSFSlZWZ2aOLDK35+JeJc8s3+bgO3vzuVCYvW0PWAbP5+xsE0qulPv6eyvanOappfuZktKqbYSpQnEedSw65dxqtjFnHfJ7MA+MNJB3LxEc282/kUtVfdnpQlnkScSy1L1mzij+9O4/M5KzmsSRb3nXkoLet5NyqpZm/aRJxzLmEa1azCS5cezoPntGX+qo2c3G8U/T6by7Ydu5IdmotDwpKIpMaShkmaKWm6pOtD+T8lzZI0RdK7krJi1rlN0jxJsyWdFFPeI5TNk3RrTHlzSWMlzZX0piS/1cO5UkgS/3dYI/5z0zF0b1OPB4fO4bTHRjFp8dpkh+YKkcgrkR3AzWZ2ENEzJldLag0MBQ4ODfNzgNsAwrxzgTZAD+AJSelhON7HgV8RPfB4XlgW4D7gITNrCawBLk/g8TjnEqxOtYo8dv5hPHdxDms3bef0x7+g72sTmPXd+mSH5nYjYUnEzJab2cQwvQGYCTQ0s0/NbEdYbAzRw4wAvYABZrbVzBYA84CO4TXPzOab2TaisU16hQchjwPeDuu/DJyeqONxzpWcE1rX49ObunLdcfszcu4qejw8kqv6T2D2dxuSHZrLo0TaRCQ1A9oDY/PMugz4OEw3BBbHzFsSynZXXhtYG5OQcsvz238fSeMljV+5cmXRD8Q5V2JqVMrgpu4HMuqWY7nm2P35fPZKejwygqtfn8icFZ5MUkXCk4ikasBA4AYzWx9T/keiKq/+uUX5rG5FKP9lodkzZpZjZjnZ2dl7Er5zLsmyqmTyu5MOZNQtx3FVtxYMn/U9Jz08gmten8hcTyZJVyGRGw8jIg4E+pvZOzHlvYFTgOPtp3uMlwCNY1ZvxE9D8uZXvgrIklQhXI3ELu+cK2NqVs3k9ye14vIu+/HcyPm89OVCPpq6nFMObcCVx+xHmwb7JDvEcimRd2cJeB6YaWYPxpT3AG4BTjOzTTGrDALOlVRRUnOgJTAO+ApoGe7EyiRqfB8Uks8w4Kywfm/g/UQdj3MuNdSqmskferRi1C3HceUxLfhs5gpO7jeKC58by+dzVno3KiUsYQ8bSuoCjASmArk3fN8O9AMqAj+EsjFmdmVY549E7SQ7iKq/Pg7lPYGHgXTgBTP7Wyjfj6ihvRbwNXChmRU41ok/bOhc2bJu03ZeH/ctL36xgO83bOXAetW54ujmnNauARUrpCc7vDLDn1gPPIk4VzZt27GLQZOX8eyI+cxesYG61StyyVHNuKBjU/apkpHs8Eo9TyKBJxHnyjYzY+TcVTw7cj4j566iSmY65+Q0pm+3FtSrUSnZ4ZVau0siCW1Yd865kiaJrgdk0/WAbGYsW89zo+bz2phFvD1hCTeeeAC9j2hKhXTv8am4+Jl0zpVZrRvU4MFz2vHZzceQ06wm93w4g1MeHcX4hauTHVqZ4UnEOVfmNa1dlRcvOZynLuzA+s3bOeup0fz+35P54ccC78NxcfAk4pwrFyTR4+B9+c/Nx3DlMS149+ulHPfA5/Qfu4hdu8pX23Bx8iTinCtXqmRW4NZfteLj64+m1b7V+eO70zjjyS+ZumRdskMrlTyJOOfKpZb1qjOgT2ce/nU7lq7ZzGmPj+LOQdPZuHVH4Su7//Ek4pwrtyRxevuGfHbzMVzUuSkvfbmQHo+M4Mt5q5IdWqnhScQ5V+7tUzmDu3sdzFu/PYIKaWmc/9xYbn93Khu2bE92aCnPk4hzzgUdm9di8HVH85ujmzNg3Lec9NAIPp/jw0cUxJOIc87FqJyZzh9Pbs3bfY+kSsUK9H5hHL//92TWbfKrkvx4EnHOuXwc1qQmH17bhauPbcE7Xy/lxIc+5z8zViQ7rJTjScQ553ajUkY6vz+pFe9ddRS1qmZyxSvj6fvaBMYvXO1dzgfeAaNzzsVh245dPDF8Hs+NXMCPW3fQat/qXNCpCae3b0j1SmW/l2DvxTfwJOKc2xsbt+5g0ORlvDZmEdOXradKZjq92jXkgk5NOLhh2R1d0ZNI4EnEOVcczIwpS9bx2phFfDBlGVu276Jd4ywu6NSEU9s2oFJG2RoQy5NI4EnEOVfc1m3azsCJS+g/dhHfrNxIVpUMrj2uJRd1bkpmhbLR9OxJJPAk4pxLFDNj7ILVPD5sHiPnrqJZ7Src1vMgureuh6Rkh7dXdpdEykaKdM65FCCJzvvV5pXLOvLipYeTkZ7Gb1+dwLnPjGHa0rLZwaMnEeecK2aSOPbAunx8/dH89fSDmff9j5z62ChuemsSy9dtTnZ4xcqTiHPOJUiF9DQu7NyUYb/vxm+7tuDDycs59l/DeXDonDLTW7C3iTjnXAlZvHoT930yiw+nLKdu9Yqce3hjTm3bgJb1qic7tEKVeJuIpMaShkmaKWm6pOtDeS1JQyXNDT9rhnJJ6idpnqQpkg6L2VbvsPxcSb1jyjtImhrW6afS3nLlnCvTGteqwmPnH8bAvkdy4L7VeWzYPE58aAQ9Hh7B48PmseiHjckOcY8l7EpEUn2gvplNlFQdmACcDlwCrDazeyXdCtQ0s1sk9QSuBXoCnYBHzKyTpFrAeCAHsLCdDma2RtI44HpgDDAY6GdmHxcUl1+JOOdSxfcbtvDx1O/4YPIyxi9aA0DbRvtwatsGnHxofervUznJEf4k6bf4SnofeCy8upnZ8pBohpvZgZKeDtNvhOVnA91yX2b221D+NDA8vIaZWatQfl7scrvjScQ5l4qWrt3MR1OW8cHk5UwNd3J1bFaLK45uTvc2+yY5ut0nkQoltPNmQHtgLFDPzJYDhERSNyzWEFgcs9qSUFZQ+ZJ8yvPbfx+gD0CTJk327mCccy4BGmZVpk/XFvTp2oIFqzby4eRlvPP1Uvq8OoEebfblrl5tqFejUrLD/IWE350lqRowELjBzNYXtGg+ZVaE8l8Wmj1jZjlmlpOdnV1YyM45l1TN61Tl2uNb8umNXbmlRyuGzf6eEx78nDfGfcuuXal1M1RCk4ikDKIE0t/M3gnFK0I1Vm67yfehfAnQOGb1RsCyQsob5VPunHNlQkZ6Gn27teCTG7rSpkENbntnKuc9O4b5K39Mdmj/k8i7swQ8D8w0swdjZg0Ccu+w6g28H1N+cbhLqzOwLlR7DQG6S6oZ7uTqDgwJ8zZI6hz2dXHMtpxzrsxoXqcqb/ymM/edeQgzlq+nxyMjeXzYPLbv3JXs0BJ6d1YXYCQwFcg90tuJ2kXeApoA3wJnm9nqkAgeA3oAm4BLzWx82NZlYV2Av5nZi6E8B3gJqAx8DFxrhRyQN6w750qz79dv4c4PpjN46ne02rc69515KG0bZyV8v0m/OytVeBJxzpUFQ6Z/x1/en8bKDVu5qHNTrj2+JXWqVUzY/rwDRuecK0NOarMvQ286hgs6NeXVMYvoev8wHvh0Nuu3bC/ROPxKxDnnSrlvVv7Ig5/O4aOpy8mqkkHfY1rQ+8hmxTowlldnBZ5EnHNl1dQl6/jnp7MZMWcl9WpU5LrjW3JOTmMy0ve+0smrs5xzrow7pNE+vHJZRwb06UzDrMr88d1pnPjg57w/aWnCni/xJOKcc2VM5/1qM7DvkTx3cQ6VMtK5fsAkTn50FCvWbyn2fZVItyfOOedKliROaF2P41rV5YMpyxg8dTnZCbh7y5OIc86VYWlpole7hvRql2/Xgnu//YRs1TnnXLngScQ551yReRJxzjlXZJ5EnHPOFZknEeecc0XmScQ551yReRJxzjlXZJ5EnHPOFVm564BR0kpgUbLjiFMdYFWygygCj7tkedwlq7zG3dTMsvMWlrskUppIGp9fr5mpzuMuWR53yfK4f86rs5xzzhWZJxHnnHNF5kkktT2T7ACKyOMuWR53yfK4Y3ibiHPOuSLzKxHnnHNF5knEOedckXkSSUGSFkqaKmmSpPHJjqcgkl6Q9L2kaTFltSQNlTQ3/KyZzBjzs5u475S0NJz3SZJ6JjPGvCQ1ljRM0kxJ0yVdH8pT+nwXEHdKn28ASZUkjZM0OcR+VyhvLmlsOOdvSspMdqyxCoj7JUkLYs55u73el7eJpB5JC4EcM0v5B5okdQV+BF4xs4ND2f3AajO7V9KtQE0zuyWZcea1m7jvBH40s38lM7bdkVQfqG9mEyVVByYApwOXkMLnu4C4zyGFzzeAJAFVzexHSRnAKOB64CbgHTMbIOkpYLKZPZnMWGMVEPeVwIdm9nZx7cuvRNxeMbMRwOo8xb2Al8P0y0QfGCllN3GnNDNbbmYTw/QGYCbQkBQ/3wXEnfIs8mN4mxFeBhwH5H4Qp+I5313cxc6TSGoy4FNJEyT1SXYwRVDPzJZD9AEC1E1yPHviGklTQnVXSlULxZLUDGgPjKUUne88cUMpON+S0iVNAr4HhgLfAGvNbEdYZAkpmBTzxm1muef8b+GcPySp4t7ux5NIajrKzA4DfgVcHapeXOI9CbQA2gHLgQeSG07+JFUDBgI3mNn6ZMcTr3ziLhXn28x2mlk7oBHQETgov8VKNqrC5Y1b0sHAbUAr4HCgFrDX1Z6eRFKQmS0LP78H3iX6wy1NVoR68Nz68O+THE9czGxF+MfbBTxLCp73UL89EOhvZu+E4pQ/3/nFXRrOdywzWwsMBzoDWZIqhFmNgGXJiqswMXH3CFWLZmZbgRcphnPuSSTFSKoaGh+RVBXoDkwreK2UMwjoHaZ7A+8nMZa45X4QB2eQYuc9NJY+D8w0swdjZqX0+d5d3Kl+vgEkZUvKCtOVgROI2nSGAWeFxVLxnOcX96yYLxsiasfZ63Pud2elGEn7EV19AFQAXjezvyUxpAJJegPoRtTN9ArgDuA94C2gCfAtcLaZpVQj9m7i7kZUtWLAQuC3uW0NqUBSF2AkMBXYFYpvJ2pfSNnzXUDc55HC5xtA0qFEDefpRF+63zKzu8P/6QCiKqGvgQvDt/uUUEDc/wWyAQGTgCtjGuCLti9PIs4554rKq7Occ84VmScR55xzReZJxDnnXJF5EnHOOVdknkScc84VmScRV65JGi4ppwT2c13oxbZ/oveVTJKyJF2V7DhcyfEk4lwRxTyxHI+rgJ5mdkGi4kkRWUTH6soJTyIu5UlqFr7FPxvGRvg0PIX7sysJSXVCN/pIukTSe5I+COMnXCPpJklfSxojqVbMLi6U9KWkaZI6hvWrhk4Bvwrr9IrZ7r8lfQB8mk+sN4XtTJN0Qyh7CtgPGCTpxjzLp0v6l6LxY6ZIujaUHx/2OzXEUTGUL5T0d0mjJY2XdJikIZK+kXRlWKabpBGS3pU0Q9JTktLCvPPCNqdJui8mjh8l/U3R+BNjJNUL5dmSBobz8JWko0L5nSGu4ZLmS7oubOpeoIWisSr+Kal+iGVS2OfRRf5DcKnJzPzlr5R+Ac2AHUC78P4toieEIeoTKCdM1wEWhulLgHlAdaIndNcRPZ0L8BBRJ4C56z8bprsC08L032P2kQXMAaqG7S4BauUTZweip7KrAtWA6UD7MG8hUCefdfoS9SlVIbyvBVQCFgMHhLJXYuJdCPSNOY4pMcf4fSjvBmwhSlzpRD3PngU0IHqiPZuoN4T/AqeHdQw4NUzfD/wpTL8OdAnTTYi6LgG4E/gSqBjO+w9E3Y03yz2HYbmbgT+G6XSgerL/nvxVvK89uRx3LpkWmNmkMD2B6MOqMMMsGr9ig6R1wAehfCpwaMxyb0A0xoikGqHPoe7AaZJ+F5apRPQhClG32vl1K9IFeNfMNgJIegc4mqhbjN05AXjKQrfiZrZaUttwvHPCMi8DVwMPh/eDYo6jWswxbsntLwkYZ2bzQxxvhNi2A8PNbGUo70+UON8DtgEfhnUnACfGxNc66moJgBoKfbsBH1nU1cdWSd8D9fI5vq+AFxR1wPhezO/QlRGeRFxpEdsv0U6gcpjewU/VspUKWGdXzPtd/PxvP2/fP0bUt9CZZjY7doakTsDG3cSo3ZQXRPnsv7DtxB5H3mPMPa7dHdPubDez3HV2xmwnDTjCzDb/LMAoqeT9nfzi8yQk5q7AycCrkv5pZq8UEIcrZbxNxJV2C4mqkeCnXlX31K/hfx0FrjOzdcAQ4NrQ2ymS2sexnRHA6ZKqKOqB+QyijgcL8ilwZW4jfWirmQU0k7R/WOYi4PM9PKaOisYBTyM6vlFEHTUeE9qO0ok6QCxsu58C1+S+UeFjcm8gql7LXb4pUTXbs0Q9+R62h8fhUpxfibjS7l/AW5IuIqrjL4o1kr4EagCXhbJ7iKqPpoREshA4paCNWDSG+EvAuFD0nJkVVJUF8BxwQNjPdqL2mcckXQr8OySXr4Cn9vCYRhM1ch9ClNzeNbNdkm4j6sZcwGAzK6wL8+uAxyVNIfq8GEE0Tne+zOwHSV9ImgZ8TNTV+O/Dsf0IXLyHx+FSnPfi61wZI6kb8DszKzDpOVccvDrLOedckfmViHPOuSLzKxHnnHNF5knEOedckXkScc45V2SeRJxzzhWZJxHnnHNF9v/WXHwyxnqK5wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(2, 35+1), error)\n",
    "plt.xlabel('number of components')\n",
    "plt.ylabel('reconstruction error')\n",
    "plt.title('Non-Negative Matrix Factorization')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the remainder of this section, lets fix the number of so that we can use these results to illustrate the power of NMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "W, H, err = NMF(TDM, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we apply NMF to a Term-Document matrix, as defined above, the resulting factors prove useful for topic detection. In particular, the W matrix provides us with the definition of each topic as a weighted distribution over all the words in the corpus:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at a couple of examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "thanks 24.718993410495905\n",
      "non 12.968158669645955\n",
      "dumb 12.968158669645955\n",
      "factual 12.968158669645955\n",
      "minute 12.968158669645955\n",
      "drops 12.968158669645955\n",
      "followers 12.968158669645955\n",
      "@applenws 12.968158669645955\n",
      "#rumors 12.968158669645955\n",
      "points 12.968158669645955\n"
     ]
    }
   ],
   "source": [
    "topic_words = top_words(W[:,2], word_list)\n",
    "for word, weight in topic_words[:10]:\n",
    "    print(word, weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "big 16.18686839488753\n",
      "defend 14.853961506200685\n",
      "court 14.853961506200685\n",
      "ip 14.853961506200685\n",
      "@samsung 14.853961506200685\n",
      "portfolio 14.853961506200685\n",
      "@robconeybeer 14.853961506200685\n",
      "battles 14.853961506200685\n",
      "look 14.751137247223213\n",
      "companies 14.089748986666253\n"
     ]
    }
   ],
   "source": [
    "topic_words = top_words(W[:,5], word_list)\n",
    "for word, weight in topic_words[:10]:\n",
    "    print(word, weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"width: 100%; overflow: hidden;\">\n",
    "     <img src=\"data/D4Sci_logo_full.png\" alt=\"Data For Science, Inc\" align=\"center\" border=\"0\" width=300px> \n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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
}