Labs/Lab04/Exercises04-Solutions.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:fafd7c3ae26dce284ffd15c8f4b6ed2040315e4c67a285ea58b7b5ad2779a79f"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import sys\n",
      "import glob\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "import pandas as pd\n",
      "%matplotlib inline\n",
      "%precision 4\n",
      "plt.style.use('ggplot')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.linalg as la\n",
      "import scipy.stats as st"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cPickle"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.sparse import csc_matrix \n",
      "from sparsesvd import sparsesvd "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.set_printoptions(formatter={'float': '{: 0.3f}'.format})"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Latent Semantic Analysis (LSA) is a method for reducing the dimnesionality of documents treated as a bag of words. It is used for document classification, clustering and retrieval. For example, LSA can be used to search for prior art given a new patent application. In this homework, we will implement a small library for simple latent semantic analysis as a practical example of the application of SVD. The ideas are very similar to PCA."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 1 (10 points)**.  Calculating pairwise distance matrices.\n",
      "\n",
      "Suppose we want to construct a distance matrix between the rows of a matrix. For example, given the matrix \n",
      "\n",
      "```python\n",
      "M = np.array([[1,2,3],[4,5,6]])\n",
      "```\n",
      "\n",
      "we want to find the new matrix\n",
      "\n",
      "````python\n",
      "D = np.array([[distance([1,2,3], [1,2,3]), distance([1,2,3], [4,5,6])],\n",
      "              [distance([4,5,6], [1,2,3]), distance([4,5,6], [4,5,6])]])\n",
      "```\n",
      "\n",
      "where `distance` is some appropriate function of two vectors (e.g. squared Euclidean).\n",
      "\n",
      "Write a function to calculate the pairwise-distance matrix given the matrix $M$ and some abritrary distance function. Your functions should have the following signature:\n",
      "```\n",
      "def func_name(M, distance_func):\n",
      "    pass\n",
      "```\n",
      "\n",
      "0. Write a distance function for the Euclidean, squared Euclidean and cosine measures.\n",
      "1. Write the function using looping for M as a collection of row vecotrs.\n",
      "2. Write the function using looping for M as a collection of column vectors.\n",
      "3. Wrtie the function using broadcasting for M as a colleciton of row vectors.\n",
      "4. Write the function using broadcasting for M as a colleciton of column vecotrs. \n",
      "\n",
      "For 3 and 4, try to avoid using transposition. Check that all four funcitons give the same result when applied to the given matrix $M$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Your code here\n",
      "\n",
      "def squared_euclidean_norm(u, axis=-1):\n",
      "    return (u**2).sum(axis)\n",
      "\n",
      "def euclidean_norm(u, axis=-1):\n",
      "    return np.sqrt(squared_euclidean_norm(u, axis))\n",
      "\n",
      "def squared_euclidean_dist(u, v, axis=-1):\n",
      "    \"\"\"Returns squared Euclidean distance between two vectors.\"\"\"\n",
      "    return squared_euclidean_norm(u-v, axis)\n",
      "\n",
      "def euclidean_dist(u, v, axis=-1):\n",
      "    \"\"\"Return Euclidean distacne between two vectors.\"\"\"\n",
      "    return np.sqrt(squared_euclidean_dist(u, v, axis))\n",
      "    \n",
      "def cosine_dist(u, v, axis=-1):\n",
      "    \"\"\"Returns cosine of angle betwwen two vectors.\"\"\"\n",
      "    # return 1 - np.dot(u, v)/(la.norm(u)*la.norm(v))\n",
      "    return 1 - (u * v).sum(axis)/(euclidean_norm(u, axis) * euclidean_norm(v, axis))\n",
      "\n",
      "def loop_row_pdist(M, f):\n",
      "    \"\"\"REturns pairwise-distance matrix assuming M consists of row vectors..\"\"\"\n",
      "    nrows, ncols = M.shape\n",
      "    return np.array([[f(M[u,:], M[v,:]) for u in range(nrows)] \n",
      "                                        for v in range(nrows)])\n",
      "\n",
      "def loop_col_pdist(M, f):\n",
      "    \"\"\"REturns pairwise-distance matrix assuming M consists of column vectors..\"\"\"\n",
      "    nrows, ncols = M.shape\n",
      "    return np.array([[f(M[:,u], M[:,v]) for u in range(ncols)] \n",
      "                                        for v in range(ncols)])\n",
      "\n",
      "def broadcast_row_pdist(M, f):\n",
      "    \"\"\"REturns pairwise-distance matrix assuming M consists of row vectors..\"\"\"\n",
      "    return f(M[None,:,:], M[:,None,:])\n",
      "\n",
      "def broadcast_col_pdist(M, f):\n",
      "    \"\"\"REturns pairwise-distance matrix assuming M consists of column vectors..\"\"\"\n",
      "    return f(M[:,None,:], M[:,:,None], axis=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = np.array([[1,2,3],[4,5,6]])\n",
      "\n",
      "# dist = euclidean_dist\n",
      "dist = cosine_dist\n",
      "\n",
      "print loop_row_pdist(M, dist), '\\n'\n",
      "print broadcast_row_pdist(M, dist), '\\n'\n",
      "print loop_col_pdist(M, dist), '\\n'\n",
      "print broadcast_col_pdist(M, dist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.000  0.025]\n",
        " [ 0.025  0.000]] \n",
        "\n",
        "[[ 0.000  0.025]\n",
        " [ 0.025  0.000]] \n",
        "\n",
        "[[ 0.000  0.009  0.024]\n",
        " [ 0.009 -0.000  0.003]\n",
        " [ 0.024  0.003  0.000]] \n",
        "\n",
        "[[ 0.000  0.009  0.024]\n",
        " [ 0.009 -0.000  0.003]\n",
        " [ 0.024  0.003  0.000]]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 2 (10 points)**. Write 3 functions to calculate the term frequency (tf), the inverse document frequency (idf) and the product (tf-idf). Each fucntion should take a single argument `docs`, which is a dictionary of (key=identifier, value=dcoument text) pairs, and return appropriately sized array. Convert '-' to ' ' (space), remove punctuation, covert text to lowercase and split on whitepsace to genrate a collection of terms from the dcoument text.\n",
      "\n",
      "- tf = the number of occurrences of term i in document j\n",
      "- idf = $\\log \\frac{n}{1 + \\text{df}_i}$ where $n$ is the total number of documents and $\\text{df}_i$ is the number of documents in which term i occurs.\n",
      "\n",
      "Print the table of tf-idf values for the follwoing document collection\n",
      "\n",
      "```\n",
      "s1 = \"The quick brown fox\"\n",
      "s2 = \"Brown fox jumps over the jumps jumps jumps\"\n",
      "s3 = \"The the the lazy dog elephant.\"\n",
      "s4 = \"The the the the the dog peacock lion tiger elephant\"\n",
      "\n",
      "docs = {'s1': s1, 's2': s2, 's3': s3, 's4': s4}\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def tf(doc):\n",
      "    \"\"\"Returns the number of times each term occurs in a dcoument.\n",
      "    We preprocess the document to strip punctuation and convert to lowercase.\n",
      "    Terms are found by splitting on whitespace.\"\"\"\n",
      "    from collections import Counter\n",
      "    from string import punctuation\n",
      "\n",
      "    terms = doc.lower().replace('-', ' ').translate(None, punctuation).split()\n",
      "    return Counter(terms)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def tfs(docs):\n",
      "    \"\"\"Create a term freqeuncy dataframe from a dictionary of documents.\"\"\"\n",
      "    from operator import add\n",
      "\n",
      "    df = pd.DataFrame({k: tf(v) for k, v in docs.iteritems()}).fillna(0)\n",
      "    return df"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def idf(docs):\n",
      "    \"\"\"Find inverse document frequecny series from a dictionry of doucmnets.\"\"\"\n",
      "    term_freq = tfs(docs)\n",
      "    num_docs = len(docs)\n",
      "    doc_freq = (term_freq > 0).sum(axis=1)\n",
      "    return np.log(num_docs/(1 + doc_freq))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def tf_idf(docs):\n",
      "    \"\"\"Return the product of the term-frequency and inverse document freqeucny.\"\"\"\n",
      "    return tfs(docs).mul(idf(docs), axis=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s1 = \"The quick brown fox\"\n",
      "s2 = \"Brown fox jumps over the jumps jumps jumps\"\n",
      "s3 = \"The the the lazy dog elephant.\"\n",
      "s4 = \"The the the the the dog peacock lion tiger elephant\"\n",
      "\n",
      "docs = {'s1': s1, 's2': s2, 's3': s3, 's4': s4}\n",
      "\n",
      "tf_idf(docs)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>s1</th>\n",
        "      <th>s2</th>\n",
        "      <th>s3</th>\n",
        "      <th>s4</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>brown</th>\n",
        "      <td> 0.287682</td>\n",
        "      <td> 0.287682</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>dog</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.287682</td>\n",
        "      <td> 0.287682</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>elephant</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.287682</td>\n",
        "      <td> 0.287682</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>fox</th>\n",
        "      <td> 0.287682</td>\n",
        "      <td> 0.287682</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>jumps</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 2.772589</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>lazy</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.693147</td>\n",
        "      <td> 0.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>lion</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.693147</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>over</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.693147</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>peacock</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.693147</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>quick</th>\n",
        "      <td> 0.693147</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>the</th>\n",
        "      <td>-0.223144</td>\n",
        "      <td>-0.223144</td>\n",
        "      <td>-0.669431</td>\n",
        "      <td>-1.115718</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>tiger</th>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.000000</td>\n",
        "      <td> 0.693147</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "                s1        s2        s3        s4\n",
        "brown     0.287682  0.287682  0.000000  0.000000\n",
        "dog       0.000000  0.000000  0.287682  0.287682\n",
        "elephant  0.000000  0.000000  0.287682  0.287682\n",
        "fox       0.287682  0.287682  0.000000  0.000000\n",
        "jumps     0.000000  2.772589  0.000000  0.000000\n",
        "lazy      0.000000  0.000000  0.693147  0.000000\n",
        "lion      0.000000  0.000000  0.000000  0.693147\n",
        "over      0.000000  0.693147  0.000000  0.000000\n",
        "peacock   0.000000  0.000000  0.000000  0.693147\n",
        "quick     0.693147  0.000000  0.000000  0.000000\n",
        "the      -0.223144 -0.223144 -0.669431 -1.115718\n",
        "tiger     0.000000  0.000000  0.000000  0.693147"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 3 (10 points)**. \n",
      "\n",
      "1. Write a function that takes a matrix $M$ and an integer $k$ as arguments, and reconstructs a reduced matrix using only the $k$ largest singular values. Use the `scipy.linagl.svd` function to perform the decomposition.\n",
      "\n",
      "2. Apply the function you just wrote to the following term-frequency matrix for a set of 9 documents using k=2 and print the reconstructed matrix $M'$.\n",
      "```\n",
      "M = np.array([[1, 0, 0, 1, 0, 0, 0, 0, 0],\n",
      "       [1, 0, 1, 0, 0, 0, 0, 0, 0],\n",
      "       [1, 1, 0, 0, 0, 0, 0, 0, 0],\n",
      "       [0, 1, 1, 0, 1, 0, 0, 0, 0],\n",
      "       [0, 1, 1, 2, 0, 0, 0, 0, 0],\n",
      "       [0, 1, 0, 0, 1, 0, 0, 0, 0],\n",
      "       [0, 1, 0, 0, 1, 0, 0, 0, 0],\n",
      "       [0, 0, 1, 1, 0, 0, 0, 0, 0],\n",
      "       [0, 1, 0, 0, 0, 0, 0, 0, 1],\n",
      "       [0, 0, 0, 0, 0, 1, 1, 1, 0],\n",
      "       [0, 0, 0, 0, 0, 0, 1, 1, 1],\n",
      "       [0, 0, 0, 0, 0, 0, 0, 1, 1]])\n",
      "```\n",
      "\n",
      "3. Calculate the pairwise correlation matrix for the original matrix M and the reconstructed matrix using $k=2$ singular values (you may use [scipy.stats.spearmanr](http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.spearmanr.html) to do teh calculations). Consider the fist 5 sets of documents as one group $G1$ and the last 4 as another group $G2$ (i.e. first 5 and last 4 columns). What is the average within group correlation for $G1$, $G2$ and the average cross-group correlaiton for G1-G2 using either $M$ or $M'$. (Do not include self-correlation in the within-group calculations.)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Your code here\n",
      "\n",
      "def svd_projection(M, k):\n",
      "    \"\"\"Returns the matrix M reconstructed using only k singluar values\"\"\"\n",
      "    U, s, V = la.svd(M, full_matrices=False)\n",
      "    s[k:] = 0\n",
      "    M_ = U.dot(np.diag(s).dot(V))\n",
      "    \n",
      "    try:\n",
      "        return pd.DataFrame(M_, index=M.index, columns=M.columns)\n",
      "    except AttributeError:\n",
      "        return M_"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = np.array([[1, 0, 0, 1, 0, 0, 0, 0, 0],\n",
      "    [1, 0, 1, 0, 0, 0, 0, 0, 0],\n",
      "    [1, 1, 0, 0, 0, 0, 0, 0, 0],\n",
      "    [0, 1, 1, 0, 1, 0, 0, 0, 0],\n",
      "    [0, 1, 1, 2, 0, 0, 0, 0, 0],\n",
      "    [0, 1, 0, 0, 1, 0, 0, 0, 0],\n",
      "    [0, 1, 0, 0, 1, 0, 0, 0, 0],\n",
      "    [0, 0, 1, 1, 0, 0, 0, 0, 0],\n",
      "    [0, 1, 0, 0, 0, 0, 0, 0, 1],\n",
      "    [0, 0, 0, 0, 0, 1, 1, 1, 0],\n",
      "    [0, 0, 0, 0, 0, 0, 1, 1, 1],\n",
      "    [0, 0, 0, 0, 0, 0, 0, 1, 1]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Md = svd_projection(M, 2)\n",
      "Md"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "array([[ 0.162,  0.400,  0.379,  0.468,  0.176, -0.053, -0.115, -0.159,\n",
        "        -0.092],\n",
        "       [ 0.141,  0.370,  0.329,  0.400,  0.165, -0.033, -0.071, -0.097,\n",
        "        -0.043],\n",
        "       [ 0.152,  0.505,  0.358,  0.410,  0.236,  0.024,  0.060,  0.087,\n",
        "         0.124],\n",
        "       [ 0.258,  0.841,  0.606,  0.697,  0.392,  0.033,  0.083,  0.122,\n",
        "         0.187],\n",
        "       [ 0.449,  1.234,  1.051,  1.266,  0.556, -0.074, -0.155, -0.210,\n",
        "        -0.049],\n",
        "       [ 0.160,  0.582,  0.375,  0.417,  0.277,  0.056,  0.132,  0.189,\n",
        "         0.217],\n",
        "       [ 0.160,  0.582,  0.375,  0.417,  0.277,  0.056,  0.132,  0.189,\n",
        "         0.217],\n",
        "       [ 0.218,  0.550,  0.511,  0.628,  0.243, -0.065, -0.143, -0.197,\n",
        "        -0.108],\n",
        "       [ 0.097,  0.532,  0.230,  0.212,  0.267,  0.137,  0.315,  0.444,\n",
        "         0.425],\n",
        "       [-0.061,  0.232, -0.139, -0.266,  0.145,  0.240,  0.546,  0.767,\n",
        "         0.664],\n",
        "       [-0.065,  0.335, -0.146, -0.301,  0.203,  0.306,  0.695,  0.977,\n",
        "         0.849],\n",
        "       [-0.043,  0.254, -0.097, -0.208,  0.152,  0.221,  0.503,  0.707,\n",
        "         0.616]])"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rho, pval = st.spearmanr(Md)\n",
      "rho"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "array([[ 1.000,  0.846,  1.000,  1.000,  0.719, -0.837, -0.837, -0.837,\n",
        "        -0.804],\n",
        "       [ 0.846,  1.000,  0.846,  0.846,  0.972, -0.557, -0.557, -0.557,\n",
        "        -0.481],\n",
        "       [ 1.000,  0.846,  1.000,  1.000,  0.719, -0.837, -0.837, -0.837,\n",
        "        -0.804],\n",
        "       [ 1.000,  0.846,  1.000,  1.000,  0.719, -0.837, -0.837, -0.837,\n",
        "        -0.804],\n",
        "       [ 0.719,  0.972,  0.719,  0.719,  1.000, -0.389, -0.389, -0.389,\n",
        "        -0.298],\n",
        "       [-0.837, -0.557, -0.837, -0.837, -0.389,  1.000,  1.000,  1.000,\n",
        "         0.977],\n",
        "       [-0.837, -0.557, -0.837, -0.837, -0.389,  1.000,  1.000,  1.000,\n",
        "         0.977],\n",
        "       [-0.837, -0.557, -0.837, -0.837, -0.389,  1.000,  1.000,  1.000,\n",
        "         0.977],\n",
        "       [-0.804, -0.481, -0.804, -0.804, -0.298,  0.977,  0.977,  0.977,\n",
        "         1.000]])"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# G1/G1, G2/G2 and G1/G2 average correlation\n",
      "\n",
      "np.mean(rho[:5, :5][np.tril_indices_from(rho[:5, :5], 1)]), \\\n",
      "np.mean(rho[5:, 5:][np.tril_indices_from(rho[5:, 5:], 1)]), \\\n",
      "rho[5:, :5].mean() "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "(0.8988, 0.9930, -0.6780)"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 4 (20 points)**. Clustering with LSA\n",
      "\n",
      "1. Begin by loading a pubmed database of selected articles using 'cPickle'. With the following:\n",
      "\n",
      "```\n",
      "import cPickle\n",
      "\n",
      "docs = cPickle.load(open('pubmed.pic'))```\n",
      "\n",
      "Create a tf-idf matrix for every term that appears at least once in any of the documents. What is the shape of the tf-idf matrix? \n",
      "\n",
      "2. Perform SVD on the tf-idf matrix to obtain $U \\Sigma V^T$ (often written as $T \\Sigma D^T$ in this context with $T$ representing the termas and $D$ representing the documents). If we set all but the top $k$ singular values to 0, the reconstructed matrix is essentially $U_k \\Sigma_k V_k^T$, where $U_k$ is $m \\times k$, $\\Sigma_k$ is $k \\times k$ and $V_k^T$ is $k \\times n$. Terms in this reduced space are represented by $U_k \\Sigma_k$ and documents by $\\Sigma_k V^T_k$. \n",
      "\n",
      "Reconstruct the matrix using the first $10$ singular values.\n",
      "\n",
      "**Note**: The SVD from scipy.linalg performs a full decomposition, which is too slow for such a large matrix. We only need to decompose unitl we get the first k singluar values. Please use the `sparsesvd` function from the [sparsesvd](https://pypi.python.org/pypi/sparsesvd/) package to perform SVD instead only for $k = 2$. You can install in the usual way with \n",
      "```\n",
      "!pip install sparsesvd\n",
      "```\n",
      "\n",
      "3. Use agglomerative hierachical clustering with complete linkage to plot a dendrogram and comment on the likely number of  document clusters with $k = 100$. Use the dendrogram function from [SciPy ](https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.cluster.hierarchy.dendrogram.html) and `sparsesvd` as before.\n",
      "\n",
      "4. Determine how similar each of the origianl 16 documents in to the new document `mystery.txt`. Since $A = U \\Sigma V^T$, we also have $V = A^T U S^{-1}$ using orthogonality and the rule for transposing matrix products. This suggest that in order to map the new doucment to the same concept space, first find the tf-idf vector $v$ for the new document - this must contain all (and only) the terms present in the existing tf-idx matrix. Then the query vector $q$ is given by $v^T U_k \\Sigma_k^{-1}$. Rank the 16 document vectors with respect to  $q$ using the cosine distance.\n",
      "\n",
      "5. Does it matter that the front and back matter of each document is essentially identical for either LSA-based clustering (part 3) or information retrieval (part 4)? Why or why not?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Your code here\n",
      "import cPickle\n",
      "\n",
      "docs = cPickle.load(open('pubmed.pic'))\n",
      "df = tf_idf(docs)\n",
      "df.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "(6488, 178)"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# %time\n",
      "# T, s, D = la.svd(df)\n",
      "# T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%time\n",
      "k = 10\n",
      "T, s, D = sparsesvd(csc_matrix(df), k=10)\n",
      "x = np.diag(s).dot(D)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CPU times: user 4 \u00b5s, sys: 1 \u00b5s, total: 5 \u00b5s\n",
        "Wall time: 8.11 \u00b5s\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print df.shape, T.shape, s.shape, D.shape\n",
      "T.T.dot(np.diag(s).dot(D)).shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(6488, 178) (100, 6488) (100,) (100, 178)\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 33,
       "text": [
        "(6488, 178)"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "doc = {'mystery': open('mystery.txt').read()}\n",
      "terms = tf_idf(doc)\n",
      "query_terms = df.join(terms).fillna(0)['mystery']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "k = 10\n",
      "T, s, D = sparsesvd(csc_matrix(df), k=100)\n",
      "q = query_terms.T.dot(T.T.dot(np.diag(1.0/s)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ranked_docs = df.columns[np.argsort(cosine_dist(q, x))][::-1]\n",
      "print \"Query article:\", \n",
      "print ' '.join(line.strip() for line in doc['mystery'].splitlines()[:2])\n",
      "print\n",
      "print \"Most similar\"\n",
      "print '='*80\n",
      "for i, title in enumerate(ranked_docs[:10]):\n",
      "    print '%03d' % i, title\n",
      "\n",
      "\n",
      "print\n",
      "print \"Most dissimilar\"\n",
      "print '='*80\n",
      "for i, title in enumerate(ranked_docs[-10:]):\n",
      "    print '%03d' % (len(docs) - i), title"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "operands could not be broadcast together with shapes (100,) (10,178) ",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-23-6c88d46bb375>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mranked_docs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margsort\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcosine_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mprint\u001b[0m \u001b[0;34m\"Query article:\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mprint\u001b[0m \u001b[0;34m' '\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstrip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdoc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mystery'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msplitlines\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mprint\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;32mprint\u001b[0m \u001b[0;34m\"Most similar\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m<ipython-input-6-7dba4e22bd75>\u001b[0m in \u001b[0;36mcosine_dist\u001b[0;34m(u, v, axis)\u001b[0m\n\u001b[1;32m     18\u001b[0m     \u001b[0;34m\"\"\"Returns cosine of angle betwwen two vectors.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0;31m# return 1 - np.dot(u, v)/(la.norm(u)*la.norm(v))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 20\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0meuclidean_norm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0meuclidean_norm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     21\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     22\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mloop_row_pdist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (100,) (10,178) "
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.cluster.hierarchy import linkage, dendrogram\n",
      "from scipy.spatial.distance import pdist, squareform\n",
      "\n",
      "plt.figure(figsize=(16,36))\n",
      "k = 150\n",
      "T, s, D = sparsesvd(csc_matrix(df), k=100)\n",
      "x = np.diag(s).dot(D).T\n",
      "data_dist = pdist(x, metric='cosine') # computing the distance\n",
      "data_link = linkage(data_dist) # computing the linkage\n",
      "labels = [c[:40] for c in df.columns[:]]\n",
      "dendrogram(data_link, orientation='right', labels=labels)\n",
      "plt.xlabel('Samples')\n",
      "plt.ylabel('Distance')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}