Reading line-structured files (FASTA) with Python

Many files contain data in a line-structured format, where a line prefix indicates a record or fields in a record. An example are FASTA files, where ';' or '>' indicate a header or comment lines, followed by multiple lines with sequence data:

  
;LCBO - Prolactin precursor
; a sample sequence in FASTA format
MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVSLRDLFDRAVMVSHYIHDLSS
EMFNEFDKRYAQGKGFITMALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHL
VTEVRGMKGAPDAILSRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDED
ARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC*

>MCHU - Calmodulin 
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
DIDGDGQVNYEEFVQMMTAK*

>gi|5524211|gb|AAD44166.1
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY

Files in these formats are surprisingly unpleasant/difficult to read if they are large and cannot be loaded completely into memory. The reason is that we need to use an Iterator to read the file line by line and cannot go back.

Before looking at the solutions below, give it a try yourself! The task is to read the file above (without loading it in memory via open().read()) and return tuples of the form (header, sequence), where header is only the first line of multiple header lines, if there are any, and sequence is the concatenation of all sequence data. In other words, we want a generator or iterator that produces the following tuples:

  
(';LCBO - Prolactin precursor', 'MDSKGSS...NNC*')

('>MCHU - Calmodulin', 'ADQLTEE...TAK*')

(>gi|5524211|gb|AAD44166.1', 'LCLYTHIG...ENY')

The sequence strings above are shortened (...) for readability. Note that newlines within the sequence data need to be removed and that there can be multiple blank lines. Here my solution in plain Python:

def read_seqs():
    header, seqs = None, []
    isheader, prev_isheader = False, False
    for line in open('fasta.txt'):
        line = line.strip()
        if not line: 
            continue
        prev_isheader = isheader
        isheader = line[0] in {'>', ';'}
        if isheader and not prev_isheader:
            if header:
                yield header, ''.join(seqs)
                header, seqs = line, []
        elif not isheader:
            seqs.append(line)
    yield header, ''.join(seqs)

Now we can read and print the sequences one-by-one:

for h,s in read_seqs():
    print '{}\n{}\n'.format(h,s)

However, the code is ugly, complex and it took me embarrassing long to write it. A simple task such as reading a line-structured file should not be that difficult and we will have a look at a much better solution later. But first a quick discussion of the code above.

We read the file 'fasta.txt' with open(), which gives us an iterator over the lines in the file. Newlines are removed and empty lines are skipped. isheader tells us if a line starts with a header prefix and prev_isheader contains this state for the previous line. If we find a header prefix and the previous line was not a header line, we start a new record. If it is not a header line we collect the sequence data. Note that we have to call yield twice. Once when we start a new record we return the record we have completed reading. And the second time when we are finished reading the file. Two flags, two yields, multiple nested if statements, and more than 10 lines to read a simple file is shockingly complex.

Luckily, with itertools the same problem can be solve much more elegantly. In just four lines (not counting imports), to be precise:

from itertools import groupby, izip

lines = (l.strip() for l in open('fasta.txt') if l.strip())
isheader = lambda line: line[0] in {'>', ';'}
chunks= (list(g) for _, g in groupby(lines, isheader))
seqs = ((h[0], ''.join(s)) for h,s in izip(chunks, chunks))

Now this is beautiful! Much shorter, more readable and easy to implement. The key is the groupby function that simplifies our life dramatically here. But let's start from the top.

The first line is a a generator that iterates over the lines in the file, strips newlines and skips empty lines. The second line defines a lambda function isheader that recognizes header lines. In the third line it gets interesting. groupby groups a stream of data (here lines) according to a key function. If you are not familiar with groupby, here a small example to understand it works:

>>> from itertools import groupby
>>> numbers = [1, 2, 3, 4, 5, 1]
>>> gt2 = lambda x : x > 2
>>> [(key, list(group)) for key, group in groupby(numbers, gt2)]
[(False, [1, 2]), (True, [3, 4, 5]), (False, [1])]

We group a list of numbers depending on if they are greater than 2 or not. groupby reads the numbers one-by-one and whenever the key function tt>gt2 changes its value it creates a new group. groupby returns an iterator over the keys and groups it produces, and each group itself is an iterator over its elements. In the example above we use list to convert the group iterator to a list and then use a list comprehension to print the keys and groups.

Note that groupby almost does what we need to read the records in our FASTA file. We can group the lines and start a new group whenever we encounter a header line. Since we don't need the result of the key function we can simply write:

chunks = (list(g) for _, g in groupby(lines, isheader))

Each chunk is a list of lines, containing either all header lines or all sequence lines. For instance, if we print the chunks

for chunk in chunks:
    print chunk

we get the following output (sequences are shortened for readability):

[';LCBO - Prolactin precursor', '; a sample sequence in FASTA format']
['MDS...LSS', 'EMF...YHL', 'VTE...DED', 'ARY...NNC*']
['>MCHU - Calmodulin']
['ADQ...GTID', 'FPE...REA', 'DID...TAK*']
['>gi|5524211|gb|AAD44166.1']
['LCL...NLV', 'EWI...DFLG', 'LLI...IVIL', 'GLM...AGX', 'IENY']

What is left to do is to bring the header lines and their sequence lines together, and there is a very neat trick to achieve this. Remember this one! As you know, zip zips iterables:

>>> zip([1, 2, 3], 'abc')
[(1, 'a'), (2, 'b'), (3, 'c')]

Under Python 2.x zip return a list and under Python 3.x it returns an iterator. izip under Python 2.x does the same as zip under Python 3.x and returns an iterator. More importantly, however, zip and izip operate on iterators and here comes the trick! If we want to group elements of an iterable let's say in tuples of size 2 we can pass the same iterator twice and zip will do the job:

>>> numbers = iter([1, 2, 3, 4, 5, 6])
>>> zip(numbers, numbers)
[(1, 2), (3, 4), (5, 6)]

Note that we wrap numbers into an iterator by calling iter(). Without that, the result would be different:

>>> numbers = [1, 2, 3, 4, 5, 6]
>>> zip(numbers, numbers)
[(1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6)]

What's going on here? zip reads the an element from numbers and reads again to build the first tuple. Since numbers is an iterator the element pointer progresses and we get (1, 2). We could equally easily group numbers in triples:

>>> numbers = iter([1, 2, 3, 4, 5, 6])
>>> zip(numbers, numbers, numbers)
[(1, 2, 3), (4, 5, 6)]

Obviously, this works also for our FASTA file reading. We can group the chunks of header or sequence lines in tuples of size 2 to bring header and sequence together. We use izip because we don't want to collect all sequences in memory. So the fourth line of

seqs = ((h[0], ''.join(s)) for h,s in izip(chunks, chunks))

the itertools solution above means that we group header lists and sequence lists in tuples h, s. But we need only the first header line and the concatenation of the sequence lines and therefore return (h[0], ''.join(s)) in the generator expression.

There you go! A much simpler implementation to read FASTA files and other line-structured files. Can we do even better? Using a library for data flows such as nuts-flow leads to even shorter and more readable code:

from nutsflow import *
from nutsflow import _

isheader = lambda line: line[0] in {'>', ';'}
seqs = (open('fasta.txt') >> Map(str.strip) >> Filter(_) >> 
        ChunkBy(isheader) >> Chunk(2) >> MapCol(0,next) >> MapCol(1,''.join))

The >> operator in nuts-flow indicates the chaining of iterators and defines the data flow. Here we open the file, map the string stripping function on each line to remove white spaces and then filter out empty lines. The underscore in Filter(_) functions as an anonymous variable here. Similar to Python's groupby, ChunkBy groups the lines according to the isheader predicate. We then call Chunk(2) to bring header and sequence lines together in tuples of the form (header, sequence) -- no need for the zip-trick we used before. Finally, we map next() on column 0 and join on column 1 of these tuples to extract the first header line and to concatenate the sequence lines.

That's it. Can't get better than that ;)

