Sets and Dictionaries

Let's start with something simpler than our actual inventory problem. Suppose we have a list of all the atoms in the warehouse, and we want to know which different kinds we have—not how many, but just their types. We could solve this problem using a list to store the unique atomic symbols we have seen. Here's a function to add a new atom to the list:

def another_atom(seen, atom):
    for i in range(len(seen)):
        if seen[i] == atom:
            return # atom is already present, so do not re-add
    seen.append(atom)

another_atom's arguments are a list of the unique atoms we've already seen, and the symbol of the atom we're adding. Inside the function, we loop over the atoms that are already in the list. If we find the one we're trying to add, we exit the function immediately: we aren't supposed to have duplicates in our list, so there's nothing to add. If we reach the end of the list without finding this symbol, though, we append it. This is a common design pattern: either we find pre-existing data in a loop and return right away, or take some default action if we finish the loop without finding a match.

Let's watch this function in action. We start with an empty list. If the first atomic symbol is 'Na', we find no match (since the list is empty), so we add it. The next symbol is 'Fe'; it doesn't match 'Na', so we add it as well. Our third symbol is 'Na' again. It matches the first entry in the list, so we exit the function immediately.

This code works, but it is inefficient. Suppose there are V distinct atomic symbols in our data, and N symbols in total. Each time we add an observation to our list, we have to look through an average of V/2 entries. The total running time for our program is therefore approximately NV/2. If V is small, this is only a few times larger than N, but what happens if we're keeping track of something like patient records rather than atoms? In that case, most values are distinct, so V is approximately the same as N, which means that our running time is proportional to N²/2. That's bad news: if we double the size of our data set, our program runs four times slower, and if we double it again, our program will have slowed down by a factor of 16.

