konniam_chan_p2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project 2: Topic Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project, you'll work with text data from newsgroup postings on a variety of topics. You'll train classifiers to distinguish between the topics based on the text of the posts. Whereas with digit classification, the input is relatively dense: a 28x28 matrix of pixels, many of which are non-zero, here we'll represent each document with a \"bag-of-words\" model. As you'll see, this makes the feature representation quite sparse -- only a few words of the total vocabulary are active in any given document. The bag-of-words assumption here is that the label depends only on the words; their order is not important.\n",
    "\n",
    "The SK-learn documentation on feature extraction will prove useful:\n",
    "http://scikit-learn.org/stable/modules/feature_extraction.html\n",
    "\n",
    "Each problem can be addressed succinctly with the included packages -- please don't add any more. Grading will be based on writing clean, commented code, along with a few short answers.\n",
    "\n",
    "As always, you're welcome to work on the project in groups and discuss ideas on the course wall, but please prepare your own write-up and write your own code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This tells matplotlib not to try opening a new window for each plot.\n",
    "%matplotlib inline\n",
    "\n",
    "# General libraries.\n",
    "import re\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# SK-learn libraries for learning.\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.naive_bayes import BernoulliNB\n",
    "from sklearn.naive_bayes import MultinomialNB\n",
    "from sklearn.grid_search import GridSearchCV\n",
    "\n",
    "# SK-learn libraries for evaluation.\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import classification_report\n",
    "\n",
    "# SK-learn library for importing the newsgroup data.\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "\n",
    "# SK-learn libraries for feature extraction from text.\n",
    "from sklearn.feature_extraction.text import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the data, stripping out metadata so that we learn classifiers that only use textual features. By default, newsgroups data is split into train and test sets. We further split the test so we have a dev set. Note that we specify 4 categories to use for this project. If you remove the categories argument from the fetch function, you'll get all 20 categories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training label shape: (2034,)\n",
      "test label shape: (677,)\n",
      "dev label shape: (676,)\n",
      "labels names: ['alt.atheism', 'comp.graphics', 'sci.space', 'talk.religion.misc']\n"
     ]
    }
   ],
   "source": [
    "categories = ['alt.atheism', 'talk.religion.misc', 'comp.graphics', 'sci.space']\n",
    "newsgroups_train = fetch_20newsgroups(subset='train',\n",
    "                                      remove=('headers', 'footers', 'quotes'),\n",
    "                                      categories=categories)\n",
    "newsgroups_test = fetch_20newsgroups(subset='test',\n",
    "                                     remove=('headers', 'footers', 'quotes'),\n",
    "                                     categories=categories)\n",
    "\n",
    "num_test = len(newsgroups_test.target)\n",
    "test_data, test_labels = newsgroups_test.data[num_test/2:], newsgroups_test.target[num_test/2:]\n",
    "dev_data, dev_labels = newsgroups_test.data[:num_test/2], newsgroups_test.target[:num_test/2]\n",
    "train_data, train_labels = newsgroups_train.data, newsgroups_train.target\n",
    "\n",
    "print 'training label shape:', train_labels.shape\n",
    "print 'test label shape:', test_labels.shape\n",
    "print 'dev label shape:', dev_labels.shape\n",
    "print 'labels names:', newsgroups_train.target_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(1) For each of the first 5 training examples, print the text of the message along with the label."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Category: 1\n",
      "Category name: comp.graphics\n",
      "Message below:\n",
      "Hi,\n",
      "\n",
      "I've noticed that if you only save a model (with all your mapping planes\n",
      "positioned carefully) to a .3DS file that when you reload it after restarting\n",
      "3DS, they are given a default position and orientation.  But if you save\n",
      "to a .PRJ file their positions/orientation are preserved.  Does anyone\n",
      "know why this information is not stored in the .3DS file?  Nothing is\n",
      "explicitly said in the manual about saving texture rules in the .PRJ file. \n",
      "I'd like to be able to read the texture rule information, does anyone have \n",
      "the format for the .PRJ file?\n",
      "\n",
      "Is the .CEL file format available from somewhere?\n",
      "\n",
      "Rych\n",
      "\n",
      "Category: 3\n",
      "Category name: talk.religion.misc\n",
      "Message below:\n",
      "\n",
      "\n",
      "Seems to be, barring evidence to the contrary, that Koresh was simply\n",
      "another deranged fanatic who thought it neccessary to take a whole bunch of\n",
      "folks with him, children and all, to satisfy his delusional mania. Jim\n",
      "Jones, circa 1993.\n",
      "\n",
      "\n",
      "Nope - fruitcakes like Koresh have been demonstrating such evil corruption\n",
      "for centuries.\n",
      "\n",
      "Category: 2\n",
      "Category name: sci.space\n",
      "Message below:\n",
      "\n",
      " >In article <1993Apr19.020359.26996@sq.sq.com>, msb@sq.sq.com (Mark Brader) \n",
      "\n",
      "MB>                                                             So the\n",
      "MB> 1970 figure seems unlikely to actually be anything but a perijove.\n",
      "\n",
      "JG>Sorry, _perijoves_...I'm not used to talking this language.\n",
      "\n",
      "Couldn't we just say periapsis or apoapsis?\n",
      "\n",
      " \n",
      "\n",
      "Category: 0\n",
      "Category name: alt.atheism\n",
      "Message below:\n",
      "I have a request for those who would like to see Charley Wingate\n",
      "respond to the \"Charley Challenges\" (and judging from my e-mail, there\n",
      "appear to be quite a few of you.)  \n",
      "\n",
      "It is clear that Mr. Wingate intends to continue to post tangential or\n",
      "unrelated articles while ingoring the Challenges themselves.  Between\n",
      "the last two re-postings of the Challenges, I noted perhaps a dozen or\n",
      "more posts by Mr. Wingate, none of which answered a single Challenge.  \n",
      "\n",
      "It seems unmistakable to me that Mr. Wingate hopes that the questions\n",
      "will just go away, and he is doing his level best to change the\n",
      "subject.  Given that this seems a rather common net.theist tactic, I\n",
      "would like to suggest that we impress upon him our desire for answers,\n",
      "in the following manner:\n",
      "\n",
      "1. Ignore any future articles by Mr. Wingate that do not address the\n",
      "Challenges, until he answers them or explictly announces that he\n",
      "refuses to do so.\n",
      "\n",
      "--or--\n",
      "\n",
      "2. If you must respond to one of his articles, include within it\n",
      "something similar to the following:\n",
      "\n",
      "    \"Please answer the questions posed to you in the Charley Challenges.\"\n",
      "\n",
      "Really, I'm not looking to humiliate anyone here, I just want some\n",
      "honest answers.  You wouldn't think that honesty would be too much to\n",
      "ask from a devout Christian, would you?  \n",
      "\n",
      "Nevermind, that was a rhetorical question.\n",
      "\n",
      "Category: 2\n",
      "Category name: sci.space\n",
      "Message below:\n",
      "AW&ST  had a brief blurb on a Manned Lunar Exploration confernce\n",
      "May 7th  at Crystal City Virginia, under the auspices of AIAA.\n",
      "\n",
      "Does anyone know more about this?  How much, to attend????\n",
      "\n",
      "Anyone want to go?\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def P1(num_examples=5):\n",
    "### STUDENT START ###\n",
    "    for message in range(num_examples):\n",
    "        print (\"Category: {}\".format(train_labels[message]))\n",
    "        print (\"Category name: {}\".format(newsgroups_train.target_names[train_labels[message]]))\n",
    "        print (\"Message below:\")\n",
    "        print (train_data[message] + \"\\n\")\n",
    "        \n",
    "### STUDENT END ###\n",
    "P1()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(2) Use CountVectorizer to turn the raw training text into feature vectors. You should use the fit_transform function, which makes 2 passes through the data: first it computes the vocabulary (\"fit\"), second it converts the raw text into feature vectors using the vocabulary (\"transform\").\n",
    "\n",
    "The vectorizer has a lot of options. To get familiar with some of them, write code to answer these questions:\n",
    "\n",
    "a. The output of the transform (also of fit_transform) is a sparse matrix: http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.csr_matrix.html. What is the size of the vocabulary? What is the average number of non-zero features per example? What fraction of the entries in the matrix are non-zero? Hint: use \"nnz\" and \"shape\" attributes.\n",
    "\n",
    "b. What are the 0th and last feature strings (in alphabetical order)? Hint: use the vectorizer's get_feature_names function.\n",
    "\n",
    "c. Specify your own vocabulary with 4 words: [\"atheism\", \"graphics\", \"space\", \"religion\"]. Confirm the training vectors are appropriately shaped. Now what's the average number of non-zero features per example?\n",
    "\n",
    "d. Instead of extracting unigram word features, use \"analyzer\" and \"ngram_range\" to extract bigram and trigram character features. What size vocabulary does this yield?\n",
    "\n",
    "e. Use the \"min_df\" argument to prune words that appear in fewer than 10 documents. What size vocabulary does this yield?\n",
    "\n",
    "f. Using the standard CountVectorizer, what fraction of the words in the dev data are missing from the vocabulary? Hint: build a vocabulary for both train and dev and look at the size of the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Part a:\n",
      "  The size of the vocabulary is 26879.\n",
      "  The average number of nonzero features is 96.71.\n",
      "  The fraction of nonzero entries is 0.00360.\n",
      "Part b:\n",
      "  The first feature is 00.\n",
      "  The last feature is zyxel.\n",
      "Part c:\n",
      "  The shape of the training matrix is (2034, 4).\n",
      "  The average number of nonzero features is now 0.27.\n",
      "Part d:\n",
      "  The size of the vocabulary is 35478.\n",
      "Part e:\n",
      "  The size of the vocabulary is 3064.\n",
      "Part f:\n",
      "  The fraction of words present in dev but not in train is 0.248.\n"
     ]
    }
   ],
   "source": [
    "def P2():\n",
    "### STUDENT START ###\n",
    "# Part a\n",
    "    vectorizer = CountVectorizer()\n",
    "    # Fit and transform with training data\n",
    "    doc_matrix = vectorizer.fit_transform(train_data)\n",
    "    print (\"Part a:\")\n",
    "    print (\"  The size of the vocabulary is {}.\".format(doc_matrix.shape[1]))\n",
    "    # Average number of nonzero features per sample (get number of nonzero entries per example, then average)\n",
    "    print (\"  The average number of nonzero features is {:.2f}.\".format(np.mean(doc_matrix.getnnz(axis=1))))\n",
    "    # What fraction of the entries are non-zero?\n",
    "    print (\"  The fraction of nonzero entries is {:.5f}.\".format(1.0*doc_matrix.nnz/(doc_matrix.shape[0]*doc_matrix.shape[1])))\n",
    "\n",
    "# Part b\n",
    "    print (\"Part b:\")\n",
    "    print (\"  The first feature is {}.\".format(vectorizer.get_feature_names()[0]))\n",
    "    print (\"  The last feature is {}.\".format(vectorizer.get_feature_names()[-1]))\n",
    "\n",
    "# Part c\n",
    "    vocab = [\"atheism\", \"graphics\", \"space\", \"religion\"]\n",
    "    vectorizer = CountVectorizer(vocabulary=vocab)\n",
    "    doc_matrix = vectorizer.fit_transform(train_data)\n",
    "    print (\"Part c:\")\n",
    "    print (\"  The shape of the training matrix is {}.\".format(doc_matrix.shape))\n",
    "    print (\"  The average number of nonzero features is now {:.2f}.\".format(np.mean(doc_matrix.getnnz(axis=1))))\n",
    "    \n",
    "# Part d: use bigram and trigram character features\n",
    "    vectorizer = CountVectorizer(ngram_range=(2,3), analyzer=\"char\")\n",
    "    doc_matrix = vectorizer.fit_transform(train_data)\n",
    "    print (\"Part d:\")\n",
    "    print (\"  The size of the vocabulary is {}.\".format(doc_matrix.shape[1]))\n",
    "    \n",
    "# Part e: mininum document frequency of 10\n",
    "    vectorizer = CountVectorizer(min_df=10)\n",
    "    doc_matrix = vectorizer.fit_transform(train_data)\n",
    "    print (\"Part e:\")\n",
    "    print (\"  The size of the vocabulary is {}.\".format(doc_matrix.shape[1]))\n",
    "\n",
    "# Part f\n",
    "    # Train data vocab\n",
    "    vectorizer_train = CountVectorizer()\n",
    "    vectorizer_train.fit(train_data)\n",
    "    features_train = vectorizer_train.get_feature_names()\n",
    "    # Dev data vocab\n",
    "    vectorizer_dev = CountVectorizer()\n",
    "    vectorizer_dev.fit(dev_data)\n",
    "    features_dev = vectorizer_dev.get_feature_names()\n",
    "    # Intersection of sets\n",
    "    common_words = set(features_dev) & set(features_train)\n",
    "    print (\"Part f:\")\n",
    "    print (\"  The fraction of words present in dev but not in train is {:.3f}.\".format(1.0 - 1.0*len(common_words)/len(features_dev)))\n",
    "    \n",
    "### STUDENT END ###\n",
    "P2()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(3) Use the default CountVectorizer options and report the f1 score (use metrics.f1_score) for a k nearest neighbors classifier; find the optimal value for k. Also fit a Multinomial Naive Bayes model and find the optimal value for alpha. Finally, fit a logistic regression model and find the optimal value for the regularization strength C using l2 regularization. A few questions:\n",
    "\n",
    "a. Why doesn't nearest neighbors work well for this problem?\n",
    "\n",
    "b. Any ideas why logistic regression doesn't work as well as Naive Bayes?\n",
    "\n",
    "c. Logistic regression estimates a weight vector for each class, which you can access with the coef\\_ attribute. Output the sum of the squared weight values for each class for each setting of the C parameter. Briefly explain the relationship between the sum and the value of C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kNN Model:\n",
      "  The optimal k is 96.\n",
      "  Accuray for the optimal k is 0.444.\n",
      "  The macro f1 score on the dev data is 0.434.\n",
      "Multinomial NB Model:\n",
      "  The optimal alpha is 0.02.\n",
      "  Accuray for the optimal alpha is 0.829.\n",
      "  The macro f1 score on the dev data is 0.765.\n",
      "Logistic Regression Model:\n",
      "  The optimal C is 0.17.\n",
      "  Accuray for the optimal C is 0.777.\n",
      "  The macro f1 score on the dev data is 0.689.\n",
      "C: 0.1\n",
      "  Sum of Squared Weights: [ 27.13  24.66  27.46  23.02]\n",
      "C: 0.2\n",
      "  Sum of Squared Weights: [ 49.75  42.74  49.32  42.67]\n",
      "C: 0.3\n",
      "  Sum of Squared Weights: [ 69.28  57.89  67.91  59.78]\n",
      "C: 0.4\n",
      "  Sum of Squared Weights: [ 86.74  71.11  84.26  75.05]\n",
      "C: 0.5\n",
      "  Sum of Squared Weights: [ 102.6    83.09   99.02   89.01]\n",
      "C: 0.6\n",
      "  Sum of Squared Weights: [ 117.23   94.08  112.53  101.84]\n",
      "C: 0.7\n",
      "  Sum of Squared Weights: [ 130.93  104.18  125.02  113.83]\n",
      "C: 0.8\n",
      "  Sum of Squared Weights: [ 143.61  113.56  136.65  125.06]\n",
      "C: 0.9\n",
      "  Sum of Squared Weights: [ 155.58  122.48  147.63  135.58]\n",
      "C: 1.0\n",
      "  Sum of Squared Weights: [ 166.97  130.88  158.    145.72]\n"
     ]
    }
   ],
   "source": [
    "def P3():\n",
    "### STUDENT START ###\n",
    "# Use default CountVectorizer options\n",
    "    vectorizer = CountVectorizer()\n",
    "    doc_matrix_train = vectorizer.fit_transform(train_data)\n",
    "    doc_matrix_dev = vectorizer.transform(dev_data)\n",
    "\n",
    "# kNN model\n",
    "    # Find optimal k first, with the default accuracy scoring option\n",
    "    ks = {'n_neighbors':[a for a in range(1,101)]}\n",
    "    knn_grid = GridSearchCV(KNeighborsClassifier(), ks)\n",
    "    knn_grid.fit(doc_matrix_train, train_labels)\n",
    "    print (\"kNN Model:\")\n",
    "    print (\"  The optimal k is {}.\".format(knn_grid.best_params_['n_neighbors']))\n",
    "    print (\"  Accuray for the optimal k is {:.3f}.\".format(knn_grid.best_score_))\n",
    "    # Compute macro f1 score on the dev data\n",
    "    knn_f1score = metrics.f1_score(dev_labels, knn_grid.predict(doc_matrix_dev), average=\"macro\")\n",
    "    print (\"  The macro f1 score on the dev data is {:.3f}.\".format(knn_f1score))\n",
    "\n",
    "# NB model\n",
    "    # Find optimal alpha first, with the default accuracy scoring option\n",
    "    alphas = {\"alpha\": [a/100. for a in range(1, 101)]}\n",
    "    nb_grid = GridSearchCV(MultinomialNB(), alphas)    \n",
    "    nb_grid.fit(doc_matrix_train, train_labels)\n",
    "    print (\"Multinomial NB Model:\")\n",
    "    print (\"  The optimal alpha is {}.\".format(nb_grid.best_params_['alpha']))\n",
    "    print (\"  Accuray for the optimal alpha is {:.3f}.\".format(nb_grid.best_score_))\n",
    "    # Compute macro f1 score on the dev data\n",
    "    nb_f1score = metrics.f1_score(dev_labels, nb_grid.predict(doc_matrix_dev), average=\"macro\")\n",
    "    print (\"  The macro f1 score on the dev data is {:.3f}.\".format(nb_f1score))\n",
    "\n",
    "# Logistic regression model\n",
    "    # Find optimal C first, with the default accuracy scoring option\n",
    "    Cs = {\"C\": [a/100. for a in range(1, 101)]}\n",
    "    lr_grid = GridSearchCV(LogisticRegression(), Cs)    \n",
    "    lr_grid.fit(doc_matrix_train, train_labels)\n",
    "    print (\"Logistic Regression Model:\")\n",
    "    print (\"  The optimal C is {}.\".format(lr_grid.best_params_['C']))\n",
    "    print (\"  Accuray for the optimal C is {:.3f}.\".format(lr_grid.best_score_))\n",
    "    # Compute macro f1 score on the dev data\n",
    "    lr_f1score = metrics.f1_score(dev_labels, lr_grid.predict(doc_matrix_dev), average=\"macro\")\n",
    "    print (\"  The macro f1 score on the dev data is {:.3f}.\".format(lr_f1score))\n",
    "\n",
    "# Weights for each class, given C\n",
    "    # Use 10 C values, from 0.1 to 1.0\n",
    "    Cs = [a/10. for a in range(1,11)]\n",
    "    # Iterate and calculate sum of weights for each C value\n",
    "    for val in Cs:\n",
    "        lr_model = LogisticRegression(C=val)\n",
    "        lr_model.fit(doc_matrix_train, train_labels)        \n",
    "        print (\"C: {}\".format(val))\n",
    "        # Square coefficients, then sum for each class. Round to 2 decimal points.\n",
    "        print (\"  Sum of Squared Weights: {}\".format(np.round(np.square(lr_model.coef_).sum(axis=1),2)))\n",
    "                \n",
    "### STUDENT END ###\n",
    "P3()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ANSWER:  \n",
    "Part a:  \n",
    "kNN model does not work well here because of the high dimension of the feature space (~27,000) and the sparse nature of our features. In general, kNN does not work well in high dimensions because all points are far apart and the training set does not give good sampling of the feature space. Because this matrix is also very sparse (the vocabulary is huge, but each individual document does not contain most words), kNN has difficulty distinguishing one document from another.  \n",
    "\n",
    "Part b:  \n",
    "Logistic regression assumes that log(odds(class k)) is linear in the features for all k. It is possible that linearity is not a good assumption for document classification. In the Naive Bayes case, we assume conditional independence of the features given the class. This assumption might be more reasonable, leading to better results from the NB model.\n",
    "\n",
    "Part c:  \n",
    "As C becomes smaller (stronger regularization), the squared weights also become smaller. Regularization penalizes large coefficients and pushes them towards zero, in order to reduce overfitting and increase parsimony."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(4) Train a logistic regression model. Find the 5 features with the largest weights for each label -- 20 features in total. Create a table with 20 rows and 4 columns that shows the weight for each of these features for each of the labels. Create the table again with bigram features. Any surprising features in this table?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unigram results:\n",
      "  Size of vocab: 26879\n",
      "  Accuracy: 0.701183431953\n",
      "[u'atheism' u'religion' u'bobby' u'atheists' u'deletion' u'computer' u'3d'\n",
      " u'file' u'image' u'graphics' u'spacecraft' u'launch' u'nasa' u'orbit'\n",
      " u'space' u'order' u'fbi' u'blood' u'christian' u'christians']\n",
      "[[ 0.94 -0.41 -0.45 -0.43]\n",
      " [ 0.95 -0.62 -0.79 -0.06]\n",
      " [ 0.99 -0.22 -0.34 -0.46]\n",
      " [ 1.03 -0.1  -0.32 -0.83]\n",
      " [ 1.13 -0.4  -0.42 -0.4 ]\n",
      " [ 0.14  0.98 -0.68 -0.49]\n",
      " [-0.36  1.13 -0.7  -0.38]\n",
      " [-0.33  1.27 -0.81 -0.63]\n",
      " [-0.58  1.35 -0.83 -0.47]\n",
      " [-0.76  1.94 -1.34 -0.76]\n",
      " [-0.36 -0.39  0.92 -0.38]\n",
      " [-0.47 -0.47  0.94 -0.33]\n",
      " [-0.57 -0.48  1.01 -0.47]\n",
      " [-0.41 -0.67  1.23 -0.63]\n",
      " [-1.26 -1.32  2.16 -1.17]\n",
      " [-0.79 -0.08 -0.15  0.91]\n",
      " [-0.31 -0.27 -0.45  0.91]\n",
      " [-0.53 -0.11 -0.32  1.05]\n",
      " [-0.61 -0.42 -0.27  1.12]\n",
      " [-0.74 -0.41 -0.53  1.15]]\n",
      "\n",
      "\n",
      "Bigram results:\n",
      "  Size of vocab: 194891\n",
      "  Accuracy: 0.627218934911\n",
      "[u'you are' u'look up' u'cheers kent' u'was just' u'claim that'\n",
      " u'out there' u'is there' u'in advance' u'comp graphics' u'looking for'\n",
      " u'sherzer methodology' u'and such' u'sci space' u'the moon' u'the space'\n",
      " u'is strength' u'the fbi' u'compuserve com' u'such lunacy' u'ignorance is']\n",
      "[[ 0.57 -0.32 -0.58  0.  ]\n",
      " [ 0.63 -0.24 -0.2  -0.17]\n",
      " [ 0.65 -0.88 -0.82  0.6 ]\n",
      " [ 0.68 -0.19 -0.2  -0.3 ]\n",
      " [ 0.77 -0.26 -0.35 -0.2 ]\n",
      " [-0.32  0.9  -0.58 -0.33]\n",
      " [-0.43  0.91 -0.57 -0.34]\n",
      " [-0.55  0.97 -0.53 -0.51]\n",
      " [-0.38  1.04 -0.47 -0.4 ]\n",
      " [-0.76  1.32 -0.61 -0.7 ]\n",
      " [-0.16 -0.23  0.68 -0.16]\n",
      " [-0.24 -0.4   0.69 -0.26]\n",
      " [-0.32 -0.39  0.74 -0.27]\n",
      " [-0.41 -0.58  0.95 -0.24]\n",
      " [-0.31 -0.65  1.03 -0.32]\n",
      " [-0.18 -0.22 -0.2   0.65]\n",
      " [-0.17 -0.26 -0.36  0.65]\n",
      " [-0.13 -0.21 -0.2   0.7 ]\n",
      " [-0.15 -0.22 -0.19  0.71]\n",
      " [-0.23 -0.25 -0.21  0.75]]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def P4():\n",
    "### STUDENT START ###\n",
    "    ngram_labels = [\"Unigram\", \"Bigram\"]\n",
    "    for n in range(2):\n",
    "        # Use either unigram or bigram\n",
    "        vectorizer = CountVectorizer(ngram_range=(n+1, n+1))\n",
    "        doc_matrix_train = vectorizer.fit_transform(train_data)\n",
    "        doc_matrix_dev = vectorizer.transform(dev_data)\n",
    "        # Fit logistic regression model\n",
    "        lr_model = LogisticRegression()\n",
    "        lr_model.fit(doc_matrix_train, train_labels)\n",
    "        # Determine 5 features with biggest weights for each of the 4 classes\n",
    "        weights_sorted = np.argsort(lr_model.coef_, axis=1)\n",
    "        # Obtain index of maximum weights for each class, then flatten to create array of 20 features\n",
    "        weights_maxindex = weights_sorted[:, -5:].flatten()    \n",
    "        # Create table\n",
    "        full_vocab = np.array(vectorizer.get_feature_names())\n",
    "        print (ngram_labels[n] + \" results:\")\n",
    "        print (\"  Size of vocab: {}\".format(len(full_vocab)))\n",
    "        print (\"  Accuracy: {}\".format(lr_model.score(doc_matrix_dev, dev_labels)))\n",
    "        print (full_vocab[weights_maxindex])\n",
    "        print (np.transpose(np.round(lr_model.coef_[:, weights_maxindex],2)))\n",
    "        print (\"\\n\")\n",
    "    \n",
    "### STUDENT END ###\n",
    "P4()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ANSWER:  \n",
    "The unigram model has descriptive terms that one would expect in each category. The bigram model seems to have more generic terms, such as \"you are\" and \"is there\". The accuracy of the bigram model is much worse than the unigram model, perhaps due to the prevalence of these generic terms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(5) Try to improve the logistic regression classifier by passing a custom preprocessor to CountVectorizer. The preprocessing function runs on the raw text, before it is split into words by the tokenizer. Your preprocessor should try to normalize the input in various ways to improve generalization. For example, try lowercasing everything, replacing sequences of numbers with a single token, removing various other non-letter characters, and shortening long words. If you're not already familiar with regular expressions for manipulating strings, see https://docs.python.org/2/library/re.html, and re.sub() in particular. With your new preprocessor, how much did you reduce the size of the dictionary?\n",
    "\n",
    "For reference, I was able to improve dev F1 by 2 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "empty_preprocessor\n",
      "  The size of the vocabulary is 26879.\n",
      "  F1 score: 0.681318340101\n",
      "better_preprocessor\n",
      "  The size of the vocabulary is 22090.\n",
      "  F1 score: 0.70382337505\n",
      "The better preprocessor reduced the vocabulary size by 4789, or 17.82%.\n"
     ]
    }
   ],
   "source": [
    "def empty_preprocessor(s):\n",
    "    return s\n",
    "\n",
    "def better_preprocessor(s):\n",
    "### STUDENT START ###\n",
    "# String operations by preproessor\n",
    "    # Turn into all lowercases first\n",
    "    s = s.lower()\n",
    "    # Remove non-letter characters\n",
    "    s = re.sub(r\"[_%#:+]\", r\"\", s)\n",
    "    # Shortening long words: stem \"s\", \"ed\", and \"y\"\n",
    "    s = re.sub(r\"s\\b|ed\\b|y\\b\", r\"\", s)\n",
    "    # Replace numbers 4 digits or longer with single token 99\"\n",
    "    s = re.sub(r\"[0-9]{4,}\", \"99\", s)\n",
    "    return s\n",
    "### STUDENT END ###\n",
    "\n",
    "def P5():\n",
    "### STUDENT START ###\n",
    "    vocab_size = []\n",
    "    # Iterate through the an empty preprocessor and the above better preprocessor\n",
    "    for preprocessor in [empty_preprocessor, better_preprocessor]:    \n",
    "        # Process raw strings in both train data and the dev data\n",
    "        train_data_parsed = [preprocessor(post) for post in train_data]\n",
    "        dev_data_parsed = [preprocessor(post) for post in dev_data]\n",
    "        # Vectorize raw strings using default CountVectorizer options\n",
    "        vectorizer = CountVectorizer()\n",
    "        doc_matrix_train = vectorizer.fit_transform(train_data_parsed)\n",
    "        doc_matrix_dev = vectorizer.transform(dev_data_parsed)\n",
    "        vocab_size.append(len(vectorizer.get_feature_names()))\n",
    "        print (preprocessor.__name__)\n",
    "        print (\"  The size of the vocabulary is {}.\".format(len(vectorizer.get_feature_names())))\n",
    "        # Fit logistic regression model and obtain F1 score\n",
    "        lr_model = LogisticRegression()\n",
    "        lr_model.fit(doc_matrix_train, train_labels)\n",
    "        print (\"  F1 score: {}\".format(metrics.f1_score(dev_labels, lr_model.predict(doc_matrix_dev), average=\"macro\")))\n",
    "    print (\"The better preprocessor reduced the vocabulary size by {}, or {:.2f}%.\"\n",
    "           .format(vocab_size[0]-vocab_size[1], 100.*(vocab_size[0]-vocab_size[1])/vocab_size[0]))\n",
    "    \n",
    "### STUDENT END ###\n",
    "P5()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(6) The idea of regularization is to avoid learning very large weights (which are likely to fit the training data, but not generalize well) by adding a penalty to the total size of the learned weights. That is, logistic regression seeks the set of weights that minimizes errors in the training data AND has a small size. The default regularization, L2, computes this size as the sum of the squared weights (see P3, above). L1 regularization computes this size as the sum of the absolute values of the weights. The result is that whereas L2 regularization makes all the weights relatively small, L1 regularization drives lots of the weights to 0, effectively removing unimportant features.\n",
    "\n",
    "Train a logistic regression model using a \"l1\" penalty. Output the number of learned weights that are not equal to zero. How does this compare to the number of non-zero weights you get with \"l2\"? Now, reduce the size of the vocabulary by keeping only those features that have at least one non-zero weight and retrain a model using \"l2\".\n",
    "\n",
    "Make a plot showing accuracy of the re-trained model vs. the vocabulary size you get when pruning unused features by adjusting the C parameter.\n",
    "\n",
    "Note: The gradient descent code that trains the logistic regression model sometimes has trouble converging with extreme settings of the C parameter. Relax the convergence criteria by setting tol=.01 (the default is .0001)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of nonzero weights with L2 regularization is 107516.\n",
      "  Accuracy is 0.701183431953.\n",
      "The number of nonzero weights with L1 regularization is 1636.\n",
      "  Accuracy is 0.690828402367.\n",
      "The number of nonzero weights with L2 regularization and reduced vocab is 3908.\n",
      "  It is much smaller after features are pruned by the L1 model.\n",
      "  Accuracy is 0.681952662722.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAEZCAYAAAA68jcWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYHVWd//H3JwkQAmEJMKxhkUVBZBk0iCATgWEbNIIo\nmwgoy+iwjMrIMKMScX4sbiigDIPsICibBofdISwiS4BAIgkQIGQBwbCZQFBCvr8/zrl09c293X07\nXX1v3/68nqefrr3OqeXUt06dW6WIwMzMzKyvDWl2AszMzKw9OcgwMzOzUjjIMDMzs1I4yDAzM7NS\nOMgwMzOzUjjIMDMzs1I4yOghSTtKelrSfEmf6qd1nifpm72Yb/2cTpWRrlYl6SZJh7bq+iVdIum7\n/ZmmZlra/Ja1PyVNlbRzXy+3GSQtL+lGSa9L+mU/r7sp21HSxZJelXR/f6+7FkkTJX2pl/OOl3R5\nX6epF+ko7ZrR4yAjb8hXJS3b14kYIE4Fzo6IkRExoXqkpJmSdu3LFUbElyPiv7qbLq97l8J8s3I6\nG3oJiqTDJb2bD7Y3JD0uad/epL0ZImLviGjaCVtcf96W91RPkv+6JWlDSYslLXGOSjpM0qS8j2ZL\nOlPS0C6WtVjSgrxf50o6W9KwRvLWSz3Ob82Z+2B/1gp0ImLLiLh7aZbbQvYH/g4YFREHlLWSVtmO\nkj4O7AasExEfrRr30Xycr1BjvkclfaWkZC3Ncd5vL6qStJ6k6yT9OQelUyQdBr2/ZvREj4IMSRsC\nY4CXgX65iy+suz8Kw55YH3iii/FLVaAupQD6KgL9fUSMBFYBzgV+IWnVPlr2e2pdPK3HlgdOAFYD\ntgd2BU7sZp6t8n7dGdgPOLrUFHZo+LhUVkZi2tAGwFMRsbjZCeknGwAzI+Lt6hERcT8whxR4vUfS\nlsDmwFX9ksLG9Po470UZejnwPOlaNgo4FHipt+vvsYjo9g/4NjAB+E/gxqpxo4HrSQHIPOCcwrij\nSBfmvwB/BLbJwxcD7ytMdwnw3dw9lnSgfAN4EbiUdMH7bV7Hq8CNwLqF+UcBFwNz8/jr8/CpwD6F\n6ZbJady6Tj6PAp4GXgF+A6ydhz8DvAu8lfOyTI15nwN2qTF8OeDHOW1zgbOAZQvjvwG8kPN8ZHHb\nVG2X1fM2eC2n727SAXp5IW3zSRebDfNyhtTZPjfUyf/hwD2F/hF5OR8u5OUHpAP1T8B5wPAG8nIe\ncBOwANgFWAe4Lu/XZ4HjCssaA0wC3sjr+mEePhy4Iu/H14AHgTXyuInAl3K3gG8CM0kn0qXASnlc\nZft8Ieflz8B/1NkmGwGvFfovAF4q9F8OnFBcP/AB4G1gUd4nr+bxF5MCt9+SjqP7KZwHVevttA+7\nOT+/CkzoYnz1+fZL4NxC/z7A5Lw9fw98qDDu74FHc3p/leetHJOdjpfqdeX8VqZdla7P4YnAf+X1\nvwlsXLU/H8vbsvK3GNg5j7uGVFa8DtwFbJGHHw38Dfhrnuc3efhMYNfuzk86yqKv5WPoBeDwLrbz\nEXSUd88ARxfG1Tx/6yznJ8As0rE/CdipznTfyXn7W87fF4HxwOX1jqO8TU8F7s3pvBVYrTD9TsB9\nOZ2zgMNI5WJ/bsd1SNebV0jl8ZF5+JeAhXScV6fUmPdk4HdVw74HXJe7PwY8lI+VB4EdCtPVLCfp\n/ti9EzgNeCDvs18DqxbyPrsqPTPJ14oa+6vmsVynDD2RVDYOKUyzHzC5znadT7rZ6LK8AXag87n2\nNvBcnm4I8O/ADFIZ/MtKXuvuz+4KsLzgGcAhwKb5YPu7PHwo6eT/Ienuajlgxzzus/nA2i73bwys\nX6fQuxg4tbBT3gFOJwUFw/PO3zd3r0gq7G4ozP+/pCh1ZWAY8PE8/N+AqwvTjQMeq5PHXUgXm22A\nZYGzgbsK42sGEd2NJ53Q95EKmdVJhWglr3vmA2rzvP2uYMlCujLt6fkAG5r/dqy3bpYsWGpunxpp\nPZx80cjr+BfSSTUyDzuLdAKtkvfDBOC0HublEtKJs0PuXx54mBQIDCNdzJ8Bds/j/wAckrtHAGNy\n9zF5vcNJgcS2hfTdCXwxd3+RVEBtCKxACmYuq9o+55OO2a1IJ9IH6myX54Ftc/eTpPPhA4VxW9dY\n/2EseQG+hHRifjhv3yuAq7o76Xtwfv66sh/qjF8MbJy7P0Aq5L+Q+7clFfwfydvzC/l4WoZ0HjwP\nHJfTuy/pQlM5Jt87XqrWVSvI6O4cnkgqfDcnFWTDituzah1Hky7mKxbSsUJO81nAo7XKllrnC12f\nn2NJZdH4nP+9SAHQynW2897ARrl75zxt5caq7vlbYzmHkC5sQ0gX5hcp3JhUTXsK+bgu9HcXZDwN\nbJL3xZ3A6XncBqTA44CcxlF0HNv9uR3vJgXjywJbky7un6h3XlXNOzqva73cPwSYTaqBH0UKng7J\nww8klW+VgKDedaQnx+4cYAtSWXVtZR9QO8gobrfxVfvrcOofy5fQuQxdjnTzvmdhmhuAr9bZNreT\ngssDyNfi7sqbvB0mAv8v95+Q9/M6OY3/Dfyiy/KpBwXYTqTosVKQTwb+NXfvkA+AJQpCUoR8XJ1l\n1goyijUZf6XOSZWn2YaOu8O1SXfySxyweUPMp6MwuhY4sc4yLwTOKPSvQAqoKoFRb4OMGVUHwe50\nRIUXVXZe7t+Y+kHGd0gXk427Wzedo9K626fGcg4nnaCv5by/RUfQKFL0XNxvOwDP9jAvlwCXFMZv\nDzxftf6TgYty912kE3D1qmmOoOpuuzCueJH/HfDPhXGb5TwNKWyfdQrjHwAOqLNdLiPVFqwFTAfO\nIAU71bUcxfUfzpIX4IuB/yn07wVMq7PO9/ZhN/vsi6Q7zlFdTLOYdIe1IHefXRh3HktePKaTLpI7\nA3Oqxt1DL4KMrs7hwrYbX29/FobtRAqKNqmz3FVyGkbWSwOdC/muzs+xpHOgeKf4Ejng7cH5dANw\nfHfnbw+W82qt4z2PG0/ni1R1f6fjKG/T/yiM/zJwc+H8u67OevplO5KChEXACoVhpwEX1zvmaizj\nduDk3P2PpGvUUNLjgfurpr2PFLg0Uk7WOnZPK/RvTrqGiQaDjG6O5UsolKF52EnAFbl7FCl4W7OL\n5Z1OquFfRKqhrNRSdzpOCvOcR6GWlBTcF681a5PL1XrbqyfPdA4DbouI+bn/mjysckA8H7WfB65H\nujPtjT9HxN8qPZJGSDo/N3B8g3QBWjk/tx1N2uFvVC8kIl4gXZD2l7QK6W77yjrrXJt011aZ901S\ndd26vcxDxTrF5ZIuCOsU1jm7MG5Ojfkrz+y+TzqRb5P0jKSTerj+utunjvsjYlXSndQE0kEMsAYp\nSn9Y0muSXgNuJt219CQvUTVsA2CdyrLy8k4mNWKDVDW6GTBN0oOS/ikPv5wUwF6dGzGeWafdTqf9\nSdruw4A1C8P+VOh+ixRY1nIXqbD4OOku6y7gH0gX4erGnd0pPgNdSLoz6hVJnyYVwHtFxKvdTL5t\nRKxIuov5gqQN8vANgK9X7Yf1SNtvHVLVcdFseqGbc7hHy5Y0mlQ9+4WImJGHDZF0hqQZebnP5clX\nr7ecKl2dnwCvVJVvb1Fnn0naS9L9kl7J23FvUrsZaOD8lXSipCdy47zXSHfWPc1PTxSP++IxOJr0\n2LI3+mo7rkMqr96sWlYj5fClpICC/P+qiHg3L3tW1bTP5+HrUaec7MWxO4t0l9/QPpM0tJtjOVjy\nHLkS+KSkEcDngLsjomY7i4h4PSJOjogtSeXgZFLgWy89x5DKuIMLgzcEbiiUFU+QApY1l1xC0mWQ\nIWn5nPBdJL0o6UXg68DWkrbKGV6/Tsv22aQquVreIl2wKtamc6PJ6Dw5XyddcMZExMqkAl75bzYw\nStLKddZ1KfB50uOb+yLixTrTvUDagADkFsqrsWQh26hOyyU1uqks80XSiV1R7O4kIhZExIkRsTGp\n6u9rkj5RGd3F+rvbPvXW9ybpLucflH6mNo9UIG0REavmv1UiYqUG8lJM5yzSnc6qhb+VImKfvP4Z\nEXFwRKwBnAlcK2n5iFgUEadGxAdJz1f3IVXxV6u13RfRu4ZOd5ECjLGkqsN7gR1Jx+HEOvN0tU+W\nmqQ9gf8htTn6Y0/ni4hrSM+Xx+dBs0g1UMX9sGJE/JK0T6sL9/UL3W9SOI8lrVVrlfl/V+dw9bRL\nyGXRr4GzIuLWwqhDSOfDrnm5G1Vm6W6ZWa3j5IVu5qmVvuVIj+S+R3qcvCrp2bmg7vm7S43lfJz0\nmPez+fxalVQLVa+BYHX+FtC5bK21T+qZRaqB7Ml6qvXJdszzjJJUDEDWp/YNWD03AOvl8nFf0jUA\nUrm7QdW0G+ThXZWTPTl216/qfodUZlafI0NJN2y1HEzXx/ISImIOqW3XfqTrXI9+jRURr5CaOaxT\nq2F/Pg5PBcZFxILCqFmkGqtieTGii+tqtzUZnyYVzJuTno1tnbvvIRXsD5AKojNytDdc0sfyvD8H\nTpT097mx+CaSKjtiMnBIjtz2JEVLXVmRdIF7Q9Io0nNHAHLmbgZ+JmkVScuo82+3byA1XjueVO1d\nz1XAEZK2zgXGaaS7+urItyvL5m1Q+RuWl/tNSatLWp3UiPaKPP2v8jo/kCPRb1Ut772DS9I+eRuK\n9Nz0XVL1FqQLZ83CoQfbp66IeI10ITs534VcAPxY0ho5TetK2r3RvGQPAvMlfUPpt/5DJW0p6cN5\n2Z+vrIdUyAawWNInJH0on6zzSSfzuzWSfxXwVaWfgq5I2p9X16l1q5fGynaYQWqz8XlSO535pCrY\nz5ACkFpeIhV0y3S3/G4Uj6fh+VzahXQHs19ETOrFMs8ADpK0Hmmf/rOkMXnZK0j6p7zN7gPelXSs\npGGSxpHablQ8BnwwnzPD6QhcKooFcd1zuGr6ei4iPVr6QdXwFUlV06/mG4PTqsa/BLyvi+XWOj97\n87PZZfPfPNJxuhfpkQFQ9/ytddyOJJW58yQtK+nbwEo1pntv0VX9k4GdJY3OF8yTezBPxS+A3SR9\nNu/v1SRtncf1y3aMiNmk4+50Scvlm9kv0lFm9mQZb5IejV9M+iXKI3nUTcBmkg7K+TuA1EbptxHx\nJ5YsJz+e5+vu2BXweUmb57LvVOCaSM8TniKdw3vnsuCbpLYUtXR3LNfbb5eRapy3JP0IoyalWt8P\n5ryPJN1EPp3L+eJ0o0nl+aGVGsOC/wZOq1zLJa2hbt4b1V2Q8QXSM/I5EfFy/nuJ1CinUoXySVKN\nxSxSNPg5gIi4Fvh/pAP3LznzlYjphDzfa3k5N1Sttzpq/jGpoeA80gF4c9U0h5IuNtNJJ8Px7y0o\n/dTpelKUXXcHRMTvSBfG60jR9EakhkGNuIlUS1P5+zapxfwk4PH8NykPIyJuITUwvZN0MP4hL+ev\nlWQV8rkJ6VnjfNI2+GlEVC5wp5NO8Nckfa0wb0Xd7VOluL6KHwOfyCf7SaQq3/uVqvNuJ0X4jeaF\nfLHfh/R881lSo9v/oaNA3QOYKmk+qQHUgRHxV1K13DWkwOMJUk1CrcLsojz87rz8t0gNGIt5rZX/\neiYC8yJibqEf4JGaU6c2IX8E/iTp5cLyq9fR3R3iAjqOpzdJDZS/SboY3az07ov5kv63i2V0WkdE\nTAX+D/haRDxM+vXAuaRn/0+Ta4Yi4h3SHdKX6Ggw91vSM1gi4ilSgXoHqUHsPVXrKua3u3O4u21x\nAPDpQn7nS9qRVMA+T7obnUo67orLuRDYIp8btc7/uudnD9LUMVEKPI8nFc6vAgeRfqFW0dX5W3RL\n/nuK1BB2IUtW8XdaNZ3PqztIj5QeJ/2K4sYaeai5j/IN1d6kO/dXSM/st8rT9ct2zA4ildcvkMrs\nb0fE/1WntxuXkmoU3ruxzI8U9yHlbx7p1xn7FB41VpeTJ+Th3R27kddzCbmRLrmMzY9fvkK66Z5D\nOp9nV81bWVZ3x3K9vF+f83pD1Phpb8HypGvta6SmDKPp/EqKyrJ3JT22vq5wrk3J435Ceox+m6S/\n5DSO6WKd6SdU7U7St4BNI6JWtXpvlrcn6cAbCvw8Is6sGn8iqUCG1A5gc1IDxtclzaTjTuadiBhT\nmG9zYAqp0euA/t17O+XFOkh6APhZRFza7cRm1i8kPQ0cUwjGWkbbBxm5euthUtXPvX2wvKGku7bd\nSBHnQ8BBETGtzvT7kH6Ns1vuf470s95Xc/++pBqQEaToe1FE7Le06WyGdsqLJfnR2lOku7hDgJ+R\nfj1S/kt8zKxbkvYj/TJys2anpZa2fuuipKNIVY0390WAkY0BZkTEzFydfDXp/Rv1HMySb5orPls7\nmlQ1N4NUVfflPkpnM7RTXix5Px0v6voqsL8DDLPWIGkiKfD/lyYnpa62r8noa5L2B/aIiKNy/+eB\n7SPiuBrTjiA9f9s4Il7Pw54ltSd4Fzg/Ii7ot8SbmZn1o1b5LshA0khU9kng3kqAke0YES/mX07c\nLml6RDT6rgUzM7OW5yCjcXNZ8n0Q9X7DfSBVj0oqvyeOiD9LuoH0+KVTkCHJ1UtmZr0QEf64Xgtx\nkNG4ScCmSl+mfYH007qDqifKv1Hv9La0/PhkaETMz7+D3p30uuEltPOJIml8RIxvdjrK4vwNbO2c\nv3bOG/gGrRU5yGhQRCySdCzp1dZDgQsjYprSK1iJiPPzpJ8Gbo2IhYXZ1yS9khXStr8yIm7rv9Sb\nmZn1HwcZvRARN5NeyFIcdn5V/6V0vM62Muw50sunzMzM2l5b/4TVWtbEZiegZBObnYCSTWx2Ako2\nsdkJKNHEZifABhf/hLUFSYp2bpNhZlYGl52txzUZZmZmVgoHGWZmZlYKBxlmZmZWCgcZZmZmVgoH\nGWZmZlYKBxlmZmZWCgcZZmZmVgoHGWZmZlYKBxlmZmZWCgcZZmZmVgoHGWZmZlYKBxlmZmZWCgcZ\nZmZmVgoHGWZmZlYKBxlmZmZWCgcZZmZmVgoHGWZmZlYKBxlmZmZWCgcZZmZmVgoHGWZmZlaKYc1O\ngJmZWXekLfeG0cfDyOEw/22YfXbE1JuanS7rmoMMMzNraSnA2OEncMEmHUOP2ljaEgcarc2PS8zM\nrMWNPr5zgAGpf/3jmpMe6ynXZJiZlcjV/H1h5PDaw1dcvn/TYY1ykGFmVhJX8/eV+W/XHr5gYf+m\nwxrlxyVmZqVxNX/fmH02HDWj87Ajn4FZ5zQnPdZTrsnoBUl7Aj8GhgI/j4gzq8afCBySe4cBmwOr\nR8Tr3c1rZu3E1fx9IWLqTdKWwN7HpW23YCHMOse1Qa3PQUaDJA0FzgV2A+YCD0maEBHTKtNExA+A\nH+Tp9wH+NQcY3c5rZu3E1fx9JQcUDioGGD8uadwYYEZEzIyId4CrgXFdTH8wcFUv5zWzAc3V/Da4\nuSajcesCswv9c4Dta00oaQSwB/CVRuc1s4HP1fw22DnIaFw0MO0ngXsj4vVezGtmbcDV/DaYOcho\n3FxgdKF/NKlGopYD6XhU0tC8ksYXeidGxMRGE2pm1s4kjQXGNjkZ1gVF+Oa6EZKGAU8CuwIvAA8C\nB1U33pS0MvAssF5ELGxw3ogIlZ0XM7N24rKz9bgmo0ERsUjSscCtpJ+hXhgR0yQdk8efnyf9NHBr\nJcDoat7+zYGZmVn/cE1GC3I0bmbWOJedrcc/YTUzM7NSOMgwMzOzUrhNxiDir0GamVl/cpAxSPhr\nkDZQOTg2G7gcZAwa9b4Gufdx+EVBpfEFcuk4ODYb2BxkDBr+GmR/8wWyLzg4NhvI3PBz0PDXIPtf\nvQvk+sc1Jz0DkYNjs4HMQcag4a9B9j9fIJeeg2OzgcyPSwaJ/voapNsgFPkCufRmnw1Hbdy5RsjB\nsdlA4SBjECn7a5Bug1DNF8il5U+lmw1sfq14Cxqor8aV9roFbt5jyTF73xJx0179n6LmS4HX+r5A\nmvWDgVp2tjPXZFgfchuEamXXHpmZtTI3/LQ+5DYIZmbWwUGG9SH/gsXMzDq4TUYLGsjPFVutDYJ/\n7WI2eAzksrNdOchoQT5R+kadX7vMgD+c4EDDrP247Gw9flxibcxv3DQzayYHGdbG/GsXM7NmcpBh\nbcy/djEzayYHGdbG/GsXM7NmcsPPFuTGS32n1X7tYmblcdnZehxktCCfKGZmjXPZ2Xr8uMTMzMxK\n4SDDzMzMSuEgw8zMzErhIMPMzMxK4SDDzMzMSjGs2QmwpeePgJmZWStykDHA1fkI2MbSljjQMDOz\nZvLjkl6QtKek6ZKelnRSnWnGSnpU0lRJEwvDZ0p6PI97cOlT44+AmZlZa3JNRoMkDQXOBXYD5gIP\nSZoQEdMK06wC/BTYIyLmSFq9sIgAxkbEq32TIn8EzMzMWpNrMho3BpgRETMj4h3gamBc1TQHA9dF\nxByAiJhXNb4P30jnj4CZmVlrcpDRuHWB2YX+OXlY0abAKEl3Spok6dDCuADuyMOPWvrk+CNgZmbW\nmvy4pHE9+djLMsDfA7sCI4A/SLo/Ip4GdoqIFyStAdwuaXpE3FO9AEnjC70TI2JizcTE1JukLYG9\n/REwMxtUJI0FxjY5GdYFBxmNmwuMLvSPJtVmFM0G5kXEQmChpLuBrYGnI+IFgIj4s6QbSI9flggy\nImJ8TxOUAwoHFWY2qOSbr4mVfkmnNC0xVpMflzRuErCppA0lLQscAEyomuY3wE6ShkoaAWwPPCFp\nhKSRAJJWAHYHpvRj2s3MzPqNazIaFBGLJB0L3AoMBS6MiGmSjsnjz4+I6ZJuAR4HFgMXRMQTkt4H\nXC8J0ra/MiJua05OzMzMyqWInjQxsP4kKSKiD3+BYmbW/lx2th4/LjEzM7NSOMgwMzOzUjjIMDMz\ns1I4yDAzM7NSOMgwMzOzUjjIMDMzs1I4yDAzM7NSOMgwMzOzUjjIMDMzs1I4yDAzM7NSOMgwMzOz\nUjjIMDMzs1IM2iBD0qckDdr8m5mZlW0wX2QPAGZI+p6kDzQ7MWZmZu1mUH/qXdLKwEHA4UAAFwNX\nRcT8JqfLnys2M2uQy87WM5hrMoiIN4BrgV8C6wD7Ao9KOr6pCTMzM2sDgzbIkDRO0g3ARGAZ4CMR\nsRewFfC1ZqbNzMysHQxrdgKaaD/grIi4uzgwIt6SdGST0mRmZtY2Bm2bDEnvA16MiIW5f3lgzYiY\n2dSE4eeKZma94bKz9QzaxyXAr4B3C/2LSe0zzMzMrA8M5iBjWET8rdITEX8ltc0wMzOzPjCYg4x5\nksZVenL3vCamx8zMrK0M5jYZmwBXkn66CjAHODQiZjQvVYmfK5qZNc5lZ+sZtEFGhaSRQETEgman\npcIniplZ41x2tp7B/BNWJO0DbAEMl9JxGRGnNjVRZmZmbWLQtsmQdD7wOeB4QLl7g6YmyszMrI0M\n2sclkqZExIckPR4RW0laEbglInZqgbS5ys/MrEEuO1vPoK3JABbm/29JWhdYBKzVxPSYmZm1lcEc\nZNwoaVXg+8DDwEzgqp7MKGlPSdMlPS3ppDrTjJX0qKSpkiY2Mq+ZmVk7GJSPSyQNAXaIiN/n/uHA\n8Ih4vQfzDgWeBHYD5gIPAQdFxLTCNKsAvwf2iIg5klaPiHk9mTfP7yo/M7MGuexsPYOyJiMiFgM/\nLfS/3ZMAIxsDzIiImRHxDnA1MK5qmoOB6yJiTl7+vAbmNTMzawuDMsjI7pC0vyq/Xe25dYHZhf45\neVjRpsAoSXdKmiTp0AbmNTMzawuD+T0Z/wx8DXhX0tt5WETESt3M15PnS8sAfw/sCowA/iDp/h7O\nC4Ck8YXeiRExsafzmpkNBpLGAmObnAzrwqANMiJixV7OOhcYXegfTaqRKJoNzMufkV8o6W5g6zxd\nd/NW0je+l+kzMxsU8s3XxEq/pFOalhiradAGGZJ2rjU8Iu7uZtZJwKaSNgReAA4ADqqa5jfAubmh\n53LA9sCPgKd6MK+ZmVlbGLRBBvANOh5fDCc1ynwY2KWrmSJikaRjgVuBocCFETFN0jF5/PkRMV3S\nLcDjwGLggoh4AqDWvH2fNTMzs+YblD9hrUXSaOAnEbFfC6TFP8MyM2uQy87WM5h/XVJtDrB5sxNh\nZmbWLgbt4xJJ5xR6hwDbkB6XtBRpy71h9PEwcjjMfxtmnx0x9aZmp8vMzKw7gzbIIAUUlWdFi4Bf\nVN4A2ipSgLHDT+CCTTqGHrWxtCUONMzMrNUN2jYZ+aurCyPi3dw/FFguIt5qbso6nitKe90CN++x\n5BR73xJx0179nzIzs9blNhmtZzC3ybgDWL7QPyIPayEjh9cevuLytYebmZm1jsEcZAyPiAWVnoiY\nTwo0Wsj8t2sPX7Cw9nAzM7PWMZiDjDclbVfpkfRhoMUu3rPPhqNmdB525DMw65za05uZmbWOwdwm\n4yOkr6C+mAetDRwQEZOal6qk+FwxNf5c/7j0iGTBQph1jht9mpktyW0yWs+gDTIAJC0LvD/3PhkR\nf2tmeip8opiZNc5lZ+sZtI9L8uu9V4iIKRExBVhB0leanS4zM7N2MWhrMiQ9FhFbVw2bHBHbNCtN\nhXQ4Gjcza5DLztYzaGsygCGS3st/fk/GMk1Mj5mZWVsZzG/8vBW4WtL5gIBjgFuamyQzM7P2MZgf\nlwwFjgZ2Jb1e/HFg7YhoersMV/mZmTXOZWfrGbSPS/LrxB8AZgJjSMHGtGamyczMrJ0Musclkt4P\nHAQcAPwZuIZUozO2mekyMzNrN4PucYmkxcBvgWMjYlYe9lxEbNTclHVwlZ+ZWeNcdraewfi4ZD/S\n68PvlvTfknYlNfw0MzOzPjToajIq8qfex5EenXwCuAy4ISJua2rCcDRuZtYbLjtbz6ANMookjQL2\nBw6MiF1aID0+UczMGuSys/U4yGhBPlHMzBrnsrP1DMY2GWZmZtYPHGSYmZlZKRxkmJmZWSkcZJiZ\nmVkpHGSYmZlZKRxkmJmZWSkcZJiZmVkpHGT0gqQ9JU2X9LSkk2qMHyvpDUmP5r9vFcbNlPR4Hv5g\n/6bczMys/wy6r7AuLUlDgXOB3YC5wEOSJkRE9Wfi74qIT9VYRABjI+LVkpNqZmbWVK7JaNwYYEZE\nzIyId4Dku+CvAAAUwElEQVSrSd9AqdbVW+f8RjozM2t7DjIaty4wu9A/Jw8rCuBjkh6TdJOkLarG\n3SFpkqSjSk6rmZlZ0/hxSeN68rGXR4DREfGWpL2AXwOb5XE7RsSLktYAbpc0PSLuqV6ApPGF3okR\nMXEp021m1lYkjQXGNjkZ1gV/IK1Bkj4KjI+IPXP/ycDiiDizi3meA7arboch6RRgQUT8sGq4P/Jj\nZtYgl52tx49LGjcJ2FTShpKWBQ4AJhQnkLSmJOXuMaRg7lVJIySNzMNXAHYHpvRv8s3MzPqHH5c0\nKCIWSToWuBUYClwYEdMkHZPHnw/sD3xZ0iLgLeDAPPtawPU5/hgGXBkRt/V3HszMzPqDH5e0IFf5\nmZk1zmVn6/HjEjMzMyuFgwwzMzMrhYMMMzMzK4WDDDMzMyuFgwwzMzMrhYMMMzMzK4WDDDMzMyuF\ngwwzMzMrhYMMMzMzK4WDDDMzMyuFgwwzMzMrhYMMMzMzK4WDDDMzMyuFgwwzMzMrhYMMMzMzK4WD\nDDMzMyuFgwwzMzMrhYMMMzMzK4WDDDMzMyuFgwwzMzMrhYMMMzMzK4WDDDMzMyuFgwwzMzMrhYMM\nMzMzK4WDDDMzMyuFgwwzMzMrhYMMMzMzK4WDDDMzMyuFg4xekLSnpOmSnpZ0Uo3xYyW9IenR/PfN\nns5rZmbWLoY1OwEDjaShwLnAbsBc4CFJEyJiWtWkd0XEp3o5r5mZ2YDnmozGjQFmRMTMiHgHuBoY\nV2M6LcW8ZmZmA56DjMatC8wu9M/Jw4oC+JikxyTdJGmLBuY1MzNrC35c0rjowTSPAKMj4i1JewG/\nBjYrN1lmZmatxUFG4+YCowv9o0k1Eu+JiPmF7psl/UzSqDxdl/NWSBpf6J0YEROXLtlmZu1F0lhg\nbJOTYV1QRE9uzK1C0jDgSWBX4AXgQeCgYuNNSWsCL0dESBoD/CoiNuzJvHn+iIhabTrMzKwOl52t\nxzUZDYqIRZKOBW4FhgIXRsQ0Scfk8ecD+wNflrQIeAs4sKt5m5EPMzOzsrkmowU5Gjcza5zLztbj\nX5eYmZlZKRxkmJmZWSkcZJiZmVkp3PCzRUl73QIjh8P8t2H22RFTb2p2mszMzBrhhp8tSFJ0fufX\nUTPgDyc40DAzq88NP1uPH5cMCBdsAusf1+xUmJmZNcJBxoCx4vLNToGZmVkjHGQMGAsWNjsFZmZm\njXCQMSAc+QzMOqfZqTAzM2uEf13Ssva+JT0iWbAQZp3jRp9mZjbQ+NclLcgtpM3MGueys/X4cYmZ\nmZmVwkGGmZmZlcJBhpmZmZXCQYaZmZmVwkGGmZmZlcJBhpmZmZXCQYaZmZmVwkGGmZmZlcJBhpmZ\nmZXCQYaZmZmVwkGGmZmZlcJBhpmZmZXCQYaZmZmVwkGGmZmZlcJBhpmZmZXCQYaZmZmVwkGGmZmZ\nlcJBRi9I2lPSdElPSzqpi+k+ImmRpM8Uhs2U9LikRyU92D8pNjMz638OMhokaShwLrAnsAVwkKTN\n60x3JnBL1agAxkbEthExpuz0tiJJY5udhjI5fwNbO+evnfNmrclBRuPGADMiYmZEvANcDYyrMd1x\nwLXAn2uMU4npGwjGNjsBJRvb7ASUbGyzE1Cysc1OQInGNjsBNrg4yGjcusDsQv+cPOw9ktYlBR7n\n5UFRGB3AHZImSTqqzISamZk107BmJ2AAiu4n4cfAv0dESBKday52jIgXJa0B3C5pekTcU0pKzczM\nmkgRPblmWoWkjwLjI2LP3H8ysDgizixM8ywdgcXqwFvAURExoWpZpwALIuKHVcO9U8zMeiEiBvvj\n6JbimozGTQI2lbQh8AJwAHBQcYKIeF+lW9LFwI0RMUHSCGBoRMyXtAKwO/Cd6hX4JDEzs3bgIKNB\nEbFI0rHArcBQ4MKImCbpmDz+/C5mXwu4Pj1BYRhwZUTcVnaazczMmsGPS8zMzKwU/nVJi+npi75a\nWa0XjkkaJel2SU9Juk3SKoXpT875nS5p9+alvDZJF0l6SdKUwrCG8yNpO0lT8rif9Hc+6qmTv/GS\n5uR9+KikvQrjBlr+Rku6U9IfJU2VdHwePuD3YRd5a4v9J2m4pAckTZb0hKTT8/ABv+8GjYjwX4v8\nkR6/zAA2BJYBJgObNztdvcjHc8CoqmHfA76Ru08CzsjdW+R8LpPzPQMY0uw8VKX948C2wJRe5qdS\nY/ggMCZ33wTs2ey8dZG/U4Cv1Zh2IOZvLWCb3L0i8CSweTvswy7y1k77b0T+Pwy4H9ipHfbdYPlz\nTUZr6emLvgaC6sarnwIuzd2XAp/O3eOAqyLinYiYSSoUWupNqJF+Yvxa1eBG8rO9pLWBkRFReZX8\nZYV5mqpO/qD2S+MGYv7+FBGTc/cCYBrp3TYDfh92kTdon/33Vu5clnQj9hptsO8GCwcZraXbF30N\nELVeOLZmRLyUu18C1szd65DyWTFQ8txofqqHz6X183mcpMckXViojh7Q+cu/CtsWeIA224eFvN2f\nB7XF/pM0RNJk0j66MyL+SJvtu3bmIKO1tEsr3B0jYltgL+BfJH28ODJSfWVXeR1Q26EH+RmIzgM2\nArYBXgR+2PXkrU/SisB1wAkRMb84bqDvw5y3a0l5W0Ab7b+IWBwR2wDrATtL+kTV+AG979qdg4zW\nMhcYXegfTefoe0CIiBfz/z8DN5Aef7wkaS2AXHX5cp68Os/r5WGtrpH8zMnD16sa3rL5jIiXIwN+\nTscjrAGZP0nLkAKMyyPi13lwW+zDQt6uqOSt3fYfQES8AfwvsB1tsu8GAwcZreW9F31JWpb0oq8J\n3czTUiSNkDQyd1deODaFlI/D8mSHAZWCfgJwoKRlJW0EbEpqoNXqGspPRPwJ+Iuk7SUJOLQwT8vJ\nBXfFvqR9CAMwfzk9FwJPRMSPC6MG/D6sl7d22X+SVq886pG0PPCPwKO0wb4bNJrd8tR/nf9Ijxie\nJDVYOrnZ6elF+jcite6eDEyt5AEYBdwBPAXcBqxSmOc/cn6nA3s0Ow818nQV6e2ufyO1mTmiN/kh\n3YFNyePObna+usjfF0kN4x4HHiMVxmsO4PztBCzOx+Sj+W/PdtiHdfK2V7vsP+BDwCM5f48D/5aH\nD/h9N1j+/DIuMzMzK4Ufl5iZmVkpHGSYmZlZKRxkmJmZWSkcZJiZmVkpHGSYmZlZKRxkmJmZWSkc\nZAxAkv5PVZ9El/Svkn7WR8s/XNI5Dc5ziaTP9MX6q5b7fkkT8+eqn5B0fh6+XV98rlnScpLuyi/o\n6ROSti5+WruB+cZL+nrubmh75he4TakxfKykGxtNS2H+FSWdL2lG/hbNnZJa6gN2XZE0TtLmS7mM\nho81STMljaoa9gFJf5D0dmU/96d6x0gP5vt9L9f3H32xnDzvr/LLtWyAcZAxMF0FHFg17ADgF320\n/N68PKWh7wdIGtrDSc8GfhgR20bEFsA5ABHxcESc0Hgyl3AI8Nvo2xfGbAvsXWuEpGFdzFfchq3y\nPYafA/MiYpOI+DDpRWSr9+UKGjgWemNf0ue/e6w6Pb081mrtu1eA44AfNJCWppXRlWM1Inbs5SJO\nLvYsxXIALgC+uhTzW5M4yBiYrgP+qVII5K8vrhMR90o6SNLjkqZIOqMyg6Q9JT0sabKk2/OwMZLu\nk/SIpN9L2qywjtH5rvUpSd+urKd4JyTpREmnVCdO0rclPZjTcH5h+ERJZ0l6CPhPSc8W8rBS7q++\n4KxF4RsDETE1T//eHbqkm3JNx6OSXpd0qNKXG7+f0/GYpKPrbMuDgN8U0nhS3n6TJZ2eh20j6f68\nnOsLrzmeKOkMSQ9IelLSTkrfkTgVOCCn53O5huJySfcCl0raINdGPSbpDknFby2oujvfSU/MNQm3\nqOObDdvlZUwGvlInfwGsJOm3kqZLOk/JFyWdVcj3UZJ+VJxR0sakb158s7D9Z0bETdUrkbRA0o8k\nTc15Wr2w3Afz9rxW6dXQlZqa/5Z0P3CmpI/UOhaVatV+Lek2Sc9JOjYfd48o1QqsWkmrpJvzNrpb\nqQbsY8Ange/nfbFRrelqpacqb8Vjbbyki/K58Yyk4+ps9yV3RMSfI2IS8E5X0+Vt+YO8X3eQ9Pl8\njD2a0zgkT/elfNw9IOkC5dpHVdWCSVpQYx0b5vw/nP92KOT1Hkm/Ib2x9735JZ2qjvNsrqQL8/Bf\n5+05Vfmry0plz/J52surliOlc3OK0rn2ucK6J0q6RtI0SVcUkjyROoG7tbhmv3LUf737A24EPpW7\n/x34Hulzxs8DqwFDgd8B44A1gFnABnn6VfL/kcDQ3L0bcG3uPpz0mulVgeGkV/FuB2wITCmk4evA\nt3P3xcBncveqhWkuA/bJ3XcC5xbGXQSMy91HA9+vkc/DgdeBm4B/BVbOw8cCN1ZNux3p9cMj8/L+\nMw9fDngI2LBq+qHAi4X+vYDfA8OrttPjwMdz93eAswr5+X5h3ttz92EUXlsMjM/rX66w7w7N3UcA\nN+TuU4CvF7bnfsAywH3Aann4AcCFhXTtlLu/V9w3hXWPBRbmfTeE9ArmzwArkF6vXNn/vwc+WDXv\np4Dre3g8LgYOyt3fAs7J3aMK03wXODZ3X0L6zkTlrcNdHYtP5/SuDrwBHJ3H/Yj01VFIx/omuXt7\n4HfF7VhIQ73pOqWnxja8sbAv7837ZTVgXiXdVfM8V8x71bj39nMX23L/3L15Tldl2/yM9N2NdfI6\nVgGGAXeTjzkK52Lun5//b1g5RoDl6TgeNwUeKuR1AbmsKM5f6F+ZdOxtWzzf8zKnFPqr56uk4zOk\n41DA35HKrLXyul/PeRPpuN+xMP9dwOY9OR791zp/XVXdWmurPDKZQLrwfBH4CDAxIl4BkHQlsDPw\nLnB3RDwPEBGv52WsAlwmaRPSHW/xeLgtIl7Ly7me9I2EWh8UKt55V6qId5H0b8AI0jcGpgK/zeN+\nWZj+58A3SDUJhwNHVi88Ii6RdCvpWxPjgGMkbb1EItKd82XAZyNivlKblQ9J2j9PshKwCTCzMNvq\nQPGT37sCF0XE23ndr0tamRTY3JOnuRS4pjDP9fn/I6RCHNI2qd4uEyLir7n/o8Cnc/cVpAChFgHv\nBz4I3KHUbGQo8EIhXffmaS8nBTq1PBgRMwEkXUUKTK6T9H/AJyVNB5aJiD9WzdfI45rFdOzbK+jY\nLh+S9F+kC9OKwC2FZV8T+epB18finRHxJvCmpNdJQRqkC9pWSh/i+xhwjTqa1ixbmL9SI7QisEOd\n6arTU08A/xsR7wCvSHoZWJMUlPeVd0m1lZCOye2ASTnNw4E/kc71uyrnsqRrgM2WXFRdywLn5nPp\nXVKgUfFgpayoppSIK0mPMB/Ng0+QVDmeR9P9Rw53An6Rt/XLku7K+flLXvcLeV2TSedUpS3HC7l/\nWgP5tCZzkDFwTQDOkrQtMCIiHlXnanfofKGr5bukO7l9JW1AqpKsRaSLyCI6P2JbnqoLkaThwE+B\n7SJirtLjlOGFSd6sdETEfbnadizpTu2JWiuP9On4i4GLlR7XfLBqnUNJQdd3qpZxbETcXidPxbx1\n1d/d9JXA4V26Pp/e6uF6al3k/hgRH+s0c35k04PlVS9Thf6fA/9JKrQvqjHfE8DWkoZExOIull+t\nuI5LSDVuUyQdRrpbrShuk66Oxb8WuhcX+heTtvkQ4LWI2LZOeippGQK83sV01fuonr8Vurvb773x\ndlWwc2lEVDeiHFc1T3H/v3ee5kcry7Kkr5Jq8Q7N58/bhXFv1pi+YjwwKyIuzcsfSwqEPhoRb0u6\nk87ney3BksdrJb/FfV29bSvlkA0gbpMxQEXEAlJ1/cV0NPh8CPgHSavlguNAUmF9P7CzUtsNKs+x\nSXf3lTuwI6pW8Y+SVlV6hj6OdDfxMvB3kkZJWg7Yp0bSKgXMK/nO8bPdZOUy0p1RrYsckvZQaueA\nUluE1Si00cjOAB6PiF8Vht0KfEUdbT42kzSiar55pLvrituBI9TRbmDViHgDeE3STnmaQ6kfjFX8\nhVT9X899dDTcPYRU1Q21a0CeBNaQ9NGcpmUkbZHvYF+XtGNhOfWMycHcEOBzwD0AEfEgsB5wMClI\n6yQingEmkR4Rkde/oaRaz8aH0LGvD66sg7R9/5T34eepXzvS1bFYj3I65wPPVWqt8jP/rfI08/Oy\niYi/dDFdTzXyK6R60zayjN8B+0taAyCfe+vTca6vko/xz9CxbWeSaj8gPfJapsZyVyLViAB8gVRD\n1iVJnyQFFMVGsCuRAry3JX2AVEtX8Y5qN3S+h9RmaUjO186kmo/utsvapEcrNoA4yBjYriJ9Cvkq\neO+O/99JwcdkYFJE3BgR80htFK7PVZBX5/m/B5wu6RFSIVP8ZcODpCrbx0jPxx/JVcSn5nG3ke50\nO8kXvwtIj0huAR7oJg+/ILX9WOIil+0OTMnpvgU4MSJepvOvL75OCooqjdL2Id2lPwE8kms/zqPq\njjMi3gWmKjf+i4hbSTVEkyQ9mpcLqY3F9yU9BmyVt0EtlfTcCWyR0/K5qnGQfmFwRF7eIXQU2kv8\noiRv8/1JjSMrn/PeIY8+AvhpTmv1Ooppegg4N2+PZ+n82OtXwL05mKrlSGBNpZ+wTiEFtS/VmO5N\nUjAzhVRbUdlG3yIdA/eyZDV3Mb1dHYtRZ57iuEOAL+VtNJV0cYV0rP+bUuPGjbqYrnrZ1emsl56u\nPC5pdv77gaQ1Jc0h1SJ8U9KsHIjXWl/qiJhGanh7Wz5ebgPWyo8UTiOdi/eS2mf8Jc92ASkAmUy6\n6C+oseyfAYflad5fZ5rq/q+S2ks8mI/t8aRzcpikJ4DTgT8U5vufvA0uLy4nIm6g4zP0vyN9vr36\nnO607hykrhcR05fYWtbS/Kl3a6p8V/nJiDisSes/HFgzIs7sbtp2pPSriR9FxJ1LuZz5EdFV7Y31\nMUkrRMSbubbgelKD4N90N99ApNTG6p+ib362bv3INRnWNEo/uTuN9Dy+WX5B+jlwn72MayDI1exP\nAm8tbYCR+W6l/43PtVhTgGfbNcDIjgTO6nYqazmuyTAzM7NSuCbDzMzMSuEgw8zMzErhIMPMzMxK\n4SDDzMzMSuEgw8zMzErhIMPMzMxK8f8BcNRuC655zgYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1094423d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def P6():\n",
    "# Keep this random seed here to make comparison easier.\n",
    "    np.random.seed(0)\n",
    "### STUDENT START ###    \n",
    "# Use default CountVectorizer options\n",
    "    vectorizer = CountVectorizer()\n",
    "    doc_matrix_train = vectorizer.fit_transform(train_data)\n",
    "    doc_matrix_dev = vectorizer.transform(dev_data)\n",
    "    \n",
    "# Logistic regression using L2 penalty with the default value of C (1.0)\n",
    "    model = LogisticRegression(penalty='l2', C=1.0)\n",
    "    model.fit(doc_matrix_train, train_labels)\n",
    "    print (\"The number of nonzero weights with L2 regularization is {}.\".format(np.count_nonzero(model.coef_)))\n",
    "    print (\"  Accuracy is {}.\".format(model.score(doc_matrix_dev, dev_labels)))\n",
    "\n",
    "# Logistic regression using L1 penalty with the default value of C (1.0)\n",
    "    model = LogisticRegression(penalty='l1', C=1.0)\n",
    "    model.fit(doc_matrix_train, train_labels)\n",
    "    print (\"The number of nonzero weights with L1 regularization is {}.\".format(np.count_nonzero(model.coef_)))\n",
    "    print (\"  Accuracy is {}.\".format(model.score(doc_matrix_dev, dev_labels)))\n",
    "    \n",
    "# Retrain using features with nonzero weights\n",
    "    # [1] gives the index of the features ([0]: index of class), np.unique yields nonrepeating indices\n",
    "    nonzero_features = np.unique(np.nonzero(model.coef_)[1])\n",
    "    # Extract vocab with nonzero features\n",
    "    nonzero_vocab = np.array(vectorizer.get_feature_names())[nonzero_features]\n",
    "    \n",
    "# Vectorize with new vocab\n",
    "    vectorizer = CountVectorizer(vocabulary=nonzero_vocab)\n",
    "    doc_matrix_train = vectorizer.fit_transform(train_data)\n",
    "    doc_matrix_dev = vectorizer.transform(dev_data)\n",
    "\n",
    "# Fit L2 model with the reduced set of features\n",
    "    model = LogisticRegression(penalty='l2', C=1.0)\n",
    "    model.fit(doc_matrix_train, train_labels)\n",
    "    print (\"The number of nonzero weights with L2 regularization and reduced vocab is {}.\".format(np.count_nonzero(model.coef_)))\n",
    "    print (\"  It is much smaller after features are pruned by the L1 model.\")\n",
    "    print (\"  Accuracy is {}.\".format(model.score(doc_matrix_dev, dev_labels)))\n",
    "    \n",
    "# Initialize C's and empty arrays for vocab_size and accuracy\n",
    "    Cs = [0.01, 0.02, 0.05, 0.1, 0.25, 0.5, 1, 3, 5, 10]\n",
    "    vocab_size = []\n",
    "    accuracy = []\n",
    "    \n",
    "# Iterate through different C's\n",
    "    for val in Cs:\n",
    "    # First vectorize with standard vocab\n",
    "        vectorizer = CountVectorizer()\n",
    "        doc_matrix_train = vectorizer.fit_transform(train_data)\n",
    "    # Fit a L1 model with C value and tol=0.01 to ensure convergence. Lower C value prunes more features.\n",
    "        model = LogisticRegression(penalty='l1', C=val, tol=0.01)\n",
    "        model.fit(doc_matrix_train, train_labels)\n",
    "    # Obtain features with nonzero weights\n",
    "        nonzero_features = np.unique(np.nonzero(model.coef_)[1])\n",
    "        nonzero_vocab = np.array(vectorizer.get_feature_names())[nonzero_features]\n",
    "        vocab_size.append(len(nonzero_vocab))\n",
    "    # Vectorize with new vocab\n",
    "        vectorizer = CountVectorizer(vocabulary=nonzero_vocab)\n",
    "        doc_matrix_train = vectorizer.fit_transform(train_data)\n",
    "        doc_matrix_dev = vectorizer.transform(dev_data)\n",
    "    # Fit L2 model with default C=1.0\n",
    "        model = LogisticRegression(penalty='l2', C=1.0, tol=0.01)\n",
    "        model.fit(doc_matrix_train, train_labels)\n",
    "        accuracy.append(model.score(doc_matrix_dev, dev_labels))\n",
    "    \n",
    "# Create plot\n",
    "    plt.plot(vocab_size, accuracy, 'bo')\n",
    "    plt.xlabel(\"Vocabulary Size (controlled by C parameter in L1 regularization)\")\n",
    "    plt.ylabel(\"Accuracy\")\n",
    "    plt.title(\"Accuracy of Logistic Regression with L2 Regularization as a function of Vocabulary Size\")\n",
    "    \n",
    "    ### STUDENT END ###\n",
    "    \n",
    "P6()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(7) Use the TfidfVectorizer -- how is this different from the CountVectorizer? Train a logistic regression model with C=100.\n",
    "\n",
    "Make predictions on the dev data and show the top 3 documents where the ratio R is largest, where R is:\n",
    "\n",
    "maximum predicted probability / predicted probability of the correct label\n",
    "\n",
    "What kinds of mistakes is the model making? Suggest a way to address one particular issue that you see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R ratio: 287.179113053\n",
      "Vector of predicted probabilities: [ 0.00242094  0.02006527  0.28226901  0.69524478]\n",
      "Correct category: 0\n",
      "Correct category name: alt.atheism\n",
      "Document below:\n",
      "\n",
      "The 24 children were, of course, killed by a lone gunman in a second story\n",
      "window, who fired eight bullets in the space of two seconds...\n",
      "\n",
      "\n",
      "R ratio: 325.004211991\n",
      "Vector of predicted probabilities: [ 0.00291046  0.97974559  0.01432939  0.00301456]\n",
      "Correct category: 3\n",
      "Correct category name: talk.religion.misc\n",
      "Document below:\n",
      "Can anyone provide me a ftp site where I can obtain a online version\n",
      "of the Book of Mormon. Please email the internet address if possible.\n",
      "\n",
      "R ratio: 929.358963659\n",
      "Vector of predicted probabilities: [ 0.00196563  0.99372025  0.00324487  0.00106925]\n",
      "Correct category: 3\n",
      "Correct category name: talk.religion.misc\n",
      "Document below:\n",
      "I am pleased to announce that a *revised version* of _The Easy-to-Read Book\n",
      "of Mormon_ (former title: _Mormon's Book_) by Lynn Matthews Anderson is now\n",
      "available through anonymous ftp (see information below). In addition to the\n",
      "change in title, the revised ETR BOM has been shortened by several pages\n",
      "(eliminating many extraneous \"that's\" and \"of's\"), and many (minor) errors\n",
      "have been corrected. This release includes a simplified Joseph Smith Story,\n",
      "testimonies of the three and eight witnesses, and a \"Words-to-Know\"\n",
      "glossary.\n",
      "\n",
      "As with the previous announcement, readers are reminded that this is a\n",
      "not-for-profit endeavor. This is a copyrighted work, but people are welcome\n",
      "to make *verbatim* copies for personal use. People can recuperate the\n",
      "actual costs of printing (paper, copy center charges), but may not charge\n",
      "anything for their time in making copies, or in any way realize a profit\n",
      "from the use of this book. See the permissions notice in the book itself\n",
      "for the precise terms.\n",
      "\n",
      "Negotiations are currently underway with a Mormon publisher vis-a-vis the\n",
      "printing and distribution of bound books. (Sorry, I'm out of the wire-bound\n",
      "\"first editions.\") I will make another announcement about the availability\n",
      "of printed copies once everything has been worked out.\n",
      "\n",
      "FTP information: connect via anonymous ftp to carnot.itc.cmu.edu, then \"cd\n",
      "pub\" (you won't see anything at all until you do).\n",
      "\n",
      "\"The Easy-to-Read Book of Mormon\" is currently available in postscript and\n",
      "RTF (rich text format). (ASCII, LaTeX, and other versions can be made\n",
      "available; contact dba@andrew.cmu.edu for details.) You should be able to\n",
      "print the postscript file on any postscript printer (such as an Apple\n",
      "Laserwriter); let dba know if you have any difficulties. (The postscript in\n",
      "the last release had problems on some printers; this time it should work\n",
      "better.) RTF is a standard document interchange format that can be read in\n",
      "by a number of word processors, including Microsoft Word for both the\n",
      "Macintosh and Windows. If you don't have a postscript printer, you may be\n",
      "able to use the RTF file to print out a copy of the book.\n",
      "\n",
      "-r--r--r--  1 dba                   1984742 Apr 27 13:12 etrbom.ps\n",
      "-r--r--r--  1 dba                   1209071 Apr 27 13:13 etrbom.rtf\n",
      "\n",
      "For more information about how this project came about, please refer to my\n",
      "article in the current issue of _Sunstone_, entitled \"Delighting in\n",
      "Plainness: Issues Surrounding a Simple Modern English Book of Mormon.\"\n",
      "\n",
      "Send all inquiries and comments to:\n",
      "\n",
      "    Lynn Matthews Anderson\n",
      "    5806 Hampton Street\n",
      "    Pittsburgh, PA 15206\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def P7():\n",
    "### STUDENT START ###\n",
    "# Vectorize with TFIDF vectorizer\n",
    "    vectorizer = TfidfVectorizer()\n",
    "    doc_matrix_train = vectorizer.fit_transform(train_data)\n",
    "    doc_matrix_dev = vectorizer.transform(dev_data)\n",
    "# Logistic regression\n",
    "    model = LogisticRegression(C=100)\n",
    "    model.fit(doc_matrix_train, train_labels)\n",
    "    probs = model.predict_proba(doc_matrix_dev)\n",
    "# Calculate R ratio\n",
    "    R_ratio = []\n",
    "    # Iterate through each dev example\n",
    "    for index, example in enumerate(probs):\n",
    "        prob_max = np.amax(example)\n",
    "        prob_correct = example[dev_labels[index]]\n",
    "        R_ratio.append(prob_max / prob_correct)\n",
    "# Print top 3 documents with largest R ratios\n",
    "    R_ratio = np.array(R_ratio)\n",
    "    # Sort R ratios and find the largest 3\n",
    "    top_3_R = np.argsort(R_ratio)[-3:]\n",
    "    for index in top_3_R:\n",
    "        print (\"R ratio: {}\".format(R_ratio[index]))\n",
    "        print (\"Vector of predicted probabilities: {}\".format(probs[index]))\n",
    "        print (\"Correct category: {}\".format(dev_labels[index]))\n",
    "        print (\"Correct category name: {}\".format(newsgroups_train.target_names[dev_labels[index]]))\n",
    "        print (\"Document below:\")\n",
    "        print (dev_data[index] + \"\\n\")\n",
    "    \n",
    "## STUDENT END ###\n",
    "P7()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ANSWER:\n",
    "Whereas CountVectorizer converts the document into a matrix of counts, the TfidfVectorizer converts it into a matrix of tfidf values. A word that appears many times in a certain document would get a high count in CountVectorizer. But this word (say \"you\") appears in many documents, the idf part of tfidf puts a lower weight on this word, as it is not as predictive of the document's class.  \n",
    "\n",
    "From the above 3 examples with the largest R ratios (i.e. model predictions were way off), we can see that religious posts from talk.religion.misc are mischaracterized as comp.graphics. In both cases, the posts talked about the Book of Mormon, but also mentioned \"ftp\", \"online\", and \"internet\". These terms combined likely have stronger weights than the term \"Mormon\". One idea would be to increase influence of proper nouns like \"Book of Mormon\", which could give more context than terms like \"ftp\" and \"internet\". This could be done by repeating them a few times in the preprocessor to increase their contribution to the topic.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(8) EXTRA CREDIT\n",
    "\n",
    "Try implementing one of your ideas based on your error analysis. Use logistic regression as your underlying model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}