Winter Break Project: Building Minesweeper from Scratch

I don’t take a ton of vacation from my work, but luckily I don’t really need to over the winter: my employer gives us christmas through new years off, which comes out to around 10 days vacation. I had planned to do some travelling, but instead tagged along with my girlfriend to do some family stuff for new years. Since I was in town for most of the break, I spent most of my time just chilling (which, frankly, was awesome) and coding.

I’ve been programming for a few years now, both for fun and for work, but my programs usually are either scripts or backend (database) programs. I never got around to learning GUI programming, so I decided to use this opportunity to build a GUI program. I had recommended some minesweeper related challenges to a friend who’s learning programming, so i had minesweeper on the brain. I ended up coding the main game out of boredom one day, and decided to slap a GUI on top of it. 

Whenever I’ve seen GUI tutorials online, for simplicity sake they combine the main program code with the GUI code. I guess this makes sense for a quick tutorial, but my understanding is that this is pretty weak GUI programming. Generally, you want to separate the “model”, “controller” and “view” components.

For my application, my model was comprised of a board class which is a container for tile class objects I also defined. The controller was a “game” class which included functions like “click” and “flip_flag.” I wanted the game component (the model and controller) to operate just fine in isolation from the GUI, so my Game class also has a play() method for playing minesweeper in the terminal. Having built the working game components, I could now slap a GUI on top in a separate “view” component.

The most common GUI toolkits for python are wx, qt and tk. Tk is built into python as Tkinter, has a repuation for being simple to use, has plenty of documentation, and has the added benefit of taking advantage of the native OS’s aesthetics, so I decided to go with Tkinter.

Now, although there are lots of Tkinter tutorials, there doesn’t seem to be a ton of consistency in how people use Tkinter. In general, the main GUI is defined as a class called “App.” But sometimes the App class inherits from the Tk class, and sometimes it inherits from the Frame class. I started using an example that had my App inheriting from Frame, and this worked fine. I later wanted to implement a feature where the buttons would be hidden when the game was paused, and I found that this required a lot of changes if I kept it inheriting from the Frame class. I changed my app to inherit from the Tk class and it was a lot easier to add that feature. I think that inheritance makes more sense intuitively as well.

Anyway, as an exercise, I strongly recommend programming minesweeper. Next time I need to learn a new language, I plan to use programming minesweeper as a comprehensive exercise. To accomplish the task takes using several primitive types and functions, defining classes, using threads (for the timer), and of course the GUI.

 

 

Downloading images from saved reddit links

Reddit user aesptux posted to /r/learnprogramming requesting a code review of their script to download images from saved reddit links. This user made the same mistake I originally made and tried to work with the raw reddit JSON data. I hacked together a little script to show them how much more quickly they could accomplish their task using praw. It took almost no time at all, and the meat of the code is only about 20 lines. What a great API.

Here’s my code: https://gist.github.com/4225456

Idea: graphing wikipedia

I’m almost certain something like this has been done before. Anyway, here’s the idea:

  1. Download the wikipedia database dump.
  2. Ingest article texts into a database
  3. Scrape wikipedia links out of the first paragraph of each article.
  4. Create a directed graph of articles where two articles share an edge if they are linked as described in (3). Treat article categories as node attributes.
  5. Investigate community structure of wikipedia articles, particularly which categories cluster together
  6. Extra challenge: Try to find articles that won’t “get you to philosophy”

There are currently over 4M articles in the english wikipedia, so for this to be feasible I will probably need to invent some criterion for including articles in the project, probably minimum length, minimum age, or minimum edits. Alternatively, I might just focus on certain categories/subcategories.

Weekend Project: Word Ladder Solver

<< I plan to add a few graphics, code snippets, and trim the code posted at the bottom, but I just haven’t gotten around to it. I wrote the bulk of this about two weeks ago now and just haven’t gotten around to finishing it. I’ll update this post when I can, but for now, I’m just going to publish since I have no time to clean it up. I’ll remove this message when I’ve made the post all nice and pretty (probably never).>>

The other day I learned about a game Lewis Carol invented where you make “word ladders” or “word bridges.” The idea is you pick two words, and try to get from one to the other by using a chain of valid words where each word in the chain differs from the previous word by a single letter. Here’s an example:

WORD
CORD
CORED
CORDED
CODDED
RODDED
RIDDED
RIDGED
RIDGE
BRIDGE

I immediately thought making a word bridge generator would be a fun little challenge to program, and it turned out to be a slightly more difficult problem than I’d anticipated. Let’s start by formalizing this game and we’ll work our way up to what specifically makes this challenging.

Word Similarity

For our purposes, two words are similar if you can get from one word to another by “changing a letter.” This change can be an insertion (CORD->CORED), a deletion (RIDGED->RIDGE), or a substitution (WORD->CORD). This is a pretty useful similarity metric and is something of a go-to tool for a lot of natural language processing tasks: it’s called “edit distance” or “levenshtein distance” where the “distance” is the minimum number of edits to transform one word into the other, so if edit_distance(w1, w2) = 0, then w1 and w2 are the identical. Note, it’s the MINIMUM number of edits. Edit distance calculation does not need to be edits-via-valid-words like I demonstrated above, so although I was able to transform WORD into BRIDGE with 9 edits, the edit distance between these two words is actually 5:

WORD
WRD — del
BRD — sub
BRID — ins
BRIDG — ins
BRIDGE — ins