Fun with itertools in Python

Let's say we are given a file 'number1.txt' that contains a header, numbers and blank lines, and the task is to sum up all numbers within the file. Here an example file

  
MyNumbers
1
2

3
4

If the file can be loaded into memory this is easy

  
>>> lines = list(open('numbers1.txt'))
>>> sum(int(l) for l in lines[1:] if l.strip())
10

open() returns an iterator over the lines in the file, which we collect in a list. Alternatively we could have used open().readlines(). We call lines[1:] to skip the header line and test for non-empty lines with if l.strip(). A sum over a list comprehension finally returns the wanted result. Note that for this simple example we don't bother with closing the file, which you should do for robust code.

If the file is too big to be loaded into memory things are getting a bit ugly. open() returns an iterator, which is good but we need to skip the first header line and cannot apply list slicing to an iterator. However, we can advance the iterator by calling next() to solve this issue

  
>>> lines = open('numbers1.txt')
>>> next(lines)
'MyNumbers\n'
>>> sum(int(l) for l in lines if l.strip())
10

But what if there are multiple header lines and multiple files we want to sum over? Let's implement this in plain Python and then use itertools to solve this more elegantly.

from glob import glob

def filesum(fname, skip=1):   # sum numbers in file
  lines = open(fname)
  for _ in xrange(skip):
      next(lines)
  return sum(int(l) for l in lines if l.strip())

sum(filesum(fname) for fname in glob('numbers*.txt'))

As you can see things are getting a tad ugly, especially with respect to the skipping of multiple header lines. We are going to use islice, ifilter, and imap from the itertools library to clean things up. Note that the names of most functions in itertools start with an 'i' to indicate that they are returning iterators and do not collect results in memory.

Similar to list slicing, islice() takes a slice out of an iterator. Here we use it to skip one header line

>>> from itertools import islice
>>> list(islice(open('numbers1.txt'), 1, None))
['1\n', '2\n', '\n', '3\n', '4']

where islice(iterable, start, stop[, step]) takes an iterable as input and skips the first start elements of it, if called with None for the stop parameter. Since islice returns an iterator we wrap it into a list to display the result.

The next step is to filter out the blank lines and we employ ifilter() for this purpose

>>> from itertools import islice, ifilter

>>> lines = islice(open('numbers1.txt'), 1, None)
>>> list(ifilter(lambda l: l.strip(), lines))
['1\n', '2\n', '3\n', '4']

In the next step we convert the number strings into integers by mapping the int function on the lines of the file using imap()