There's a better way to solve this problem that is simpler to use and runs much faster. The trick is to use a set to store the symbols. A set is an unordered collection of distinct items. The word "collection" means that a set can hold zero or more values. The word "distinct" means that any particular value is either in the set or not: a set can't store two or more copies of the same thing. And finally, "unordered" means that values are simply "in" the set. They're not in any particular order, and there's no first value or last value. (They actually are stored in some order, but as we'll discuss in the next section, that order is as random as the computer can make it.)

To create a set, we simply write down its elements inside curly braces:

>>> primes = {3, 5, 7}

However, we have to use set() to create an empty set, because the symbol {} was already being used for something else when sets were added to Python:

>>> even_primes = set() # not '{}' as in math

We'll meet that "something else" later in this chapter.

To see what we can do with sets, let's create three holding the integers 0 through 9, the first half of that same range of numbers (0 through 4), and the odd values 1, 3, 5, 7, and 9:

>>> ten  = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
>>> lows = {0, 1, 2, 3, 4}
>>> odds = {1, 3, 5, 7, 9}

If we ask Python to display one of our sets, it shows us this:

>>> print lows
set([0, 1, 2, 3, 4])

rather than using the curly-bracket notation. I personally regard this as a design flaw, but it does remind us that we can create always create a set from a list.

Sets have methods just like strings and lists, and, like the methods of strings and lists, most of those methods create new sets instead of modifying the set they are called for. These three come straight from mathematics:

>>> print lows.union(odds)
set([0, 1, 2, 3, 4, 5, 7, 9])
>>> print lows.intersection(odds)
set([1, 3])
>>> print lows.difference(odds)
set([0, 2, 4])

Another method that creates a new set is symmetric_difference, which is sometimes called "exclusive or":

>>> print lows.symmetric_difference(odds)
set([0, 2, 4, 5, 7, 9])

It returns the values that are in one set or another, but not in both.

Not all set methods return new sets. For example, issubset returns True or False depending on whether all the elements in one set are present in another:

>>> print lows.issubset(ten)
True

A complementary method called issuperset also exists, and does the obvious thing:

>>> print lows.issuperset(odds)
False

We can count how many values are in a set using len (just as we would to find the length of a list or string), and check whether a particular value is in the set or not using in:

>>> print len(odds)
7
>>> print 6 in odds
False

Some methods modify the sets they are called for. The most commonly used is add, which adds an element to the set:

>>> lows.add(9)
>>> print lows
set([0, 1, 2, 3, 4, 9])

If the thing being added is already in the set, add has no effect, because any specific thing can appear in a set at most once:

>>> lows.add(9)
>>> print lows
set([0, 1, 2, 3, 4, 9])

This behavior is different from that of list.append, which always adds a new element to a list.

Finally, we can remove individual elements from the set:

>>> lows.remove(0)
>>> print lows
set([1, 2, 3, 4])

or clear it entirely:

>>> lows.clear()
>>> print lows
set()

Removing elements is similar to deleting things from a list, but there's an important difference. When we delete something from a list, we specify its location. When we delete something from a set, though, we must specify the value that we want to take out, because sets are not ordered. If that value isn't in the set, remove does nothing.

To help make programs easier to type and read, most of the methods we've just seen can be written using arithmetic operators as well. For example, instead of lows.issubset(ten), we can write lows <= ten, just as if we were using pen and paper. There are even a couple of operators, like the strict subset test <, that don't have long-winded equivalents.

The fact that the values in a set are distinct makes them a convenient way to get rid of duplicate values, like the "unique atoms" problem at the start of this section. Suppose we have a file containing the names of all the atoms in our warehouse, and our task is to produce a list of the their types. Here's how simple that code is:

def unique_atoms(filename):
    atoms = set()
    with open(filename, 'r') as source:
        for line in source:
            name = line.strip()
            atoms.add(name)
    return atoms

We start by creating an empty set which we will fill with atomic symbols and opening the file containing our data. As we read the lines in the file, we strip off any whitespace (such as the newline character at the end of the line) and put the resulting strings in the set. When we're done, we print the set. If our input is the file:

Na
Fe
Na
Si
Pd
Na

then our output is:

set(['Na', 'Fe', 'Si', 'Pd'])

The right atoms are there, but what are those extra square brackets for? The answer is that if we want to construct a set with values using set(), we have to pass those values in a single object, such as a list. This syntax:

set('Na', 'Fe', 'Si', 'Pd')

does not work, even though it seems more natural. On the other hand, this means that we can construct a set from almost anything that a for loop can iterate over:

>>> set('lithium')
set(['i', 'h', 'm', 'l', 'u', 't'])

But hang on: if we're adding characters to the set in the order 'l', 'i', 't', 'h', 'i', 'u', 'm', why does Python show them in the order 'i', 'h', 'm', 'l', 'u', 't'? To answer that question, we need to look at how sets are actually stored, and why they're stored that way.

Let's create a set, add the string 'lithium' to it (as a single item, not character by character), and print the result:

>>> things = set()
>>> things.add('lithium')
>>> print things
set(['lithium'])

As expected, the string is in the set. Now let's try adding a list to the same set:

>>> things.add([1, 2, 3])
TypeError: unhashable type: 'list'

Why doesn't that work? And what does that word "unhashable" mean?

When we create a set, the computer allocates a block of memory to store references to the set's elements. When we add something to the set, or try to look something up, the computer uses a hash function to figure out where to look. A hash function is any function that produces a seemingly-random number when given some data as input. For example, one way to hash a string is to add up the numerical values of its characters. If the string is "zebra", those values are 97 for lower-case 'a', 98 for lower-case 'b', and so on up to 122 for lower-case 'z'. When we add them up, we will always get the same result: in this case, 532. If our hash table has 8 slots, we can take the remainder 532%8=4 to figure out where to store a reference to our string in the hash table (Figure 2).

Now let's take a look at how a list would be stored. If the list contains the same five characters, so that its hash code is still 4, it would be stored as shown in Figure 3:

But what happens if we change the characters in the list after we've added it to the set? For example, suppose that we change the first letter in the list from 'z' to 'X'. The hash function's value is now 498 instead of 532, which means that the modified list belongs in slot 2 rather than slot 4. However, the reference to the list is still in the old location: the set doesn't know that the list's contents have changed, so it hasn't moved its reference to the right location (Figure 4):

This is bad news. If we now ask, "Is the list containing 's', 'e', 'b', 'r', and 'a' in the set?" the answer will be "no", because the reference to the list isn't stored in the location that our hash function tells us to look. It's as if someone changed their name from "Tom Riddle" to "Lord Voldemort", but we left all the personnel records filed under 'R'.

This problem arises with any mutable structure—i.e., any structure whose contents or value can be changed after its creation. Integers and strings are safe to hash because their values are fixed, but the whole point of lists is that we can grow them, shrink them, and overwrite their contents.

Different languages and libraries handle this problem in different ways. One option is to have each list keep track of the sets that it is in, and move itself whenever its values change. However, this is expensive: every time a program touched a list, it would have to see if it was in any sets, and if it was, recalculate its hash code and update all the references to it.

A second option is to shrug and say, "It's the programmer's fault." This is what most languages do, but it's also expensive: programmers can spend hours tracking down the bugs that arise from data being in the wrong place.

Python uses a third option: it only allows programmers to put immutable values in sets. After all, if something's value can't change, neither can its hash code or its location in a hash table.

But if sets can only hold immutable values, what do we do with mutable ones? In particular, how should we store things like (x,y) coordinates, which are naturally represented as lists, or people's names, which are naturally represented as lists of first, middle, and last names? Again, there are several options.

The first is to concatenate those values somehow. For example, if we want to store "Charles" and "Darwin", we'd create the string "Charles Darwin" and store that. This is simple to do, but our code will wind up being littered with string joins and string splits, which will make it slower to run and harder to read. More importantly, it's only safe to do if we can find a concatenator that can never come up in our data. (If we join "Paul Antoine" and "St. Cyr" using a space, there would be three possible ways to split it apart again.)

The second option—the right one—is to use tuples instead of lists. A tuple is an immutable list, i.e., a sequence of values that cannot be changed after its creation. Tuples are created exactly like lists, except we use parentheses instead of square brackets:

>>> full_name = ('Charles', 'Darwin')

They are indexed the same way, too, and functions like len do exactly what we'd expect:

>>> print full_name[0]
Charles
>>> print len(full_name)
2

What we cannot do is assign a new value to a tuple element, i.e., change the tuple after it has been created:

>>> full_name[0] = 'Erasmus'
TypeError: 'tuple' object does not support item assignment

This means that a tuple's hash code never changes, and that means that tuples can be put in sets:

>>> names = set()
>>> names.add(('Charles', 'Darwin'))
>>> print names
set([('Charles', 'Darwin')])

Now that we know how to find out what kinds of atoms are in our inventory, we want to find out how many of each we have. Our input is a list of several thousand atomic symbols, and the output we want is a list of names and counts.

Once again, we could use a list to store names and counts, but the right solution is to use another new data strucure called a dictionary. A dictionary is a unordered collection of key-value pairs (Fixture 5). The keys are immutable, unique, and unordered, just like the elements of a set. There are no restrictions on the values stored with those keys: they don't have to be immutable or unique. However, we can only look up entries by their keys, not by their values.

We create a new dictionary by putting key-value pairs inside curly braces with a colon between the two parts of each pair:

>>> birthdays = {'Newton' : 1642, 'Darwin' : 1809}

The dictionary's keys are the strings 'Newton' and 'Darwin'. The value associated with 'Newton' is 1642, while the value associated with 'Darwin' is 1809. We can think of this as a two-column table:

but it's important to remember that the entries aren't necessarily stored in this order (or any other specific order).

We can get the value associated with a key by putting the key in square brackets:

>>> print birthdays['Newton']
1642

This looks just like subscripting a string or list, except dictionary keys don't have to be integers—they can be strings, tuples, or any other immutable object. It's just like using a phonebook or a real dictionary: instead of looking things up by location using an integer index, we look things up by name.

If we want to add another entry to a dictionary, we just assign a value to the key, just as we create a new variable in a program by assigning it a value:

>>> birthdays['Turing'] = 1612
>>> print birthdays
{'Turing' : 1612, 'Newton' : 1642, 'Darwin' : 1809}

If the key is already in the dictionary, assignment replaces the value associated with it rather than adding another entry (since each key can appear at most once). Let's fix Turing's birthday by replacing 1612 with 1912:

>>> birthdays['Turing'] = 1912
>>> print birthdays
{'Turing' : 1912, 'Newton' : 1642, 'Darwin' : 1809}

Trying to get the value associated with a key that isn't in the dictionary is an error, just like trying to access a nonexistent variable or get an out-of-bounds element from a list. For example, let's try to look up Florence Nightingale's birthday:

>>> print birthdays['Nightingale']
KeyError: 'Nightingale'

If we're not sure whether a key is in a dictionary or not, we can test for it using in:

>>> print 'Nightingale' in birthdays
False
>>> print 'Darwin' in birthdays
True

And we can see how many entries are in the dictionary using len:

>>> print len(birthdays)
3

and loop over the keys in a dictionary using for:

>>> for name in birthdays:
...     print name, birthdays[name]
...
Turing 1912
Newton 1642
Darwin 1809

This is a little bit different from looping over a list. When we loop over a list we get the values in the list. When we loop over a dictionary, on the other hand, the loop gives us the keys, which we can use to look up the values.

We're now ready to count atoms. The main body of our program looks like this:

def main(filename):
    counts = count_atoms(filename)
    for atom in counts:
        print atom, counts[atom]

count_atoms reads atomic symbols from a file, one per line, and creates a dictionary of atomic symbols and counts. Once we have that dictionary, we use a loop like the one we just saw to print out its contents.

Here's the function that does the counting:

def count_atoms(filename):
    '''Count unique atoms, returning a dictionary.'''

    result = {}
    with open(filename, 'r') as reader:
        for line in reader:
            atom = line.strip()
            if atom not in result:
                result[atom] = 1
            else:
                result[atom] = result[atom] + 1
    return result

We start with a docstring to explain the function's purpose to whoever has to read it next. We then create an empty dictionary to fill with data, and use a loop to process the lines from the input file one by one. Notice that the empty dictionary is written {}: this is the "previous use" we referred to when explaining why an empty set had to be written set().

After stripping whitespace off the atom's symbol, we check to see if we've seen it before. If we haven't, we set its count to 1, because we've now seen that atom one time. If we have seen it before, we add one to the previous count and store that new value back in the dictionary. When the loop is done, we return the dictionary we have created.

Let's watch this function in action. Before we read any data, our dictionary is empty. After we see 'Na' for the first time, our dictionary has one entry: its key is 'Na', and its value is 1. When we see 'Fe', we add another entry to the dictionary with that string as a key and 1 as a value. Finally, when we see 'Na' for the second time, we add one to its count.

Just as we use tuples for multi-part entries in sets, we can use them for multi-part keys in dictionaries. For example, if we want to store the years in which scientists were born using their full names, we could do this:

birthdays = {
    ('Isaac', 'Newton') : 1642,
    ('Charles', 'Robert', 'Darwin') : 1809,
    ('Alan', 'Mathison', 'Turing') : 1912
}

If we do this, though, we always have to look things up by the full key: there is no way to ask for all the entries whose keys contain the word 'Darwin', because Python cannot match part of a tuple.

If we think of a dictionary as a two-column table, it is occasionally useful to get one or the other column, i.e., just the keys or just the values:

all_keys = birthdays.keys()
print all_keys
[('Isaac', 'Newton'), ('Alan', 'Mathison', 'Turing'), ('Charles', 'Robert', 'Darwin')]
all_values = birthdays.values()
print all_values
[1642, 1912, 1809]

These methods should be used sparingly: the dictionary doesn't store the keys or values in a list, these methods both actually create a new list as their result. In particular, we shouldn't loop over a dictionary's entries like this:

for key in some_dict.keys():
    ...do something with key and some_dict[key]

since "for key in some_dict" is shorter and much more efficient.

To see how useful dictionaries can be, let's switch tracks and do some birdwatching. We'll start by asking how early in the day we saw each kind of bird? Our data consists of the date and time of the observation, the bird's name, and an optional comment:

2010-07-03    05:38    loon
2010-07-03    06:02    goose
2010-07-03    06:07    loon
2010-07-04    05:09    ostrich   # hallucinating?
2010-07-04    05:29    loon
     …           …        …

Rephrasing our problem, we want the minimum of all the times associated with each bird name. If our data was stored in memory like this:

loon = ['05:38', '06:07', '05:20', ...]

the solution would simply be min(loon), and similarly for the other birds. However, we have to work with the data we have, so let's start by reading our data file and creating a list of tuples, each of which contains a date, time, and bird name as strings:

def read_observations(filename):
    '''Read data, returning [(date, time, bird)...].'''

    reader = open(filename, 'r')
    result = []

    for line in reader:
        fields = line.split('#')[0].strip().split()
        assert len(fields) == 3, 'Bad line "%s"' % line
        result.append(fields)

    return result

This function follows the pattern we've seen many times before. We set up by opening the input file and creating an empty list that we'll append records to. We then process each line of the file in turn. Splitting the line on the '#' character and taking the first part of the result gets rid of any comment that might be present; stripping off whitespace and then splitting breaks the remainder into fields.

To prevent trouble later on, we check that there actually are three fields before going on. (An industrial-strength version of this function would also check that the date and time were properly formatted, but we'll skip that for now.) Once we've done our check, we append the triple containing the date, time, and bird name to the list we're going to return.

Here's the function that turns that list of tuples into a dictionary:

def earliest_observation(data):
    '''How early did we see each bird?'''

    result = {}
    for (date, time, bird) in data:
        if bird not in result:
            result[bird] = time
        else:
            result[bird] = min(result[bird], time)

    return result

Once again, the pattern should by now be familiar. We start by creating an empty dictionary, then use a loop to inspect each tuple in turn. The loop explodes the tuple into separate variables for the date, time and bird. If the bird's name is not already a key in our dictionary, this must be the first time we've seen it, so we store the time we saw it in the dictionary. If the bird's name is already there, on the other hand, we keep the minimum of the stored time and the new time. This is almost exactly the same as our earlier counting example, but instead of either storing 1 or adding 1 to the count so far, we're either storing the time or taking the minimum of it and the least time so far.

Now, what if we want to find out which birds were seen on particular days? Once again, we are aggregating values, i.e., combining many separate values to create one new one. However, since we probably saw more than one kind of bird each day, that "new value" needs to be a collection of some kind. We're only interested in which birds we saw, so the right kind of collection is a set. Here's our function:

def birds_by_date(data):
    '''Which birds were seen on each day?'''

    result = {}
    for (date, time, bird) in data:
        if date not in result:
            result[date] = {bird}
        else:
            result[date].add(bird)

    return result

Again, we start by creating an empty dictionary, and then process each tuple in turn. Since we're recording birds by date, the keys in our dictionary are dates rather than bird names. If the current date isn't already a key in the dictionary, we create a set containing only this bird, and store it in the dictionary with the date as the key. Otherwise, we add this bird to the set associated with the date. (As always, we don't need to check whether the bird is already in that set, since the set will automatically eliminate any duplication.)

Let's watch this function in action for the first few records from our data:

For our last example, we'll figure out which bird we saw least frequently—or rather, which birds, since two or more may be tied for the low score. Forgetting that values may not be unique is a common mistake in data crunching, and often a hard one to track down.

Our first strategy is simple: figure out how many times we've seen each bird, then find the minimum of those counts and get the set of birds we've seen that many times. The function below implements this fairly directly:

def least_common_birds(data):
    '''Which bird or birds have been seen least frequently?'''
    
    counts = count_by_bird(data)
    least = min(counts.values())
    result = set()
    for bird in counts:
        if counts[bird] == least:
            result.add(bird)
    return result

least_common_birds depends on a function count_by_bird, but this is yet another example of using a dictionary to aggregate values (in this case, to sum the number of birds we have seen). Just for variety's sake, we'll use a slightly different strategy that we've used before: whenever we see a new kind of bird, we'll set its count to zero, and then always add one to the stored count:

def count_by_bird(data):
    '''How many times was each bird seen?'''
    result = {}
    for (date, time, bird) in data:
        if bird not in result:
            result[bird] = 0
        result[bird] += 1
    return result

Finally, we'll test our function:

print least_common_birds(entries)
set(['goose', 'ostrich'])

This does the job, but is somewhat inefficient: we do one pass through all the data while counting birds, then another pass through all the birds to find those that we've seen the least number of times. We can actually do the whole job with a single pass through the data, but as we'll see in the challenges, the resulting code is significantly more complex than what we have written so far. Unless we're sure that the second pass is really a performance bottleneck, we should stick with this simple implementation.

We can now solve Fan's original nanotech inventory problem. As explained in the introduction, our goal is to find out how many molecules of various kinds we can make using the atoms in our warehouse. The number of molecules of any particular type we can make is limited by the scarcest atom that molecule requires. For example, if we have five nitrogen atoms and ten hydrogen atoms, we can only make three ammonia molecules, because we need three hydrogen atoms for each.

The formulas for the molecules we know how to make are stored in a file like this:

# Molecular formula file

helium : He 1
water : H 2 O 1
hydrogen : H 2

and our inventory is stored in a file like this:

# Atom inventory file

He 1
H 4
O 3

Let's start by reading in our inventory. It consists of pairs of strings and numbers, which by now should suggest using a dictionary for storage. The keys will be atomic symbols, and the values will be the number of atoms of that kind we currently have (Figure 6). If an atom isn't listed in our inventory, we'll assume that we don't have any.

What about the formulas for the molecules we know how to make? Once again, we want to use strings—the names of molecules—as indices, which suggests a dictionary. Each of its values will be something storing atomic symbols and the number of atoms of that type in the molecule—the same structure, in fact, that we're using for our inventory. Figure 7 shows what this looks like in memory if the only molecules we know how to make are water and ammonia.

Finally, we'll store the results of our calculation in yet another dictionary, this one mapping the names of molecules to how many molecules of that kind we can make (Figure 8).

The main body of the program is straightforward: it reads in the input files, does our calculation, and prints the result:

def main(inventory_file, formula_file):
    inventory = read_inventory(inventory_file)
    formulas = read_formulas(formula_file)
    counts = calculate_counts(inventory, formulas)
    show_counts(counts)

Reading the inventory file is simple. We take each interesting line in the file, split it to get an atomic symbol and a count, and store them together in a dictionary:

def read_inventory(inventory_file):
    result = {}
    with open(inventory_file, 'r') as reader:
        for line in reader:
            name, count = line.strip().split()
            result[name] = int(count)
    return result

Let's test it:

print read_inventory('inventory-03.txt')
ValueError                                Traceback (most recent call last)
 in ()
----> 1 print read_inventory('inventory-03.txt')

 in read_inventory(inventory_file)
      3     with open(inventory_file, 'r') as reader:
      4         for line in reader:
----> 5             name, count = line.strip().split()
      6             result[name] = int(count)
      7     return result

ValueError: too many values to unpack

Our mistake was to forget that files can contain blank lines and comments. It's easy enough to modify the function to handle them, though it complicates the logic:

def read_inventory(inventory_file):
    result = {}
    with open(inventory_file, 'r') as reader:
        for line in reader:
            line = line.strip()
            if (not line) or line.startswith('#'):
                continue
            name, count = line.split()
            result[name] = int(count)
    return result

print read_inventory('inventory-03.txt')
{'H': 4, 'O': 3, 'He': 1}

The next step is to read the files containing formulas. Since the file format is more complicated, the function is as well. In fact, it's complicated enough that we'll come back later and simplify it.

def read_formulas(formula_file):
    result = {}
    with open(formula_file, 'r') as reader:
        for line in reader:
            line = line.strip()
            if (not line) or line.startswith('#'):
                continue
            name, atoms = line.split(':')
            name = name.strip()
            atoms = atoms.strip().split()
    
            formula = {}
            for i in range(0, len(atoms), 2):
                symbol = atoms[i].strip()
                count = int(atoms[i+1])
                formula[symbol] = count
            result[name] = formula

    return result

We start by creating a dictionary to hold our results. We then split each interesting line in the data file on the colon ':' to separate the molecule's name (which may contain spaces) from its formula. We then split the formulas into a list of strings. These alternate between atomic symbols and numbers, so in the inner loop, we move forward through those values two elements at a time, storing the atomic symbol and count in a dictionary. Once we're done, we store that dictionary as the value for the molecule name in the main dictionary. When we've processed all the lines, we return the final result. Here's a simple test:

print read_formulas('formulas-03.txt')
{'water': {'H': 2, 'O': 1}, 'hydrogen': {'H': 2}, 'helium': {'He': 1}}

Now that we have all our data, it's time to calculate how many molecules of each kind we can make. inventory maps atomic symbols to counts, and so does formulas[name], so let's loop over all the molecules we know how to make and "divide" the inventory by each one:

def calculate_counts(inventory, formulas):
    '''Calculate how many of each molecule can be made with inventory.'''

    counts = {}
    for name in formulas:
        counts[name] = dict_divide(inventory, formulas[name])

    return counts

We say we're "dividing" the inventory by each molecule because we're trying to find out how many of that molecule we can make without requiring more of any particular atom than we actually have. (By analogy, when we divide 11 by 3, we're trying to find out how many 3's we can make from 11.) The function that does the division is:

def dict_divide(inventory, molecule):
    number = None
    for atom in molecule:
        required = molecule[atom]
        available = inventory.get(atom, 0)
        limit = available / required
        if (number is None) or (limit < number):
            number = limit

    return number

This function loops over all the atoms in the molecule we're trying to build, see what limit the available inventory puts on us, and return the minimum of all those results. This function uses a few patterns that come up frequently in many kinds of programs:

1. The first pattern is to initialize the value we're going to return to None, then test for that value inside the loop to make sure we re-set it to a legal value the first time we have real data. In this case, we could just as easily use -1 or some other impossible value as an "uninitialized" flag for number.
2. Since we're looping over the keys of molecule, we know that we can get the value stored in molecule[atom]. However, that atom might not be a key in inventory, so we use inventory.get(atom, 0) to get either the stored value or a sensible default. In this case zero, the sensible default is 0, because if the atom's symbol isn't in the dictionary, we don't have any of it. This is our second pattern.
3. The third is using calculate, test, and store to find a single value—in this case, the minimum—from a set of calculated values. We could calculate the list of available over required values, then find the minimum of the list, but doing the minimum test as we go along saves us having to store the list of intermediate values. It's probably not a noticeable time saving in this case, but it would be with larger data sets.

The last step in building our program is to show how many molecules of each kind we can make. We could just loop over our result dictionary, printing each molecule's name and the number of times we could make it, but let's put the results in alphabetical order to make it easier to read:

def show_counts(counts):
    names = counts.keys()
    names.sort()
    for name in names:
        print name, counts[name]

It's time to test our code. Let's start by using an empty inventory and a single formula:

# inventory-00.txt

# formulas-00.txt

There's no output, which is what we expect. Let's add a formula but no atoms:

# inventory-00.txt

# formulas-01.txt
helium : He 1

helium 0

and now an atom:

# inventory-01.txt
He 1

# formulas-01.txt
helium : He 1

helium 1

That seems right as well. Let's add some hydrogen and another formula:

# inventory-02.txt
He 1
H 4

# formulas-01.txt
helium : He 1
water : H 2 O 1

helium 1
water 0

The output doesn't change, which is correct. Our final test adds some oxygen:

# inventory-03.txt
He 1
H 4
O 3

# formulas-03.txt
helium : He 1
water: H 2 O 1
hydrogen : H 2

helium 1
water 2

That's right too: we can make two water molecules (because we don't have enough hydrogen to pair with our three oxygen atoms).

def readlines(filename):
    result = []
    with open(filename, 'r') as reader:
        for line in reader:
            line = line.strip()
            if line and (not line.startswith('#')):
                result.append(line)
    return result

def read_inventory(inventory_file):
    result = {}
    for line in readlines(inventory_file):
        name, count = line.split()
        result[name] = int(count)
    return result

def read_formulas(formula_file):
    result = {}
    for line in readlines(formula_file):
        name, atoms = line.split(':')
        name = name.strip()
        atoms = atoms.strip().split()

        formula = {}
        for i in range(0, len(atoms), 2):
            symbol = atoms[i].strip()
            count = int(atoms[i+1])
            formula[symbol] = count
        result[name] = formula

    return result

def read_formulas(formula_file):
    result = {}
    for line in readlines(formula_file):
        name, atoms = line.split(':')
        name = name.strip()
        result[name] = make_formula(atoms)
    return result

def make_formula(atoms):
    formula = {}
    atoms = atoms.strip().split()
    for i in range(0, len(atoms), 2):
        symbol = atoms[i].strip()
        count = int(atoms[i+1])
        formula[symbol] = count
    return formula

The example above used two data file formats: one for storing molecular formulas, the other for storing inventory. Both formats were specific to this application, which means we needed to write, debug, document, and maintain functions to handle them. Those functions weren't particularly difficult to create, but they still took time to create, and if anyone ever wants to read our files in Java, MATLAB, or Perl, they'll have to write equivalent functions themselves.

A growing number of programs avoid these problems by using a flexible data format called JSON, which stands for "JavaScript Object Notation". Despite the name, it is a language-independent way to store nested data structures made up of strings, numbers, Booleans, lists, dictionaries, and the special value null (equivalent to Python's None)—in short, the basic data types that almost every language supports. For example, let's convert a dictionary of scientists' birthdays to a string:

>>> import json
>>> birthdays = {'Curie' : 1867, 'Hopper' : 1906, 'Franklin' : 1920}
>>> as_string = json.dumps(birthdays)
>>> print as_string
{"Curie": 1867, "Hopper": 1906, "Franklin": 1920}
>>> print type(as_string)
<type 'str'>

json.dumps doesn't seem to do much, but that's kind of the point: the textual representation of the data structure looks pretty much like what a programmer would type in to re-create it. The advantage is that this representation can be saved in a file:

>>> import json
>>>
>>> writer = open('/tmp/example.json', 'w')
>>> json.dump(birthdays, writer)
>>> writer.close()
>>>
>>> reader = open('/tmp/example.json', 'r')
>>> duplicate = json.load(reader)
>>> reader.close()
>>>
>>> print 'original:', birthdays
original: {'Curie': 1867, 'Hopper': 1906, 'Franklin': 1920}
>>> print 'duplicate:', duplicate
duplicate: {u'Curie': 1867, u'Hopper': 1906, u'Franklin': 1920}
>>> print 'original == duplicate:', birthdays == duplicate
original == duplicate: True
>>> print 'original is duplicate:', birthdays is duplicate
original is duplicate: False

As the example above shows, saving and loading data is as simple as opening a file and then calling one function. The data file holds what we'd type in to create the data in a program:

$ cat /tmp/example.json
{"Curie": 1867, "Hopper": 1906, "Franklin": 1920}

which makes it easy to edit by hand.

How is this different in practice from what we had? First, our inventory file now looks like this:

{"He" : 1,
 "H" : 4,
 "O" : 3}

while our formulas files look like:

{"helium"   : {"He" : 1},
 "water"    : {"H" : 2, "O" : 1},
 "hydrogen" : {"H" : 2}}

Those aren't as intuitive for non-programmers as the original flat text files, but they're not too bad. The worst thing is the lack of comments: unfortunately—very unfortunately—the JSON format doesn't support them. (And note that JSON requires us to use a double-quote for strings: unlike Python, we cannot substitute single quotes.)

The good news is that given files like these, we can rewrite our program as:

'''Calculate how many molecules of each type can be made with the atoms on hand.'''
import json

def main(inventory_file, formulas_file):
    '''Main driver for program.'''
    with open(inventory_file, 'r') as reader:
        inventory = json.load(reader)
    with open(formulas_file, 'r') as reader:
        formulas = json.load(reader)
    counts = calculate_counts(inventory, formulas)
    show_counts(counts)

def calculate_counts(inventory, formulas):
    ...as before...

def dict_divide(inventory, molecule):
    ...as before...

def show_counts(counts):
    ...as before...

The two functions that read formula and inventory files have been replaced with a single function that reads JSON. Nothing else has to change, because the data structures loaded from the data files are exactly what we had before. The end result is 51 lines long compared to the 80 we started with, a reduction of more than a third.

Nothing's Perfekt

JSON's greatest weakness isn't its lack of support for comments, but the fact that it doesn't recognize and manage aliases. Instead, each occurrence of an aliased structure is treated as something brand new when data is being saved. For example:

>>> inner = ['name']
>>> outer = [inner, inner] # Creating an alias.
>>> print outer
[['name'], ['name']]
>>> print outer[0] is outer[1]
True
>>> as_string = json.dumps(outer)
>>> duplicate = json.loads(as_string)
>>> print duplicate
[[u'name'], [u'name']]
>>> print duplicate[0] is duplicate[1]
False

Figure 9 shows the difference between the original data structure (referred to by outer) and what winds up in duplicate. If aliases might be present in our data, and it's important to preserve their structure, we must either record the aliasing ourself (which is tricky), or use some other format. Luckily, a lot of data either doesn't contain aliases, or the aliasing in it isn't important.

As Theodosius Dobzhansky said almost a century ago, nothing in biology makes sense except in the light of evolution. Since mutations usually occur one at a time, the more similarities there are between the DNA of two species, the more recently they had a common ancestor. We can use this idea to reconstruct the evolutionary tree for a group of organisms using a hierarchical clustering algorithm.

We don't have to look at the natural world very hard to realize that some organisms are more alike than others. For example, if we look at the appearance, anatomy, and lifecycles of the seven fish shown in Figure 10, we can see that three pairs are closely related. But where does the seventh fit in? And how do the pairs relate to each other?

The first step is to find the two species that are most similar, and construct their plausible common ancestor. We then pair two more, and two more, and start joining pairs to individuals, or pairs with other pairs. Eventually, all the organisms are connected. We can redraw those connections as a tree, using the heights of branches to show the number of differences between the species we're joining up (Figure 11).

Let's turn this into an algorithm:

S = {all organisms}
while S != {}:
  a, b = two closest entries in U
  p = common parent of {a, b}
  S = S - {a, b}
  S = S + {p}

Initially, the set S contains all the species we're interested in. Each time through the loop, we find the two that are closest, create their common parent, remove the two we just paired up from the set, and insert the newly-created parent. Since the set shrinks by one element each time (two out, one in), we can be sure this algorithm eventually terminates.

But how do we calculate the distance between an inferred parent and other species? One simple rule is to use the average distance between that other species and the two species that were combined to create that parent. Let's illustrate it by calculating a phylogenetic tree for humans, vampires, werewolves, and mermaids. The distances between each pair of species is shown in Figure 12 and in the table below. (We only show the lower triangle because it's symmetric.)

The closest entries—i.e., the pair with minimum distance—are human and werewolf. We replace this with a common ancestor, which we will call HW, then set the distance between it and each other species X to be (HX + WX)/2, i.e., the average of the human-to-X and werewolf-to-X distances. This gives us a new table:

Repeating this step, we combine HW with V:

and finally HWV with M.

We illustrated our algorithm with a triangular matrix, but the order of the rows and columns is arbitrary. The matrix is really just a lookup table mapping species to distances, and as soon as we think of lookup tables, we should think of dictionaries. The keys are species—either the ones we started with, or the ones we created—and the values are the distances between them, so our original table becomes:

{
    ('human',   'mermaid')  : 12,
    ('human',   'vampire')  : 13,
    ('human',   'werewolf') :  5,
    ('mermaid', 'vampire')  : 15,
    ('mermaid', 'werewolf') : 29,
    ('vampire', 'werewolf') :  6
}

There is one trick here. Whenever we have a distance, such as that between mermaids and vampires, we have to decide whether to use the key ('mermaid', 'vampire') or ('vampire', 'mermaid') (or to record the value twice, once under each key).

Let's start by setting up our test case and then calling a top-level function to process our data:

if __name__ == '__main__':

    species = {'human', 'mermaid', 'werewolf', 'vampire'}

    scores = {
        ('human',   'mermaid')  : 12,
        ('human',   'vampire')  : 13,
        ('human',   'werewolf') :  5,
        ('mermaid', 'vampire')  : 15,
        ('mermaid', 'werewolf') : 29,
        ('vampire', 'werewolf') :  6
    }

    order = main(species, scores)
    print order

In a real program, of course, the data would be read in from a file, and the set of actual species' names would be generated from it, but this will do for now.

Next, let's translate our algorithm into something that could be runnable Python:

def main(species, scores):
    result = []
    while len(species) > 1:
        left, right = find_min_pair(species, scores)
        result.append(make_pair(left, right))
        species -= {left, right}
        make_new_pairs(species, scores, left, right)
        species.add(make_name(left, right))
    return result

This is almost a direct translation of our starting point; the only significant difference is that we're keeping the set of "active" species in a set, and the scores in a dictionary. As species are combined, we remove their names from the set and add a made-up name for their parent. We never actually remove scores from the table; once the name of a species is out of the set species, we'll never try to look up anything associated with it in scores again.

The next step is to write find_min_pair to find the lowest score currently in the table:

def find_min_pair(species, scores):
    min_pair = None
    min_val = None
    for left in species:
        for right in species:
            if left < right:
                this_pair = make_pair(left, right)
                if (min_val is None) or (scores[this_pair] < min_val):
                    min_pair = this_pair
                    min_val = scores[this_pair]
    return min_pair

This function loops over all possible combinations of species names, but only actually uses the ones that pass our ordering test (i.e., the ones for which the first species name comes before the second species name). If this is the first score we've looked at, or if it's lower than a previously-seen score, we record the pair of species and the associated score. When we're done, we return the pair of species.

The function that makes new entries for the table is fairly straightforward as well. It just loops over all the active species, averages the distances between them and the two species being combined, and puts a new score in the table:

def make_new_pairs(species, scores, left, right):
    for current in species:
        left_score = scores[make_pair(current, left)]
        right_score = scores[make_pair(current, right)]
        new_score = (left_score + right_score) / 2.0
        scores[make_pair(current, make_name(left, right))] = new_score

Finally, the make_pair and make_name functions are simply:

def make_pair(left, right):
    if left < right:
        return (left, right)
    else:
        return (right, left)

def make_name(left, right):
    return '<%s, %s>' % make_pair(left, right)

Let's try running the program:

$ python phylogen.py
[('human', 'werewolf'), ('<human, werewolf>', 'vampire'), ('<<human, werewolf>, vampire>', 'mermaid')]

This shows that humans and werewolves were combined first, that their pairing was then combined with vampires, and that mermaids were added to the cluster last. We obviously should do a lot more testing, but so far, we seem to be on the right track.