The edit distance algorithm is a pretty neat little implementation of dynamic programming, but luckily I don’t need to build it up from scratch: python’s natural language toolkit has it built-in (be careful, it’s case sensitive):

from nltk import edit_distance
edit_distance('WORD','BRIDGE') # 5
edit_distance('WORD','bridge') # 6

GRAPH TRAVERSAL

In the context of this problem, we can think of our word list (i.e. the English dictionary) as an undirected graph where the nodes are words and similar words are connected. Using this formulation of the problem, a word bridge is a path between two nodes. There are several path finding algorithms available to us, but I’m only going to talk about two of the more basic ones: depth-first search (DFS) and breadth-first search (BFS).

DFS essentially goes as far along a particular branch as it can from a start node until it reaches a “terminal node” (can’t go any further) or a “cycle” (found itself back where it had already explored) and then it backtracks, doing the same thing along all the unexplored paths. If you think of a tree climber trying to map out a tree, it would be like them climbing up the trunk until they reached a fork, and then following that branch in the same direction all the way out until it reached a leaf, then working backwards to the last fork it passed, and then going forwards again until it hit a leaf, and so on.

Depth First Search

On the other hand, BFS looks at the graph in terms of layers. BFS first gets all the child nodes of the starting node, so now the furthest nodes are 1 edge from the start. Then it does the same thing for those children, so the furthest nodes in our searched subgraph are 2 edges away, and so on. A good visualization for this projection is a really contagious virus spreading through a population. First, people in immediate contact with patient zero get sick. Then everyone those people are in immediate contact with get sick, and so on. Replace “person” with “node” and “sick” with searched, and that’s BFS.

Breadth First Search

If we use DFS to traverse this graph, will we find the word bridge (assuming one exists, since a path does not exist between all nodes in our graph) eventually, but a) there’s no guarantee it will be the shortest path, and b) how long it takes really depends on the order of nodes in our graph, which is sort of silly. BFS will necessarily find the shortest path for us: given that we’ve searched to a depth of N edges and not found a path, we know there cannot be word bridge of length N from our start word to our target word. Therefore, if we reach our target word on the N+1 iteration of BFS, we know that the shortest path is N+1 edits. It’s possible that there are multiple “shortest paths,” but we’ll stop when we reach the first one we see.

THE BOTTLENECK

Here’s a not-so-obvious question: what part of this problem is going to cause us the biggest problems? Edit distance? The graph search? The real problem here is scale, the sheer size of the word list. Once we have everything set up, traversing the graph will be fairly easy, but building up the network will take some time.

To find a path between two arbitrary words, we need a reasonably robust dictionary. The one I had on hand for this project was 83,667 words. We need to convert this dictionary into a network. Conceptually, we do this with a similarity matrix.  We take all the words in the dictionary and put all of them on each axis. For each cell in the matrix, we populate the value with the edit distance between the two words that correspond with that cell. Since we’re only concerned with entries that had edit distance “1,” we can set every other value in the matrix to zero. This reduces our similarity matrix to an adjacency matrix, describing which nodes in the graph share an edge.

Adjacency Matrix of an undirected graph. (Image via bytehood.com)

Obviously, implementing this in an array structure would be wildly inefficient. Such a matrix would need to be 83,667 x 83,667, or 7,000,166,889 cells, and presumably most of them will be 0’s. If we allocate a byte to each cell in the array, this would take up 7 GB in memory. In actuality, elements of a list are allocated 4 bytes each with additional 36 bytes allocated for the list itself.

>>>import sys
>>>sys.getsizeof([])
36
>>>sys.getsizeof([1])
40
>>>sys.getsizeof([1,2])
44

Therefore, if our array is implemented as a list of lists (there are other options, such as a pandas dataframe or a numpy array, but this is the naive solution), the array would take up:

36 * (83,667+1) + 4 * 83,667^2 = 28,003,679,604 bytes = 28 GB

I certainly don’t have 28 GB of RAM. We’re going to need a “sparse” data structure instead, one where we get to keep the ones but can ignore the zeros. The simple solution here in python is a dictionary (which is basically just a hash table). For each word in our word list, find all similar words. In our python dictionary, the key is the searched word and the collection of matches will be our value (as a list, set, tuple…whatever).

# The adjacency matrix for the graph above converted into a sparser
# dictionary representation
## MATRIX FORM: 6 X 6 = 36 Elements
## SPARSE FORM: 18 elements
{'A':['B','C','D','E']
,'B':['A','D','E']
,'C':['A','F']
,'D':['A','B','F']
,'E':['A','B','F']
,'F':['C','D','E']
}

HEURISTICS

We’re not quite there yet. For each word, we need to look over the entire dictionary. This means our code scales with n^2 (where n is the number of words in our wordlist), which sucks. Even without that knowledge, you can try the code out with this naive search and you’ll see, it will search like 20-100 words a minute. We have 83,667 words of varying size to get through before we’ve processed the graph we need for this problem, so we’ve got to come up with some tricks to make this faster.

One thing we can do is try to limit the number of words we run edit_distance() on. This function is reasonably fast, but it’s not as fast as a letter-by-letter comparison because it needs to build a matrix and do a bunch of calculations. For two words to have edit distance = 1, they can only differ by one letter. So we can compare the two letter prefix of our search word against every word in the dictionary and skip all the words whose prefixes differ by two letters. BAM! Way faster. But on my computer, this code will still take way too long to run. So let’s do the sae thing looking at suffixes. Even faster! I can’t remember, but I think at this point the estimated total run time for my code was 4 hours. Tolerable, but I’m impatient.