>>> from itertools import islice, ifilter, imap

>>> alllines = islice(open('numbers1.txt'), 1, None)
>>> numlines = ifilter(lambda l: l.strip(), alllines)
>>> list(imap(int, numlines))
[1, 2, 3, 4]

Note that imap, in contrast to map returns an iterator in Python 2, while in Python 3, map is equivalent to imap. Having an iterator over integer numbers, we can now sum them up

>>> from itertools import islice, ifilter, imap
>>> SKIP = 1
>>> alllines = islice(open('numbers1.txt'), SKIP , None)
>>> numlines = ifilter(lambda l: l.strip(), alllines)
>>> sum(imap(int, numlines))
10

What is left to do is to aggregate over multiple files and similarly to our first implementation we define a separate function filesum for this purpose

from glob import glob
from itertools import islice, ifilter, imap, chain

def filesum(fname, skip=1):
    isnumber = lambda l: l.strip()
    lines = islice(open(fname), skip, None)
    return sum(imap(int, ifilter(isnumber, lines)))
    

sum(imap(filesum, glob('numbers*.txt')))

If you compare both implementations, I hope you will find the itertools-based code a bit more readable and elegant than the plain Python code presented at the start. But still the code seems overly complex for such a simple problem. An alternative to itertools is nuts-flow, which allows implementing data flows, such as the one above, much more easily. For instance, to sum up the numbers within the file we could simply write

>>> from nutsflow import Drop, Map, Sum, nut_filter
>>> IsNumber = nut_filter(lambda l: l.strip())

>>> open('numbers1.txt') >> Drop(1) >> IsNumber() >> Map(int) >> Sum()
10

where Drop(1) skips the header line, IsNumber filters out the numbers (non-empty lines), Map(int) converts to integers and Sum sums the numbers up. The explicit data flow and linear sequence of processing steps renders this code considerably more readable and it can easily be extended to operate over multiple files

>>> from glob import glob
>>> from nutsflow import Drop, Map, Sum, nut_filter

>>> IsNumber = nut_filter(lambda l: l.strip())
>>> filesum = lambda f: open(f) >> Drop(1) >> IsNumber() >> Map(int) >> Sum()
>>> glob('numbers*.txt') >> Map(filesum) >> Sum()
10

More about nuts-flow in a follow-up post. That's it for today ;)

lambda functions in Python

The most common (and recommended) way to implement a function in Python is via the def keyword and function name. For instance

  
def add(a, b): return a + b 

is a function that adds two numbers and can be called as

  
>>> add(1, 2)
3

We can (re)define the same function as a lambda function

  
add = lambda a, b: a + b

and it can be invoked with exactly the same syntax as before

  
>>> add(1, 2)
3

However, we can also define and invoke a lambda function anonymously (without giving it a name). Here some examples

  
>>> (lambda a, b: a + b)(1, 2)      # add
3
>>> (lambda a: a * 2)(2)            # times 2
4
>>> (lambda a: sum(a))([1, 2, 3])   # sum of list
6
>>> (lambda *a: min(a))(1, 2, 3)    # min of list
1

The syntactical differences between conventional function definitions via def and lambda functions are that lambda functions have no brackets around their arguments, don't use the return statement, do not have a name and must be implemented in a single expression.

So, what are they good for? Fundamentally there is no need for lambda functions and we could happily live without them but they make code shorter and more readable when a simple function is needed only once. The most common use case is as key function for sort:

  
>>> names = ['Fred', 'Anna', 'John']
>>> sorted(names)                        # alphabetically
['Anna', 'Fred', 'John']

>>> sorted(names, key=lambda n: n[1])    # second letter
['Anna', 'John', 'Fred']

Especially for collections of data structures key functions come in handy for sorting

  
>>> contacts = [('Fred', 18), ('Anna', 22)]

>>> sorted(contacts)                               # name then age
[('Anna', 22), ('Fred', 18)]

>>> sorted(contacts, key=lambda (name, _): name)   # name
[('Anna', 22), ('Fred', 18)]

>>> sorted(contacts, key=lambda (_, age): a)       # age
[('Fred', 18), ('Anna', 22)]

>>> sorted(contacts, key=lambda (n, a): a, n)      # age then name
[('Fred', 18), ('Anna', 22)]

Note the brackets in lambda (name, _) and similar calls that result in tuple unpacking and a very readable sorting implementation by name or age (removed in Python 3, see PEP 3113). We could have achieved the same by defining a conventional function

  
>>> def by_age((name, age)): return age
>>> sorted(contacts, key=by_age)
[('Fred', 18), ('Anna', 22)]

which is also very readable but longer and pollutes the name space with a by_age function that is needed only once. Note that lambda functions are not needed if we already have a function! For instance

  
>>> names = ['Fred', 'Anna', 'John']

>>> sorted(names, key=lambda n: len(n))    # don't do this
['Fred', 'Anna', 'John']

>>> sorted(names, key=len)                 # simply do this
['Fred', 'Anna', 'John']

An elegant example for the usage of lambda functions is the implementation of argmax, which returns the index of the largest element in a list (and the element itself)

  
>>> argmax = lambda ns: max(enumerate(ns), key=lambda (i,n): n)
>>> argmax([2, 3, 1])
(1, 3)

