Showing posts with label map-reduce. Show all posts
Showing posts with label map-reduce. Show all posts

Friday, July 18, 2014

Clustering Medical Procedure Codes with Scalding


A colleague pointed out that there exists an inverse use-case to finding outliers in medical claims. This one is to group procedure codes into clusters, or what Health Insurance companies call Episode Treatment Groups (ETG). Essentially an ETG is a way to cluster a group of services (procedures) into a medically relevant unit.

The CMS.gov dataset provides slightly under 16 million anonymized outpatient claims for Medicare/Medicaid patients. Each outpatient record can have upto 6 ICD-9 procedure codes, upto 10 ICD-9 diagnosis codes and upto 45 HCPCS codes. So just like the outlier case, we can derive a measure of similarity between a pair of codes as the average co-occurrence within claims across the dataset.

I decided to use a variant of the DBSCAN clustering algorithm. This post provides some tips on how to implement DBSCAN in a distributed manner - I used the ideas in this post to develop my implementation. The intuition behind my clustering algorithm goes something like this.

We calculate the similarity sAB between a pair of codes A and B as the number of times they co-occur in the corpus. Clustering algorithms need a distance measure, so we treat the distance dAB as the reciprocal of their similarity, ie 1/sAB. The DBSCAN clustering algorithm works by selecting other points around each point that are within a specified distance ε from each other. Candidate cluster centroids are those that have at least MinPoints codes within this distance ε. My algorithm deviates from DBSCAN at this point - instead of finding density-reachable codes I just find the Top-N densest clusters. Density is calculated as the number of codes within a circular area of the mean radius, i.e. N2 / πΣi=0..Nd2. We then calculate the top N densest code clusters - these are our derived ETGs.

The Scalding code below does just this. We simplify a bit by not calculating using some constants such as π but otherwise the code is quite faithful to the algorithm described above.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
// Source: src/main/scala/com/mycompany/cmspp/cluster/CodeCluster.scala
package com.mycompany.cmspp.clusters

import com.twitter.scalding.Job
import com.twitter.scalding.Args
import com.twitter.scalding.TextLine
import com.twitter.scalding.Tsv
import scala.io.Source

class CodeCluster(args: Args) extends Job(args) {

  def extractPairs(line: String): List[(String,String)] = {
    val cols = line.split(",").toList
    val codes = (cols.slice(22, 27)  // ICD9 procedure code cols
      .map(x => if (x.isEmpty) x else "ICD9:" + x) 
      ::: cols.slice(31, 75)         // HCPCS (CPT4) procedure cols
      .map(x => if (x.isEmpty) x else "HCPCS:" + x))          
      .filter(x => (! x.isEmpty))
    val cjoin = for {codeA <- codes; codeB <- codes} yield (codeA, codeB)
    cjoin.filter(x => x._1 < x._2)
  }

  val Epsilon = args("epsilon").toDouble
  val MinPoints = args("minpoints").toInt
  val NumClusters = args("nclusters").toInt
  
  val output = Tsv(args("output"))