The problem is that although we no longer need to look at every letter of every word during a single pass, we still need to look at every word in the dictionary. But we can use the prefix/suffix heurstic we just came up with to significantly reduce the size of the seach space during any given pass though.

DIVIDE AND CONQUER

Bubble sort of an unordered list

The classic divide and conquer algorithm is binary search. Consider a sorted list of numbers and some random number within the range of this list. You want to figure out where this number goes. The naive solution is called bubble sort: start at the top of the list, compare your number against that number, if yours is bigger, move to the next item. Do this until you find your number’s spot in the order. At worst, you have to look through every number in the list. For any given problem, bubble sort is almost always the wrong algorithm.

Binary Search

The idea behind binary search is because we know the list is already ordered, we can throw away half the list each time. Let’s say the list you’re given is 100 numbers. Jump to the middle and compare. Say your number is smaller: we don’t need to consider the top 50 numbers at all now. Now jump to the 25th number. Bigger or smaller? Throw away the other 12. And so on. In a list of 100 numbers, binary search will find the position of your random number in at most 7 steps (log_2(100)). That’s quite an improvement from bubble sort’s 100 steps.

OVERLAPPING SUBPROBLEMS

Binary search isn’t really applicable to our problem, but we can still discount a signficant amount of the search space if we appropriately pre-process our word list. As we discovered above, we can divide the word list into 3 groups relative to our target word: both two-letter affixes match, prefix matches and/or suffix matches, or neither affix matches. There are a limited number of prefixes and suffixes in our word list, so if we’re searching a word that has the same affix we’ve already searched, we’re unnecessarily repeating work. This is what’s known as an overlapping subproblem, and wherever you can identify one you can make your code faster. Oftentimes a lot faster.

We’re going to handle this affix problem by grouping the words by prefix and suffix. In other words, we’re going to build two search indexes: an index on prefix and an index on suffix. I did this using a dictionary where the affix is the key and the matching words are the values.

def build_index(words, n_letters=2, left=True, verbose=True):
    n=n_letters
    res={}
    for w in words:
        if left:
            ix = w[:n]
        else:
            ix = w[(-1)*n:]
        d = res.get(ix, [])
        d.append(w)
        res[ix] = d
    return res

>>> build_index(words)

{'AA':
['AA',
 'AAH',
 'AAHED',
 'AAHING',
 'AAHS',
 'AAL',
 'AALII',
 'AALIIS',
 'AALS',
 'AARDVARK',
 'AARDWOLF',
 'AARGH',
 'AARRGH',
 'AARRGHH',
 'AAS',
 'AASVOGEL']
,'AB':
['AB',
 'ABA',
 'ABACA',
 'ABACAS',
 'ABACI',
 'ABACK',
 'ABACUS',
 'ABACUSES',
...
 'ZYMOSES',
 'ZYMOSIS',
 'ZYMOTIC',
 'ZYMURGY',
 'ZYZZYVA',
 'ZYZZYVAS']
,'ZZ':
['ZZZ']}

This significantly speeds up our code, but if you watch the code run, it seems to speed up and slow down in bursts. This isn’t because your RAM is wonky. This is because the elements in our indexes aren’t evenly distributed. Specifically, two-letter prefixes are pretty evenly distributed, but two letter-suffixes are not.

 To remedy this we’re going to use a three-letter index. An unfortunate result of my index implementation is that words that are smaller than the index size (here two letter words) don’t get considered when looking up words of size greater than or equal to the index. So although “AH” and “AAH” are similar, they’re contained in separate indexes with no overlap, so when I search for words similar to “AA” I’ll get “AT” but not “AAH”, and similarly when I look for words similar to “AAH” I’ll get “AAHS” but not “AA” or “AH.” This isn’t pefect, but I stronglly suspect the effect on the final product is negligible. I could just use a separate index size for the two prefixes, but this was simpler to code and it only just occured to me I could use a two-letter index on the left and a three-letter index on the right. So sue me. I regret nothing. The three letter index is faster anyway.

## Most Common suffixes: two letters
# ES     6810
# ED     6792
# ER     5296
# NG     4222
# RS     3480
# TS     2332
## Most Common suffixes:three letters
# ING    4062
# ERS    2743
# IER    1142
# EST    1035
# TED    1014
# IES     985
#   NB: The top two three-letter suffixes combined comprise about
#   the same amount of the wordlist as the top two most popular
#   two-letter suffixes separately, and the rest of the suffixes
#   are pretty uniformly distributed.

If you look at the X axis of the histograms below for the three letter indices you’ll notice that two letter words are segregated into their own buckets. If we made similar histograms for indices of four or five letters, these buckets would manifest as an unusual density on the far left of the graph.

With our indexes in place, we can now generate our similarity network in a reasonable amount of time: about 90 minutes on my computer.