In Python 3 we can use the following implementation instead, which is also short but less readable

  
>>> argmax = lambda ns: max(enumerate(ns), key=lambda t: t[1])
>>> argmax([2, 3, 1])
(1, 3)

While applications for lambda functions in typical Python code are limited they occur frequently in APIs for graphical user interface, e.g. to react on buttons pressed or in code that employs a functional programming style, which uses filter, map, reduce or similar functions that take functions as arguments. For instance, the accumulative difference of positive numbers provided as strings could be computed (in functional style) as follows

>>> from operator import sub  
>>> numbers = ['10', '-5', '3', '2']  # we want to compute 10-3-2
>>> reduce(sub, filter(lambda x: x > 0, map(int, numbers)))  
5

However, list comprehension and generators are often more readable and recommended. Here the same computation using a generator expression

>>> reduce(sub, (int(n) for n in numbers if int(n) > 0)) 
5

Note, however, that we have to convert strings to integers twice and that there is no elegant way to compute the accumulative difference without using reduce, which is not available in Python 3 anymore :( So in plain, boring, imperative Python we would have to write

numbers = ['10', '-5', '3', '2']
diff = None
for nstr in numbers:
   n = int(nstr)
   if n <= 0: 
      continue
   diff = n if diff is None else diff - n

Pay special attention to the ugly edge case of the first element when accumulating the differences. We need to initialize diff with a marker value, here None, and then have to test for it in the last line. Ugly! The functional implementation is clearly more readable and shorter in this case (and similar problems).

Luckily the functools library comes to the rescue and brings back reduce. Also have a look at itertools, which provides many useful functions for functional programming. Similarily, there is nuts-flow, which allows even more elegant code. Here the accumulative difference using nuts-flow

>>> from operator import sub 
>>> from nutsflow import _, Reduce, Map, Filter

>>> numbers = ['10', '-5', '3', '2']
>>> numbers >> Map(int) >> Filter(_ > 0) >> Reduce(sub)
5

That's it. Have fun ;)

The beauty of FP (and Scala)

I recently came across a coding problem that beautifully demonstrates the benefits of functional programming and - in this case - Scala. The following (Python) code shows a simplified version of the problem.

numbers1 = [1,2,3]
numbers2 = [1,2,3]
pairs = [(n1,n2) for n1 in numbers1 for n2 in number2]
processed = [process(n1,n2) for n1,n2 in pairs]
for p in processed: print p,  

As you can see it is trivial code that pairs the numbers in both lists, processes the pairs and prints out the results of the process function. Obviously, one could simply write

for n1 in numbers1:
  for n2 in numbers2:
    print process(n1,n2) 

and it would be better code but the previous example resembles the structure of the real, more complex code that has several additional processing steps.

To come back to the previous example, there are two possible problems: 1) let's assume it is not a few numbers but many 2) and that the process function is computationally expensive. For large number-lists memory consumption becomes an issue (since all pairs are stored in a list) and the processing will be time consuming.

Let's address the first problem first. We can easily replace list comprehensions by generators and would then iterate over number-pairs instead of storing them and their results. Easy.

pairs = ((n1,n2) for n1 in numbers1 for n2 in number2)
processed = (process(n1,n2) for n1,n2 in pairs)
for p in processed: print p,  

Problem two is harder. Let's assume the order of the results is of no importance and on a multi-core computer it would be nice to distribute the processing over the cores. But as far as I am aware there is no really elegant way to achieve this with Python. I'll show you soon what I mean by "elegant". Let's switch to Scala.

  
  val numbers1 = List(1,2,3)
  val numbers2 = List(1,2,3)
  val pairs = for(n1 <- numbers1; n2 <- numbers2) yield (n1,n2)
  val processed = pairs map { case (n1,n2) => process(n1,n2)}
  processed foreach (print)

As you can see the Scala code does the same as the first Python program. To change from list comprehensions to generators/iterators in the Scala code it is sufficient to transform the first list into an iterator (using toIterator).

  
  val numbers1 = List(1,2,3).toIterator
  val numbers2 = List(1,2,3)
  val pairs = for(n1 <- numbers1; n2 <- numbers2) yield (n1,n2)
  val processed = pairs map {case (n1,n2) => process(n1,n2)}
  processed foreach (print)

In the Python version we had to replace two list comprehensions by generators. In the Scala version we need to modify only one list. pairs and processed automatically become iterators because their inputs are iterators. In languages without generators this simple change would require substantial restructuring of the code.

Let's assume memory is not the problem but we want to speed things up. We simply add .par after pair to convert the list to a parallel collection and the processing is now performed concurrently. Beautiful!

    
  val numbers1 = List(1,2,3)
  val numbers2 = List(1,2,3)
  val pairs = for(n1 <- numbers1; n2 <- numbers2) yield (n1,n2)
  val processed = pairs.par map {case (n1,n2) => process(n1,n2)}
  processed foreach (print)

Can we have both - speed and low memory consumption? Not without effort. A parallel collection cannot be an iterator and vice versa. We could use futures but it wouldn't be trivial. However, this example shows how easy it is to tune functional code or parts of code to be memory efficient or fast/concurrent.

Test if string has all unique characters