  val dists = TextLine(args("input"))
    .read
    // compute pair-wise distances between procedure codes
    .flatMapTo('line -> ('codeA, 'codeB)) { line: String => extractPairs(line) }
    .groupBy('codeA, 'codeB) { group => group.size('sim) }
    .map('sim -> 'radius) { x: Int => (1.0D / x) }
    .discard('sim)
    // group by codeA and retain only records which are within epsilon distance
    .groupBy('codeA) { group => group.sortBy('radius).reverse }
    .filter('radius) { x: Double => x < Epsilon }
    
  val codeCounts = dists
    .groupBy('codeA) { group => 
      group.sizeAveStdev('radius -> ('count, 'mean, 'std)) 
    }
    // only retain codes that have at least MinPoints points within Epsilon
    .filter('count) { x: Int => x > MinPoints }
    .discard('std)

  val densities = dists.joinWithSmaller(('codeA -> 'codeA), codeCounts)
    .map(('mean, 'count) -> 'density) { x: (Double,Int) => 
      1.0D * Math.pow(x._2, 2) / Math.pow(x._1, 2)
    }
    .discard('radius, 'count)
    
  // sort the result by density descending and find the top N clusters
  val densestCodes = densities.groupAll { group => 
      group.sortBy('density).reverse }
    .unique('codeA)
    .limit(NumClusters)
    
  // join code densities with densest codes to find final clusters
  densities.joinWithTiny(('codeA -> 'codeA), densestCodes)
    .groupBy('codeA) { group => group.mkString('codeB, ",")}
    .write(output)
}

object CodeCluster {
  def main(args: Array[String]): Unit = {
    // populate redis cache
    new CodeCluster(Args(List(
      "--local", "",
      "--epsilon", "0.3",
      "--minpoints", "10",
      "--nclusters", "10",
      "--input", "data/outpatient_claims.csv",
      "--output", "data/clusters.csv"
    ))).run
    Source.fromFile("data/clusters.csv")
      .getLines()
      .foreach(Console.println(_))
  }
}

I ran this locally with 1 million claims (out of the 16 million claims in my dataset) and got results like this:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
HCPCS:00142    HCPCS:85025,  HCPCS:36415, HCPCS:93005, HCPCS:80053, ...
HCPCS:00300    HCPCS:85025,  HCPCS:80053, HCPCS:36415, HCPCS:99284, ...
HCPCS:00400    HCPCS:85025,  HCPCS:93005, HCPCS:80053, HCPCS:36415, ...
HCPCS:00532    HCPCS:85025,  HCPCS:36415, HCPCS:80048, HCPCS:G0378, ...
HCPCS:0073T    HCPCS:77417,  HCPCS:77336, HCPCS:77427, HCPCS:77280, ...
HCPCS:00740    HCPCS:85025,  HCPCS:36415, HCPCS:93005, HCPCS:85610, ...
HCPCS:00750    HCPCS:36415,  HCPCS:85025, HCPCS:85610, HCPCS:J3010, ...
HCPCS:00790    HCPCS:36415,  HCPCS:85025, HCPCS:80048, HCPCS:80053, ...
HCPCS:00810    HCPCS:36415,  HCPCS:85025, HCPCS:93005, HCPCS:80053, ...
HCPCS:00830    HCPCS:36415,  HCPCS:85025, HCPCS:80048, HCPCS:93005, ...

[Edit: 07/22/2014: This approach does not produce clusters. Notice in the data the same HCPCS:85025 is part of the first 3 clusters, which is obviously not wanted. I will implement the last part of the DBSCAN algorithm and update this page when I am done.]

And thats all I have for today. I'd like to point out a new book on Scalding, Programming MapReduce with Scalding by Antonios Chalkiopoulos. I was quite impressed by this book, you can read my review on Amazon if you are interested.

Sunday, March 02, 2014

Cleaning UMLS data and Loading into Graph


The little UMLS ontology I am building needs to support two basic features in its user interface - findability and navigability. I now have a reasonable solution for the findability part, and I am planning to use Neo4j (a graph database) for the navigability part.

As before, the nodes are extracted from the MRCONSO table. The relationships between nodes are extracted from the MRREL table. Both SQL queries are shown below:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
mysql> select CUI, STR from MRCONSO
...    where LAT = 'ENG'
...    into outfile '/tmp/cuistr.csv'
...    fields terminated by '\t'
...    lines terminated by '\n';
Query OK, 7871075 rows affected (27.03 sec)

mysql> select CUI1, CUI2, RELA from MRREL 
...    into outfile '/tmp/cuirel.csv' 
...    fields terminated by '\t'
...    lines terminated by '\n';
Query OK, 58024739 rows affected (1 min 17.78 sec)

The Neo4j community seems to have standardized on Michael Hunger's batch-import tool for loading data into Neo4j. It takes as input tab separated files for the nodes and relationships, and writes out the graph into an embedded database. The node file(s) should specify a nodeId, and one or more properties separated by tabs. The relationship file(s) should specify the start node, end node, relationship name, and zero or more relationship properties separated by tabs.

Since my node file (cuistr.csv) was normalized (one row per synonym), I needed to transform this file to a (cui, list(str)) format. I decided to use MRJob (a Python based Map-Reduce framework from Yelp that you can use to run your jobs on Hadoop and Amazon EMR, although I just ran them locally) to write a little Map-Reduce job to do this.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
# syns_aggregator_job.py
from mrjob.job import MRJob

class SynsAggregatorJob(MRJob):
  """
  Groups unique synonyms by CUI. 
  Input format: (CUI,DESCR)
  Output format: (CUI,[DESCR,...])
  """

  def mapper(self, key, value):
    (cui, descr) = value.split("\t")
    yield cui, descr

  def reducer(self, key, values):
    uniqSyns = set()
    for value in values:
      uniqSyns.add(value)
    print "%s\t%s" % (key, list(uniqSyns))

if __name__ == "__main__":
  SynsAggregatorJob.run()

I ran this locally with the following command to aggregate the 7,871,075 records into an aggregated file cuistr_agg.csv with 2,880,385 records.

1
2
sujit@tsunami:umls$ python syns_aggregator_job.py \
    /path/to/cuistr.csv > /path/to/cuistr_agg.csv 

I also built another job to remove edges that referred to non-existent nodes. Notice that in the SQL I only retrieved English names (LAT='ENG'), and there is no corresponding filter on the MRREL query. This step is actually unnecessary because the batch-import tool checks and skips such rows, but I include it here anyway because it seems to me to be quite a nice way to remove non-existent rows without having to look up a dictionary. But if you are trying to replicate, you should skip doing this step.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
# rels_filter_job.py
from mrjob.job import MRJob
from mrjob.step import MRStep
from mrjob.compat import get_jobconf_value

class RelsFilterJob(MRJob):
  """
  Removes records from CUIREL where either node in a relation
  does not exist in CUISYN. Needs to be run twice - first run
  removes one non-existent CUI, second run removes second non
  existent CUI.
  Input format: (CUI, SYN_LIST) - from CuiSynsJob OR
                (CUI1, REL, CUI2)   - from cuirels.csv
  Output format: (CUI1, CUI2, REL)
  """

  def mapper_init(self):
    self.cui_idx = int(get_jobconf_value("cui_idx"))

  def mapper(self, key, value):
    ncols = len(value.split("\t"))
    if ncols == 2:
      # from the output of SynsAggregatorJob
      (cui, payload) = value.split("\t")
      yield (cui, "XXX")
    else:
      # from cuirels
      cols = value.split("\t")
      yield (cols[self.cui_idx], value)

  def reducer(self, key, values):
    # if one of the records in the reduced set has value XXX 
    # then all the values (except the XXX one) are good
    include = False
    vallist = []
    for value in values:
      if value == 'XXX':
        include = True
        continue
      vallist.append(value)
    if include:
      for value in vallist:
       print value
    
if __name__ == "__main__":
  RelsFilterJob.run()

I ran the above job twice, first to remove relationship rows which had non-existent source CUIs, and the second to remove ones with non-existent target CUIs. The JobConf parameter specifies which CUI to check. Here are the commands:

1
2
3
4
5
6
7
8
sujit@tsunami:umls$ python rels_filter_job.py \
    --jobconf cui_idx=0 \
    /path/to/cuistr_agg.csv /path/to/cuirel.csv > \
    /path/to/cuirel_filt_left.csv
sujit@tsunami:umls$ python rels_filter_job.py \
    --jobconf cui_idx=1 \
    /path/to/cuirel_agg.csv /path/to/cuirel_filt_left.csv > \
    /path/to/cuirel_filt_right.csv

This resulted in a much less dramatic reduction from 58,024,739 records in the source cuirel.csv file to 58,021,093 records in the cuirel_filt_right.csv target. However, as mentioned above, this step is unnecessary (and time-consuming), we can provide the cuirel.csv file to batch-import and it will do the right thing.

Batch Import did not run for me out of the box. In order to make it run, I had to parse the instructions on the README file multiple times, as well as do several searches on Neo4j's mailing list on Google Groups. I describe below what I had to do to make it run for my data.

My input files are cuistr_agg.csv (for the nodes) and cuirel_filt_right.csv (for the relationships). I needed to put headers on both of them to indicate to batch-import what the property names were and which fields I should be able to look up. This is because internally Neo4j uses longs to refer to node IDs - since my unique key for a node is the CUI (a string field), it creates a Lucene index to map the CUI to the internal node ID. Here are the first 10 rows from both files, showing the headers - the empty space between fields are tabs. The header creates a Lucene index called "concepts" that maps the string field "CUI" to the internal Neo4j nodeID.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# cuistr_agg.csv
cui:string:concepts  syns
C0000005             ['(131)I-MAA', '(131)I-Macroaggregated Albumin']
C0000039             ['Dipalmitoylglycerophosphocholine', ...]
C0000052             ['1,4-alpha-D-Glucan:1,4-alpha-D-glucan ...]
C0000074             ['1 Alkyl 2 Acylphosphatidates', ...]
C0000084             ['1 Carboxyglutamic Acid', ...]
C0000096             ['Isobutyltheophylline', ...]
C0000097             ['Methylphenyltetrahydropyridine (substance)', ...]
C0000098             ['1 Methyl 4 phenylpyridine', ...]
C0000102             ['a- Naphthylamine', '1 Naphthylamine', ...]
...

# cuirel_filt_right.csv
cui:string:concepts  cui:string:concepts  rela
C0000039             C0000039             entry_version_of
C0000039             C0000039             has_entry_version
C0000039             C0000039             has_permuted_term
C0000039             C0000039             has_permuted_term
C0000039             C0000039             has_permuted_term
C0000039             C0000039             has_permuted_term
C0000039             C0000039             has_sort_version
C0000039             C0000039             has_sort_version
C0000039             C0000039             has_translation
...

To download and compile batch-import, run the following commands:

1
2
3
sujit@tsunami:Downloads$ git clone https://github.com/jexp/batch-import.git
sujit@tsunami:Downloads$ cd batch-import
sujit@tsunami:batch-import$ mvn clean compile assembly:single

My first attempt to run the importer just hung. I needed to add the following two properties to the batch.properties file supplied with batch-import.

1
2
3
4
5
6
7
# create lucene index "concepts" for exact lookup
batch_import.node_index.concepts=exact

# input CSVs don't have quoted fields. Apparently this speeds
# things up considerably since it allows use of a simpler CSV
# parser.
batch_import.csv.quotes=false

Finally, batch-import sets an upper limit on the length of a property value (possibly for performance) in Chunker.BUFSIZE. Since I was using JSON-ified lists for synonyms, this field can be very long and the import would fail until I set Chunker.BUFSIZE from 32*1024 to 128*1024. I had to rebuild the JAR (mvn assembly:single) after this change. The following command created my Neo4j database in target/db and loaded my two files into it.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
sujit@tsunami:batch-import$ java -server -Dfile.encoding=UTF-8 -Xmx4G \
    -jar target/batch-import-jar-with-dependencies.jar \
    target/db /path/to/cuistr_agg.csv /path/to/cuirel_filt_right.csv

Using Existing Configuration File
............................
Importing 2880384 Nodes took 117 seconds 
..................................................... 58903 ms for 10000000
..................................................... 194346 ms for 10000000
..................................................... 419848 ms for 10000000
..................................................... 274616 ms for 10000000
..................................................... 507095 ms for 10000000
..................................................... 
Importing 58021093 Relationships took 1764 seconds 

Total import time: 1901 seconds 

In order to verify that the database was built correctly, I copied the contents of target/db over to /var/lib/neo4j/data/graph.db, the data directory of a Neo4j installation I had installed using apt-get install. Unfortunately, there is a version mismatch, so the database was unreadable. To get the correct version of Neo4j, I looked at the POM file of batch-import (neo4j.version was set to 1.9) and found a tarball download of the same version here. Installation consisted of exploding the tarball and starting the server with bin/neo4j start).

1
2
3
sujit@tsunami:opt$ sudo tar xvzf neo4j-community-1.9.6-unix.tar.gz
sujit@tsunami:opt$ cd neo4j-community-1.9.6
sujit@tsunami:neo4j-community-1.9.6$ bin/neo4j start

The server exposes a Web Admin client (similar to the Solr Admin client) at port 7474 (http://localhost:7474/webadmin/). The dashboard shows 2,880,385 nodes, 3,375,083 properties, 58,021,093 relationships and 653 relationship types, which matches with what we put in.

Thats all I have for today. Next week I hope to learn more about Neo4j, specifically its Cypher Query Language, and see if I can model some common use-cases using it.

Saturday, June 01, 2013

MapReduce with Python and mrjob on Amazon EMR


I've been doing the Introduction to Data Science course on Coursera, and one of the assignments involved writing and running some Pig scripts on Amazon Elastic Map Reduce (EMR). I've used EMR in the past, but have avoided it ever since I got burned pretty badly for leaving it on. Being required to use it was a good thing, since I got over the inertia and also saw how much nicer the user interface had become since I last saw it.

I was doing another (this time Python based) project for the same class, and figured it would be educational to figure out how to run Python code on EMR. From a quick search on the Internet, mrjob from Yelp appeared to be the one to use on EMR, so I wrote my code using mrjob.

The code reads an input file of sentences, and builds up trigram, bigram and unigram counts of the words in the sentences. It also normalizes the text, lowercasing, replacing numbers and stopwords with placeholder tokens, and Porter stemming the remaining words. Heres the code, as you can see, its fairly straightforward:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
from __future__ import division
from mrjob.job import MRJob
import nltk
import string

class NGramCountingJob(MRJob):

  def mapper_init(self):
#    self.stopwords = nltk.corpus.stopwords.words("english")
    self.stopwords = set(['i', 'me', 'my', 'myself', 'we',
      'our', 'ours', 'ourselves', 'you', 'your', 'yours',
      'yourself', 'yourselves', 'he', 'him', 'his', 'himself',
      'she', 'her', 'hers', 'herself', 'it', 'its', 'itself',
      'they', 'them', 'their', 'theirs', 'themselves', 'what',
      'which', 'who', 'whom', 'this', 'that', 'these', 'those',
      'am', 'is', 'are', 'was', 'were', 'be', 'been', 'being',
      'have', 'has', 'had', 'having', 'do', 'does', 'did',
      'doing', 'a', 'an', 'the', 'and', 'but', 'if', 'or',
      'because', 'as', 'until', 'while', 'of', 'at', 'by',
      'for', 'with', 'about', 'against', 'between', 'into',
      'through', 'during', 'before', 'after', 'above', 'below',
      'to', 'from', 'up', 'down', 'in', 'out', 'on', 'off',
      'over', 'under', 'again', 'further', 'then', 'once',
      'here', 'there', 'when', 'where', 'why', 'how', 'all',
      'any', 'both', 'each', 'few', 'more', 'most', 'other',
      'some', 'such', 'no', 'nor', 'not', 'only', 'own', 'same',
      'so', 'than', 'too', 'very', 's', 't', 'can', 'will',
      'just', 'don', 'should', 'now'])
    self.porter = nltk.PorterStemmer()

  def mapper(self, key, value):

    def normalize_numeric(x):
      xc = x.translate(string.maketrans("", ""), string.punctuation)
      return "_NNN_" if xc.isdigit() else x

    def normalize_stopword(x):
      return "_SSS_" if str(x) in self.stopwords else x

    cols = value.split("|")
    words = nltk.word_tokenize(cols[1])
    # normalize number and stopwords and stem remaining words
    words = [word.lower() for word in words]
    words = [normalize_numeric(word) for word in words]
    words = [normalize_stopword(word) for word in words]
    words = [self.porter.stem(word) for word in words]
    trigrams = nltk.trigrams(words)
    for trigram in trigrams:
      yield (trigram, 1)
      bigram = trigram[1:]
      yield (bigram, 1)
      unigram = bigram[1:]
      yield (unigram, 1)

  def reducer(self, key, values):
    yield (key, sum([value for value in values]))

if __name__ == "__main__":
  NGramCountingJob.run()

The class must extend MRJob and call its run() method when invoked from the shell. The MRJob class implements a sequence of methods that will be called (which can be overriden) - so we just override the appropriate methods. The mrjob framework runs over Hadoop Streaming, but offers many convenience features.

I first tried using the EMR console to create a job flow with mrjob, but the closest I found was "Streaming Jobs". Streaming jobs require the mapper and reducer scripts and the input files to reside on Amazon's S3 storage. Output is also written to S3. However, I was not able to make this setup work with the mrjob code above.

Reading some more, it turns out that mrjob jobs can be started from your local shell with the "-r emr" switch, and it will copy your input and scripts to S3, create a job flow, run your job, write output to S3, and then copy the output back to STDOUT of your local shell, where you can capture it. The first thing that is needed is an .mrjob.conf file. Mine looks like this (with the secret bits appropriately sanitized).

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# Source: $HOME/.mrjob.conf
runners:
  emr:
    aws_access_key_id: 53CR3T53CR3T53CR3T53
    aws_region: us-west-1
    aws_secret_access_key: SuperSecretAccessKeyIfITellYouGottaKillU
    ec2_key_pair: EMR
    ec2_key_pair_file: /path/to/pem/file.pem
    bootstrap_cmds:
    - sudo easy_install http://nltk.googlecode.com/files/nltk-2.0b5-py2.6.egg
    ec2_instance_type: m1.small
    num_ec2_core_instances: 4
    cmdenv:
      TZ: America/Los_Angeles
  local:
    base_tmp_dir: /tmp/$USER

The configuration will start up a EC2 m1.small master node with 4 slave nodes of the same type. The bootstrap_cmds installs NLTK on all the worker nodes, since my code is using it and because it doesn't come standard with Python installs. I also had a call to nltk.corpus to read English stopwords, but I just changed the code to declare the list explicitly since I didn't want to install the full corpus.

You can run the code locally (for testing, generally on a subset of the data) as follows.

1
sujit@cyclone:parw$ python ngram_counting_job.py input.txt > output.txt

Or you can run on EMR by adding a "-r emr" switch, or running on your own Hadoop cluster by adding a "-r hadoop" switch to your command. The EMR version is shown below.

1
sujit@cyclone:parw$ python ngram_counting_job.py input.txt -r emr > output.txt

Of course, you can monitor your job from the EMR console as it is running. This is all I've done with mrjob so far, but I hope to do much more with it.

Friday, December 28, 2012

Analyzing the Enron Data: Frequency Distribution, Page Rank and Document Clustering


I've been using the Enron Dataset for a couple of projects now, and I figured that it would be interesting to see if I could glean some information out of the data. One can of course simply read the Wikipedia article, but that would be too easy and not as much fun :-).

My focus on this analysis is on the "what" and the "who", ie, what are the important ideas in this corpus and who are the principal players. For that I did the following:

  • Extracted the words from Lucene's inverted index into (term, docID, freq) triples. Using this, I construct a frequency distribution of words in the corpus. Looking at the most frequent words gives us an idea of what is being discussed.
  • Extract the email (from, {to, cc, bcc}) pairs from MongoDB. Using this, I piggyback on Scalding's PageRank implementation to produce a list of emails by page rank. This gives us an idea of the "important" players.
  • Using the triples extracted from Lucene, construct tuples of (docID, termvector), then cluster the documents using KMeans. This gives us an idea of the spread of ideas in the corpus. Originally, the idea was to use Mahout for the clustering, but I ended up using Weka instead.

I also wanted to get more familiar with Scalding beyond the basic stuff I did before, so I used that where I would have used Hadoop previously. The rest of the code is in Scala as usual.

Unfortunately, my Scalding-fu was not strong enough, because I couldn't figure out how to run Scalding jobs in non-local mode, and I was running out of memory for some of the jobs in local mode. Also, Mahout expects its document vectors to be in SequenceFile format, but Scalding does not allow you to write SequenceFiles in local mode. As a result, I converted the document vector file to ARFF format and used Weka for the clustering instead. I have asked about this on the Cascading-Users mailing list.

Frequency Distribution


Here is the the code to extract the (term, docID, freq) tuples from the Lucene index. As you can see the generate method takes three cutoff parameters, the minimum document frequency, the minimum total term frequency and minimum term frequency. I added these cutoffs in because the unfiltered output contained about 700 million triples and was causing the next job (FreqDist) to throw an OutOfMemoryException.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
// Source: src/main/scala/com/mycompany/solr4extras/corpus/Lucene4TermFreq.scala
package com.mycompany.solr4extras.corpus

import java.io.{PrintWriter, FileWriter, File}

import org.apache.lucene.index.{MultiFields, IndexReader, DocsEnum}
import org.apache.lucene.search.DocIdSetIterator
import org.apache.lucene.store.NIOFSDirectory
import org.apache.lucene.util.BytesRef

/**
 * Reads a Lucene4 index (new API) and writes out a
 * text file as (term, docID, frequency_of_term_in_doc).
 * @param indexDir the location of the Lucene index.
 * @param outputFile the output file name.
 * @param minDocFreq terms which are present in fewer
 *        documents than minDocFreq will be ignored. 
 * @param minTTF the minimum Total Term Frequency a 
 *        term must have to be considered for inclusion.
 * @param minTermFreq the minimum term frequency within
 *        a document so the term is included.
 */
class Lucene4TermFreq(indexDir: String) {

  def generate(outputFile: String, minDocs: Int,
      minTTF: Int, minTermFreq: Int): Unit = {
    val reader = IndexReader.open(
      new NIOFSDirectory(new File(indexDir), null))
    val writer = new PrintWriter(new FileWriter(outputFile), true)
    val terms = MultiFields.getTerms(reader, "body").iterator(null)
    var term: BytesRef = null
    var docs: DocsEnum = null
    do {
      term = terms.next
      if (term != null) {
        val docFreq = terms.docFreq
        val ttf = terms.totalTermFreq
        if (docFreq > minDocs && ttf > minTTF) {
          docs = terms.docs(null, docs)
          var docID: Int = -1
          do {
            docID = docs.nextDoc
            if (docID != DocIdSetIterator.NO_MORE_DOCS) {
              val termFreq = docs.freq
              if (termFreq > minTermFreq)
                writer.println("%s\t%d\t%d".format(
                  term.utf8ToString, docID, docs.freq))
            }
          } while (docID != DocIdSetIterator.NO_MORE_DOCS)
        }
      }
    } while (term != null)
    writer.flush
    writer.close
    reader.close
  }
}

To decide the cutoffs, I first plotted the sorted total term frequencies of all words in the corpus, that results in the Zipf distribution shown below, with the last 20,000 terms contributing to most of the occurrence count. Cutting the TTFs off at about 1,000 occurrences, and setting the minimum document frequency to 5 (term must exist in at least 5 documents) and minimum term frequency to 10 (term must exist at least 10 times in a document to be counted) resulted in a more manageable number of about 1.27 million triples.


Here is the code to compute the Frequency Distribution of the terms in the corpus.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
// Source: src/main/scala/com/mycompany/solr4extras/corpus/FreqDist.scala
package com.mycompany.solr4extras.corpus

import com.twitter.scalding.{Tsv, Job, Args}

import cascading.pipe.joiner.LeftJoin

/**
 * Reads input of the form (term docID freq), removes stopword
 * terms based on a stop word list, sums up the term frequency
 * across docs and outputs the term frequency counts sorted by
 * count descending as (term count).
 * NOTE: this can also be done directly from Lucene using 
 * totalTermFreq.
 */
class FreqDist(args: Args) extends Job(args) {

  val stopwords = Tsv(args("stopwords"), ('stopword)).read
  val input = Tsv(args("input"), ('term, 'docID, 'freq))
  val output = Tsv(args("output"))
  input.read.
    joinWithSmaller('term -> 'stopword, stopwords, joiner = new LeftJoin).
    filter('stopword) { stopword: String => 
      (stopword == null || stopword.isEmpty) 
    }.
    groupBy('term) { _.sum('freq) }.
    groupAll { _.sortBy('freq).reverse }.
    write(output)
}

The top 100 words (and their raw frequencies) from the resulting frequency distribution are shown below:

enron (1349349), ect (1133513), hou (578328), subject (446732), pm (388238), http (325908), power (309063), cc (303380), enron.com (290452), energy (286475), corp (262331), message (245042), mail (234409), time (217656), gas (216405), company (189878), market (181420), information (180539), ees (176279), original (170590), call (153876), california (152559), business (148100), forwarded (145700), day (138167), na (132233), td (132135), font (132130), price (131493), week (131246), state (130245), year (127376), email (124309), attached (121637), houston (119962), image (118982), john (113327), meeting (111627), agreement (111621), mark (111518), deal (108930), make (106030), group (105643), trading (105130), questions (103417), enron_development (102359), contact (97935), date (96189), back (95138), million (93875), services (93825), work (92014), jeff (89695), today (89682), report (89632), electricity (88541), service (88489), monday (87089), prices (85148), free (83838), friday (83750), credit (83487), contract (82981), system (80955), financial (80927), good (80489), review (78720), fax (78219), management (76054), companies (75371), david (74984), news (74666), number (73917), file (73378), jones (72288), thursday (72182), order (71953), list (71845), send (71824), forward (71651), tuesday (71515), office (70850), october (70826), based (70518), enronxgate (69905), wednesday (69339), risk (69091), change (68911), received (68688), mike (68609), issues (68449), team (67302), bill (67289), click (66993), plan (66953), customers (66047), communications (65620), november (65478), phone (65434), provide (65347)

Page Rank


The next thing I wanted to find was to somehow rank the players (the email authors) in importance. In any organization, people tend to email at their own level, but include their immediate bosses on the CC or BCC as a form of CYA. So top management would get relatively few (but highly ranked) emails from middle management, and middle management would get many (low ranked) emails from their underlings. Sorting the email authors by descending order of page rank should tell us about he major players.

The input data comes from the MongoDB database I populated for my previous project. I had to go through a few hoops because the data is encrypted, but the net result is that I end up with two files, the first containing the (from_email, from_id) tuples and the second containing (from_id, {to_id, cc_id, bcc_id}) tuples. The reason I generated two files is because the Scalding distribution comes with a PageRank implementation which expects the input data to be numeric. Here is the extraction code:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
// Source: src/main/scala/com/mycompany/solr4extras/corpus/MongoEmailPairs.scala
package com.mycompany.solr4extras.corpus

import java.io.{PrintWriter, FileWriter, File}
import java.util.concurrent.atomic.AtomicInteger

import org.apache.commons.codec.binary.Hex
import org.apache.lucene.index.IndexReader
import org.apache.lucene.store.NIOFSDirectory

import com.mongodb.casbah.Imports.{wrapDBObj, wrapDBList, MongoDBObject, MongoConnection, BasicDBList}
import com.mycompany.solr4extras.secure.CryptUtils

class MongoEmailPairs(host: String, port: Int, db: String, 
    indexDir: String) {

  val conn = MongoConnection(host, port)
  val emails = conn(db)("emails")
  val users = conn(db)("users")
  val reader = IndexReader.open(
    new NIOFSDirectory(new File(indexDir), null))

  def generate(refFile: String, outputFile: String): Unit = {
    val counter = new AtomicInteger(0)
    val userKeys = users.find().map(user => 
      user.as[String]("email") -> 
      (Hex.decodeHex(user.as[String]("key").toCharArray), 
      Hex.decodeHex(user.as[String]("initvector").toCharArray),
      counter.incrementAndGet)).toMap
    // write out dictionary file for reference
    val refWriter = new PrintWriter(new FileWriter(new File(refFile)), true)
    userKeys.map(user =>
      refWriter.println("%s\t%d".format(user._1, user._2._3))
    )
    refWriter.flush
    refWriter.close
    // write out main file as required by PageRank
    val dataWriter = new PrintWriter(new FileWriter(new File(outputFile)), true)
    val numdocs = reader.numDocs
    var i = 0
    while (i < numdocs) {
      val doc = reader.document(i)
      val messageID = doc.get("message_id").asInstanceOf[String]
      val author = doc.get("from").asInstanceOf[String]
      val mongoQuery = MongoDBObject("message_id" -> messageID)
      val cur = emails.find(mongoQuery)
      emails.findOne(mongoQuery) match {
        case Some(email) => {
          try {
            val from = CryptUtils.decrypt(
              Hex.decodeHex(email.as[String]("from").toCharArray), 
              userKeys(author)._1, userKeys(author)._2)
            val fromId = userKeys(from)._3
            val targets = 
              (try {
                email.as[BasicDBList]("to").toList  
              } catch {
                case e: NoSuchElementException => List()
              }) ++
              (try {
                email.as[BasicDBList]("cc").toList
              } catch {
                case e: NoSuchElementException => List()
              }) ++
              (try {
                email.as[BasicDBList]("bcc").toList
              } catch {
                case e: NoSuchElementException => List()
              })
            targets.map(target => {
              val targetEmail = CryptUtils.decrypt(Hex.decodeHex(
                target.asInstanceOf[String].toCharArray), 
                userKeys(author)._1, userKeys(author)._2).trim
              val targetEmailId = userKeys(targetEmail)._3
              dataWriter.println("%d\t%d".format(fromId, targetEmailId))
            })
          } catch {
            // TODO: BadPaddingException, likely caused by, 
            // problems during population. Fix, but skip for now
            case e: Exception => println("error, skipping")
          }
        }
        case None => // skip
      }
      i = i + 1
    }
    dataWriter.flush
    dataWriter.close
    reader.close
    conn.close
  }
}

We then create our own subclass of PageRank and override the initialize method that produces a Source tap for PageRank from the email pairs we just generated from MongoDB. The work that we do here is to group by the from_id and aggregate the to_ids into a comma-separated list, then add a column with the initial value of the page rank (1.0), and rename the columns so it is usable by the parent class. Once the job has finished, we post-process the output to produce a list of from email addresses and their associated page rank in descending rank order. Here is the code for these two classes:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
// Source: src/main/scala/com/mycompany/solr4extras/corpus/MailRank.scala
package com.mycompany.solr4extras.corpus

import com.twitter.scalding.examples.PageRank
import com.twitter.scalding.{Tsv, Job, Args}

import cascading.pipe.joiner.LeftJoin
import cascading.pipe.Pipe

/**
 * Converts data generated by MongoEmailPairs of the 
 * form: (from_id, {to_id|cc_id|bcc_id}) to the format
 * required by PageRank, ie (from_id, List(to_ids), pagerank)
 */
class MailRank(args: Args) extends PageRank(args) {

  override def initialize(nodeCol: Symbol, neighCol: Symbol, 
      pageRank: Symbol): Pipe = {
    val input = Tsv(args("input"), ('from, 'to))
    input.read.
      groupBy('from) { _.toList[String]('to -> 'tos) }.
      map('tos -> 'tosf) { tos: List[String] => 
        tos.foldLeft("")(_ + "," + _).substring(1) }.
      map('from -> ('from, 'prob)) { from: String => 
        (from, 1.0) 
      }.project('from, 'tosf, 'prob).
      mapTo((0, 1, 2) -> (nodeCol, neighCol, pageRank)) {
        input : (Long, String, Double) => input
      }
  }
}

/**
 * Converts the format returned by PageRank, ie:
 * (from_id, List(to_id), final_pagerank) to 
 * (from_email, final_pagerank) sorted by pagerank 
 * descending.
 */
class MailRankPostProcessor(args: Args) extends Job(args) {
  
  val input = Tsv(args("input"), ('from, 'tos, 'rank))
  val output = Tsv(args("output"))

  val reference = Tsv(args("reference"), ('ref_email, 'ref_from)).read
  input.read.
    project('from, 'rank).
    joinWithSmaller('from -> 'ref_from, reference, joiner = new LeftJoin).
    project('ref_email, 'rank).
    groupAll { _.sortBy('rank).reverse }.
    write(output)
}

Here is the list of the top 15 email addresses which had the highest page rank. If you followed the Enron trial, then you may recognize a few names here:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
vince.kaminski@enron.com    129.6456344199462
sara.shackleton@enron.com   123.65630393847148
louise.kitchen@enron.com    117.4938083777274
jeff.dasovich@enron.com     116.83443565190146
tana.jones@enron.com        113.12379367608482
mark.taylor@enron.com       109.67593395984632
sally.beck@enron.com        108.4905670421611
ebass@enron.com             103.54919812912469
jeff.skilling@enron.com     100.57545577565519
steven.kean@enron.com        99.69270847011283
john.lavorato@enron.com      90.16447940199485
gerald.nemec@enron.com       88.4001157213643
kenneth.lay@enron.com        88.1467699737448
richard.shapiro@enron.com    82.10524578705625
kay.mann@enron.com           69.19222780384432

Document Clustering


Finally, I decided to cluster the documents in order to find if there were multiple topics that were being discussed in the corpus. The input to this process is the (term, docID, freq) triples that we generated from Lucene, and the output is a set of (docID, {term frequency vector}) tuples. The sets up partial vectors and aggregates them together to form document vectors using Mahout's SequentialAccessSparseVector class. I found post on Software Anatomy very useful when writing this code - much of the DocVector class below is based on the code shown there.

1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
// Source: src/main/scala/com/mycompany/solr4extras/corpus/DocVector.scala
package com.mycompany.solr4extras.corpus

import java.io.{PrintWriter, FileWriter, File}

import scala.collection.mutable.ListBuffer
import scala.io.Source

import org.apache.mahout.math.{VectorWritable, SequentialAccessSparseVector}

import com.twitter.scalding.{Tsv, TextLine, Job, Args}

import cascading.pipe.joiner.LeftJoin

/**
 * Reads input of the form (term freq) of most frequent terms,
 * and builds a dictionary file. Using this file, creates a 
 * collection of docID to sparse doc vector mappings of the 
 * form (docID, {termID:freq,...}).
 */
class DocVector(args: Args) extends Job(args) {

  val input = Tsv(args("input"), ('term, 'docID, 'freq))
  val termcounts = TextLine(args("termcounts"))
  
  val dictOutput = Tsv(args("dictionary"))
  val output = Tsv(args("docvector"))
  
  // (term freq) => (term num)
  val dictionary = termcounts.read.
    project('num, 'line).
    map('line -> 'word) { line: String => line.split('\t')(0) }.
    project('word, 'num)

  // input: (term, docID, freq)
  // join with dictionary ond write document as (docId, docvector) 
  input.read.
    joinWithSmaller('term -> 'word, dictionary, joiner = new LeftJoin).
    filter('word) { word: String => (!(word == null || word.isEmpty)) }.
    project('docID, 'num, 'freq).
    map(('docID, 'num, 'freq) -> ('docId, 'pvec)) { 
      doc: (String, Int, Int) =>
        val pvec = new SequentialAccessSparseVector(
          args("vocabsize").toInt)
        pvec.set(doc._2, doc._3)
      (doc._1, new VectorWritable(pvec))
    }.
    groupBy('docId) { 
      group => group.reduce('pvec -> 'vec) {
        (left: VectorWritable, right: VectorWritable) => 
          new VectorWritable(left.get.plus(right.get).normalize)
    }}.
    write(output)
    
    // save the dictionary as (term, idx)    
    dictionary.write(dictOutput)
}

/**
 * Converts the Document Vector file to an ARFF file for 
 * consumption by Weka.
 */
class DocVectorToArff {
  
  def generate(input: String, output: String, 
      numDimensions: Int): Unit = {
    val writer = new PrintWriter(new FileWriter(new File(output)), true)
    // header
    writer.println("@relation docvector\n")
    (1 to numDimensions).map(n => 
      writer.println("@attribute vec" + n + " numeric"))
    writer.println("\n@data\n")
    // body
    Source.fromFile(new File(input)).getLines.foreach(line => { 
      writer.println(line.split('\t')(1).replaceAll(":", " "))
    })
    writer.flush
    writer.close
  }
}

/**
 * Reads output from Weka Explorer SimpleKMeans run (slightly
 * modified to remove header information) to produce a list
 * of top N words from each cluster.
 */
class WekaClusterDumper {
  
  def dump(input: String, dictionary: String, 
      output: String, topN: Int): Unit = {
    
    // build up map of terms from dictionary
    val dict = Source.fromFile(new File(dictionary)).getLines.
      map(line => { 
        val cols = line.split("\t")
        cols(1).toInt -> cols(0)
    }).toMap
    // build up elements list from weka output
    var clusterScores = new Array[ListBuffer[(Int,Double)]](5)
    Source.fromFile(new File(input)).getLines.
      foreach(line => {
        val cols = line.split("\\s+")
        val idx = cols(0).toInt - 1
        val scores = cols.slice(2, 7) 
        (0 to 4).foreach(i => 
          if (scores(i).toDouble > 0.0D) {
            if (clusterScores(i) == null)
              clusterScores(i) = new ListBuffer[(Int,Double)]
            clusterScores(i) += Tuple(idx, scores(i).toDouble)
        })
    })
    // sort each clusterScore by score descending and get the
    // corresponding words from the dictionary by idx
    val writer = new PrintWriter(new FileWriter(new File(output)), true)
    var i = 0
    clusterScores.foreach(clusterScore => {
      writer.println("Cluster #" + i)
      clusterScore.toList.sort(_._2 > _._2).
        slice(0, topN).map(tuple => {
          val word = dict(tuple._1)
          writer.println("  " + word + " (" + tuple._2 + ")")
        })
      i = i + 1
    })
    writer.flush
    writer.close
  }
}

As I mentioned above, the output of DocVector was supposed to be passed into Mahout's KMeans driver (and Canopy driver for the initial centroids) but I could not generate Sequence files which Mahout expects, so I converted the DocVector output to an ARFF file so I could pass it into Weka.

I then ran SimpleKMeans requesting 5 clusters with 10 iterations, and the first attempt resulted in an OutOfMemoryException. Since I have a pretty bad-ass MacBook Pro with 8GB RAM (okay, it was bad-ass when I bought it 3 years ago), I upped the default -Xmx for Weka from 256MB to 2GB (this is in the Info.plist file under the /Applications/weka directory for those of you who use Macs), and life was good again.

In any case, after a while, Weka dumped out the results of the computation to the Explorer console, from which I copy-pasted it and passed it through my WekaClusterDumper class. This class parses the Weka clustering output (minus the headers which I removed manually) to print the top N words in each cluster. Since Weka does not know the actual terms being clustered (only their position in the dictionary file generated by DocVector), the Cluster dumper uses this file to look up the actual terms in each cluster. Here is the output of the dump.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
Cluster #0
  ect (13.3978)
  image (4.639)
  enron.com (3.7695)
Cluster #1
  enron (20.7614)
Cluster #2
  http (0.3507)
  enron.com (0.309)
  enron_development (0.2871)
  mail (0.1964)
  ees (0.1615)
  image (0.1558)
  hou (0.1097)
  message (0.1039)
  ect (0.0909)
  energy (0.0892)
Cluster #3
  mail (56.5979)
Cluster #4
  study (12.2941)

As you can see, according the clustering, there does not seem to be too much variety in the discussion going on in the Enron dataset.

Finally, I used a single object to call the classes described above from sbt (using "sbt run"). The whole process was very interactive, so don't run the code as is, it will fail. I ran each step individually, some multiple times, and commented out the previous blocks as I went forward.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
// Source: src/main/scala/com/mycompany/solr4extras/corpus/CorpusAnalyzer.scala
package com.mycompany.solr4extras.corpus

import com.twitter.scalding.Args

object CorpusAnalyzer extends App {

  //////////// document clustering ///////////
  
  (new Lucene4TermFreq("/Users/sujit/Downloads/apache-solr-4.0.0/example/solr/collection1/data/index")).
    generate("data/input/corpus_freqs.txt", 50, 1000, 10)

  (new FreqDist(Args(List(
    "--local", "", 
    "--input", "data/input/corpus_freqs.txt", 
    "--output", "data/output/freq_dist.txt",
    "--stopwords", "/Users/sujit/Downloads/apache-solr-4.0.0/example/solr/collection1/conf/stopwords.txt")))).
    run

  //////// mail rank ////////
  
  (new MongoEmailPairs("localhost", 27017, "solr4secure", 
    "/Users/sujit/Downloads/apache-solr-4.0.0/example/solr/collection1/data/index")).
    generate("data/input/email_refs.txt",
    "data/input/email_pairs.txt")
  
  (new MailRank(Args(List(
    "--local", "",
    "--input", "data/input/email_pairs.txt",
    "--output", "data/output/mailrank.txt",
    "--iterations", "10")))).
    run

  (new MailRankPostProcessor(Args(List(
    "--local", "",
    "--input", "data/output/mailrank.txt",
    "--reference", "data/input/email_refs.txt",
    "--output", "data/output/mailrank_final.txt")))).run

  ////////////// clustering terms ///////////////
    
  (new DocVector(Args(List(
    "--local", "",
    "--input", "data/input/corpus_freqs.txt",
    "--termcounts", "data/input/freq_words.txt",
    "--vocabsize", "1326", // cat freq_words | cut -f1 | sort | uniq | wc
    "--dictionary", "data/output/dictionary.txt",
    "--docvector", "data/output/docvector.txt")))).
    run

  (new DocVectorToArff()).generate(
    "/Users/sujit/Projects/solr4-extras/data/output/docvector.txt", 
    "data/output/docvector.arff", 1326)

  (new WekaClusterDumper()).dump(
    "/Users/sujit/Projects/solr4-extras/data/output/weka_cluster_output.txt",
    "/Users/sujit/Projects/solr4-extras/data/output/dictionary.txt",
    "data/output/cluster.dump",10)
    
}

Well, this is all I have for today. The source code for this stuff is also available on my GitHub project page if you want to play around with it on your own. Hope you enjoyed it. If I don't get a chance to post again before next year (unlikely given that there are just 3 more days left and my mean time between posts is now about 14 days), be safe and have a very Happy New Year. Heres hoping for a lot of fun in 2013 (the International Year of Statistics).