def build_wordnet(wordlist, max_dist=1, index_size=None):
    """
    Returns a word similarity network (as a dict) where similarity is
    measured as edit_distance(word1,word2) <= max_dist.
    This implementation doesn't properly account for words smaller than
    the index size.

    index_size defaults to max_dist+1, but you should check the
    distribution  of prefixes/affixes in your wordlist to judge.
    In an 80K dictionary, I found a 3 letter index to be suitable,
    although this did result in the isolation of two letter words from
    the major component(s) of the network.
    """
    if index_size is None:
        index_size = max_dist + 1
    if verbose:
         print "Building right index..."
    R_ix = build_index(wordlist, n_letters=index_size, left=False)
    if verbose:
         print "Building left index..."
    L_ix = build_index(wordlist, n_letters=index_size, left=True)

    def check_index(j, affix, index, left=True):
        right = not left
        for k in index[affix]:
            if abs(len(j)-len(k))>max_dist:
                continue
            if right and edit_distance(j[-1*index_size:], k[-1*index_size:]) > max_dist:
                continue
            if left and edit_distance(j[:index_size], k[:index_size]) > max_dist:
                continue
            if j == k:
                continue
            if edit_distance(j,k) <= max_dist:
                similarity_network[j].update([k])

    similarity_network = {}
    n=0
    start = time.time()
    now = 0
    last = 0
    for j in wordlist:
        similarity_network[j]=set()
        prefix = j[:index_size]
        suffix = j[-1*index_size:]
        check_index(j,suffix, R_ix, left=False)
        check_index(j,prefix, L_ix, left=True)

        if verbose:
            n+=1
            now = int(time.time() - start)/60
            if now%5 == 0 and now == last+1:
                print n, int(time.time() - start)/60
            last = now
    return similarity_network

At this point you should really think about saving your result so you don’t have to wait 90 minutes every time you want to play with this tool we’re building (i.e. we should only ever have to run these calculations once). I recommend pickling the dictionary to a file.

import cPickle

datafile = 'similarity_matrix.dat'

f=open(datafile,'wb')
p = cPickle.Pickler(f)
g = build_wordnet(words)
p.dump(g)
f.close()

NETWORKS

Now that we’ve got the structure in place, we’re ready to find some word bridges! When I first started this project, I tweaked a recipe I found online for my BFS function. It gave me results, but it was messy so I’ll spare you my code. Because we’re working with a graph, we can take advantage of special graph libraries. There are several options in python, one of which is igraph which I’ve played with a little in R and it’s fine. After some googling around, I got the impression that the preference for python is generally the networkx library, so I decided to go with that.

All we have left to do is convert our graph into a networkx.Graph object, and call the networkx path search algorithm.

import networkx as nx
net = nx.Graph()
net.add_nodes_from(words)
for k,v in g.iteritems():
    for v_i in v:
        net.add_edge(k,vi)

w1 = raw_input("Start Word: ")
w2 = raw_input("End Word:   ")
nx.shortest_path(g, w1, w2)   # BAM! Don't need to reinvent the motherfucking wheel.

Really, we could have just used this data structure from the start when we were building the network up, but the dictionary worked just fine, the conversion step is still fast, I strongly suspect the native dict pickles to a smaller file than an nx.Graph would (feel free to prove it to yourself), and it was also a better pedagogical tool. Moreover, networkx graphs are dictionary-like, so the idea is still the same.

MAKING IT BETTER: OOP

A feature we’re lacking here is the ability to add missing words. After playing with this a little I’ve found that the word list I started with is short a few words, and doesn’t contain any proper nouns (preventing me from using my own name in word bridges, for example). In the current implementation, to add a single word I would need to rebuild the whole network. I could build in a new function that adds a word to the graph, but I’d want to add this word into indexes too. This is getting complicated. But if I rebuild my tools using OOP principles, this feature will be pretty simple to add and I can reuse most of my existing code (by wrapping several things that currently appear in loops as methods). There are two kinds of entities we’re going to want to represent: indexes, and word networks.


class Index(object):
    '''
    Container class for wordbridge index, which gets stored as a dict
    in the index attribute. affix
    '''
    def __init__(self, index_size=2, left=True):
        self.index_size = index_size
        self.left = left
        self.index = {}
        self._check_if_word_in_index = False
    def in_index(self, word):
        try:
            for w in self.index[self.get_affix(word)]:
                if w==word:
                    return True
        except:
            pass
        return False
    def get_affix(self, word):
        n = self.index_size
        if self.left:
            ix = word[:n]
        else:
            ix = word[(-1)*n:]
        return ix
    def add_word(self, word):
        if self._check_if_word_in_index:
            if in_index(word):
                return
        ix = self.get_affix(word)
        d = self.index.get(ix, [])
        d.append(word)
        self.index[ix] = d
        #self.words.append(word)
    def add_words(self, words):
        if type(words) in [type(()),type([]), type(set())]:
            for w in words:
                self.add_word(w)
        elif type(words) == type(''):
            self.add_word(words)
        else:
            # This should really be an error
            return "Invalid input type"