Here comes another typical question from coding interviews: "Test whether all characters within a string are unique". An inefficient solution with O(n^2) complexity is to loop over all characters of the string and test in a second, inner loop whether an equal character can be found:

  
def all_unique(str):
    for ch1 in str:
        for ch2 in str:
            if ch1==ch2:
                 return False
    return True 

This is obviously not what the interviewer is after. A simple and efficient solution is to convert the string to a set of characters and to compare the number of elements in the set with the length of the string:

 
def all_unique(str): 
    return len(set(str)) == len(str)   

Too easy and therefore the interviewer doesn't like it either ;-). What is wanted is a mapping of characters to a boolean/bit array that indicates which character has been found in the string so far:

  
def all_unique(str):
    test = [False]*256
    for ch in str:
        idx = ord(ch)
        if test[idx]: return False
        test[idx] = True
    return True 

This approach has linear time complexity but requires an array with the length of the size of the character set.

Implement a queue using two stacks

Implementing a queue using two stacks is a typical question for coding interviews. It is a bit artificial, at least I haven't come across any real need, but hey, maybe you need to implement a queue in FORTH. Apart from that it is a nice puzzle ;)

The key insight is that a stack reverses order (while a queue doesn't). A sequence of elements pushed on a stack comes back in reversed order when popped. Consequently, two stacks chained together will return elements in the same order, since reversed order reversed again is original order.

How to do it? We use an in-stack that we fill when an element is enqueued and the dequeue operation takes elements from an out-stack. If the out-stack is empty we pop all elements from the in-stack and push them onto the out-stack. Here is the implementation in Python:

  
class Queue(object):
    def __init__(self):
        self.instack = []
        self.outstack = []
    
    def enqueue(self,element):
        self.instack.append(element)
   
    def dequeue(self):
        if not self.outstack:
            while self.instack:
                self.outstack.append(self.instack.pop())
        return self.outstack.pop()    

Let's try it:

  
q = Queue()
for i in xrange(5):
    q.enqueue(i)
for i in xrange(5):
    print q.dequeue()

0 1 2 3 4

The complexity is amortized O(1) for enqueue and dequeue.

Tree iterator in Python

Sometimes you have a tree data structure and want to iterate over the nodes of the tree in breadth-first or depth-first order. This also seems to be a question that is occasionally posed in coding interviews. Since it is fun and especially easy to implement in Python let's do it. First things first. We design a Node class that represents the nodes of the tree.

  
class Node(object):

    def __init__(self, value, children=[]):
        self.value = value
        self.children = children

The constructor __init__ takes the value of the node and a list of child nodes that is empty by default. We can now construct a tree such as the following

  
root = Node(1, [Node(2, [Node(4), Node(5)]), Node(3)])

To be able to display the tree we add two method to the Node class. __str__ returns a string representation of a node and showTree performs a recursive, depth-first traversal of the tree and prints the node values indented by their hierarchical level.

  
class Node(object):

    def __init__(self, value, children=[]):
        self.value = value
        self.children = children

    def __str__(self):
        return str(self.value)

    def showTree(self, level=0):
        print "%s%s" % ("."*level, self)
        for node in self.children:
            node.showTree(level+2)

Running root.showTree() for the tree defined above produces the following output. 1 is the root node with two children 2 and 3, where node 2 has 4 and 5 as children.

1
..2
....4
....5
..3

Now, the iterator. There is two ways (actually there are more) to iterate over the nodes. Level by level (= breadth first) or the same way showTree works (depth first). Both are equally easy. The only difference is one implementation uses a stack while the other uses a queue. Let's start with breadth first. The following code implements an iterator class. As such it has to provide two methods: __iter__, which returns the class object itself and next, which returns the next node every time it is called.

  
class BreadthFirst(object):

    def __init__(self, node):
        self.queue = [node]

    def __iter__(self):
        return self

    def next(self):
        if not self.queue: raise StopIteration
        node = self.queue.pop(0)
        self.queue += node.children
        return node

The constructor takes a node and puts it into a list that we will treat as a queue. Therefore we name it "queue" and not "nodes" or even worse "list". What happens in the next method? When the queue is empty (not self.queue) a StopIteration exception is raised that stops the iterator. If the queue is not empty a node is dequeued (terrible word) via queue.pop(0) from the beginning of the list and the children of the node are added to the end of the list (enqueued). The iterator is then called as follows and returns the tree nodes in breadth-first order.

  
for node in BreadthFirst(root):
    print node,

1 3 2 5 4

Using lists as queues is not very efficient and there is a Queue class in Python, which is likely to be faster but the beauty lies in the fact that the implementation of the breadth-first iterator is mirror image of the depth-first iterator. As mentioned before, the only difference is that the list is now used as a stack:

  
class DepthFirst(object):

    def __init__(self, node):
        self.stack = [node]

    def __iter__(self):
        return self

    def next(self):
        if not self.stack: raise StopIteration
        node = self.stack.pop()
        self.stack += node.children
        return node

C'mon, admit it! This IS beautiful. Finally, using the iterator generates the following output

  
for node in DepthFirst(root):
    print node,

1 2 3 4 5

That's it.

Fun with Hopfield and Numpy

Hopfield networks are fun to play with and are very easily implemented in Python using the Numpy library. Simple as they are, they are the basis of modern machine learning techniques such as Deep Learning and programming models for quantum computers such as Adiabatic quantum computation.

We are going to use a Hopfield network for optical character recognition. This is obviously a toy example and there are much more better machine learning techniques for this kind of problem (e.g. deep learning) but it is a nice demonstration.

Among other things Hopfield networks work as autoassociative memories, allowing to store and to recall patterns, e.g. images of letters. We want to store two letters and then try to recognize noisy versions of those letters. In the following we first create images of the letters "A" and "Z" from a string representation:

  
A = """
.XXX.
X...X
XXXXX
X...X
X...X
"""

Z = """
XXXXX
...X.
..X..
.X...
XXXXX
"""

Hopfield networks operate on binary vectors and we therefore define a function to_pattern() that takes a letter image and transforms it into a bipolar vector by gluing all image rows together, and converting "X" to +1 and "." to -1.

  
def to_pattern(letter):
    from numpy import array
    return array([+1 if c=='X' else -1 for c in letter.replace('\n','')])

For instance, calling to_pattern(A) returns then

[-1  1  1  1 -1  1 -1 -1 -1  1  1  1  1  1  1  1 -1 -1 -1  1  1 -1 -1 -1  1]

To visualize a pattern as an image we write a function display(pattern) that reshapes the pattern vector to a matrix and displays it as an image using the pylab function imshow().

def display(pattern):
    from pylab import imshow, cm, show
    imshow(pattern.reshape((5,5)),cmap=cm.binary, interpolation='nearest')
    show()

Calling display(to_pattern(A)); display(to_pattern(Z)) produces the following pictures:

Having all the basic elements in place let's get serious. First thing to do is to create the Hopfield network from the given patterns. Training a Hopfield network is extremely simple and requires only computation of the outer products of all pattern vectors with each others and summation of the resulting matrices.

$$W = ∑↙{k}↖n \bo p_{k} \bo p^T_{k} $$

The resulting matrix W specifies the weights of the edges between the nodes of the network and the nodes represent the features of the patterns (= pixels in an image). There are variants but typically the weight matrix is normalized (division by the number of patterns) and diagonal elements are set to zero (subtraction of the identity matrix I), which leads to the following expression:

$$W = 1/n ∑↙{k}↖n \bo p_{k} \bo p^T_{k} - I$$

Let us organize all our patterns in a matrix (rows are patterns)

 
patterns = array([to_pattern(A), to_pattern(Z)])

and the implementation of the training formula is straight forward:

  
def train(patterns):
    from numpy import zeros, outer, diag_indices 
    r,c = patterns.shape
    W = zeros((c,c))
    for p in patterns:
        W = W + outer(p,p)
    W[diag_indices(c)] = 0
    return W/r

To recall a stored pattern from a Hopfield network the pattern itself (or a partial or nosy version of it) the dot product of the pattern vector and the weight matrix is computed. This produces a new pattern vector that is binarized and then fed back into the system:

$$\bo p^{(t+1)} = sgn( \bo W \bo p^{(t)})$$

The operation is repeated until the pattern vector becomes stable or after a specified number of iteration has been performed. We keep the implementation simple and stop after 5 steps, which is typically sufficient.

  
def recall(W, patterns, steps=5):
    from numpy import vectorize, dot
    sgn = vectorize(lambda x: -1 if x<0 else +1)
    for _ in xrange(steps):        
        patterns = sgn(dot(patterns,W))
    return patterns

Note that in contrast to the recall formula the code performs the recall on all pattern within the pattern matrix at once, while the formula describes the recall for an individual pattern vector. Also the pattern matrix contains the patterns a rows and therefore the dot product W*p in the formula becomes dot(patterns,W) in the code (Since W is symmetric, transposing W is not required).

To test the network I created noisy version of the training patterns by randomly flipping 10 bits. The image sequence below shows the resulting input patterns and the output patterns after the second and third iteration. As you see the original patterns have been recalled perfectly after three iterations (from left to right).

Patterns are stored as low energy states (attractors) within the system and the recall procedure describes the trajectory of an input pattern to its attractor. Each step reduces the energy and the energy for a pattern is defined as follows:

$$E(\bo p) = \bo p \bo W \bo p^T$$

Here is the corresponding code to compute the energy for the pattern in the pattern matrix.

  
def hopfield_energy(W, patterns):
    from numpy import array, dot
    return array([-0.5*dot(dot(p.T,W),p) for p in patterns])

A few final comments. The Hopfield network does not only store the original patterns but also their inverse versions and possibly linear combinations of patterns. Furthermore the capacity of the network is limited. For random patterns the network can store roughly 0.14*N pattern without errors (N = Number of network nodes). There are many variations of the Hopfield network with different training and recall procedures that give better performance but the beauty of the original Hopfield network is in its simple implementation but complex behavior.

Last but not least cudos to the developers of jqmath, the JavaScript module used here to display the mathematical expressions.

more reading

Hopfield network
Hopfield network
Weave a neural net with Python
Associative recall of memory without errors
Minimum Probability Flow Learning
Finite memory loading in hairy neurons
Error tolerant associative memory

Testing for a path through a maze

This is another coding interview question. Find a path through a labyrinth. Let's assume the labyrinth is given as a string of the following form:

maze = """
XXXXXXXXXXX
X         X
X XXXXX   X
X  XX   X X
X   XX    X
XXXXXXXXXXX
"""

The first thing to do is to convert this string into a more suitable data structure. A boolean matrix springs to mind, where true indicates an empty cell and false a stone. It is a one-liner in Python to convert the string given in maze into the boolean matrix:

maze = [[cell!='X' for cell in row] for row in maze.split('\n')[1:]]

The function to test for a path from a starting point start given as tuple (row,col) to an end point end can then be implemented as follows:

def find(maze,start,end):
    stack = [start]
    while stack:
        i,j = stack.pop()
        if not maze[i][j]: continue
        if (i,j) == end: return True
        maze[i][j] = False        
        if maze[i][j+1]: stack.append((i,j+1))
        if maze[i][j-1]: stack.append((i,j-1))
        if maze[i+1][j]: stack.append((i+1,j))
        if maze[i-1][j]: stack.append((i-1,j))
    return False

Instead of a recursive solution we use a stack to explore all possible paths with backtracking. We initialize the stack with the starting point and then loop as long as there are points on the stack or the end point has been found. In each iteration we take the top point from the stack (i,j = stack.pop()) and test if the corresponding cell of the maze is empty. If not we jump to the beginning of the loop. If the cell is empty we mark it as visited (maze[i][j] = False) and leave the loop if it is the end point. If we have not reached the end point we push the coordinates of all neighboring empty cells on the stack and keep iterating.

Let's give it a go:

start = (1,1)
end   = (4,9)      
print find(maze, start, end)

True

One last thing. The code to test for empty neighboring cells is simple and clear but repetitive. Another solution is to replace those four lines by a loop. It isn't much shorter (three instead of four lines) and the code becomes more complex though:

def find(maze,start,end):
    stack = [start]
    while stack:
        i,j = stack.pop()
        if (i,j) == end: return True
        maze[i][j] = False
        if (i,j) == end: return True
        for io,jo in [(0,+1),(0,-1),(+1,0),(-1,0)]:
            ni,nj = i+io,j+jo
            if maze[ni][nj]: stack.append((ni,nj))
    return False

Counting nucleotides

Rosalind is a nice web site I recently discovered. Like 99 scala problems, Topcoder or similar sites it poses little programming challenges but with a special focus on bioinformatics.

The first and simplest programming challenges is to count the nucleotides within a DNA sequence (see here). For instance, given a sequence "AGCTTTTCATTCTGACTGC", what are the frequencies of the four nucleotides A,C,T,G?

It is a trivial problem, I used to solve with a dictionary but reading the problem description, which requires an ordered output of the frequencies, I realized that there is an even more elegant solution.

An efficient but ugly solution would be to count the individual nucleotide occurrences within a loop, e.g.

s = "AGCTTTTCATTCTGACTGC"
A,C,T,G = 0,0,0,0
for e in s:
    if e=='A': A +=1
    if e=='C': C +=1
    if e=='G': T +=1
    if e=='T': G +=1            

print A,C,T,G

The repetition of the four ifs is not nice and could be replaced as follows

s = "AGCTTTTCATTCTGACTGC"
idx = {'A':0, 'C':1, 'G':2, 'T':3}
counts = [0]*4
for e in s:
    counts[idx[e]] += 1           
print counts 

However, it remains a rather ugly solution and I find the following line much more elegant

print map(s.count, "ACTG") 

The good news is, it is even pretty efficient. The sequence gets parsed 4 times and the complexity is therefore O(4n) = O(n), which is not worse than the complexity of the dictionary approach:

counts = dict()
for e in s:
    counts[e] = counts.get(e,0)+1
print counts  

print map(s.count, "ACTG") 

Just for fun and comparison the same in Scala. First the dictionary approach

def counts = s.foldLeft(Map[Char,Int]() withDefaultValue 0){
  (m,c) => m + ((c,m(c)+1))
}
println(counts)

which is rather convoluted. Less efficient but much more readable is the version that uses "groupBy":

val s = "AGCTTTTCATTCTGACTGC"
def counts = s.groupBy(identity).mapValues(_.size)
println(counts)

and mapping the count function onto the nucleotides again leads to the most elegant code

val s = "AGCTTTTCATTCTGACTGC"
val counts = "ACTG".map(ch => s.count(_==ch))
println(counts)

A slightly different version I saw on Rosalind

val s = "AGCTTTTCATTCTGACTGC"
println( "ACGT" map (s count _.==) mkString " " )

Permutations in Python and Scala

A permutations of a collection of objects describes all possible, ordered arrangements of those objects. For instance, all permutations of (1,2,3) are (1,2,3), (1,3,2), (2,1,3), (2,3,1), (3,1,2), and (3,2,1). The following Python code generates all permutation of the elements in the given list:

def permutations(ls):
    if len(ls) == 1: yield ls
    for i in xrange(len(ls)):
        for p in permutations(ls[:i]+ls[i+1:]):
            yield [ls[i]]+p

Here a usage example:

ls = [1,2,3]            
for p in permutations(ls):
    print p

[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]

Now, let's compute permutations in Scala (2.9):

def permutations[T](ls: List[T]): List[List[T]] = ls match {
  case List() => List(List())
  case _ => for(e <- ls; r <- permutations(ls filterNot(_==e))) yield e::r
}

In contrast to the Python version the Scala version above is not a generator and creates all permutations in memory. While the "yield" keyword in Python defines a generator in Scala it is a way to describe a list comprehension. To create a generator version Iterators can be used:

def permutations[T](ls: List[T]): Iterator[List[T]] = ls match {
  case List() => Iterator(List())
  case _ => for(e <- ls.iterator; r <- permutations(ls.filterNot(_==e))) yield e::r
}

And a more compact but slightly less readable version:

def perms[T](ls:List[T]): Iterator[List[T]] =
    if(ls.isEmpty) Iterator(sl) else for (h <- ls.iterator; t <- perms(ls.filterNot(_==h))) yield h::t

The time complexity is n factorial (O(n!)), which is easy to see. Let's say there is one element, then there is only one permutations. If there are two elements, the there is two possible permutations. With three elements let's hold one, that there is two left, which gives us two permutations times three. Obviously the sequence is 1*2*3*...*n, which is n!. BTW, this can be beautifully written in Scala:

def fac(n:Int) = (1 to n) reduceLeft (_ * _)

or more cryptically-compact using left fold instead of reduce:

def fac(n:Int) = (1/:(2 to n))(_*_)

Fun with cartesian products

Cartesian product is a fancy term for all possible, ordered combinations of things taken from 2 or more collections of things (There is an exact mathematical definition for it but I don't bother you with it.) The Cartesian product is frequently used to define a grid of coordinates. For instance the following code (Python 2.7) prints out the Cartesian product of the two lists xs and ys:

xs = [1,2]
ys = [1,2,3]
for x in xs:
    for y in ys:
        print x,y

And here is the output:

1 1
1 2
1 3
2 1
2 2
2 3

Simple enough. Just two nested loops for two dimensions. This can be written very nicely as a list comprehension

product = [(x,y) for x in xs for y in ys]

which gives you

[(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3)]

However, it gets more challenging if you want to generalize to arbitrary numbers of dimensions. An iterative implementation is rather ugly but the recursive version is beautiful (though hard to remember):

def product(ls):
    if not ls: return [[]]
    return [[e]+p for e in ls[0] for p in product(ls[1:])]

And here a usage example:

ls = [[1,2], ['a','b','c'], ['A','B']]
for p in product(ls): 
    print p

[1, 'a', 'A']
[1, 'a', 'B']
...
[2, 'c', 'A']
[2, 'c', 'B']

How does it work? If the input list ls is empty a list of empty lists is returned ([[]]). Otherwise a list comprehension is performed, where the elements ([e]+p) are the elements (e) of first list (ls[0]) within the input list, prepended to the cartesian products (p) of the remainder (ls[1:]) of the input list. While really nice, it has the disadvantage that the entire Cartesian product is created in memory. Frequently, only the individual products are needed and it would be much more efficient, if one could iterate over them. The following implementation show the generator version:

def product(ls):
    if not ls: return iter( [[]] ) 
    return ([e]+p for e in ls[0] for p in product(ls[1:]))

Often only a more constraint version of this very general implementation is need. Here, for instance, a binary counter with n digits:

def bin_counter(n):
    if not n: return iter( [[]] ) 
    return ([e]+p for e in [0,1] for p in bin_counter(n-1)

or a bit more general, a counter with arbitrary digit ranges:

def counter(ls):
    if not ls: return iter( [[]] ) 
    return ([e]+p for e in xrange(ls[0]) for p in counter(ls[1:]))

ls = [2,3]                
for p in counter(ls): print p

[0, 0]
[0, 1]
[0, 2]
[1, 0]
[1, 1]
[1, 2]

Finally, let's do the same in Scala (2.9). First the general version of the Cartesian Product:

def product[T](ls: List[List[T]]): List[List[T]] = ls match {
  case Nil => List(List())  
  case h::t => for(e <- h; r <- product(t)) yield e::r
}  
  
val ls = List(List(1,2),List(1,2,3))
println( product(ls) )

The type declarations make it a bit more wordy but type-safe and it still looks elegant. The main difference to the Python version is, that the dimensions need to be of the same type (e.g. all Ints or Char but not mixed). The counter implementation is very similar. Apart from the different type (List[Int]) the main difference is the usage of List.range(0,h). Alternatively, e <-(0 to h).toList, could have been used but the required conversion to a list would make it ugly.

def counter(ls: List[Int]): List[List[Int]] = ls match {
  case Nil => List(List())  
  case h::t => for(e <- List.range(0,h); r <- counter(t)) yield e::r
}

val ls = List(2,3)
println( counter(ls) )

None of the two Scala implementations above is a generator though. Don't get deceived by the "yield " keyword, which defines a list comprehension in Scala but not a generator as in Python. However, creating a version that generates products or counts on demand using an iterator is easy; here for the counter:

def counter(ls: List[Int]): Iterator[List[Int]] = ls match {
  case Nil => Iterator(List())  
  case h::t => for(e <- Iterator.range(0,h); r <- counter(t)) yield e::r
}

val ls = List(2,3)
counter(ls).foreach(println)

That's it.