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.

Reading FASTA files with Scala

FASTA is a simple file format to store nucleotide or peptide sequences. A sequence in FASTA format starts with a greater-than sign followed by the sequence name (and possibly other information) and then several lines of nucleoide or amino acid letters. A multi-FASTA file contains multiple FASTA sequences. Here an example taken from wikipedia:

>SEQUENCE_1
MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG
LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK
IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL
MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL
>SEQUENCE_2
SATVSEINSETDFVAKNDQFIALTKDTTAHIQSNSLQSVEELHSSTINGVKFEEYLKSQI
ATIGENLVVRRFATLKAGANGVVNGYIHTNGRVGVVIAAACDSAEVASKSRDLLRQICMH

It is a common task to read FASTA sequences from a file and I was aiming for an elegant/short implementation in Scala. Here is a recursive version I like:

  
  def readFasta(filename: String) = {
     import scala.io.Source
     def parse(lines: Iterator[String]): List[(String,String)] = {
       if(lines.isEmpty) return List()
       val name = lines.next.drop(1)
       val (seq,rest) = lines.span(_(0)!='>')
       (name, seq.mkString)::parse(rest)
     }
     parse(Source.fromFile(filename).getLines())
  }

It defines an internal function parse that takes an iterator over the lines in the file (generated by Source.fromFile(filename).getLines() and returns a list of tuples, where each tuple contains the sequence name and the string of nucleotide or amino acid letters.

The recursion stops when the iterator is empty (lines.isEmpty) otherwise a line is read, which must be the name and the first letter is dropped to get rid of the greater-than sign (lines.next.drop(1)). After the line with the name span(_(0)!='>') is used to take the following lines until the next sequence name is found and store them in seq. _(0)!='>' tests whether the first letter of the line is the greater sign. The remaining lines are stored in rest.

Actually, "storing" is the wrong word, since Scala beautifully returns two new iterators when span is called on an iterator. It is worth noting, that nowhere a list of the read lines is stored. It is turtles/iterators all the way through except for the last line ((name, seq.mkString)::parse(rest)), where the return list is created. (name, seq.mkString) builds a tuple with the sequence name and the the sequence letters - the latter is constructed by gluing all sequence lines together with mkString. This tuple is the head of the return list and the tail is build recursively by calling parse for the remaining lines.

Let's add three minor improvements. It will make the code less elegant but more robust. Firstly, lines could have leading or trailing white spaces. They shouldn't but hey, you know how it is. To play it save we add a to the part where the lines are read from the file (last line in code). Secondly, there could be empty lines that we filter out via filterNot(_.isEmpty()).

Thirdly, the line with the sequence name frequently contains additional annotation that is typically separated by the pipe "|" symbol. But there is really no solid standard and it could be in any format - and it will. We therefore define a regular expression to extract the name from the header line and tweak it if necessary. The expression """>(.+?)(\\|.+)?""".r has two groups. The first captures the name the second captures any annotation starting with the pipe "|" symbol, which needs to be escaped (because it functions as OR in Regex language). The final code looks like this:

  
  def readFasta(filename:String) = { 
     import scala.io.Source
     val header = """>(.+)(\\|.+)?""".r
     def parse(lines: Iterator[String]): List[(String,String)] = {
       if(lines.isEmpty) return List()
       val header(name,annotation) = lines.next
       val (seq,rest) = lines.span(_(0)!='>')
       (name, seq.mkString)::parse(rest)
     }
     parse(Source.fromFile(filename).getLines().map(_.trim).filterNot(_.isEmpty()))
  }    

While elegant the above solution has the problem of being recursive (specifically not even tail recursive), which will cause a stack overflow for files with large numbers of sequences. The following iterative implementation is ugly but I haven't found any really nice alternative. Take it as a temporary implementation ;)

  
  def readFasta(filename:String) = {
    val header = """>(.+)(\\|.+)?""".r
    var lines = Source.fromFile(filename).getLines().filterNot(_.isEmpty())
    var sequences = List[(String,String)]()
    while(lines.hasNext) {
      val header(name,annotation) = lines.next
      val (seq,rest) = lines.span(_(0)!='>')
      sequences = (name, seq.mkString)::sequences
      lines = rest
    }
    sequences.reverse
  }  

Also interesting is this solution that uses parser combinators to read FASTA sequences.

Alexandre Masselot has posted an interesting problem (see below). What if you want an iterator over the sequences instead of a list? We can easily modify the code above to achieve this by creating and returning an iterator when calling readFasta. An iterator needs a hasNext method, which essentially replaces the condition of the while loop used above and a next method, which contains the body of the while loop:

  
  def readFasta(filename: String) = 
    new Iterator[(String, Iterator[Char])] {
    val header = """>(.+)(\|.+)?""".r
    var lines = Source.fromFile(filename).getLines().filterNot(_.isEmpty())
    def hasNext = lines.hasNext
    def next = {
      val header(name, annotation) = lines.next
      val (seq, rest) = lines.span(_(0) != '>')
      lines = rest
      (name, seq.flatMap(_.iterator))
    }
  }

If you look at the code closely, you will notice on more change. I have replaced (name, seq.mkString) by (name, seq.flatMap(_.iterator)). This is where Scala shines. Instead of concatenating the lines to a sequence via mkString we can simply call flatMap(_.iterator) to convert lines to iterators over line characters and then concatenate the iterators via flatMap.

In practice, you will find an iterator over sequences frequently useful, but an iterator over sequence characters less so. Many computational evaluations of sequence data will require random access. If memory consumption is an issue, compressing the sequence is probably a better way. Also in many cases a set of sequence data will be related and storing only the differences to a reference sequence (e.g. the human reference genome) can be a very efficient representation of large sequence data.

Java to Scala: Motivation

If you are a Java programmer and haven't heard about Scala already you might wonder if it is worth the effort learning. Here I show two examples that motivated me to switch from Java to Scala.

This first one I stole from somewhere ― I just can't remember where from, sigh ― and modified it slightly. It is a simple implementation of a class that stores the contact information (name and phone number) of a person.

public class Contact implements Serializable {
    private final String name;
    private final String phone;

    public Contact(String name, String phone) {
        this.name = name;
        this.phone = phone;
    }

    public String getName() {
        return name;
    }

    public String getPhone() {
        return phone;
    }

    public Contact withName(String name) {
        return new Contact(name, phone);
    }

    public Contact withPhone(String phone) {
        return new Contact(name, phone);
    }

    public boolean equals(Object o) {
        if (this == o)  return true;
        if (o == null || getClass() != o.getClass()) return false;
        
        Contact p = (Contact) o;
        if (name != null ? !name.equals(p.name) : p.name != null) {
            return false;
        }
        if (phone != null ? !phone.equals(p.phone) : p.phone != null) {
            return false;
        }
        return true;
    }

    public int hashCode() {
        int result = name != null ? name.hashCode() : 0;
        result = 31 * result + (phone != null ? phone.hashCode() : 0);
        return result;
    }

    public String toString() {
        return "Contact(" + name + "," + phone + ")";
    }
} 

There are no surprises or difficulties in it. Just a lot of boilerplate code. Any Java programmer will have written hundreds of classes similar to this one - most of them a lot more complex. With the Contact class one can create new contacts from scratch or use one of the copy constructors as shown below:

Contact stefanOffice = new Contact("Stefan", "62606");
Contact johnOffice   = new Contact("John", stefan.getPhone());
Contact stefanHome   = stefanOffice.withPhone("401251497");

These 20 plus lines of the the Contact class in Java boil down to a single line in Scala with essentially identical functionality:

case class Contact(name: String, phone: String)


val stefanOffice = Contact("Stefan", "62606")
val johnOffice   = Contact("John", stefan.phone)
val stefanHome   = stefanOffice.copy(phone = "401251497")

A case class in Scala provides default implementations for setters, getters, equality tests and copy constructors. In many cases they are sufficient and if not they can easily be overridden.

While this example is a bit extreme with respect to the amount of code reduction (typically the factor is between two and five) it is not an unrealistic example. But I also wanted to demonstrate the benefits of Scala on a slightly more complex problem.

We are going to look at some code that generates an index of over keywords. For instance, let's say we have the keywords "Apple", "Ananas", "Mango", "Banana", "Blackberry" and we would like to have an alphabetically sorted index that lists the keywords (again alphabetically sorted) starting with the same letter:

A: Ananas, Apple
B: Banana, Blackberry
M: Mango

Not very challenging; should be easy. Well, here comes the Java implementation and as you can see it is far from elegant:

import java.util.*;

public class Index {    
  public static void main(String[] args) {
    List<String> keywords = 
       Arrays.asList("Apple", "Ananas", "Mango", "Banana", "Blackberry"); 
    Map<Character, List<String>> ch2words = new TreeMap<>(); 
    for(String keyword : keywords) {   
      char firstChar = keyword.charAt(0);     
      if(!ch2words.containsKey(firstChar))      
        ch2words.put(firstChar, new ArrayList<String>());       
      ch2words.get(firstChar).add(keyword); 
    } 
    
    for(List<String> list : ch2words.values()) 
      Collections.sort(list); 
      
    StringBuilder str = new StringBuilder();  
    for(Character firstChar : ch2words.keySet()) {
      str.append(firstChar+": ");
      Iterator<String> words = ch2words.get(firstChar).iterator(); 
      while(words.hasNext()) {
        str.append(words.next());  
        str.append(words.hasNext() ? ", " : "\n");  
      }
    }
    System.out.print(str);      
  }
}

We take the keywords and put them into a TreeMap to map the first letter to the list of words. I chose a TreeMap to have the map entries automatically sorted according to the first character. We then need to sort the keywords in each of the entries of the TreeMap and finally use a StringBuilder to write the formatted index to a String.

Now I am going to write some bad and truly disgusting Scala code that essentially is a one-to-one mapping from Java to Scala. While far from optimal it shows what already be gained without using any novel concepts such as functional programming or pattern matching.

import scala.collection.mutable.{Map,StringBuilder}

val keywords = List("Apple", "Ananas", "Mango", "Banana", "Blackberry")
val ch2words = Map[Char, List[String]]()
for(keyword <- keywords) {
  val firstChar = keyword(0)
  if(!ch2words.contains(firstChar))      
    ch2words(firstChar) = List()       
  ch2words(firstChar) = keyword::ch2words(firstChar)
}
val str = new StringBuilder()
for(ch <- ch2words.keys.toSeq.sorted) {
  str ++= ch2words(ch).sorted.mkString(ch+": ", ", ", "\n")      
}  
println(str)

Largely due to type inference in Scala the code becomes less cluttered and the logical flow is easier to recognize. Also the confusing distinction between primitive data types such as char and their corresponding objects (Character) is gone.

Now let's do it again and see how a solution that takes advantage of Scala features could look like:

val keywords = List("Apple", "Ananas", "Mango", "Banana", "Blackberry")
val ch2words = keywords.groupBy(_(0)).mapValues(_.sorted).toList.sortBy(_._1)
val lines = ch2words.map{case (ch,words) => ch+": "+words.mkString(", ")}
print(lines.mkString("\n"))

First of all it is much shorter. But in addition, it is also more readable provided an understanding of some basic concepts in Scala has been acquired. Let's have a closer lock at the code.

The first line is trivial. It creates a list of all the keywords and that's it. The second line is interesting. First the keywords get grouped by their first character (keywords.groupBy(_(0))). That produces a map of letters to keyword lists:

Map(M -> List(Mango), A -> List(Apple, Ananas), B -> List(Banana, Blackberry))

As you can see neither the letters nor the lists of keyword are sorted. We take the values of the map (those are the lists of keywords) and sort them. mapValues(_.sorted) does the job.

Map(M -> List(Mango), A -> List(Ananas, Apple), B -> List(Banana, Blackberry))

To sort the map according to letters we convert the map to a list, which returns a list of tuples, containing letters and their keywords, and sort the list according to the first tuple component ― the letter (toList.sortBy(_._1)).

List((A,List(Ananas, Apple)), (B,List(Banana, Blackberry)), (M,List(Mango)))

Now the data structure is sorted and we simply need to create a formatted string. Pattern matching is used to extract the letters and the corresponding keyword lists (case (ch,words)) and the lines of the index are constructed by mapping (ch,words) tuples to strings. Finally, all lines are glued together via mkString("\n") and printed out.

A: Ananas, Apple
B: Banana, Blackberry
M: Mango

To conclude, in almost all cases Scala allows one to write shorter code than in Java; and it is frequently even more readable. But it also makes it possible to write code that causes a PERL coder to blush. To quote Spider-man: "With great powers come great responsibilities!".

An example class in Scala

This is just a little example class to remind me of the syntax and some of the nifty features in Scala that I tend to forget. I essentially take the example of the Rational class from the excellent book "Programming in Scala" by Odersky et al. and modify it slightly.

Let's start very simple. We call the class "Frac" for fraction instead of Rational and its constructor takes the numerator and denominator of the fraction. The val specifier ensures that the class is immutable. We also override the toString method to get a nice print out.

class Frac(val num:Int, val den:Int) {
  override def toString = "%d/%d" format (num,den)
}

With this class at hand we can now create a fraction and print it and its components:

val a = new Frac(1,2)
println(a)
println(a.num)

1/2
1

Now let's add a companion object Frac that allows us to omit the new keyword when creating a fraction. We also want to ensure that the denominator is greater than zero and Scala provides the handy require function to test pre-conditions. Furthermore, We add a toDouble variable that contains the floating point value of the fraction. Note that we could also have said def toDouble = .. instead of val toDouble = ... . It is a trade off between speed and memory consumption. The nice thing about Scala is that the call interface remains the same and the implementation could be changed without any dramas.

class Frac(val num:Int, val den:Int) {
  require(den > 0)
  val toDouble = num/den.toDouble 
  override def toString = "%d/%d" format (num,den)
}

object Frac {
  def apply(num:Int, den:Int) = new Frac(num,den)
}  

println(Frac(1,2))
println(Frac(1,2).toDouble)

1/2
0.5

So far so good. But one problem with the current implementation is that comparisons such as Frac(1,2) == Frac(1,2) return false. We need to override the equals() method and consequently the method hashCode should also be overridden:

class Frac(val num:Int, val den:Int) {
  ...
  override def equals(other:Any) = other match {
    case that:Frac => (this eq that) || 
       (that.num == this.num && this.den == that.den)
    case _ => false
  }

  override def hashCode = 13*(num+13*den)
  ...
}

Apart from testing for equality it would be convenient if fractions could be ordered. The current class implementation does not allow to say Frac(1,3) < Frac(1,2), for instance. But implementing the compare method of the Ordered trait does the trick.

class Frac(val num:Int, val den:Int) extends Ordered[Frac] {
  ...
  def compare(that:Frac):Int = 
    (this.num*that.den) - (that.num*this.den)     
  ...
}

Operator overloading is another feature of Scala and useful when applied appropriately. For a fraction class it certainly makes sense. To keep things simple the following code is limited to the multiplication of fractions with fractions and fractions with scalars.

class Frac(val num:Int, val den:Int)  {
  ...
  def *(that:Frac):Frac = Frac(this.num*that.num, this.den*that.den)
  def *(c:Int):Frac = Frac(num*c, den)   
  ...
}

This enables us to compute Frac(1,3) * Frac(1,2) and Frac(1,3) * 2 but we still cannot calculate 2 * Frac(1,3) because there is no multiplication method for Int that takes a fraction. For that we need an implicit function, preferably defined within the companion object:

object Frac {
  def apply(num:Int, den:Int) = new Frac(num,den)
  implicit def int2Frac(num:Int):Frac = Frac(num,1)
}  

Depending on scope it might be necessary to import the implicit function of the companion object but then things work as expected:

import Frac.int2Frac
println(Frac(1,3) * Frac(1,2))    // prints 1/6
println(Frac(1,3) * 2)            // prints 2/3
println(2 * Frac(1,3))            // prints 2/3

When multiplying a fraction with an integer scalar we would like to get a fraction back. However, when multiplying with a floating point number we need a floating point number as a result. The implementation is similar; just the direction is different.

class Frac(val num:Int, val den:Int) {
  ...
  val toDouble = num/den.toDouble
  def *(c:Int):Frac = Frac(num*c, den)   // Frac * Int => Frac
  def *(c:Double):Double = toDouble*c    // Frac * Double => Double
  ...
}

object Frac {
  def apply(num:Int, den:Int) = new Frac(num,den)
  implicit def int2Frac(num:Int):Frac = Frac(num,1)
  implicit def frac2double(frac:Frac):Double = frac.toDouble
}  

import Frac._
println(Frac(1,2) * 2.5)   // prints 1.25
println(2.5 * Frac(1,2))   // prints 1.25

Alright, that it's. Of course, there is some functionality missing such as addition, subtraction and division operations. Furthermore, fractions should be normalized to the greatest common divisor as shown in Odersky's code. But I wanted to keep this example simple. To conclude let us put together what have now and see what we can do:

class Frac(val num:Int, val den:Int) extends Ordered[Frac] {
  require(den > 0)
  val toDouble = num/den.toDouble
  def *(that:Frac):Frac = Frac(this.num*that.num, this.den*that.den)
  def *(c:Int):Frac = Frac(num*c, den)   
  def *(c:Double):Double = toDouble*c   
  def compare(that:Frac):Int = 
    (this.num*that.den) - (that.num*this.den)     
  override def equals(other:Any) = other match {
    case that:Frac => (this eq that) || 
       (that.num == this.num && this.den == that.den)
    case _ => false
  }
  override def hashCode = 13*(num+13*den)
  override def toString = "%d/%d" format (num,den)
}

object Frac {
  def apply(num:Int, den:Int) = new Frac(num,den)
  implicit def int2Frac(num:Int):Frac = Frac(num,1)
  implicit def frac2double(frac:Frac):Double = frac.toDouble
}  

import Frac._
println(Frac(1,3) == Frac(1,2))
println(Frac(1,3) < Frac(1,2))
println(Frac(1,3) * Frac(1,2))
println(Frac(1,2) * 2.5)
println(2.5 * Frac(1,2))
println(Frac(1,2) * 2)
println(2 * Frac(1,2))

Newick tree format parser in Scala

Parser combinators in Scala are nice for small, non-time-critical parsers. While typically slower than parser generators such as ANTLR or hand-written recursive descent parser (see also here), they are easier to implement and do not require external syntax or output files. If speed is an issue there is also Parboiled, a parser generator very similar in style to the parser combinators in Scala but apparently faster. I haven't tried it, however.

I needed a parser to read phylogentic tree data and Scala's parser combinators seemed just right. Without going into any details concerning phylogenetic trees, the following figure shows an example tree with a root node R, two internal nodes X and Y, and the leaf nodes A,B,C,D. Each branch has a length annotation, e.g. 0.1.

There are various formats to describe phylogentic trees but to keep it simple let's start with a format defined by the following EBNF:

  tree    ::= identifier length [subtree]
  subtree ::= "(" tree {"," tree} ")"
  length  ::= ":" floatingPointNumber 

Here a tree has a node identifier followed by the branch length (distance to parent), followed by an optional subtree. A subtree is a sequence of trees enclosed in brackets and separated by commas. The branch length is given by a colon followed by a floating point number.

Assuming branch lengths of 0.1 for all nodes (except the root node R), the tree above would be written as

R:0.0 (X:0.1 (A:0.1, B:0.1), Y:0.1 (A:0.1, B:0.1))

The given EBNF can easily be translated into a parser combinator that can read this tree description:

import scala.util.parsing.combinator._

class SimpleParser extends JavaTokenParsers {
  def tree                = ident ~ length ~ opt(subtree)
  def subtree:Parser[Any] = "(" ~> repsep(tree, ",") <~ ")"
  def length              = ":" ~> floatingPointNumber
}

Note that [subtree] translates to opt(subtree) and tree {"," tree} to repsep(tree, ","). Furthermore, subtree requires a return type, since it is a recursive function. I am lazy here and say Parser[Any], since we are going to improve on this code anyway. Let's add a method read() to read a tree description and return the parsing result:

class SimpleParser extends JavaTokenParsers {
  def tree = ident ~ length ~ opt(subtree)
  def subtree:Parser[Any] = "(" ~> repsep(tree, ",") <~ ")"
  def length = ":" ~> floatingPointNumber

  def read(text:String) = parseAll(tree, text)
}

val parser = new SimpleParser()
println( parser.read("R:0.0(X:0.1(A:0.1,B:0.1),Y:0.1(C:0.1,D:0.1))") )

The output is a bit convoluted but the structure reflects the tree structure. Note that opt() returns an Option, which is the reason for the Some and None in the output:

parsed: ((R~0.0)~Some(List(((X~0.1)~Some(List(((A~0.1)~None), ((B~0.1)~None)))), ((Y~0.1)~Some(List(((C~0.1)~None), ((D~0.1)~None)))))))

So far so good, but the current parser output is not very useful. What we want to get back from the parser is a proper data structure of nodes and their children. In the next steps we are going to refactor and extend the current code a bit. We start with a Node class that has a display() method to print a tree:

case class Node(name:String, length:Double, descendants:List[Node]) {
  def display(level:Int=0) {
    printf("%s%s:%.2f\n","  "*level,name,length)
    descendants.foreach(_.display(level+1))
  }
}

To keep the grammar separated from the rest and also to allow comments we factor out the read() method into an abstract super class that takes the node type as type parameter T:

abstract class TreeParser[T] extends JavaTokenParsers {
  val comment = """\\[.+?\\]"""

  def tree:Parser[T]

  def read(file:File):T = read(Source.fromFile(file).getLines().mkString("\n"))

  def read(text:String):T = parseAll(tree, text.replaceAll(comment,"")) match {
    case Success(result, next) => result
    case failure => throw new Exception(failure.toString)
  }
}

Comments are everything between two rectangular brackets (defined by a non-greedy, regular expression) and removed from the input text before running the parser. We also overload the read() method to read a tree from a file and tree is an abstract method to be implement by a subclass. Note that the read() method returns a type T (a node), which will be the root node of the tree.

parseAll() returns a ParseResult with sub classes Success, Failure or Error. We are only interested in Success and match against it to retrieve the parsing result (the root node of type T). In all other cases an exception is thrown, showing the parsing error.

Now we can rewrite the SimpleParser class as follows:

class SimpleParser[T](nf: (String,Double,List[T]) => T) extends TreeParser[T] {
  def tree = (ident ~ length ~ opt(subtree)) ^^ {case n~l~d => nf(n,l,d.getOrElse(Nil))}
  def subtree:Parser[List[T]] = "(" ~> repsep(tree, ",") <~ ")"
  def length = ":" ~> floatingPointNumber ^^ { _.toDouble }
}

There are two main difference to the previous implementation. First, a node factory (nf) is provided that takes a node name (String), the branch length (Double) and a list of child nodes (List[T]) to create a node of type T. This allows us to use the same parser to create trees of different node types, which is especially handy for testing. Note that the function subtree() now has a tight return type Parser[List[T]] and is not Parser[Any] anymore.

The second change is that actions have been added to the parsers. { _.toDouble } returns a Double for the parsed floating point number of the branch length and {case n~l~d => nf(n,l,d.getOrElse(Nil))} calls the node factory to create a node. If there is no subtree the list of descendants is Nil (d.getOrElse(Nil)).

If we run the parser now, we get back a tree structure of nodes to play with. I am a bit sneaky with the node factory and simply provide the apply() method of the node singleton that is automatically generated for the case class.

val parser = new SimpleParser(Node.apply)
val root = parser.read("R:0.0(X:0.1(A:0.1,B:0.1),Y:0.1(C:0.1,D:0.1))")
root.display()

R:0.00
  X:0.10
    A:0.10
    B:0.10
  Y:0.10
    C:0.10
    D:0.10

Now we are ready for the real deal: A parser for phylogenetic trees in Newick format. I took the following grammar from here and adapted it bit.

  tree        ::= subtree ";"
  descendants ::= "(" subtree {subtree ","} ")"
  subtree     ::= descendants name length | leaf
  leaf        ::= name length
  name        ::= [quoted | unquoted]
  unquoted    ::= identifier
  quoted      ::= "'" { not "'" | "''"} "'"
  length      ::= [":" floatingPointNumber]

From that we can derive the corresponding parser:

class NewickParser[T](nf: (String,Double,List[T]) => T) extends TreeParser[T] {
  def tree = subtree <~ ";"
  def descendants:Parser[List[T]] = "(" ~> repsep(subtree, ",") <~ ")"
  def subtree = descendants~name~length ^^ {case t~n~l => nf(n,l,t)} | leaf
  def leaf = name~length ^^ {case n~l => nf(n,l,Nil)}
  def name = opt(quoted | unquoted) ^^ { _.getOrElse("") }
  def unquoted = ident
  def quoted = """'([^']|'')*'""".r  ^^ { _.drop(1).dropRight(1).replace("''","'") }
  def length = opt(":" ~> floatingPointNumber) ^^ { _.getOrElse("0").toDouble }
}

The Newick format allows to leave out node names or branch length information, which can lead to rather wicked tree definitions such as

"(A,(B,C)E)F;"
(((,),));
(,B)C:0.3;

With the lengthy introduction in mind the NewickParser class is hopefully not beyond comprehension. The only difficult and admittedly ugly part is the action { _.drop(1).dropRight(1).replace("''","'") }, which simply removes the enclosing quotes from a quoted identifier and replaces pairs of single quotes ('') by single quotes. The quoted parser could be implemented more nicely as:

def quoted = "'" ~> """([^']|'')*""".r <~ "'"  ^^ { _.replace("''","'") }

but in this case leading and trailing white spaces of quoted identifiers are lost, which might be acceptable. The astute reader will notice two other limitation. Firstly, nested comments, e.g. [ a comment [ another comment]] are not permitted and identifiers cannot contain brackets. I have not found a really nice way of supporting this functionality. Either one extends the grammar to allow comments where ever white spaces can appear, which makes the grammar unwieldy. Or one overrides the white space parser by a custom parser that supports comments. Or the last option that I chose is to implement a separate comment parser to remove comments from the input text instead of the regular expression employed in TreeParser.

class CommentParser extends JavaTokenParsers {
  override val skipWhitespace = false

  /** removes comments from given text */
  def remove(text:String) = parseAll(rest, text)

  def rest = rep(quoted | comment | any) ^^ { _.mkString }
  def any = """.""".r
  def quoted = """'([^']|'')*'""".r
  def comment: Parser[String] = ("["~rep(not("]")~(comment | any))~"]") ^^ {_ => ""}
}

I found the increase in complexity not worth the additional functionality and stayed with the limited parser that does not allow nested comments.

Finally, some links I came across recently that I found useful for parser developement:
Code examples for combinator parsing
Ten grammars for simple integer arithmetic expressions

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.