class Wordbridge(object):
    def __init__(self, words=None, wordslist_filepath=None):
        self.g = nx.Graph()
        self.indexes = []
        self.words = set()
        if type(wordslist_filepath) == type(''):
            self.seed_words_from_file(wordslist_filepath)
        if words:
            self.seed_words(words)
    def graph_from_dict(self, net_dict):
        for k,v in net_dict.iteritems():
            for v_i in v:
                self.g.add_edge(k,v_i)
        self.words.update(net_dict.keys())
    def update_indexes(self, word):
        for ix in self.indexes:
            ix.add_word(word)
    def add_word(self, word):
        w = word.upper()
        try:
            _ = self.g.node(w) # check if node is in graph
        except:
            self.update_indexes(w)
            self.g.add_node(w)
            self.check_indexes(w)
    def add_index(self, index_size, left=True):
        '''
        Creates a new index, stored in Wordnet.indexes list"
        '''
        ix = Index(index_size, left)
        ix.add_words(self.words)
        self.indexes.append(ix)
    def seed_words(self, words):
        self.words.update(words)
        self.g.add_nodes_from(self.words)
    def seed_words_from_file(self,path):
        with open(path) as f:
            word_list = f.read().split()
        self.seed_words(word_list)
    def check_indexes(self, j, max_dist=1):
        '''
        Adds a word to the graph
        '''
        candidates = set()
        n = max_dist+2
        for ix in self.indexes:
            right = not ix.left
            left = ix.left
            affix = ix.get_affix(j)
            candidates.update(ix.index[affix])
        for k in candidates:
            if abs(len(j)-len(k))>max_dist:
                continue
            # this is sort of redundant because of the index design, but probably helps a tick.
            if edit_distance(j[-1*n:], k[-1*n:]) + edit_distance(j[:n], k[:n]) > max_dist:
                continue
            if j == k:
                continue
            if edit_distance(j,k) <= max_dist:
                #similarity_network[j].update([k])
                self.g.add_edge(j,k)
    def build_network(self, verbose=False):
        n=0
        start = time.time()
        now = 0
        last = 0
        for j in wordlist:
            check_indexes(j)
            if verbose:
                n+=1
                now = int(time.time() - start)/60
                if now%5 == 0 and now == last+1:
                    print n, int(time.time() - start)/60
                last = now
    def wordbridge(self, w1, w2):
        try:
            return nx.shortest_path(self.g, w1.upper(), w2.upper())
        except: # NetworkXNoPath error
            return "No path between %s and %s" % (w1, w2)

Idea: Correlate Reddit Usage With Weather

Clearly, this blog has largely become  a project blog for my analyses of reddit (particular /r/warhingtondc), so I’m going to stop starting posts by going into the background of the project.

Normally, I’ve just been posting results, but I want to get more into the habit of commenting on my whole process, so this is the only “brainstorming” post so far but will probably be the first of many.

Hypothesis

Do people use reddit more when it’s shitty out. Can we show this analytically? Are there some very active users who are likely to stop redditing when it’s really nice out? Can we identify homebodies (people who’s redditing behavior is affected very little by how nice the weather is)?

Potential Methodology

1: Aggregate reddit usage BY DAY for r/dc. Might not have enough data to do something meaningful here with just r/dc, but we can try. Consider determining average usage for each user by day of week and normalizing relative to that mean. Then derive an average based on available users (otherwise we’ll always see more usage later in the dataset since there will be more users). Consider limiting dataset to users who have commented within a certain threshold time frame towards the end of the dataset to ensure we don’t grab any “dead” users.

2: download wunderground data for the appropriate time series. Try to draw predictions based on precipitation, temperature, etc.

One option: work backwards. Identify days during which it rained all day, during the afternoon only, etc (not just any day it rained). See how these buckets affect the analysis.

Getting Pandas

No, I didn’t go to the store and buy a panda. In fact, the national zoo here in DC lost a panda fairly recently 😦

The title of this post is meant to contrast my last posts frustration: I’m starting to get the whole pandas thing.

When I last played with my reddit data, I produced an interesting graphic that showed the reddit activity of DC area redditors. In my analysis, I treated the data as though it were descriptive of the redditors as individuals, as though it were an average of each individual’s usage pattern, but really it was aggregate data. I’m interested in trying to characterize different “kinds” of redditors from their usage patterns, and hope to also build a classifier that can predict a user’s timezone. The obvious next step in my analysis was to get more granular and work with the individual users. As I mentioned in my last post I’m challenging myself to do this analysis in python.

My dataset resides in a SQLite database. I have extensive experience with SQL at work, where I basically live in an oracle database, but I’m lucky in that I can do most of my analysis directly in the database. I’m still getting the hang of how applications I program should interface with the database. So, like any noob, I started out by making a pretty stupid, obvious mistake. Let me show you what my code looked like when I started this analysis:

"""
NB: This code is for sample purposes only and represents
what NOT to do. The better solution is posted directly below.

Many of my posts describe good code in the context of bad code:
if you're going to cut and paste from anything I post here, please
make sure to read the whole post to make sure you aren't using the
bad solution to your problem.
"""
import database as db
from datamodel import *
import sqlalchemy as sa

s = db.Session()

def get_users(min_comments=900):
    usertuples = s.query(Comment.author) \
    .group_by(Comment.author) \
    .having(sa.func.count(Comment.comment_id)>min_comments) \
    .all()
    usernames = []
    for u in usertuples:
        usernames.append(u[0])
    return usernames

def get_comment_activity(username):
    """
    Gets comments from database, develops a "profile" for the user by
    summarizing their % activity for each hour in a 24 horus clock.
    Could make this more robust by summarizing each day separately
    """
    comments = s.query(Comment.created_utc)\
            .filter(Comment.author == username)\
            .all()
    timestamps = []
    for c in comments:
        timestamps.append(c[0])
    return timestamps

