unit12_cluster/demo1_doc_cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Demo:  Document Clustering and Latent Semantic Analysis\n",
    "\n",
    "An important application of clustering is for sorting documents into groups.  In this demo, we will illustrate how to use the k-means algorithms for this task.  This example is taken mostly from one of the [sklearn examples](http://scikit-learn.org/stable/auto_examples/text/document_clustering.html).\n",
    "\n",
    "Through the demo, you will learn how to:\n",
    "* Represent a corpus as a set of strings\n",
    "* Build a vocabulary from a corpus\n",
    "* Compute the TF-IDF scores for the documents in the corpus based on the vocabulary\n",
    "* Run k-means to automatically discover document clusters\n",
    "* Display key terms in each document cluster\n",
    "* Perform an LSA on a corpus with a sparse SVD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the Data\n",
    "\n",
    "First, we load the standard packages along with a number of `sklearn` sub-packages for text processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_20newsgroups\n",
    "from sklearn.decomposition import TruncatedSVD\n",
    "from sklearn.feature_extraction.text import TfidfVectorizer\n",
    "from sklearn.feature_extraction.text import HashingVectorizer\n",
    "from sklearn.feature_extraction.text import TfidfTransformer\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import Normalizer\n",
    "from sklearn import metrics\n",
    "\n",
    "from sklearn.cluster import KMeans, MiniBatchKMeans\n",
    "\n",
    "import sys\n",
    "from time import time\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [UseNet newsgroups](https://en.wikipedia.org/wiki/Usenet_newsgroup) were popular 20 years ago as online forums for discussing various issues.  Although they are not used much today for topic discussions, the posts from that era are still widely-used in machine learning classes for demonstrating various text processing methods.  Due to their wide use, the `sklearn` package has a built-in routine `fetch_20newsgroups` for extracting the newsgroup examples.  The code below can be used to get a subset of categories or all categories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading 20 newsgroups dataset for selected categories\n"
     ]
    }
   ],
   "source": [
    "use_all_categories = True\n",
    "if use_all_categories:\n",
    "    categories = None\n",
    "else:\n",
    "    categories = [\n",
    "        'alt.atheism',\n",
    "        'talk.religion.misc',\n",
    "        'comp.graphics',\n",
    "        'sci.space',\n",
    "    ]\n",
    "\n",
    "print(\"Loading 20 newsgroups dataset for selected categories\")\n",
    "dataset = fetch_20newsgroups(subset='all', categories=categories,\n",
    "                             shuffle=True, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print out the categories and the number of documents in each category.  We see the `dataset` has two key fields:\n",
    "\n",
    "* `dataset.target_names`:  The names of the categories\n",
    "* `dataset.target`:  The index of the category for each post"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alt.atheism                   : 799 \n",
      "comp.graphics                 : 973 \n",
      "comp.os.ms-windows.misc       : 985 \n",
      "comp.sys.ibm.pc.hardware      : 982 \n",
      "comp.sys.mac.hardware         : 963 \n",
      "comp.windows.x                : 988 \n",
      "misc.forsale                  : 975 \n",
      "rec.autos                     : 990 \n",
      "rec.motorcycles               : 996 \n",
      "rec.sport.baseball            : 994 \n",
      "rec.sport.hockey              : 999 \n",
      "sci.crypt                     : 991 \n",
      "sci.electronics               : 984 \n",
      "sci.med                       : 990 \n",
      "sci.space                     : 987 \n",
      "soc.religion.christian        : 997 \n",
      "talk.politics.guns            : 910 \n",
      "talk.politics.mideast         : 940 \n",
      "talk.politics.misc            : 775 \n",
      "talk.religion.misc            : 628 \n"
     ]
    }
   ],
   "source": [
    "categories = dataset.target_names\n",
    "labels = dataset.target\n",
    "ncat = len(categories)\n",
    "for i in range(ncat):\n",
    "    cnt = np.sum(labels == i)\n",
    "    print('%-30s: %d ' % (categories[i], cnt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data is in the `data` field of `dataset`.  Each entry `dataset.data[i]` is a string corresponding to the post to the newsgroup.  We can print an example as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Post from rec.autos\n",
      "\n",
      "From: jackw@boi.hp.com (jack wood)\n",
      "Subject: Re: Chevy/GMC 4x4 Fullsize Pickups, Opinions?\n",
      "Distribution: na\n",
      "Organization: Hewlett-Packard / Boise, Idaho\n",
      "X-Newsreader: TIN [version 1.1.4 PL6]\n",
      "Lines: 41\n",
      "\n",
      "Dick Grady (grady@world.std.com) wrote:\n",
      ": \n",
      ": I am considering buying a 1993 Chevy or GMC 4x4 full-size pickup with\n",
      ": the extended cab.  Any opinions about these vehicles?  Have there been\n",
      ": any significant problems?\n",
      ": \n",
      ": -- \n",
      ": Dick Grady           Salem, NH,  USA            grady@world.std.com\n",
      ": So many newsgroups, so little time!\n",
      "\n",
      "\n",
      "I bought a brand new 1992 Chevrolet K2500 HD 4x4 extended cab last\n",
      "May.  It has had many, many problems.  See my earler post that describes\n",
      "the situation.  I went to BBB arbitration, and they ruled that Chevrolet\n",
      "must buy it back from me.  If you do get one, stay away from the 5 speed\n",
      "manual with the deep low first gear.  They have put three of them in my\n",
      "truck so far.  After about 1,500 miles, overdrive either starts\n",
      "rattling or hissing loudly.  There is no way to fix them.  Chevrolet \n",
      "says that the noise is \"a characteristic of the transmission.\"\n",
      "\n",
      "Also, if you are planning to use your truck to tow, the\n",
      "gear ratios in that tranny suck.  On a steep hill, you get up to about\n",
      "55 MPH in second gear at 4,000 RPM (yellow line).  If you shift to third,\n",
      "the RPM drop to only 2,500, and you begin to loose speed.  I should\n",
      "point out that the 350 V8 they put in the HD (8600 GVW) trucks is a\n",
      "detuned motor compared to the one they put in the light duty ones.  They\n",
      "dropped the compression ratio, supposedly for \"engine longevity\"\n",
      "reasons.  So the light duty 350 may pull better than my truck does.\n",
      "Other things that have gone wrong include the ventilation fan (3 times \n",
      "so far), paint (had specs of rust embedded in the paint from being\n",
      "shipped by rail with no covering), and suspension parts (link between\n",
      "stabilizer and control arm fell off).\n",
      "\n",
      "Any company can make a bad individual car, Chevrolet included.  What\n",
      "really bothered me was the way they reacted.  They made no attempt\n",
      "to deal with me except to tell me to take it back to the dealer for \n",
      "them to attempt to fix it one more time.  So I bought a brand new\n",
      "Ford F250 HD Super Cab with a 460 and an automatic.  I will never\n",
      "buy another Chevrolet.\n",
      "\n",
      "jackw@hpdmd48.boi.hp.com\n",
      "\n"
     ]
    }
   ],
   "source": [
    "doc_ind = 20  # Index of an example document\n",
    "data_ex = dataset.data[doc_ind]\n",
    "cat_ex  = dataset.target_names[labels[doc_ind]]\n",
    "print('Post from {0:s}'.format(cat_ex))\n",
    "print()\n",
    "print(data_ex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing TF-IDF features\n",
    "\n",
    "Documents are natively text.  To apply most machine learning algorithms, we need to convert the documents to vectors.  One popular method is the so-called TF-IDF score.  First, we select a set of words in the corpus.  Each word is sometimes called a *token*.  For each token `n` and document `i`, we then compute the data matrix:\n",
    "      \n",
    "    X[n,i] = TF-IDF score of word i in document n\n",
    "           = term freq[n,i] * inverse doc frequency[i]\n",
    "           \n",
    "where\n",
    "\n",
    "    term freq[n,i]  = (#occurances of word i in doc n)/(#words in doc n)  \n",
    "    inverse doc freq[i] = log(#docs in corpus/#docs with word i)\n",
    "        \n",
    "In the data matrix `X`, each document `n` is represented by a vector `X[n,:]`.\n",
    "\n",
    "The data matrix `X` can be computed by a *vectorizer*.  Writing an efficient vectorizer is somewhat time-consuming.  Luckily, `sklearn` has very good routines to compute the TF-IDF representations of a corpus.  We first create a `TfidfVectorizer` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorizer = TfidfVectorizer( \n",
    "                max_df=0.5, # max doc freq (as a fraction) of any word to include in the vocabulary\n",
    "                min_df=2,   # min doc freq (as doc counts) of any word to include in the vocabulary\n",
    "                max_features=10000,           # max number of words in the vocabulary\n",
    "                stop_words='english',         # remove English stopwords\n",
    "                use_idf=True )        # use IDF scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we create the data matrix from the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracting features from the training dataset using a sparse vectorizer\n",
      "done in 3.737130s\n",
      "n_samples: 18846, n_features: 10000\n"
     ]
    }
   ],
   "source": [
    "print(\"Extracting features from the training dataset using a sparse vectorizer\")\n",
    "t0 = time()\n",
    "X = vectorizer.fit_transform(dataset.data)\n",
    "print(\"done in %fs\" % (time() - t0))\n",
    "print(\"n_samples: %d, n_features: %d\" % X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display the terms with the highest TF-IDF scores in a post as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "grady                0.352524 \n",
      "truck                0.217618 \n",
      "gear                 0.216398 \n",
      "hd                   0.192578 \n",
      "boi                  0.169038 \n",
      "chevy                0.166363 \n",
      "rpm                  0.163060 \n",
      "350                  0.147700 \n",
      "dick                 0.143673 \n",
      "duty                 0.141442 \n",
      "paint                0.137826 \n",
      "extended             0.134956 \n",
      "brand                0.129220 \n",
      "std                  0.127186 \n",
      "fix                  0.126137 \n",
      "500                  0.115667 \n",
      "attempt              0.115309 \n",
      "bought               0.110747 \n",
      "light                0.105704 \n",
      "hp                   0.101480 \n",
      "speed                0.100635 \n",
      "buy                  0.094978 \n",
      "com                  0.088686 \n",
      "problems             0.088631 \n",
      "salem                0.088521 \n",
      "opinions             0.088382 \n",
      "rail                 0.086027 \n",
      "embedded             0.084239 \n",
      "characteristic       0.083437 \n",
      "far                  0.082767 \n"
     ]
    }
   ],
   "source": [
    "doc_ind = 20  # Index of an example document\n",
    "xi = X[doc_ind,:].todense()\n",
    "term_ind = np.fliplr(xi.argsort())\n",
    "xi_sort = xi[0,term_ind]\n",
    "terms = vectorizer.get_feature_names()\n",
    "\n",
    "for i in range(30):\n",
    "    term = terms[term_ind[0,i]]\n",
    "    tfidf = xi[0,term_ind[0,i]]\n",
    "    print('{0:20s} {1:f} '.format(term, tfidf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run k-Means\n",
    "\n",
    "We now run k-means on the TF-IDF vectors to try  to automatically detect clusters.  First, we construct a `kMeans` object to perform the clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_clusters = 4\n",
    "km = KMeans(n_clusters=n_clusters, init='k-means++', max_iter=100, n_init=1,\n",
    "                verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Then, we run the k-means clustering.  This will run through several iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Clustering sparse data with KMeans(max_iter=100, n_clusters=4, n_init=1, verbose=True)\n",
      "Initialization complete\n",
      "Iteration 0, inertia 36308.55050753375\n",
      "Iteration 1, inertia 18483.78491841037\n",
      "Iteration 2, inertia 18449.502064743752\n",
      "Iteration 3, inertia 18435.418734316725\n",
      "Iteration 4, inertia 18425.497261784905\n",
      "Iteration 5, inertia 18421.03470869525\n",
      "Iteration 6, inertia 18418.94116476523\n",
      "Iteration 7, inertia 18417.230428206916\n",
      "Iteration 8, inertia 18415.384950171392\n",
      "Iteration 9, inertia 18413.987962558203\n",
      "Iteration 10, inertia 18412.99329233192\n",
      "Iteration 11, inertia 18412.171536205227\n",
      "Iteration 12, inertia 18411.60837661232\n",
      "Iteration 13, inertia 18411.225307704204\n",
      "Iteration 14, inertia 18410.892197952897\n",
      "Iteration 15, inertia 18410.599707570695\n",
      "Iteration 16, inertia 18410.367291957944\n",
      "Iteration 17, inertia 18410.20232294416\n",
      "Iteration 18, inertia 18410.102683963843\n",
      "Iteration 19, inertia 18410.022130359208\n",
      "Iteration 20, inertia 18409.953056320624\n",
      "Iteration 21, inertia 18409.89323704692\n",
      "Iteration 22, inertia 18409.84093257174\n",
      "Iteration 23, inertia 18409.81243694899\n",
      "Iteration 24, inertia 18409.796024149007\n",
      "Iteration 25, inertia 18409.781423535125\n",
      "Iteration 26, inertia 18409.768889107483\n",
      "Iteration 27, inertia 18409.75872799079\n",
      "Iteration 28, inertia 18409.749737643684\n",
      "Iteration 29, inertia 18409.742443900308\n",
      "Iteration 30, inertia 18409.73016934399\n",
      "Iteration 31, inertia 18409.719185963066\n",
      "Iteration 32, inertia 18409.710309146118\n",
      "Iteration 33, inertia 18409.70443808642\n",
      "Iteration 34, inertia 18409.699969241658\n",
      "Iteration 35, inertia 18409.696832537327\n",
      "Iteration 36, inertia 18409.69508185692\n",
      "Iteration 37, inertia 18409.694703918136\n",
      "Converged at iteration 37: center shift 1.5772211300955778e-33 within tolerance 9.844308442311186e-09\n",
      "done in 5.209s\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"Clustering sparse data with %s\" % km)\n",
    "t0 = time()\n",
    "km.fit(X)\n",
    "print(\"done in %0.3fs\" % (time() - t0))\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get an idea of the clusters, we print out the terms corresponding to the 10 largest components of the centroid in each cluster.  You can clearly see the clustering of topics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cluster 0: com netcom article key hp posting clipper nntp host access\n",
      "Cluster 1: windows drive card dos file thanks university use pc com\n",
      "Cluster 2: god people don jesus think com say just believe israel\n",
      "Cluster 3: university ca posting article host nntp cs just like know\n"
     ]
    }
   ],
   "source": [
    "order_centroids = np.fliplr(np.argsort(km.cluster_centers_,axis=1))\n",
    "for i in range(n_clusters):\n",
    "    print(\"Cluster %d:\" % i, end='')\n",
    "    for ind in order_centroids[i, :10]:\n",
    "        print(' %s' % terms[ind], end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next print the percentage each category appears in one of the four clusters.  We see that certain categories are mostly contained in specific clusters found by k-means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alt.atheism                              0.168 0.004 0.631 0.198 \n",
      "comp.graphics                            0.113 0.584 0.003 0.300 \n",
      "comp.os.ms-windows.misc                  0.074 0.819 0.004 0.103 \n",
      "comp.sys.ibm.pc.hardware                 0.118 0.738 0.000 0.144 \n",
      "comp.sys.mac.hardware                    0.077 0.591 0.001 0.331 \n",
      "comp.windows.x                           0.133 0.671 0.000 0.196 \n",
      "misc.forsale                             0.123 0.248 0.001 0.628 \n",
      "rec.autos                                0.334 0.005 0.011 0.649 \n",
      "rec.motorcycles                          0.427 0.007 0.004 0.562 \n",
      "rec.sport.baseball                       0.149 0.000 0.007 0.844 \n",
      "rec.sport.hockey                         0.042 0.000 0.003 0.955 \n",
      "sci.crypt                                0.812 0.046 0.035 0.106 \n",
      "sci.electronics                          0.205 0.177 0.000 0.618 \n",
      "sci.med                                  0.246 0.007 0.040 0.706 \n",
      "sci.space                                0.203 0.014 0.016 0.767 \n",
      "soc.religion.christian                   0.067 0.007 0.772 0.153 \n",
      "talk.politics.guns                       0.376 0.001 0.386 0.237 \n",
      "talk.politics.mideast                    0.073 0.000 0.701 0.226 \n",
      "talk.politics.misc                       0.261 0.000 0.383 0.356 \n",
      "talk.religion.misc                       0.161 0.002 0.588 0.250 \n"
     ]
    }
   ],
   "source": [
    "cnt = np.zeros((ncat,n_clusters), dtype=int)\n",
    "for i, cat in enumerate(categories):\n",
    "    print('%-40s ' % cat, end='')\n",
    "    for j in range(n_clusters):\n",
    "        cnt[i,j] = np.sum((labelkm==j) & (labels==i))\n",
    "    cntp = cnt[i,:] / np.sum(cnt[i,:])\n",
    "    for j in range(n_clusters):\n",
    "        print('%5.3f ' % cntp[j], end='')\n",
    "    print()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LSA Analysis\n",
    "\n",
    "The TF-IDF matrix is typically low rank since there are usually a small number of underlying features.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse.linalg\n",
    "U1,S1,V1 = scipy.sparse.linalg.svds(X,k=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1acb18a4400>]"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAfhUlEQVR4nO3deXxU5d3+8c83CQlbWBPWJBDCvqkQIriBVStuRdEiotatIra2ahfrU61d9KFqqxZ/bg8q4o67oqIobqjIEvY1IQQCYUvYAgGy388fGZ5fmiYkwCRnluv9euWVmTOHmev2JJcn95xzxpxziIhI8IvwOoCIiPiHCl1EJESo0EVEQoQKXUQkRKjQRURCRJRXLxwXF+e6d+/u1cuLiASlxYsX73LOxdf0mGeF3r17d9LT0716eRGRoGRmObU9pikXEZEQoUIXEQkRKnQRkRChQhcRCREqdBGREFFnoZvZNDPLM7NVtTxuZva4mWWZ2QozG+L/mCIiUpf67KFPB0Yf5fELgF6+r4nA0yceS0REjlWdhe6cmwvsOcoqY4CXXKX5QBsz6+yvgCIioSJvfxFT527ghw27G+T5/XFiUVdgS5X7ub5l26uvaGYTqdyLJykpyQ8vLSIS2IrLyvlibR5vL87lm8x8yisct45KYURKe7+/lj8K3WpYVuOnZjjnpgJTAVJTU/XJGiISkpxzrMgt4O3Fucxcvo2Cw6V0atWUW87qweVDE0iJb9kgr+uPQs8FEqvcTwC2+eF5RUSCSt7+It5bupW3F+eyPq+QmKgIzh/QiSuGJnB6zzgiI2ra//UffxT6TOA2M5sBnAoUOOf+Y7pFRCQUFZUemVLZwjeZ+VQ4GNqtLX8fO4iLBnemVdMmjZalzkI3s9eBUUCcmeUCfwaaADjnngFmARcCWcAh4IaGCisiEghqmlLp3Lopt45KYeyQhptSqUudhe6cu6qOxx3wS78lEhEJUDVNqYweWDmlclpKw0+p1MWzy+eKiASD0vIKvlibx5vpW/g6I8/TKZW6qNBFRGqQlVfIm+lbeHdJLrsKS+gQG8OkkSlcMTSBHh5NqdRFhS4i4nOwuIyPV2znjfQtLM7ZS1SEcU6/DoxLTWRk73iiIgP78lcqdBEJa845lmzeyxuLtvDRiu0cKiknJb4Ff7ywL5edkkB8bIzXEetNhS4iYWlXYTHvLsnlzfRcsvIKaR4dycWDO3PlsESGJLXFzNs3OI+HCl1EwkZZeQVz1+fzxqItfLE2j7IKx5CkNjx0+SAuGtyFljHBXYnBnV5EpB427TrIW4u38PbiXHbuLyauZTQ3npHMuNQEenaI9Tqe36jQRSQkFZeVM3v1Tl5fsJkfsncTYTCqTwf++pNEzunXgSYB/gbn8VChi0hIyc4vZMaiyr3xPQdLSGzXjN+f34fLhyTQqXVTr+M1KBW6iAS96nvjURHGuf06MuHUJM7oGUeEx2dwNhYVuogErdr2xn86NIEOrUJ7b7wmKnQRCSo17Y2f178jV6WF1954TVToIhIUat0bT02gQ2z47Y3XRIUuIgFLe+PHRoUuIgEnZ/dBXluwmbe0N35MVOgiEhDKKxxfZ+Tx8vwcvsnMJ9LC80iVE6FCFxFP7Sos5s30Lbw6fzNb9x2mY6sYbj+nF1elJdExDI9UOREqdBFpdEeucPjyDznMWrmDkvIKTktpz70X9ePc/h1D8izOxqBCF5FGc7C4jA+WbePl+Tms3b6f2JgoJpyaxDXDk0LqmipeUaGLSIPLyivklfk5vLM4lwPFZfTr3Iq/jx3EmJO70DxaNeQv+i8pIg2itLyCOWt28vL8HOZt2E10ZAQXDurEtSO6Be31xgOdCl1E/GpXYTGvL9jMKwty2Lm/mK5tmnHX6D6MS00krmXwfPpPMFKhi4hfrMwtYPq8TXy4fBsl5RWc1TueyZd1Y1SfDkTqkMNGoUIXkeNWWl7B7NU7mP79JtJz9tIiOpKr0hL52WndSYlv6XW8sKNCF5FjtruwmNcXbuaV+ZvZsb+Ibu2bc9/F/bkiNYFWTZt4HS9sqdBFpN5Wba2cVpm5fBslZRWc2SuOyWMHMqp3B53JGQBU6CJyVKXlFXy2eifT521k0aa9NI+O5MrURK47rZuOHQ8wKnQRqdHegyW8tnAzr8zPYXtBEUntmvOni/tzxdAEWjfTtEogUqGLyL/Jyitk2vcbeXdJLkWlldMqD1w6UEerBAEVuojgnOOHDbt57ruNfLkuj+ioCMae0pUbz0imd0dNqwQLFbpIGCsuK+fD5dt57tts1u04QFzLaO48tzdXD0/SSUBBSIUuEob2HCzhtQU5vPhDDvkHiundsSUPXz6Yn5zchaZNIr2OJ8dJhS4SRqrPj5/VO55HfprMmb3idG2VEFCvQjez0cAUIBJ4zjn3YLXHWwOvAEm+5/ync+4FP2cVkeOg+fHwUWehm1kk8CRwHpALLDKzmc65NVVW+yWwxjl3iZnFAxlm9qpzrqRBUotInUrKKvhw+Tae+24ja7fvp32LaO44txfXDO+m+fEQVZ899DQgyzmXDWBmM4AxQNVCd0CsVf7N1hLYA5T5OauI1ENhcRkzFm7m+e82sr2giF4dWvLQ5YMYc3JXzY+HuPoUeldgS5X7ucCp1dZ5ApgJbANigSudcxXVn8jMJgITAZKSko4nr4jUIu9AES/O28TLP+Swv6iMU5PbMfmyQYzqE6/58TBRn0Kv6SfBVbt/PrAM+BGQAnxuZt865/b/2z9ybiowFSA1NbX6c4jIccjOL+TZbzfyzpJcSssrGD2gExPP6sEpSW29jiaNrD6FngskVrmfQOWeeFU3AA865xyQZWYbgb7AQr+kFJH/sHTzXv7nm2xmr9lBk8gIrhiawM1n9iA5roXX0cQj9Sn0RUAvM0sGtgLjgQnV1tkMnAN8a2YdgT5Atj+DighUVDi+zszjmW+yWbhxD62aRvHLUT257rTuxMfqjc5wV2ehO+fKzOw2YDaVhy1Oc86tNrNJvsefAe4HppvZSiqnaP7gnNvVgLlFwkpJWQUzl29j6twNZO4spHPrptx7UT/GpyXRMkank0ilev0kOOdmAbOqLXumyu1twI/9G01EDhaX8XqVI1b6dIzl0XEncclJXWgSGeF1PAkw+l+7SADad6iEF+fl8MK8jew7VFp5xMrYQYzqrSNWpHYqdJEAknegiOe/28grP+RwsKScc/t14Bdn92SIjliRelChiwSA3L2HmDo3mzcWbaG0vIKLBnfhF6NS6Ne5ldfRJIio0EU8tCG/kKe/3sD7S7diBmNPSWDSqBQdeijHRYUu4oFVWwt4+usNzFq1nZioCK4Z3o2JZ/WgS5tmXkeTIKZCF2lE6Zv28ORXWXyVkU9sTBS3jkzhxjOSdbEs8QsVukgjmJ+9mylz1vND9m7aNm/C737cm2tHdNeHLYtfqdBFGohzjh98Rb5g4x7iWsZw70X9mHBqEs2j9asn/qefKhE/c84xb0NlkS/ctIf42Bj+dHF/JqQl0Sxal6+VhqNCF/ET5xzfZe1iypz1pOfspWOrGP5ySX/GpyXpOuTSKFToIifIOcfc9buYMieTJZv30alVU/42ZgDjUhNV5NKoVOgix8k5x9eZ+UyZs55lW/bRpXVT7r90IONSE4iJUpFL41Ohixwj5xxfZeQxZc56lucW0LVNM/77soFcMVRFLt5SoYvU05GplUc/y2B5bgEJbZvx4NhBjB2SQHSUrnwo3lOhi9TD/OzdPPJZBos27aVrm2Y8dHllkesSthJIVOgiR7E4Zy+Pfp7B91m76dgqhvvHDGDcsERNrUhAUqGL1GDV1gIe+SyDrzLyad8imnsv6sc1w7vpqBUJaCp0kSoydhzgsc8z+XT1Dlo3a8Jdo/tw3YjutNDHvEkQ0E+pCJCdX8i/5qznwxXbaBkdxR3n9uLGM5Jp1VTXWpHgoUKXsLZlzyGmfLGed5fkEhMVya0jU5h4Vg/aNI/2OprIMVOhS1jaVVjME19m8eqCHMyMG05P5tZRKbqMrQQ1FbqElcLiMp6dm81z32ZTVFbBuNREbj+nF51aN/U6msgJU6FLWCguK+e1BZt54sssdh8s4cJBnfjtj/uQEt/S62gifqNCl5BWUeH4YPlWHvksk9y9hzktpT1/GN2XkxLbeB1NxO9U6BKSnHN8nZHPQ5+uY92OAwzo0orJlw3izF5xmJnX8UQahApdQs7inL089Ok6Fm7cQ7f2zXn8qlO4eFBnIiJU5BLaVOgSMrLzC3nwk3V8tmYncS0rT9O/cliSLpwlYUOFLkFvz8ESpszJ5NUFm4mJiuC35/XmxjOSdXanhB39xEvQKiotZ/q8TTz5ZRYHS8q4Ki2JO87tTXysjiWX8KRCl6DjnOPDFdt56JN1bN13mB/17cB/XdCXXh1jvY4m4ikVugSVRZv28MDHa1m+ZR/9Orfi4SsGc3rPOK9jiQQEFboEhU27DvLgJ+v4dPUOOraK4R9XDGbskAQideSKyP9RoUtA219UyhNfZvHC9xtpEhnBb87rzc/PTKZ5tH50Raqr12+FmY0GpgCRwHPOuQdrWGcU8C+gCbDLOTfSjzklzJRXON5K38I/P8tg98ESrhiSwO/P70OHVrrmikht6ix0M4sEngTOA3KBRWY20zm3pso6bYCngNHOuc1m1qGhAkvoW7hxD3/9cDWrt+0ntVtbXrg+jUEJrb2OJRLw6rOHngZkOeeyAcxsBjAGWFNlnQnAu865zQDOuTx/B5XQt3XfYf4+ay0frdhO59ZNefyqU7hkcGedqi9ST/Up9K7Alir3c4FTq63TG2hiZl8DscAU59xL1Z/IzCYCEwGSkpKOJ6+EoMMl5TzzzQae+WYDALef04tJI1NoFq3P7xQ5FvUp9Jp2j1wNzzMUOAdoBvxgZvOdc5n/9o+cmwpMBUhNTa3+HBJmjhxP/vdZa9leUMQlJ3Xh7gv60rVNM6+jiQSl+hR6LpBY5X4CsK2GdXY55w4CB81sLnASkIlIDTJ2HOBPH6xi4cY9DOjSiinjTyEtuZ3XsUSCWn0KfRHQy8ySga3AeCrnzKv6AHjCzKKAaCqnZB7zZ1AJDYXFZUyZk8m07zcR2zSKyZcN4sphiTqeXMQP6ix051yZmd0GzKbysMVpzrnVZjbJ9/gzzrm1ZvYpsAKooPLQxlUNGVyCy5Hplf/+eA15B4oZPyyRu87vS9sW+jBmEX8x57yZyk5NTXXp6emevLY0rvU7D3DfB6v5IXs3g7q25m9jBnBKUluvY4kEJTNb7JxLrekxnW4nDeZgcRmPf7me57/dSIuYKB64dCBXpSVpekWkgajQxe+cc8xevYO/friG7QVFjEtN4A+j+9K+pS5rK9KQVOjiV9v2Hea+D1YxZ20e/Tq34okJQxjaTdMrIo1BhS5+UV7hmD5vE498loFz8McL+3Lj6clERerj30QaiwpdTtiqrQX817srWbm1gFF94rl/zEAS2zX3OpZI2FGhy3E7WFzGY59nMu37jbRrEcMTE07hokG69oqIV1Tocly+XLeTP72/mq37DjPh1CT+MLovrZs18TqWSFhTocsx2V1YzJ9nruajFdvp1aElb08aQWp3nbIvEghU6FIvR870/MvM1RQWlfGb83ozaWQK0VF601MkUKjQpU55+4u49/1VfLZmJycltuEfVwymd8dYr2OJSDUqdKmVc453l2zlbx+toai0XIciigQ4FbrUaHvBYf747kq+ysgntVtbHr5iMD3iW3odS0SOQoUu/8Y5x4xFW5j88VrKKhx/vqQ/143oToSuvyIS8FTo8n+2FxzmrrdX8O36XYzo0Z6HLh9MUnudICQSLFTognOO95dt5b4PVlNe4Xjg0oFMSEvSXrlIkFGhh7k9B0u4572VfLJqB6nd2vLIuJPo1r6F17FE5Dio0MPYnDU7ufvdlew/XMrdF/Tl5jN76FrlIkFMhR6GDhSV8sBHa3kjfQt9O8Xy8k1p9OvcyutYInKCVOhhZn72bn731nK27TvML0alcPu5vYiJivQ6loj4gQo9TJSWV/Do55k8880Gkto1561JIxjaTddgEQklKvQwsGnXQW6fsZTluQVcmZrIfZf0p0WMNr1IqNFvdQhzzvH24lz+PHM1URHGU1cP4cJBnb2OJSINRIUeogoOlfLH91fy8YrtnJrcjseuPJkubZp5HUtEGpAKPQQt3LiHO99Yxs79Rfz+/D5MGpmiwxFFwoAKPYRUVDie/CqLx+ZkktiuOW/fehonJ7bxOpaINBIVeojYVVjMnW8s49v1u/jJSV2YPHYQLfXGp0hY0W98CJifvZtfv76UgsOl/H3sIMYPS9QHNYuEIRV6EKuocDz1dRaPfp5J9/YtePFGnfEpEs5U6EFKUywiUp0aIAgdmWLZpykWEalChR5Eqk+xTL8hjf5dNMUiIpVU6EGi4FApd765jC/X5WmKRURqpEYIAmu372fSK4vZtu8w948ZwDXDu2mKRUT+Q0R9VjKz0WaWYWZZZnb3UdYbZmblZnaF/yKGtw+WbeWyp77ncEk5MyYO59oR3VXmIlKjOvfQzSwSeBI4D8gFFpnZTOfcmhrWewiY3RBBw01peQWTZ63lhe83kda9HU9cfQodYpt6HUtEAlh9plzSgCznXDaAmc0AxgBrqq33K+AdYJhfE4ahvANF3PbqUhZu2sMNp3fnjxf2o0lkvf6YEpEwVp9C7wpsqXI/Fzi16gpm1hW4DPgRRyl0M5sITARISko61qxhYenmvUx6ZTEFh0uZMv5kxpzc1etIIhIk6rPbV9OErat2/1/AH5xz5Ud7IufcVOdcqnMuNT4+vr4Zw8Z7S3O5cup8oqMieO8Xp6vMReSY1GcPPRdIrHI/AdhWbZ1UYIbvzbo44EIzK3POve+XlCGuosLx8OwMnvlmA8N7tOOpq4fSrkW017FEJMjUp9AXAb3MLBnYCowHJlRdwTmXfOS2mU0HPlKZ109hcRl3zFjKnLV5TDg1ib/+ZIDmy0XkuNRZ6M65MjO7jcqjVyKBac651WY2yff4Mw2cMWRt2XOIn7+YTlZ+IX8bM4BrdXy5iJyAep1Y5JybBcyqtqzGInfOXX/isULfok17uOXlxZSVV/DiDWmc0SvO60giEuR0pqgHPly+jd++tZyENs14/vphJMe18DqSiIQAFXojcs7xP3OzefCTdQzr3pZnf5ZKm+Z681NE/EOF3kjKyiv4y4ereWX+Zi4e3Jl//vQkmjaJ9DqWiIQQFXojOFhcxq9eX8qX6/KYNDKFu87vQ0SE3vwUEf9SoTewvANF3DQ9ndXbCnjg0oFcM7yb15FEJESp0BvQ5t2HuOb5BewqLOa561L5Ud+OXkcSkRCmQm8gGTsOcO3zCygpr+C1m4dzcmIbryOJSIhToTeAJZv3csMLi2jaJII3bxlB746xXkcSkTCgQvezb9fnM/GlxXRoFcMrN51KYrvmXkcSkTChQvejT1dt51evLyUlviUv3ZSmD6QQkUalQveTD5Zt5TdvLuekhNa8cEMarZs18TqSiIQZFbofvLM4l9+/vZzU7u144fphtIjRf1YRaXy6TusJemPRZn739nJGpLRn+g0qcxHxjtrnBLy6IId73lvFWb3jmXrtUJ3KLyKe0h76cXp7cS73vLeKs/uozEUkMKjQj8OHy7dx19vLOaNnHE9fozIXkcCgQj9Gn6/ZyZ1vLCO1Wzum/kxlLiKBQ4V+DNI37eG215YwoEsrnr8+lebRegtCRAKHCr2eNuQX8vOX0unSphkv3JBGbFMdZy4igUWFXg95B4q4btpCoiKMF29Io10LfcqQiAQezRnUobC4jBunL2J3YQlv3DKcpPa6NouIBCYV+lGUllfwi1eXsHb7AZ77WSqDE3QJXBEJXJpyqYVzjnveW8nczHwmXzaQs/t28DqSiMhRqdBr8fx3G3kzPZdf/6gnVw5L8jqOiEidVOg1+Cojj8mz1nLBwE7ccW5vr+OIiNSLCr2arLxCfv3aUvp2asUj404iIsK8jiQiUi8q9CoKDpVy80vpxDSJ4NnrdOKQiAQXNZZPeYXjtteXkLv3EK/fPJyubZp5HUlE5Jio0H2mzMnk2/W7eOjyQaR2b+d1HBGRY6YpFyrfBH38yyzGpSboiBYRCVphX+hb9x3mzjeW0bdTLH8bM9DrOCIixy2sC72krPJM0LJyp+uai0jQC+s59Mmz1rJ8yz6evnoIyXEtvI4jInJC6rWHbmajzSzDzLLM7O4aHr/azFb4vuaZ2Un+j+pfs1fvYPq8Tdx4ejIXDOrsdRwRkRNWZ6GbWSTwJHAB0B+4ysz6V1ttIzDSOTcYuB+Y6u+g/rRzfxF3v7OCAV1acfcFfb2OIyLiF/XZQ08Dspxz2c65EmAGMKbqCs65ec65vb6784EE/8b0n4oKx2/fXM7h0nKmjD+F6KiwfhtBREJIfdqsK7Clyv1c37La3AR8UtMDZjbRzNLNLD0/P7/+Kf3o+e828l3WLu67eAA9O7T0JIOISEOoT6HXdDETV+OKZmdTWeh/qOlx59xU51yqcy41Pj6+/in9ZNXWAh6evY7zB3TkqrTERn99EZGGVJ+jXHKBqu2XAGyrvpKZDQaeAy5wzu32Tzz/OVxSzu0zltKuRTQPjh2MmS66JSKhpT576IuAXmaWbGbRwHhgZtUVzCwJeBe41jmX6f+YJ+4fszPYkH+QR8edTFt9JqiIhKA699Cdc2VmdhswG4gEpjnnVpvZJN/jzwD3Ae2Bp3x7vmXOudSGi31sftiwm2nfb+S6Ed04vWec13FERBqEOVfjdHiDS01Ndenp6Q3+OoXFZYz+11yiIoxZt5+pS+KKSFAzs8W17TCHfLtNnrWWrfsO89YtI1TmIhLSQvog7HkbdvHags3cfGYPXRJXREJeyBZ6cVk5976/iqR2zfnNefpcUBEJfSE7B/Hs3Gyy8w8y/YZhuoqiiISFkNxDz9l9kP/3ZRYXDerMqD4dvI4jItIoQq7QnXPc98FqmkRG8KeLq19DTEQkdIVcoX+2ZiffZObzm/N606l1U6/jiIg0mpAq9OKycibPWkuvDi352YhuXscREWlUIVXoL83LIWf3Ie69uD9RkSE1NBGROoVM6+0uLObxL9Zzdp94RvZu/Cs5ioh4LWQK/bE5mRwqLeeei/p5HUVExBMhUejrdx7gtQWbuXZ4N3p2iPU6joiIJ0Ki0B+bk0nz6Ch+fU4vr6OIiHgm6At91dYCZq3cwY1nJNNO1zkXkTAW9IX+2OeZtGoaxU1nJHsdRUTEU0Fd6Cty9/HFujxuGZlC62ZNvI4jIuKpoC70ad9tpGVMlE4iEhEhiAs9b38RH6/czk9TE4htqr1zEZGgLfRXF2ymrMJx3YjuXkcREQkIQVnoxWXlvLogh7P7dKB7XAuv44iIBISgLPQv1uaxq7CE60/r7nUUEZGAEZSF/vmanbRrEc3pPeO8jiIiEjCCrtDLyiv4KiOPUX3iiYwwr+OIiASMoCv0xTl72XeolPP6dfQ6iohIQAm6Qo+MMEb2judMXSJXROTfRHkd4Fildm/HizemeR1DRCTgBN0euoiI1EyFLiISIlToIiIhQoUuIhIiVOgiIiFChS4iEiJU6CIiIUKFLiISIsw5580Lm+UDOcf5z+OAXX6M4yWNJTBpLIFJY4FuzrkaT5X3rNBPhJmlO+dSvc7hDxpLYNJYApPGcnSachERCREqdBGREBGshT7V6wB+pLEEJo0lMGksRxGUc+giIvKfgnUPXUREqlGhi4iEiKArdDMbbWYZZpZlZnd7nedYmdkmM1tpZsvMLN23rJ2ZfW5m633f23qdsyZmNs3M8sxsVZVltWY3s//ybacMMzvfm9Q1q2UsfzGzrb5ts8zMLqzyWECOxcwSzewrM1trZqvN7Hbf8qDbLkcZSzBul6ZmttDMlvvG8lff8obdLs65oPkCIoENQA8gGlgO9Pc61zGOYRMQV23Zw8Ddvtt3Aw95nbOW7GcBQ4BVdWUH+vu2TwyQ7NtukV6PoY6x/AX4XQ3rBuxYgM7AEN/tWCDTlzfotstRxhKM28WAlr7bTYAFwPCG3i7BtoeeBmQ557KdcyXADGCMx5n8YQzwou/2i8ClHmaplXNuLrCn2uLaso8BZjjnip1zG4EsKrdfQKhlLLUJ2LE457Y755b4bh8A1gJdCcLtcpSx1CaQx+Kcc4W+u018X44G3i7BVuhdgS1V7udy9A0eiBzwmZktNrOJvmUdnXPbofKHGujgWbpjV1v2YN1Wt5nZCt+UzJE/h4NiLGbWHTiFyr3BoN4u1cYCQbhdzCzSzJYBecDnzrkG3y7BVuhWw7JgO+7ydOfcEOAC4JdmdpbXgRpIMG6rp4EU4GRgO/CIb3nAj8XMWgLvAHc45/YfbdUalgX6WIJyuzjnyp1zJwMJQJqZDTzK6n4ZS7AVei6QWOV+ArDNoyzHxTm3zfc9D3iPyj+rdppZZwDf9zzvEh6z2rIH3bZyzu30/RJWAM/y///kDeixmFkTKgvwVefcu77FQbldahpLsG6XI5xz+4CvgdE08HYJtkJfBPQys2QziwbGAzM9zlRvZtbCzGKP3AZ+DKyicgzX+Va7DvjAm4THpbbsM4HxZhZjZslAL2ChB/nq7cgvms9lVG4bCOCxmJkBzwNrnXOPVnko6LZLbWMJ0u0Sb2ZtfLebAecC62jo7eL1u8HH8e7xhVS++70BuMfrPMeYvQeV72QvB1YfyQ+0B74A1vu+t/M6ay35X6fyT95SKvcobjpaduAe33bKAC7wOn89xvIysBJY4fsF6xzoYwHOoPJP8xXAMt/XhcG4XY4ylmDcLoOBpb7Mq4D7fMsbdLvo1H8RkRARbFMuIiJSCxW6iEiIUKGLiIQIFbqISIhQoYuIhAgVuohIiFChi4iEiP8F5UvWMzYN4u0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s = np.flipud(S1)\n",
    "pov = np.cumsum(s**2) / np.sum(s**2)\n",
    "plt.plot(pov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import TruncatedSVD\n",
    "pca = TruncatedSVD(n_components=100)\n",
    "X_tr = pca.fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1acbd422a00>]"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxV1bn/8c9DIMwzQTAQEiBMCgiGSdTigAWt4lAronUu1ypVW63Xqbej1vZXa7FVuVRREBUFRWm1TjgrQhKZCUOAQAYgzAgBMpzn98c59qYxyAESTrLzfb9eecE+a+99nmXC15111t7L3B0REQmuerEuQEREqpeCXkQk4BT0IiIBp6AXEQk4Bb2ISMDVj3UBlWnXrp0nJyfHugwRkVojMzNzm7snVNZWI4M+OTmZjIyMWJchIlJrmNmGQ7Vp6EZEJOAU9CIiAaegFxEJOAW9iEjAKehFRAJOQS8iEnAKehGRgFPQi4jEWHFpiLeWbWbSR2ur5fw18oYpEZGgc3eW5u/mlcw85iwuYGdRCYmtGnPD8BTi61ftNbiCXkTkOCrcc4DZC/OZlZnHmsK9NKxfj/NO6sClAxM5o3s76sdV/UCLgl5EpJodLC1jblYhMzNy+Wj1VkIOp3Zpze8v7cv5fTvSsnGDan1/Bb2ISDVZXrCbmRl5vLYon11FJXRo0Ygfj+jGZQM70TWh2XGrQ0EvIlKFdheV8PrifF5Kz2V5wR7i4+px3kkncHlaZ07v3o64enbca4oq6M1sFDARiAOecveHK7T3Ap4BBgL3u/ufKrTHARlAvrt/ryoKFxGpKUIh54t123kpI5d/LdtMcWmIPh1b8KsL+3DxgERaNYmPaX2HDfpISD8OjATygHQzm+PuK8rttgO4Dbj4EKe5HcgCWhxbuSIiNceWPQeYlZnHS+m5bNxRRItG9Rk7qDM/SOvMyYktY13ev0VzRT8YyHb3dQBmNgMYA/w76N29ECg0swsqHmxmnYALgAeBn1VF0SIisVJaFuLDVVuZkb6RD1ZtpSzkDO3ahp+N7MGokzvQqEFcrEv8hmiCPhHILbedBww5gvf4C3A30PzbdjKz8cB4gKSkpCM4vYhI9cvdUcRL6bnMzMxly56DJDRvyPgzu/KDtM6ktGsa6/K+VTRBX9knBx7Nyc3se0Chu2ea2Yhv29fdJwOTAdLS0qI6v4hIdSopCzE3awvPz9/Ip9nbMOA7PRL4zZgkzu7VngbVMOe9OkQT9HlA53LbnYCCKM8/HLjIzM4HGgEtzGy6u199ZGWKiBw/uTuKmJG+kZcz8tj61UE6tmzEbWencsWgzpzYqnGsyzti0QR9OpBqZilAPjAWGBfNyd39XuBegMgV/V0KeRGpiUrLQsxdWcjz8zfyyZqtGHBWz/aMG5LEiJ7tYzItsqocNujdvdTMJgBvE55eOcXdl5vZzZH2SWbWgfD0yRZAyMzuAPq4+55qrF1E5Jht2r2fGQtymZG+kS17DtKhRe2+eq+Mude84fC0tDTPyMiIdRkiElChkPPZ2m1M/2ID72UVEnLnzNQErhoSHnuvjufNVDczy3T3tMradGesiNQZu4qKmZWZx/PzN7J+2z7aNI3npjNSuGpwF5LaNol1edVGQS8igbc4dxfPfbGBfywu4GBpiFO7tOa2c7pzft+ONKxf8+a9VzUFvYgE0oGSMv6xuIDpX2xgcd5umsTHcdmpnfjh0C707li3btJX0ItIoOTuKGL6/A28lJ7LrqISUts34zdjTuKSAYk0b1S9jwOuqRT0IlLrhULOp9nbmDYvh7krC6lnxnl9TuCaYckM7doGs9o7NbIqKOhFpNbac6CEWRl5PPfFBtZv20e7ZvFMOKs744Yk0bFlMKZGVgUFvYjUOqu3fMXUz3OYvTCfouIyBia14vYrTmF03w514sPVI6WgF5FaobQsxHtZhUz9PId567YTX78eF/U/kWuHJdO3U815JHBNpKAXkRptV1ExM9JzeW7eBvJ37efElo24e1RPxg5Kok3T2C7oUVso6EWkRlq5ec+/h2cOlIQY2rUNv/heb87tfUKtvHM1lhT0IlJjlIWcuVlbeOaz8PBMw/r1uHRgIteelkyvDnVr7ntVUtCLSMztOVDCy+m5TJ2XQ+6O8PDMf4/qxdhBnWmt4ZljpqAXkZjJ2baPZz/PYWZGLvuKyxic3IZ7R/fmvD4anqlKCnoROa7cnfnrd/DUJ+uZu3IL9esZF/Y7keuHp2j2TDVR0IvIcVFcGuLNpZt46tN1LMvfQ+smDbh1RHd+OKwLJ7RoFOvyAk1BLyLVandRCS8s2MjUz3PYvOcA3RKa8tAlfbl0YCKNGujmpuNBQS8i1SJn2z6e+Ww9L2fksb+kjNO7t+P3l/XlO6kJ1KvFy/LVRgp6Eaky7k7mhp38/ZN1vLMiPP4+5pREbjw9pc49GrgmUdCLyDErCzlvL9/M3z9Zx8KNu2gVGX+/ZlgX2mv8PeaiCnozGwVMJLw4+FPu/nCF9l7AM8BA4H53/1Pk9c7ANKADEAImu/vEqitfRGKpqLiUmRl5PP3pejbuKKJL2yb8ZsxJfP/UTjSJ13VkTXHY74SZxQGPAyOBPCDdzOa4+4pyu+0AbgMurnB4KXCnu39pZs2BTDN7t8KxIlLLbN97kKnzNjBtXg67ikoYkNSK+87vxcg+HYjT+HuNE83/cgcD2e6+DsDMZgBjgH+HtbsXAoVmdkH5A919E7Ap8vevzCwLSCx/rIjUHhu27+Pvn6xjZkYeB0tDjOxzAv91ZlfSktvEujT5FtEEfSKQW247DxhypG9kZsnAAGD+IdrHA+MBkpKSjvT0IlKNluXv5smP1vKvpZuoX68elwxI5EdndqV7+2axLk2iEE3QV/Z7mB/Jm5hZM+AV4A5331PZPu4+GZgMkJaWdkTnF5Gq5+58lr2dSR+t5dPsbTRvWJ/xZ3bj+uHJusGplokm6POAzuW2OwEF0b6BmTUgHPLPu/urR1aeiBxvZSHnrWWbmfTRWpbm76Z984bcM7oX44Yk0aKOLq5d20UT9OlAqpmlAPnAWGBcNCe38Iq8TwNZ7v7no65SRKrdwdIyXsnMZ/LHa8nZXkRKu6b8/tLwHaxanq92O2zQu3upmU0A3iY8vXKKuy83s5sj7ZPMrAOQAbQAQmZ2B9AH6Af8EFhqZosip7zP3d+shr6IyFHYe7CUF+Zv4KlP1lP41UH6JrbkyasGct5JmkETFFFNdI0E85sVXptU7u+bCQ/pVPQplY/xi0iM7dxXzLOf5/Ds5zns3l/Cad3a8ucfnMLw7m0J/zIuQaE7GkTqmMI9B/j7J+t4fv5GiorLOK/PCfx4RDcGJLWOdWlSTRT0InVE7o4i/vfjtbyckUdpWYgL+5/ILSO607ND81iXJtVMQS8ScNmFe3niw2xeX1RAnBmXndqJm7/TlS5tm8a6NDlOFPQiAbW8YDePf5DNv5ZtpmH9elw7LJnxZ3alQ0vNga9rFPQiAfPlxp387f1s3l9ZSPOG9bllRDduGJ5C22YNY12axIiCXiQAvl6H9a/vr+Gz7O20btKAu87rwQ+HJdOysW5yqusU9CK1mLvzafY2Hpu7hvScnbRr1pD7zu/FVUO60LSh/nlLmH4SRGohd+fD1Vt5bO4aFm7cRceWjfj1RSdxxaDOWodVvkFBL1KLuDtzswp57P01LMnbTWKrxjx0SV8uO1WPKZBDU9CL1ALuzjsrtvDY3DUsL9hDUpsm/OGyvlw6sBMN4urFujyp4RT0IjVYKBQO+Ilz15C1aQ/JbZvwp8v7M+aUExXwEjUFvUgNVDHgU9o15c8/6M9F/U+kvgJejpCCXqQG+XqIZuJ7a1gRCfhHr+jPhf0U8HL0FPQiNYC78/7KQh59bzXL8sNDNLqCl6qioBeJIXfno9VbefTd1SzO201Sm/AY/MWnKOCl6ijoRWLg6/VY//zuKr7cuItOrRtrFo1UGwW9yHE2f912Hnl3NQvW76Bjy0Y8eMnJXH5qZ+LrK+CleijoRY6ThRt38sg7q/k0exvtmzfUnaxy3ER1CWFmo8xslZllm9k9lbT3MrN5ZnbQzO46kmNFgi5r0x5umprOJU98TtamPTxwQW8+vvssrj0tWSEvx8Vhr+jNLA54HBgJ5AHpZjbH3VeU220HcBtw8VEcKxJI67bu5dH31vCPxQU0b1Sfu87rwfXDU/SwMTnuovmJGwxku/s6ADObAYwB/h3W7l4IFJrZBUd6rEjQFOzaz8T31jDryzwa1q/HrWd1Y/wZ3WjZRI8LltiIJugTgdxy23nAkCjPH/WxZjYeGA+QlJQU5elFao7tew/yxIdreW7eBgCuGdaFW0Z0J6G5FvyQ2Iom6K2S1zzK80d9rLtPBiYDpKWlRXt+kZjbe7CUpz5Zx98/Xsf+kjK+f2onbj+3B4mtGse6NBEguqDPAzqX2+4EFER5/mM5VqRGO1haxvNfbORvH2SzY18xo0/uwJ3n9aB7++axLk3kP0QT9OlAqpmlAPnAWGBclOc/lmNFaqSykDN7YT6Pvrua/F37Gd69LXd/txf9O7eKdWkilTps0Lt7qZlNAN4G4oAp7r7czG6OtE8ysw5ABtACCJnZHUAfd99T2bHV1RmR6vT182j+8NZKVm/ZS9/Eljx8WV/OSE2IdWki38rca95weFpammdkZMS6DJF/W5S7i4fezGLB+h2ktGvKz7/bk9End8Csso+hRI4/M8t097TK2jShV+RbbNi+jz++vYo3lmyiXbN4fjvmJMYOTtLzaKRWUdCLVGLnvmL++n42z32RQ/169bjt7O6M/043mulmJ6mF9FMrUs6BkjKmfp7D3z7IZt/BUi4/tTM/O68HJ7RoFOvSRI6agl6E8NJ9/1hSwB/fWkX+rv2M6JnAvaN707ODpkpK7aeglzpvwfodPPjGChbn7aZPxxb88fv9GN69XazLEqkyCnqps9Zv28fD/8ri7eVb6NCiEY9c3p9LBiRSr55m0kiwKOilztlVVMzEuWt4bt4GGtavx13n9eDG07vSOF6PDJZgUtBLnVFcGuK5Lzbw2Nw1fHWghCsGJfHTkam0b64PWiXYFPQSeO7O3KxCHnwzi/Xb9nFGajseuKCPPmiVOkNBL4G2estX/PafK/hkzTa6JTTlmesGMaJngu5olTpFQS+BtKuomEffXc30+RtpGh/HLy/sw9VDu+iOVqmTFPQSKKVlIV5csJFH3l3Nnv0lXDWkCz8b2YPWTeNjXZpIzCjoJTAWrN/B/7y+jJWbv2JY17b88qI+9OrQItZlicScgl5qvS17DvD7N7N4bVEBia0a88RVA/VkSZFyFPRSa5WUhXjms/VMfG8NJSHntrO78+MR3TUfXqQCBb3USp9nb+N/5iwnu3Av5/Rqz/9c2IcubZvGuiyRGklBL7XKlj0H+N0bWfxjcQFJbZrw9LVpnNP7hFiXJVKjKeilVigtCzF13gYefXc1xWUhbj8nlR+P6EajBhqmETmcqILezEYBEwmv+/qUuz9cod0i7ecDRcB17v5lpO2nwE2AA0uB6939QJX1QAIvc8NOHnhtGVmb9vCdHgn8+qKTSG6nYRqRaB026M0sDngcGAnkAelmNsfdV5TbbTSQGvkaAjwJDDGzROA2wguF7zezl4GxwLNV2gsJpJ37ivnDWyuZkZ5LhxaNePKqgYzSbBqRIxbNFf1gINvd1wGY2QxgDFA+6McA0zy80vgXZtbKzDqWe4/GZlYCNAEKqqx6CSR355Uv83nozSx27y/hR2ekcPu5PbSMn8hRiuZfTiKQW247j/BV++H2SXT3DDP7E7AR2A+84+7vHEO9EnDZhXu5f/ZS5q/fwcCkVjx4SV96d9RNTyLHIpqgr+z3ZI9mHzNrTfhqPwXYBcw0s6vdffo33sRsPDAeICkpKYqyJEgOlJTxxAfZPPnRWho3iOOhS/oydlBnLQIiUgWiCfo8oHO57U58c/jlUPucC6x3960AZvYqcBrwjaB398nAZIC0tLSK/yORAPssexv3z15KzvYiLhmQyH3n9yahecNYlyUSGNEEfTqQamYpQD7hD1PHVdhnDjAhMn4/BNjt7pvMbCMw1MyaEB66OQfIqLLqpVbbvvcgD76RxasL80lu24TnbxqitVpFqsFhg97dS81sAvA24emVU9x9uZndHGmfBLxJeGplNuHplddH2uab2SzgS6AUWEjkql3qLndn9sJ8fvvPFew9WMpPzu7OrWd115x4kWpi4YkyNUtaWppnZOjCP4gKdu3nvtlL+XDVVgYmteIPl/Uj9QSt9CRyrMws093TKmvTfDU5LkIh58X0jfz+zZWUhZxfXtiHa4YlE6cPW0WqnYJeqt3G7UX89ytLmLduO6d1a8vDl/YjqW2TWJclUmco6KXahELOtHk5/OGtVdSvZ/z+0vCUSd3ZKnJ8KeilWuRs28fds5awIGcHZ/VM4KFL+9KxZeNYlyVSJynopUp9fRX/8FsraRBXjz9d3p/LBibqKl4khhT0UmU27d7Pz2cu4dPsbYzomcDDl/ajQ8tGsS5LpM5T0EuV+MfiAu6fvZSSMufBS05m3OAkXcWL1BAKejkmu4tK+MXry5izuIABSa348w9OIUXPihepURT0ctQ+y97GnS8vZtveg9w5sgc/HtGN+nH1Yl2WiFSgoJcjVlwa4pF3VzH543WktGvKq9ecRr9OrWJdlogcgoJejkjOtn3cNmMhS/J2M25IEr+4oA+N4/WMGpGaTEEvUXttYT73z15K/bh6TLp6IKNO7nj4g0Qk5hT0clhFxaX88vXlzMzMY1ByayaOHcCJrXTzk0htoaCXb7Vq81fc+sKXrN26l5+c3Z3bz0nVB64itYyCXirl7ryUnssv5yyneaMGTL9Ri4KI1FYKevmGvQdLuX/2Ul5fVMDp3dvx6BWnaGk/kVpMQS//IWvTHm59/ktytu/jzpE9uOWs7npmvEgtp6AX4D+Halo0bsDzNw1lWLe2sS5LRKqAgl4oKi7lgdnLeHVhvoZqRAIoqukTZjbKzFaZWbaZ3VNJu5nZY5H2JWY2sFxbKzObZWYrzSzLzIZVZQfk2GQX7uXixz9j9qJ87jg3lak3DFbIiwTMYa/ozSwOeBwYCeQB6WY2x91XlNttNJAa+RoCPBn5E2Ai8Ja7f9/M4gGtIVdDvLFkE3fPWkzDBnFMu2EwZ6QmxLokEakG0QzdDAay3X0dgJnNAMYA5YN+DDDN3R34InIV3xHYB5wJXAfg7sVAcdWVL0ejpCzEw/9aydOfrmdAUiueuGqgVn8SCbBogj4RyC23ncf/Xa1/2z6JQCmwFXjGzPoDmcDt7r6v4puY2XhgPEBSUlK09csRKvzqABOeX8iCnB1cd1oy953fm/j6ugFKJMii+Rde2dw6j3Kf+sBA4El3H0D4Cv8bY/wA7j7Z3dPcPS0hQUMI1SFzw04u/OunLMnfxcSxp/Cri05SyIvUAdFc0ecBncttdwIKotzHgTx3nx95fRaHCHqpPu7OCws28qs5y+nYsjGzbxlM744tYl2WiBwn0VzOpQOpZpYS+TB1LDCnwj5zgGsis2+GArvdfZO7bwZyzaxnZL9z+M+xfalmB0vLuPfVpdw/exmndWvHnAnDFfIidcxhr+jdvdTMJgBvA3HAFHdfbmY3R9onAW8C5wPZQBFwfblT/AR4PvI/iXUV2qQabd59gJunZ7Iodxe3ntWNn43sqbtcReogC0+UqVnS0tI8IyMj1mXUapkbdnLz9Ez2HSzlkcv7M7qvnh0vEmRmlunuaZW16c7YAHo5PZcHXltGx1aNmH7jEHp2aB7rkkQkhhT0AVJaFuJ3b2Tx7Oc5nN69HX8bN4BWTeJjXZaIxJiCPiB27y9hwgtf8smabdwwPIX7zu+lBUJEBFDQB8KG7fu44dl0Nmwv4g+X9eWKQbrhTET+j4K+lkvP2cH4aRk48NyNQ/RoYRH5BgV9LTZ7YR7/PWspnVo35unrBpHSrmmsSxKRGkhBXwu5O4++u5rH3s9mWNe2PHn1QH3oKiKHpKCvZQ6WlnH3rCW8vqiAH6R14ncX99XzakTkWynoa5Gd+4r5r+cyWZCzg59/tye3jOiGme50FZFvp6CvJTZuL+K6ZxaQt3M/j105gIv6nxjrkkSkllDQ1wKLc3dxw7PplIac6TcNYXBKm1iXJCK1iIK+hpubtYUJLyykXfN4nr1+MN0SmsW6JBGpZRT0NdgL8zfywGtLOTmxJU9fO0iLdovIUVHQ10Duzl/eW8PEuWs4u1d7/jZuAE3i9a0SkaOj9KhhykLOA68t48UFG/lBWiceuqSvnlkjIsdEQV+DHCgp47YXF/LOii3celY37jqvp6ZPisgxU9DXELv3l/CjaRmk5+zgVxf24brhKbEuSUQCQkFfAxTuOcA1UxawduteJo7VHHkRqVpRDf6a2SgzW2Vm2WZ2TyXtZmaPRdqXmNnACu1xZrbQzP5ZVYUHRc62fVw26XM27ihiynWDFPIiUuUOG/RmFgc8DowG+gBXmlmfCruNBlIjX+OBJyu03w5kHXO1AbOiYA/fnzSPvQdKeeFHQzkjNSHWJYlIAEVzRT8YyHb3de5eDMwAxlTYZwwwzcO+AFqZWUcAM+sEXAA8VYV113rpOTu4YvI8GsQZM28eximdW8W6JBEJqGiCPhHILbedF3kt2n3+AtwNhL7tTcxsvJllmFnG1q1boyir9np/5Raufmo+Cc0bMuvHp9G9vRbvFpHqE03QVza/z6PZx8y+BxS6e+bh3sTdJ7t7mrunJSQEdwjj9UX5jJ+WSeoJzZj5X8NIbNU41iWJSMBFE/R5QOdy252Agij3GQ5cZGY5hId8zjaz6UddbS333Lwc7nhpEad2ac2LPxpK22Z6pIGIVL9ogj4dSDWzFDOLB8YCcyrsMwe4JjL7Ziiw2903ufu97t7J3ZMjx73v7ldXZQdqiyc+zOYXry/nnF7tmXrDYJo3ahDrkkSkjjjsPHp3LzWzCcDbQBwwxd2Xm9nNkfZJwJvA+UA2UARcX30l1y7uzv97exVPfLiWMaecyJ8u708DPdJARI4jc6843B57aWlpnpGREesyjlko5Pzmnyt49vMcrhycxO8uPpm4enqkgYhUPTPLdPe0ytp0Z2w1CYWc+yMPJ7vp9BTuv6C3nlsjIjGhoK8GZSHn3leX8HJGnh5OJiIxp6CvYmUh5+ezFvPql/ncdk4qPz03VSEvIjGloK9CZSHnrpmLmb0wn5+e24Pbz02NdUkiIgr6qlI+5O8c2YOfnKOQF5GaQUFfBcqH/F3n9WDC2Qp5Eak5NKH7GIVCzj2vLFHIi0iNpaA/Bu7hefIzM/O47ZxUhbyI1EgK+mPw/95exbOf53DT6Sn8VB+8ikgNpaA/SpM+WssTH65l3JAk3QwlIjWagv4ovJyey8P/WsmF/U/kt2NOVsiLSI2moD9Cby/fzD2vLuGM1HY8cnl/PbtGRGo8Bf0RmL9uOz95cSF9O7Vi0tWnEl9f//lEpOZTUkVp5eY93DQtg86tG/PMdYNo2lC3IIhI7aCgj0L+rv1cNyWdJvFxTL1hMG2axse6JBGRqOmy9DB2FRVz7ZQF7CsuZebNw+jUukmsSxIROSK6ov8WB0rKuGlqBhu3F/H3a9Lo1aFFrEsSETliuqI/hFDI+elLi8jYsJPHxw1kaNe2sS5JROSo6Ir+EH73Rhb/WraZBy7ozQX9Osa6HBGRoxZV0JvZKDNbZWbZZnZPJe1mZo9F2peY2cDI653N7AMzyzKz5WZ2e1V3oDo8/el6pny2nhuGp3DTGV1jXY6IyDE5bNCbWRzwODAa6ANcaWZ9Kuw2GkiNfI0Hnoy8Xgrc6e69gaHArZUcW6O8s3wzv3tjBaNO6sADF/SOdTkiIscsmiv6wUC2u69z92JgBjCmwj5jgGke9gXQysw6uvsmd/8SwN2/ArKAxCqsv0otydvF7TMW0a9TKx694hTq6a5XEQmAaII+Ecgtt53HN8P6sPuYWTIwAJhf2ZuY2XgzyzCzjK1bt0ZRVtXK37WfG6dm0LZZPE9dk0bj+LjjXoOISHWIJugru6z1I9nHzJoBrwB3uPueyt7E3Se7e5q7pyUkJERRVtXZd7CUG59N50BJGc9cN4iE5g2P6/uLiFSnaII+D+hcbrsTUBDtPmbWgHDIP+/urx59qdXj62mUq7d8xePjBpJ6QvNYlyQiUqWiCfp0INXMUswsHhgLzKmwzxzgmsjsm6HAbnffZOHn9z4NZLn7n6u08iryyLureGfFFn7xvT6c2eP4/iYhInI8HPaGKXcvNbMJwNtAHDDF3Zeb2c2R9knAm8D5QDZQBFwfOXw48ENgqZktirx2n7u/WbXdODqvL8rn8Q/WcuXgzlx3WnKsyxERqRbmXnG4PfbS0tI8IyOjWt9j9ZavGPO3z+ib2JLpNw3RI4dFpFYzs0x3T6usrU6m296Dpdw8PZOmDevzt3EDFPIiEmh1LuHcnXtfXUrOtn389coBtG/RKNYliYhUqzoX9M99sYF/LC7gzvN6MqybHlQmIsFXp4I+u3AvD76RxYieCfz4O91iXY6IyHFRZ4K+tCzEnTMX0zg+jj9e1k+PNxCROqPOPI/+fz9ex+LcXRqXF5E6p05c0a8o2MNf3lvNBf06cmH/E2NdjojIcRX4oC8uDQ/ZtGwcz2/HnBzrckREjrvAD91M+mgtWZv28L8/PJU2TeNjXY6IyHEX6Cv61Vu+4q/vr+F7/Try3ZM6xLocEZGYCGzQl5aF+PnMxTRv1IBfX3RSrMsREYmZwA7dTPlsPYvzdvPYlQNo20zPlxeRuiuQV/S7ioqZ+N4azunVngv7dYx1OSIiMRXIoJ82bwP7isv4+aiehB+JLyJSdwUu6IuKS3nms/Wc3as9vTq0iHU5IiIxF7igfyk9l51FJdwyQs+yERGBgAV9cWmIv3+8jkHJrUlLbhPrckREaoRABf2cxQUU7D7ALSO6x7oUEZEaI6qgN7NRZrbKzLLN7J5K2s3MHou0LzGzgdEeW1VCIWfSR2vp1aE5I3pqkW8Rka8ddh69mcUBjwMjgTwg3czmuPuKcruNBlIjX1l9Y4UAAASmSURBVEOAJ4EhUR5bJYpKykjr0pozUhM000ZEpJxobpgaDGS7+zoAM5sBjAHKh/UYYJqHVxr/wsxamVlHIDmKY6tEs4b1efiyflV9WhGRWi+aoZtEILfcdl7ktWj2ieZYAMxsvJllmFnG1q1boyhLRESiEU3QVzYO4lHuE82x4RfdJ7t7mrunJSRojF1EpKpEM3STB3Qut90JKIhyn/gojhURkWoUzRV9OpBqZilmFg+MBeZU2GcOcE1k9s1QYLe7b4ryWBERqUaHvaJ391IzmwC8DcQBU9x9uZndHGmfBLwJnA9kA0XA9d92bLX0REREKmXhiTI1S1pammdkZMS6DBGRWsPMMt09rbK2QN0ZKyIi36SgFxEJuBo5dGNmW4ENR3l4O2BbFZZTG9TFPkPd7Hdd7DPUzX4faZ+7uHulc9NrZNAfCzPLONQ4VVDVxT5D3ex3Xewz1M1+V2WfNXQjIhJwCnoRkYALYtBPjnUBMVAX+wx1s991sc9QN/tdZX0O3Bi9iIj8pyBe0YuISDkKehGRgAtM0B+vJQtjzcw6m9kHZpZlZsvN7PbI623M7F0zWxP5s3Wsa61qZhZnZgvN7J+R7brQ51ZmNsvMVka+58OC3m8z+2nkZ3uZmb1oZo2C2Gczm2JmhWa2rNxrh+ynmd0bybdVZvbdI3mvQAR9uSULRwN9gCvNrE9sq6o2pcCd7t4bGArcGunrPcBcd08F5ka2g+Z2IKvcdl3o80TgLXfvBfQn3P/A9tvMEoHbgDR3P5nwwxDHEsw+PwuMqvBapf2M/BsfC5wUOeaJSO5FJRBBT7nlDt29GPh6ycLAcfdN7v5l5O9fEf6Hn0i4v1Mju00FLo5NhdXDzDoBFwBPlXs56H1uAZwJPA3g7sXuvouA95vwU3Ubm1l9oAnhNSwC12d3/xjYUeHlQ/VzDDDD3Q+6+3rCTwoeHO17BSXoo16yMEjMLBkYAMwHToisAUDkz/axq6xa/AW4GwiVey3ofe4KbAWeiQxZPWVmTQlwv909H/gTsBHYRHhti3cIcJ8rOFQ/jynjghL0US9ZGBRm1gx4BbjD3ffEup7qZGbfAwrdPTPWtRxn9YGBwJPuPgDYRzCGLA4pMiY9BkgBTgSamtnVsa2qRjimjAtK0Eez3GFgmFkDwiH/vLu/Gnl5i5l1jLR3BApjVV81GA5cZGY5hIflzjaz6QS7zxD+uc5z9/mR7VmEgz/I/T4XWO/uW929BHgVOI1g97m8Q/XzmDIuKEFfZ5YsNDMjPGab5e5/Ltc0B7g28vdrgdePd23Vxd3vdfdO7p5M+Hv7vrtfTYD7DODum4FcM+sZeekcYAXB7vdGYKiZNYn8rJ9D+HOoIPe5vEP1cw4w1swamlkKkAosiPqs7h6IL8JLGa4G1gL3x7qeauzn6YR/ZVsCLIp8nQ+0Jfwp/ZrIn21iXWs19X8E8M/I3wPfZ+AUICPy/X4NaB30fgO/BlYCy4DngIZB7DPwIuHPIUoIX7Hf+G39BO6P5NsqYPSRvJcegSAiEnBBGboREZFDUNCLiAScgl5EJOAU9CIiAaegFxEJOAW9iEjAKehFRALu/wNYsXON/gIRrwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = np.cumsum(pca.explained_variance_ratio_)\n",
    "plt.plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialization complete\n",
      "Iteration 0, inertia 210.81517939508086\n",
      "Iteration 1, inertia 174.7319485592479\n",
      "Iteration 2, inertia 172.1255842873564\n",
      "Iteration 3, inertia 170.9544371946473\n",
      "Iteration 4, inertia 170.47439248151383\n",
      "Iteration 5, inertia 170.19920016091964\n",
      "Iteration 6, inertia 170.03584414279382\n",
      "Iteration 7, inertia 169.9388364297175\n",
      "Iteration 8, inertia 169.87629164183969\n",
      "Iteration 9, inertia 169.8420998619317\n",
      "Iteration 10, inertia 169.82513929958827\n",
      "Iteration 11, inertia 169.81388007000578\n",
      "Iteration 12, inertia 169.80926262653355\n",
      "Converged at iteration 12: center shift 2.861467613817101e-07 within tolerance 3.139215698600419e-07\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "KMeans(max_iter=100, n_clusters=4, n_init=1, verbose=True)"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_clusters = 4\n",
    "km = KMeans(n_clusters=n_clusters, init='k-means++', max_iter=100, n_init=1,\n",
    "                verbose=True)\n",
    "km.fit(X_tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cluster 0: com article netcom stratus hp sun posting don nntp ibm\n",
      "Cluster 1: university windows ca posting host nntp cs know like thanks\n",
      "Cluster 2: people government key don just com article clipper think like\n",
      "Cluster 3: god jesus people bible christian believe don christians say think\n"
     ]
    }
   ],
   "source": [
    "V = pca.inverse_transform(km.cluster_centers_)\n",
    "order_centroids = np.fliplr(np.argsort(V, axis=1))\n",
    "if 1:\n",
    "    for i in range(n_clusters):\n",
    "        print(\"Cluster %d:\" % i, end='')\n",
    "        for ind in order_centroids[i, :10]:\n",
    "            print(' %s' % terms[ind], end='')\n",
    "        print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alt.atheism                              0.000 0.121 0.707 0.171 \n",
      "comp.graphics                            0.000 0.883 0.004 0.113 \n",
      "comp.os.ms-windows.misc                  0.000 0.886 0.004 0.110 \n",
      "comp.sys.ibm.pc.hardware                 0.003 0.852 0.000 0.145 \n",
      "comp.sys.mac.hardware                    0.000 0.916 0.001 0.083 \n",
      "comp.windows.x                           0.000 0.858 0.001 0.141 \n",
      "misc.forsale                             0.017 0.841 0.002 0.139 \n",
      "rec.autos                                0.000 0.518 0.020 0.462 \n",
      "rec.motorcycles                          0.001 0.425 0.012 0.562 \n",
      "rec.sport.baseball                       0.658 0.192 0.005 0.145 \n",
      "rec.sport.hockey                         0.893 0.067 0.002 0.038 \n",
      "sci.crypt                                0.000 0.156 0.052 0.791 \n",
      "sci.electronics                          0.012 0.753 0.000 0.235 \n",
      "sci.med                                  0.084 0.531 0.070 0.315 \n",
      "sci.space                                0.059 0.640 0.036 0.264 \n",
      "soc.religion.christian                   0.000 0.118 0.813 0.068 \n",
      "talk.politics.guns                       0.000 0.127 0.498 0.375 \n",
      "talk.politics.mideast                    0.000 0.146 0.778 0.077 \n",
      "talk.politics.misc                       0.005 0.228 0.381 0.386 \n",
      "talk.religion.misc                       0.002 0.162 0.662 0.174 \n"
     ]
    }
   ],
   "source": [
    "cnt = np.zeros((ncat,n_clusters), dtype=int)\n",
    "for i, cat in enumerate(categories):\n",
    "    print('%-40s ' % cat, end='')\n",
    "    for j in range(n_clusters):\n",
    "        cnt[i,j] = np.sum((labelkm==j) & (labels==i))\n",
    "    cntp = cnt[i,:] / np.sum(cnt[i,:])\n",
    "    for j in range(n_clusters):\n",
    "        print('%5.3f ' % cntp[j], end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of clusters to newsgroup categories\n",
    "\n",
    "The clusters found by k-means were not based on the newsgroup category in which the post came from.  To compare the two, we create a sort of confusion matrix where:\n",
    "\n",
    "`C[i,j] = ` fraction of cluster `j` came from newsgroup `i`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.0015873  0.07790627 0.55480984 0.79090909]\n",
      " [0.01269841 0.58673159 0.00111857 0.        ]\n",
      " [0.97936508 0.22337188 0.0033557  0.        ]\n",
      " [0.00634921 0.11199026 0.44071588 0.20909091]]\n"
     ]
    }
   ],
   "source": [
    "labelkm = km.labels_\n",
    "from sklearn.metrics import confusion_matrix\n",
    "C = confusion_matrix(labels,labelkm)\n",
    "\n",
    "Csum = np.sum(C,axis=0)\n",
    "Cnorm = C / Csum[None,:]\n",
    "print(Cnorm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To interpret this confusion matrix, let's print out the newsgroup names:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['alt.atheism', 'comp.graphics', 'sci.space', 'talk.religion.misc']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset.target_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, returning to the confusion matrix, we see that some discovered clusters lie almost entirely within one of the newsgroup subjects.  This is especially true for `comp.graphics` and `sci.space`.  However, some discovered clusters tend to have entries of both `alt.atheism` and `talk.religon.misc`, whose topics are likely to have a lot of overlap."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print an example of a post that came from a newsgroup that is different from the most common newsgroup in that cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Actual newsgroup: talk.religion.misc\n",
      "Most common newsgroup in cluster:  alt.atheism\n",
      "\n",
      "From: nrp@st-andrews.ac.uk (Norman R. Paterson)\n",
      "Subject: Re: After 2000 years, can we say that Christian Morality is\n",
      "Organization: St. Andrews University, Scotland.\n",
      "Lines: 16\n",
      "\n",
      "In article <1r59na$e81@fido.asd.sgi.com> livesey@solntze.wpd.sgi.com (Jon Livesey) writes:\n",
      ">In article <1993Apr21.141259.12012@st-andrews.ac.uk>, nrp@st-andrews.ac.uk (Norman R. Paterson) writes:\n",
      ">|> In article <1r2m21$8mo@fido.asd.sgi.com> livesey@solntze.wpd.sgi.com (Jon Livesey) writes:\n",
      "...\n",
      ">> Ok, so you don't claim to have an absolute moral system.  Do you claim\n",
      ">> to have an objective one?  I'll assume your answer is \"yes,\" apologies\n",
      ">> if not.\n",
      ">\n",
      ">I've just spent two solid months arguing that no such thing as an\n",
      ">objective moral system exists.\n",
      ">\n",
      ">jon.\n",
      "\n",
      "Apologies, I've not been paying attention.\n",
      "\n",
      "-Norman\n",
      "\n"
     ]
    }
   ],
   "source": [
    "I = np.where((labels==3) & (labelkm == 3))[0]\n",
    "doc_ind = I[3]\n",
    "ind_cluster = labelkm[doc_ind]\n",
    "km_cat = dataset.target_names[np.argmax(Cnorm[:,ind_cluster])]\n",
    "\n",
    "data_ex = dataset.data[doc_ind]\n",
    "true_cat  = dataset.target_names[labels[doc_ind]]\n",
    "print('Actual newsgroup: {0:s}'.format(true_cat))\n",
    "print('Most common newsgroup in cluster:  {0:s}'.format(km_cat))\n",
    "print()\n",
    "print(data_ex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Latent Semantic Analysis\n",
    "\n",
    "Another important tool in document analysis is [latent semantic analysis (LSA)](https://en.wikipedia.org/wiki/Latent_semantic_analysis).  In LSA, we simply compute an SVD of the TF-IDF matrix,\n",
    "\n",
    "    X = U diag(S) V\n",
    "    \n",
    "This is equivalent to performing a PCA on `X`.  If we let `A = U diag(S)` then `X = AV`.  First, we compute the PCs of `X`.  Since `X` is a sparse matrix, it is preferable to use the sparse `svds` method in the `scipy` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse.linalg\n",
    "U1,S1,V1 = scipy.sparse.linalg.svds(X,k=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can next plot the singular values. We see that the first few singular values are significantly larger than the remaining singular values suggesting that the term-document matrix `X` has a low rank structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x18b0992fc18>]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFfJJREFUeJzt3X2QXXV9x/H3Nw9ENgkkgUAICQmIEAkiYMCqoLeDUoby\nYO1YtcNYrTJ22gpFpqJ1OtmZtk47nVqtU9v6hKU1daytPIwMitbLgyCgJIRAIlgaQgiEhDyRRDbZ\n8Osfv3vZzbKbvZc99949d9+vmTP37Lkn5/5+czaf+9vveYqUEpKkcprU6QZIkl49Q1ySSswQl6QS\nM8QlqcQMcUkqMUNckkqsoRCPiKsj4uHadFWrGyVJasyoIR4RS4GPAMuAM4FLIuKkVjdMkjS6Rkbi\nrwfuSyn1pZQOAHcC72ltsyRJjWgkxNcA50fE7IjoAS4GFra2WZKkRkwZbYWU0rqI+BvgdmA3sBI4\n0OqGSZJGF83eOyUi/gp4KqX0z0OWexMWSWpSSinG8u8bPTtlbu31BOC3gBUjNKYrp+XLl3e8DfbP\n/tm/7puKMGo5pea/ImIOsB/4w5TSrkI+XZI0Jg2FeErp7a1uiCSpeV6x2YBKpdLpJrSU/Ss3+zex\nNX1gc8QNRaSitiVJE0FEkNpxYFOSND4Z4pJUYoa4JJWYIS5JJWaIS1KJGeKSVGKGuCSVmCEuSSVm\niEtSiRniklRihrgklZghLkklZohLUokZ4pJUYoa4JJWYIS5JJWaIS1KJNfq0+2siYk1ErI6Ib0bE\nYa1umCRpdKOGeETMBz4OnJ1SOoP8cOX3t7phkqTRNfS0e2AyMD0iXgJ6gE2ta5IkqVGjjsRTSpuA\nvwM2AE8DO1JKP2x1wyRJo2uknDILuBxYBMwHZkTE77a6YZKk0TVSTnkn8ERKaRtARPw38FZgxdAV\ne3t7X56vVCpUKpVCGilJ3aBarVKtVgvdZqSUDr1CxLnA14BzgD7geuCBlNI/DlkvjbYtSdKAiCCl\nFGPZRiM18fuB7wArgYeAAL48lg+VJBVj1JF4wxtyJC5JTWnLSFySNH4Z4pJUYoa4JJWYIS5JJWaI\nS1KJGeKSVGKGuCSVmCEuSSVmiEtSiRniklRihrgklZghLkklZohLUokZ4pJUYoa4JJWYIS5JJWaI\nS1KJGeKSVGKGuCSV2KghHhGnRMTKiHiw9rozIq5qR+MkSYfW1IOSI2ISsBF4c0rpqSHv+aBkSWpC\nJx6U/E7gf4cGuCSpM5oN8fcB/9GKhkiSmjel0RUjYipwGfCpkdbp7e19eb5SqVCpVMbQNEnqLtVq\nlWq1Wug2G66JR8RlwB+mlC4a4X1r4pLUhHbXxD+ApRRJGlcaGolHRA/wJHBSSumFEdZxJC5JTShi\nJN7UKYajNMYQl6QmdOIUQ0nSOGKIS1KJGeKSVGKGuCSVmCEuSSVmiEtSiRniklRihrgklZghLkkl\nZohLUokZ4pJUYoa4JJWYIS5JJWaIS1KJGeKSVGKGuCSVmCEuSSVmiEtSiRniklRiDYV4RBwZEf8Z\nEWsj4pGIeHOrGyZJGt2UBtf7AnBrSum9ETEF6GlhmyRJDRr1afcRcQSwMqX02lHW82n3ktSEdj3t\n/kRga0RcHxEPRsSXI+Lw4VY0wyWpvRopp0wBzgb+KKX0s4j4PPApYPnQFZcv72VS7WuhUqlQqVSK\na6kklVy1WqVarRa6zUbKKccC96aUTqr9fB5wXUrp0iHrpb6+xGGHFdo+SepabSmnpJQ2A09FxCm1\nRRcAjw63bn//WJoiSWpWo2enXAV8MyKmAk8AHx5uJUNcktpr1HJKwxuKSNu2JWbPLmRzktT12nV2\nSsMciUtSexniklRihYb4gQNFbk2SNBpH4pJUYoa4JJWY5RRJKjFH4pJUYo7EJanECg3xPXuK3Jok\naTSFhvj27UVuTZI0mkJDfNu2IrcmSRqNIS5JJWaIS1KJGeKSVGKGuCSVmCEuSSVmiEtSiRniklRi\nhrgklVhDz9iMiPXATuAlYH9K6dxh1kmTJyf27YNJhX41SFJ3KuIZm40+7f4loJJSOuSF9dOnw65d\nMGvWWJokSWpUo2PmaGTdOXMsqUhSOzUa4gm4PSIeiIgrR1rJEJek9mq0nPK2lNIzETGXHOZrU0p3\nD11p+/ZePv95OPlkqFQqVCqVQhsrSWVWrVapVquFbrOhA5sH/YOI5cALKaXPDVmePvaxxKmnwjXX\nFNlESepORRzYHLWcEhE9ETGjNj8duBBYM9y6558Pd901luZIkpox6kg8Ik4Evkuui08BvplS+uth\n1ksbNiTe9CbYvBliTN8tktT9ihiJN11OOURjUkqJxYvhtttgyZJCNitJXast5ZRmnX8+3P2KQ56S\npFYoPMTPOgtWry56q5Kk4RQe4qefDmuGPewpSSqaIS5JJVZ4iB93HPT3w3PPFb1lSdJQhYd4hKNx\nSWqXltw01hCXpPZoSYifcQasWtWKLUuSBmtJiC9bBj/7WSu2LEkarPArNgH6+mD2bNi6FXp6Ctm8\nJHWdcXnFJsC0abB0qSUVSWq1lj0Nc9kyeOCBVm1dkgQtDvGf/7xVW5ckQQtD/LTTYN26Vm1dkgQt\nDPHXvQ4efxwKOm4qSRpGy0L8qKPy1Ztbt7bqEyRJLQvxiDwaf+yxVn2CJKllIQ5wyim5pCJJao2W\nh7gjcUlqnYZDPCImRcSDEXFzo/+mfnBTktQazYzErwYebWbjS5fCQw811yBJUuMaCvGIWABcDHy1\nmY2fdhps3uwZKpLUKo2OxP8e+FOgqbO+J0+GN78Z7r236XZJkhowZbQVIuI3gc0ppVURUQFGvONW\nb2/vy/OVSoVKpcJb3wr33AOXXlpAayWpxKrVKtVqtdBtjnor2oj4LHAF0A8cDswE/jul9MEh66Xh\ntvWDH8Bf/iXceWdhbZakrlDErWibup94RLwDuDaldNkw7w0b4i+8AAsX5vuozJs3lqZKUncZt/cT\nH2zmTPjt34ZvfKPVnyRJE09Lnuwz1P33wwc+kM8Zn9Tyrw1JKodSjMQBzjkHjjgCfvSjdnyaJE0c\nbQnxCPjYx+Bf/qUdnyZJE0dbyikAu3bB4sXw9a/Du99dyEdKUqmVppwCuZxyyy1w7bWOyCWpKG0b\nidf94hdw3nlwxx35snxJmqhKNRKvO/VUuPJKWLGi3Z8sSd2nIyf8XXgh3H57Jz5ZkrpL28spAH19\nMHcuPPkkzJ5dyMdLUumUspwCMG1arovfemsnPl2SuseodzFslT/5E7jiCnjxRfjIRzrVCkkqt46U\nU+pWrYJLLoH162FKx75OJKkzSltOqTvzTDjhBPje9zrZCkkqr47fjurqq+FDH8o3yNqypdOtkaRy\n6XiIv+99+WHKixblR7nt3NnpFklSeXS0Jj7UlVfCa14DX/xiIU2SpHGt7U/2GaUxYw7xbdvgjW/M\nZ6584hP57oeS1K2KCPFxdU7InDnwk5/ARRfB5Mk5zCVJIxtXIQ75bJVbb4W3vCX//PGP50CXJL3S\nuCqnDPbYY/kioPnz882yDHJJ3aYt54lHxLSIuC8iVkbEwxGxfCwf2KhTTsk3ydq6FSoV+Pa3Yd++\ndnyyJJVHQyPxiOhJKe2NiMnAT4CrUkr3D1mn0JF43b59cOON8E//BM8+m5/TOX9+4R8jSW3Xtis2\nU0p7a7PTyHX04tN6BIcdBr/zO/DjH8MHPwi/9mtw/fXQgu8LSSqdRkfik4CfA68F/jGl9Olh1mnJ\nSHyoO++Ea66BJUvgq1+Fww9v+UdKUku07RTDlNJLwFkRcQRwY0ScllJ6dOh6vb29L89XKhUqlcpY\n2jast78d7r4bfv/3c638u9+1vCKpHKrVKtVqtdBtNn12SkT8ObAnpfS5IcvbMhKvSwk++1n40pfg\nK1+Biy9u20dLUiHadXbK0RFxZG3+cOBdwLqxfGgRIuAzn4EbbsgXBV12Gfzf/3W6VZLUXo0c2DwO\n+HFErALuA76fUho3z+S54AJ4+OF8cdA55+SLg1asgOef73TLJKn1xu3FPq/G+vXwrW/BT38Kd9yR\nb3H70Y/C0qUdbZYkDavrboBVpGeegc99Dv7933OQf/KTMHNmp1slSQNK/2SfVjruOPjbv833Kr/3\nXpg3L59j/u1vQ39/p1snScXo2pH4UAcOwM0359H5pk1w7bX5YOiCBZ1umaSJynLKq3THHfDlL8P3\nvw9XXAGXXgrnnmu5RVJ7GeJj9Oyz8IUvwF13wapV8N735vr5W94Ck7q20CRpvDDEC7RjR34s3IoV\neVT+jW/4ZCFJrWWIt8DevXDeeXDMMXDWWfnhzZdd5shcUvEM8RbZti3fy/yxx+Cmm+BXv4Jf/3W4\n7jpYuLDTrZPULQzxNnjppXzx0Pe+l+/TcvzxuWZ+8cXwznd6MFTSq2eIt9nOnfn+LHfckZ8Des89\ncOqpuYb+B3+QR+mzZ3e6lZLKwhDvsN274ZFH4JZb8uX+W7bAe94Db3gDvO51cPrpsGiR9XRJwzPE\nx5nnn89PHXrqqVxPX7MGtm+Hd7wjh/vSpfCmN8HUqZ1uqaTxwBAvgR078pWiP/whrF4Njz6aSy4X\nXACXXAJnnw0nnZQfQydpYjHES2j//nxzrttuy3X1NWtg48Z8wHTZMrjwQjjzTDjtNB89J3U7Q7xL\n7NuXD5jedVc+aLp6dS7HnHginHFGnpYuzc8VPekkyzFStzDEu9i+fbBuXQ70hx6CtWvzzxs35nBf\nsgRe//oc7vWAf81rOt1qSc0wxCegF1+Exx/Pgb52bT47Zs0aeOKJHOpLluSrTM87L5/yePTRnh0j\njVeGuF7W1wcPPpgDvlrN8xs3wq5dMHcuHHtsvsf6smVw/vn5/PZjjoFp0zrdcmniakuIR8QC4Abg\nWOAl4CsppX8YZj1DfBx68cV8/vrmzfD00/nq07vvzjX4557L4f7GN8JRR+Vp4UI44YR8fvsJJ+Qv\nAG8EJrVGu0J8HjAvpbQqImYAPwcuTymtG7KeIV4yKcEvf5lPe9y+HbZuzee4P/kkbNiQpz17cpjX\np3q41+cXLHA0L71aHSmnRMSNwBdTSj8astwQ70K7d+dg37Dh4HCvz2/aBHPmHBzu8+ZBTw/Mn59P\nnTz++Dzinzy5072Rxpe2h3hELAaqwOkppd1D3jPEJ6ADB/LDNQYH/LPP5hH8M8/kEs7TT+c7Qx5z\nTK7Fn3xyDv4TT8w3EDvyyFzSmT/f0o0mlraGeK2UUgX+IqV00zDvG+Ia0f79edT+yCM58J9/Htav\nzyP9bdvyk5V27cpBXh/Bv/a1+f4zs2bBEUcMTEcdBTNmGPgqvyJCfEqDHzQF+A7wb8MFeF1vb+/L\n85VKhUqlMpa2qYtMnZpLLosWjbzO3r056OvTunX5fu67dsELL+TXnTvzF0B/fz7bZtGiHPhz5+bp\n6KNfOT9njqdZanyoVqtUq9VCt9nQSDwibgC2ppQ+cYh1HImrbfbuzeWaDRtyuWbLljxt3TowX592\n7cr3q5k7N9frFy/OPx95ZJ5mzcrTnDl5mj07v3rxlFqtXWenvA24E3gYSLXpz1JKtw1ZzxDXuNTf\nn0fvW7bk4F+/Pt+YbOfOPO3Ykaft23Nppz5NmTIQ6IPDfbhlg+ePPNKRvxrjxT5Si6SUR/vbtg2E\n++CQH26+/rp7d67dNxv+jv4nHkNcGof6+/MIv5ngr0+TJjUW/rNmHVwKmjUrfwF4sLdcDHGpi6SU\nH8rdyOi/XgbauTMv37kzPw92cKgPNw0N/sHT4Yf7JdBuhrikl7344sE1/qHTod7bsSP/BTFrVj53\nf/r0/HrEEQMHgOvzM2fmUzzrU09P/gLo6Rn4opg50+MCjTDEJRWmry+H+e7dedq1a2CqHwTeuXPg\n/T17Bl5/9at8DKH+RbFnTw79adPyRV4LF+agHzz19Ax8EcycOTAN/rk+362lIkNc0rjU35/Dv68v\nnxG0aVMO+qHT7t35GoAXXjh4fuiy/v5Xhnv958Gvh1o2+IKxnp7x8aVgiEuaEPbvPzjU638NNLqs\nfsFY/aKxvr5cMhpcFqqXkQZPM2a8cn7w69BlzR5XaNsVm5LUSVOnDpyhU4T+/oFy0NAvgL1783uD\np23bBuaHKyfVl+3bd3CZaKTwnz49P8ClCIa4pAlnypSBA7ZFOnBg5LAfGvxF3dXTcookdUgR5RRP\nApKkEjPEJanEDHFJKjFDXJJKzBCXpBIzxCWpxAxxSSoxQ1ySSswQl6QSM8QlqcRGDfGI+FpEbI6I\n1e1okCSpcY2MxK8HfqPVDRnPqtVqp5vQUvav3OzfxDZqiKeU7ga2t6Et41a3/xLZv3KzfxObNXFJ\nKjFDXJJKrKH7iUfEIuCWlNIZh1jHm4lLUpPa9Xi2qE0ta4gkqXmNnGK4ArgHOCUiNkTEh1vfLElS\nIwp7PJskqf3GfGAzIi6KiHUR8VhEXFdEozotItZHxEMRsTIi7q8tmx0RP4iIX0TE9yOi4Eests5w\nF2wdqj8R8emIeDwi1kbEhZ1pdeNG6N/yiNgYEQ/WposGvVea/kXEgoj4n4h4JCIejoirasu7Yv8N\n07+P15Z3y/6bFhH31bLk4YhYXlte3P5LKb3qifwl8EtgETAVWAUsGcs2x8MEPAHMHrLsb4BP1uav\nA/660+1soj/nAWcCq0frD3AasJJ8vGRxbf9Gp/vwKvq3HPjEMOu+vkz9A+YBZ9bmZwC/AJZ0y/47\nRP+6Yv/V2txTe50M/BQ4t8j9N9aR+LnA4ymlJ1NK+4FvAZePcZvjQfDKv1IuB/61Nv+vwLvb2qIx\nSMNfsDVSfy4DvpVS6k8prQceJ+/ncWuE/sHwB+Mvp0T9Syk9m1JaVZvfDawFFtAl+2+E/h1fe7v0\n+w8gpbS3NjuNHM6JAvffWEP8eOCpQT9vZGAHlFkCbo+IByLio7Vlx6aUNkP+xQOO6VjrinHMCP0Z\nuk+fprz79I8jYlVEfHXQn6ul7V9ELCb/xfFTRv597Ib+3Vdb1BX7LyImRcRK4Fng9pTSAxS4/7zY\nZ3hvSymdDVwM/FFEnE8O9sG67Yhwt/XnS8BJKaUzyf95/q7D7RmTiJgBfAe4ujZi7arfx2H61zX7\nL6X0UkrpLPJfUOdGxFIK3H9jDfGngRMG/bygtqzUUkrP1F63ADeS/5zZHBHHAkTEPOC5zrWwECP1\n52lg4aD1SrlPU0pbUq3ICHyFgT9JS9e/iJhCDrh/SyndVFvcNftvuP510/6rSyntAqrARRS4/8Ya\n4g8AJ0fEoog4DHg/cPMYt9lREdFTGxUQEdOBC4GHyf36UG213wNuGnYD49fQC7ZG6s/NwPsj4rCI\nOBE4Gbi/XY0cg4P6V/uPUfceYE1tvoz9+zrwaErpC4OWddP+e0X/umX/RcTR9VJQRBwOvItc9y9u\n/xVw5PUi8hHlx4FPdfpIcAH9OZF8ls1Kcnh/qrZ8DvDDWl9/AMzqdFub6NMKYBPQB2wAPgzMHqk/\nwKfJR8XXAhd2uv2vsn83AKtr+/JGcg2ydP0D3gYcGPQ7+WDt/9yIv49d0r9u2X9vqPVpVa0/n6kt\nL2z/ebGPJJWYBzYlqcQMcUkqMUNckkrMEJekEjPEJanEDHFJKjFDXJJKzBCXpBL7fw5Hz7lpq1Ix\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x18b09895748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(S1[::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the low rank structure of `X`, we can write \n",
    "\n",
    "     X[n,i] = \\sum_k A[n,k] V[k,i]\n",
    "     \n",
    "where the sum is over some relatively small number of components.  There are two uses for this representation:\n",
    "* Word and document embeddings:  A[n,:] provides a low-dimensional vector representation of each document.  This is useful pre-processing step in many natural processing (NLP) methods.  This type of representation is closely related to an important topic of *word embeddings* and *document embeddings*.\n",
    "* Topic modeling:  One interpretation of the PCA is that each PC `k` represents some common *topic* in the corpus.  Then, `A[n,k] =` the component of topic `k` in document `n` and `V[k,i]` represents the occurance of word `i` in topic `k`.\n",
    "\n",
    "To get an idea of the words within each PC, we print the words for the largest components in the first 5 PCs.  On a small corpus like 20 newsgroups, the PCs in this case are not very useful.  But, the technique can yield more useful results in larger corpi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PC 0: muenchen does dan bockamp targa colour shall ch uk man\n",
      "PC 1: phigs government uci 3d khomeini colorado uni rh rayshade screen\n",
      "PC 2: clarke asimov wesleyan values vga physics fl ed pluto tyre\n",
      "PC 3: ericsson color wesleyan point program convenient boeing scott targa jpeg\n",
      "PC 4: thanks muenchen format mac earth pluto color uci true ether\n"
     ]
    }
   ],
   "source": [
    "V1sort = np.abs(V1).argsort()[:, ::-1]\n",
    "for i in range(5):\n",
    "    print(\"PC %d:\" % i, end='')\n",
    "    for ind in V1sort[i, :10]:\n",
    "        print(' %s' % terms[ind], end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3387, 300)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialization complete\n",
      "Iteration 0, inertia 86.27353297947184\n",
      "Iteration 1, inertia 53.17605492696268\n",
      "Iteration 2, inertia 52.017132158220114\n",
      "Iteration 3, inertia 51.32589878421361\n",
      "Iteration 4, inertia 50.99974695298535\n",
      "Iteration 5, inertia 50.85708605986888\n",
      "Iteration 6, inertia 50.80220190835376\n",
      "Iteration 7, inertia 50.758690599261435\n",
      "Iteration 8, inertia 50.73585498021275\n",
      "Iteration 9, inertia 50.71972442121352\n",
      "Iteration 10, inertia 50.713337634106075\n",
      "Iteration 11, inertia 50.71068661044915\n",
      "Iteration 12, inertia 50.706919021802136\n",
      "Iteration 13, inertia 50.70432669115986\n",
      "Converged at iteration 13: center shift 2.474934615666404e-07 within tolerance 6.349559876836782e-07\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "KMeans(max_iter=100, n_clusters=4, n_init=1, verbose=True)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "km = KMeans(n_clusters=true_k, init='k-means++', max_iter=100, n_init=1,\n",
    "                verbose=True)\n",
    "km.fit(X_tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "V_tr =  km.cluster_centers_\n",
    "V_tr.shape\n",
    "V = pca.inverse_transform(V_tr)\n",
    "\n",
    "order_centroids = np.fliplr(np.argsort(V, axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cluster 0: henry space toronto access digex nasa zoo pat spencer alaska\n",
      "Cluster 1: god com don people space graphics university know just like\n",
      "Cluster 2: sandvik kent apple newton com alink ksand jesus net cookamunga\n",
      "Cluster 3: sgi livesey keith wpd solntze jon caltech com morality moral\n"
     ]
    }
   ],
   "source": [
    "for i in range(true_k):\n",
    "    print(\"Cluster %d:\" % i, end='')\n",
    "    for ind in order_centroids[i, :10]:\n",
    "        print(' %s' % terms[ind], end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}