def profile_user_hourly_activity(data):
    """
    Returns profile as a total count of activity.
    Should ultimately normalize this by dividing by total user's activity
    to get percentages. Will need to do that with numpy or some other
    numeric library.
    """
    pivot = {}
    for d in data:
        timestamp = time.gmtime(d)
        t = timestamp.tm_hour
        pivot[t] = pivot.get(t,0) +1
    return pivot

users = get_users()
profiles = {}
times = []
for user in users:
    comments = get_comment_activity(user)
    profiles[user] = profile_user_hourly_activity(comments)

I’ve simplified my code slightly for putting it up on this blog, but it used to contain some print statements inside the for loop at the bottom so I could see how long each iteration was taking. On average, it took about 2.2 seconds to process a single user. That means to process all 2651 users in my database, it would take about an hour. HOLY CRAP! Instead of fixing my problem, I decided to trim the dataset down to only those users for whom I had been able to download >900 comments (hence the min_comments parameter to get_users), reducing the dataset to 564 users. I actually ran the code and it took 21 minutes. Absolutely unacceptable. What’s the problem?

I/O IS EXPENSIVE.

This is a problem I’ve made before in the past, and it always leads to significant performance gains if it can be eliminated. Like, code that once took 18 hours finished in about 4 minutes running on the same data set because I minimized the I/O (in that instance it was file read/writes).

Instead of a separate database call for each user, let’s hit the database once and slice the data as necessary from inside python.

import database as db
import pandas as pd

s = db.Session()
data = s.query(Comment.author, Comment.comment_id, Comment.created_utc).all()

author        = [a for a,b,c in data]
comment_id    = [b for a,b,c in data]
timestamps    = [time.gmtime(c) for a,b,c in data]
created_dates = [dt.datetime(t.tm_year, t.tm_mon, t.tm_mday, t.tm_hour, t.tm_min, t.tm_sec) for t in timestamps]

df = pd.DataFrame({'author':author, 'comment_id':comment_id}, index=created_dates)

Much better. Also, putting it in the data frame object makes the rest of the analysis much much easier. Well… easier in terms of how fast the code runs and how little I need to write. I’m still learning pandas so doing it the “hard” way actually produces results faster for me, but I’m not on a deadline of anything so let’s do this the right way.

grouped = df.groupby(['author', lambda x: x.hour]).agg(len)
grid    = grouped.unstack()

Holy shit that was easy.Since that happened so fast you might have missed it, those two lines of code are summarizing my data set by user to give a profile of their sum activity by hour. Did I mention this code takes about a minute to run (compare this with the code at the top of the article that would have taken a full HOUR).

Here’s how to plot this and what we get:

import pylab as pl
pl.plot(grid.T)
pl.show()

NB: Times here are in UTC (EST=UTC-5). Also, although I didn’t go to the trouble of labeling axis, you can tell pretty quickly: the x-axis is hour, the y-axis is count.

You can see that there’s a lot of diversity in these users, but there’s a very clear pattern here which is clearly consistent with the graphic I produced in my first pass at this. Let’s take a look at that pattern in a broad way by taking the average of the elements composing this graph:

pl.plot(grid.mean())
pl.show()

Man…pandas makes everything so…simple! And it integrates crazy nice with matplotlib (imported here as part of pylab).

As I mentioned at the top, I suspect that there are distinct user categories that I can identify from how people use the site, and my immediate hypothesis is that I can find these categories based on the times people use the site (e.g. people who browse predominantly at work, people who browse predominantly after dinner, people who browse predominantly on the weekends, etc.). I was having a ton of trouble installing scikit-learn so I cheated and did some analysis in R. Here’s what PCA of this data looks like (just to add more pretty pictures):

Obviously the clusters just aren’t there yet. Maybe there are two: the horizontal ellipse on the left and the vertical ellipse on the right (maybe), but more likely instead of distinct user types, it looks like there may be a predominant user type (the right blob), and then a fully populated spectrum of other usage styles. I can’t remember if I generated this using the full data set, or just the users for whom I had >900 data points, so I could be wrong. I guess I’m going to need some more features.

Anyway, enough blogging. I’ve got a mid term this week and am taking a rare night off from my books to play around with this data, and instead of coding I’ve been writing up analyses I’ve already done. Time to play!

Playing With Pandas: DataFrustration

After putting up the post a few weeks ago where I began to analyze the /r/washingtondc community , I received some feedback regarding my methodology. The comment that most resonated with me was “Why are you trying to do your analysis in R if you’ve already got the data in python? Why not just do the analysis and develop your visualizations there?” I had good reasons to use R, mainly: I saw an opportunity to learn plyr, and I don’t really know my way around the numerical programming and plotting tools in python (primarily numpy, scipy, pylab,matplotlib, scikit-learn, and pandas).

As I’ve mentioned previously, I recently started grad school and I’m going to need R specifically in my program, so my need to learn R is much more urgent than my need to develop expertise in python. But frankly, I just enjoy using python more, and I’d rather do everything in one environment if I can so learning the python data analytics libraries is a very attractive prospect to me. Wes McKinney, the author of the pandas library, is publishing a book through O’Reilly called Python for Data Analysis that I’m super interested in. The problem is right before I started my graduate program, I splurged a little and bought some textbooks I had been drooling over but have barely had opportunity to crack into (Drew Conway and John Myles White‘s Machine Learning for Hackers and Christopher M. Bishop’s Pattern Recognition and Machine Learning) so I’m apprehensive to buy more books until I can bite into the books I’ve already bought (while I’m dreaming of a world where I have time, I also really want to pick up Nate Silver’s new book The Signal And The Noise).

Instead of buying a new book (I expect I”ll buy it later) for the moment let’s do this the ol’ fashioned way and hit the documentation. The documentation for pandas looks to be pretty thorough: I recommend starting with the Intro to Data Structures.

Now, it might be because I’m coming into this with very little NumPy experience, but right off the bat the DataFrame class, the heart of pandas, seems awkward to me. To illustrate my problem, let me compare how to take a particular slice in R vs. pandas.

What I want to do is isolate those columns of the dataframe whose first row is below some threshold value:

A <- runif(10)
B <- runif(10)
C <- runif(10)
D <- runif(10)
E <- runif(10)

df <- data.frame(A,B,C,D,E)
sliced_df <- df[ , df[1,]<.5 ]

This is pretty straightforward. The first 5 lines I create labeled vectors of 10 random numbers selected from the uniform distribution (in the range (0,1)). Then I form a dataframe from these vectors. Finally, from the dataframe I select all rows, but only those columns from the where the value of the first row is less than 0.5.

> df
            A           B         C           D          E
1  0.45274205 0.543755858 0.2225730 0.643710467 0.44527644
2  0.55692168 0.687034039 0.8480953 0.494917616 0.98080695
3  0.19127556 0.419290813 0.2744206 0.005422064 0.58559636
4  0.58410947 0.094669003 0.3284746 0.891122109 0.05962251
5  0.94561895 0.022608545 0.9431832 0.951050056 0.38312492
6  0.72316858 0.003073411 0.3336150 0.201627465 0.89433597
7  0.02145899 0.685167549 0.5754166 0.371717998 0.06746820
8  0.47334489 0.143967454 0.4463423 0.959687645 0.64947595
9  0.75215197 0.068791088 0.0343898 0.117595073 0.28861395
10 0.78567118 0.398529395 0.6467450 0.883467028 0.86369047

> sliced_df
            A         C          E
1  0.45274205 0.2225730 0.44527644
2  0.55692168 0.8480953 0.98080695
3  0.19127556 0.2744206 0.58559636
4  0.58410947 0.3284746 0.05962251
5  0.94561895 0.9431832 0.38312492
6  0.72316858 0.3336150 0.89433597
7  0.02145899 0.5754166 0.06746820
8  0.47334489 0.4463423 0.64947595
9  0.75215197 0.0343898 0.28861395
10 0.78567118 0.6467450 0.86369047

Here’s how we do the same thing in python, using the pandas.DataFrame datatype (this might not look so bad, but DataFrame.T is the transpose method): <<Edit: Actually, since finishing this post I figured out a better way which I discuss towards the end o fthis post>>

import pandas as pd
import random

def runif(n):
    res = []
    for i in range(n):
        res.append(random.random())
    return res

A = runif(10)
B = runif(10)
C = runif(10)
D = runif(10)
E = runif(10)

df = pd.DataFrame({'A':A,'B':B,'C':C,'D':D,'E':E})
sliced_df = df.T[ df.T[0]>0.5 ].T

To start with, I have to create a function to build sequences of random numbers. I strongly suspect that some such function exists in NumPy or somewhere, but I don’t know what it is so I made my own (which has the added benefit of making my code align better with the R example). <<EDIT: Since typing this blogpost, I’ve learned that the function I needed was numpy.random.randn.>>

Next, populating the dataframe is already extremely awkward. From the documentation, it seems like the pandas.DataFrame is primarily designed to be generated by dictionaries. What particularly annoys me about this is that, considering the DataFrame is supposed to be a container for the pandas.Series class, I can’t just feed pandas.DataFrame() a list of pandas.Series objects, they have to be in a dictionary. I can create a DataFrame from a single Series like so:

>>> df = pd.DataFrame(A)
>>> df

0
0  0.446686
1  0.696637
2  0.265383
3  0.165647
4  0.861040
5  0.347737
6  0.280177
7  0.980586
8  0.334544
9  0.645143

Note how the column name was changed to “0”. I would have expected that the variable name becomes the column name like in R, but nope, no such luck. The same thing happens if we feed DataFrame() two series, which makes it hard to distinguish them for slicing.

>>> df = pd.DataFrame(A,B)
>>> df

0
0.689893  0.446686
0.373250  0.696637
0.619527  0.265383
0.775759  0.165647
0.819883  0.861040
0.099763  0.347737
0.143239  0.280177
0.303460  0.980586
0.975462  0.334544
0.672082  0.645143

>>> df[0]

0.689893    0.446686
0.373250    0.696637
0.619527    0.265383
0.775759    0.165647
0.819883    0.861040
0.099763    0.347737
0.143239    0.280177
0.303460    0.980586
0.975462    0.334544
0.672082    0.645143
Name: 0

If I try to create a dataframe from three or more series, I just get an error.

>>> df = pd.DataFrame(A,B,C)

Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "C:\Python27\lib\site-packages\pandas\core\frame.py", line 412, in __init__
copy=copy)
File "C:\Python27\lib\site-packages\pandas\core\frame.py", line 506, in _init_ndarray
block = make_block(values.T, columns, columns)
File "C:\Python27\lib\site-packages\pandas\core\internals.py", line 461, in make_block
do_integrity_check=do_integrity_check)
File "C:\Python27\lib\site-packages\pandas\core\internals.py", line 26, in __init__
assert(len(items) == len(values))
AssertionError

Frankly, I’m not sure why two series even worked. Anyway, the DataFrame object wants a dictionary so it knows explicitly how to name each data set. You’d think the variable names would be sufficient, and that’s the intuitive solution we get in R, but no such luck in pandas. For the time being, this is a minor problem but we’ll see what happens when I start working with larger datasets.

What really irks me is the slicing operations are a little funny, and things that should be able to handle booleans can’t. Basically, slicing dataframes is much less intuitive in pandas than in R or octave (matlab) and this concerns me, since this sort of thing is what enables vectorizing operations (i.e. instead of using slower for loops).

The main problem is that the df[x] operator doesn’t take both rows and columns at the same time like in R (e.g. df[rows, columns] ). But weirdly, it can take one or the other separately. In pandas, df[x] takes a couple of different kinds of ‘x’: a column name, a list of column names, indices, or a conditional statement. If I give df[x] a conditional statement or indices (e.g. df[0:3]), they apply to the rows. Did you catch that? If I feed it labels that match column names, I get columns back. If I feed it slice indices I get rows back. What the fuck, pandas. I can’t be the only person who thinks this is spectacularly confusing. A particular slice method (e.g. [ ] ) should give me back either columns or rows, not both.

>>> df['A']

0    0.446686
1    0.696637
2    0.265383
3    0.165647
4    0.861040
5    0.347737
6    0.280177
7    0.980586
8    0.334544
9    0.645143
Name: A

>>> df[['A','B']]

A         B
0  0.446686  0.689893
1  0.696637  0.373250
2  0.265383  0.619527
3  0.165647  0.775759
4  0.861040  0.819883
5  0.347737  0.099763
6  0.280177  0.143239
7  0.980586  0.303460
8  0.334544  0.975462
9  0.645143  0.672082

>>> df[ df['A']>0.5 ]

A         B         C         D         E
1  0.696637  0.373250  0.778385  0.918453  0.580366
4  0.861040  0.819883  0.397820  0.735031  0.110483
7  0.980586  0.303460  0.117398  0.033969  0.731914
9  0.645143  0.672082  0.138268  0.730780  0.969545

>>> df[1:3]

A         B         C         D         E
1  0.696637  0.373250  0.778385  0.918453  0.580366
2  0.265383  0.619527  0.633576  0.553009  0.599043

This means if I want just those columns that meet a particular row condition, I have to first take the transpose of the dataframe to make the rows into columns, apply my conditional statement to the appropriate column (which requires taking the transpose again), and take the transpose one last time to flip the returned data set back into the original orientation. Hence:

sliced_df = df.T[ df.T[0]>0.5 ].T

Having written out this entire rant, I finally figured out the better way of doing this. I’m still annoyed by all the stuff I described above, but here’s how you accomplish this slice without using transpose at all: you use the DataFrame.ix method. This allows for more intuitive R-style indexing of the form df.ix[rows, columns]. It accepts indices and booleans. Do I sound dejected? This took me a while to figure out.

>>> df.ix[:,:]

A         B         C         D         E
0  0.242808  0.829024  0.027734  0.510985  0.430466
1  0.301553  0.834208  0.600806  0.773148  0.119008
2  0.968252  0.098827  0.290203  0.555629  0.652359
3  0.351365  0.391068  0.352370  0.531282  0.478862
4  0.513526  0.138082  0.538826  0.252554  0.486603
5  0.705628  0.362105  0.800225  0.977828  0.454140
6  0.097671  0.613972  0.712334  0.473130  0.886449
7  0.386206  0.520115  0.589156  0.722709  0.293428
8  0.337381  0.102242  0.296870  0.725426  0.475001
9  0.076314  0.894782  0.115159  0.592838  0.402849

>>> df.ix[2:3,'B':'E']

B         C         D         E
2  0.098827  0.290203  0.555629  0.652359
3  0.391068  0.352370  0.531282  0.478862

>>> df.ix[0,:] < 0.5

A     True
B    False
C     True
D    False
E     True
Name: 0

>>> df.ix[ :, df.ix[0,:] < 0.5 ]

A         C         E
0  0.242808  0.027734  0.430466
1  0.301553  0.600806  0.119008
2  0.968252  0.290203  0.652359
3  0.351365  0.352370  0.478862
4  0.513526  0.538826  0.486603
5  0.705628  0.800225  0.454140
6  0.097671  0.712334  0.886449
7  0.386206  0.589156  0.293428
8  0.337381  0.296870  0.475001
9  0.076314  0.115159  0.402849

One thing that is much less confusing about these dataframes is if I want to plot something, all I have to do is:

from pylab import plot, show

plot(df)
show()

And I get a single plot where each line is a separate column of data. I would have had to have used a for loop otherwise, so that at least is nice and intuitive.

These are normal growing pains, as with learning any new tool. If you’ve got a few hours on your hands and don’t feel like mucking through the documentation, here’s a pycon tutorial Wes gave earilier this year that’s pretty interesting (I’ve only made my way through about an hour’s worth myself, but it looks like a thorough introduction). NB: The first 15 minutes or so he’s just setting up IPython, so I recommend you skip ahead a bit.