Edit distance in Scala

How handy.

I was just looking for a way to highlight differences between 2 strings for my Behavior-Driven Development library, specs. And reading the excellent Algorithm Design Manual. Intuitively, I was thinking that there was a better way to show string differences than the one used in jUnit:

Expected kit<t>en but was kit<ch>en
Expected <skate> but was <kite>

The first message shows that <t> has been replaced by <ch>, while the second message says that everything is different! There must be a way to find the minimal set of operations necessary to transform one string to another, right?

The Levenshtein distance

There is indeed. The "Edit" distance, also called "Levenshtein" distance, computes exactly this, the minimal number of insertions, deletions and substitutions required to transform "World" into "Peace" (5, they're very far apart,...).

The computation is usually done by:

1. examining the first 2 or last 2 letters of each string
2. assuming the eventual transformation which may occur for those 2 letters:

• they're equal: "h..." and "h...." --> distance = 0
  
• the first one has been removed: "ah..." and "h..." --> distance = 1
• the second one has been inserted: "h..." and "ah..." --> distance = 1
  
• they've been substituted: "ha..." and "ga..." --> distance = 1
  

3. saying that the final distance is the minimum of the distance when choosing one the

4 possibilities + the distance resulting from the consequences of that choice

You can notice that many variations are possible: the distance could be different depending on the letters being added or removed, there could be more allowed operations such as the swapping of 2 letters (distance 1 instead of 2 substitutions of distance 2),...

Then the final result can either be constructed incrementally or recursively. Incrementally, for "skate" and "kite" you will compute costs from any small prefix to another prefix:

• "s" to "k" --> distance 1
  
• "s" to "ki" --> distance 2
  
  
• "sk" to "ki" --> distance 2
• "ska" to "kit" --> distance 3
  

Then you increase the prefix size and you reuse the distance numbers found for smaller prefixes. The different results can be stored in the following matrix:

s k a t e

k 1 1 2 3 4
i 2 2 2 3 4
t 3 3 3 2 3
e 4 4 4 4 2

Recursively, you do the inverse and you establish that the distance between 2 strings can be computed from knowing the distance between smaller prefixes and you travel the matrix to its upper left corner.

This is a very good example of Dynamic Programming, where the optimal quantity (the distance at [i, j] in the matrix) you're looking for only depends on:

• the optimal quantity for a subset of the data (the distance at [i-1, j-1], [i-1, j], [i, j-1])
  
• the resulting state (the substrings defined by a position [i-x, j-x] in the matrix)
• not how you got to that resulting state

The EditDistance trait

For the record, I will write here the code I'm using for the implementation, which is constructing the solution incrementally, as opposed to the recursive solution presented on Tony's blog, here.

The few interesting things to notice about this code are:

• the algorithm itself ("initializing the matrix"), which is very straightforward
• the minimum function which is a bit strange because it is not comparing all alternatives. The cost for an insertion is only computed if the cost for a suppression is not better than the cost for a substitution. I would like to work out a formal proof on why this is correct by my intuition tells me that a substitution is generally a "good" operation. It has the same cost as an insertion or a suppression but processes 2 letters at the time. So if we find a better cost using a suppression, it is going to be the best
  
• the matrix traversal to retrieve one of the possible paths showing the letters which needs to be transformed for each string: not so fun code with a lot of corner cases
• the "separators" feature which allows to select the separators of your choice to display the string differences


trait EditDistance {
/**
* Class encapsulating the functions related to the edit distance of 2 strings
*/
case class EditMatrix(s1: String, s2: String) {
/* matrix containing the edit distance for any prefix of s1 and s2:
matrix(i)(j) = edit distance(s1[0..i], s[0..j])*/
val matrix = new Array[Array[int]](s1.length + 1, s2.length + 1)

/* initializing the matrix */
for (i <- 0 to s1.length;
j <- 0 to s2.length) {
if (i == 0) matrix(i)(j) = j // j insertions
else if (j == 0) matrix(i)(j) = i // i suppressions
else matrix(i)(j) = min(matrix(i - 1)(j) + 1, // suppression
matrix(i - 1)(j - 1) + (if (s1(i - 1) == s2(j - 1)) 0
else 1), // substitution
matrix(i)(j - 1) + 1) // insertion

}
/** @return the edit distance between 2 strings */
def distance = matrix(s1.length)(s2.length)

/** prints the edit matrix of 2 strings */
def print = {
for (i <- 0 to s1.length) {
def row = for (j <- 0 to s2.length) yield matrix(i)(j)
println(row.mkString("|"))
}
this
}

/** @return a (String, String) displaying the differences between each
input strings. The used separators are parenthesis: '(' and ')'*/
def showDistance: (String, String) = showDistance("()")

/**
* @param sep separators used to hightlight differences. If sep is empty,
* then no separator is used. If sep contains one character, it is
* taken as the unique separator. If sep contains 2 or more characters,
* the first 2 characters are taken as opening separator and closing
* separator.
*
* @return a (String, String) displaying the differences between each
* input strings. The used separators are specified by the caller.
*/
def showDistance(sep: String) = {
val (firstSeparator, secondSeparator) = separators(sep)
def modify(s: String, c: Char): String = modifyString(s, c.toString)
def modifyString(s: String, mod: String): String = {
(firstSeparator + mod + secondSeparator + s).
removeAll(secondSeparator + firstSeparator)
}
def findOperations(dist: Int, i: Int, j:Int, s1mod: String, s2mod: String): (String, String) = {
if (i == 0 && j == 0) {
("", "")
}
else if (i == 1 && j == 1) {
if (dist == 0) (s1(0) + s1mod, s2(0) + s2mod)
else (modify(s1mod, s1(0)), modify(s2mod, s2(0)))
}
else if (j < 1) (modifyString(s1mod, s1.slice(0, i)), s2mod)
else if (i < 1) (s1mod, modifyString(s2mod, s2.slice(0, j)))
else {
val (suppr, subst, ins) = (matrix(i - 1)(j), matrix(i - 1)(j - 1),
matrix(i)(j - 1))
if (suppr < subst)
findOperations(suppr, i - 1, j, modify(s1mod, s1(i - 1)), s2mod)
else if (ins < subst)
findOperations(ins, i, j - 1, s1mod, modify(s2mod, s2(j - 1)))
else if (subst < dist)
findOperations(subst, i - 1, j - 1, modify(s1mod, s1(i - 1)),
modify(s2mod, s2(j - 1)))
else
findOperations(subst, i - 1, j - 1, s1(i - 1) + s1mod, s2(j - 1) + s2mod)
}
}
findOperations(distance, s1.length, s2.length, "", "")
}
def min(suppr: Int, subst: Int, ins: =>Int) = {
if(suppr < subst) suppr
else if (ins < subst) ins
else subst
}
}
def editDistance(s1: String, s2: String): Int = EditMatrix(s1, s2).distance
def showMatrix(s1: String, s2: String) = EditMatrix(s1, s2).print
def showDistance(s1: String, s2: String) = EditMatrix(s1, s2).showDistance
def showDistance(s1: String, s2: String, sep: String) = EditMatrix(s1, s2).showDistance(sep)

private def separators(s: String) = (firstSeparator(s), secondSeparator(s))
private def firstSeparator(s: String) = if (s.isEmpty) "" else s(0).toString
private def secondSeparator(s: String) = {
if (s.size < 2) firstSeparator(s) else s(1).toString
}
}

Conclusion

What's the conclusion? Computer science is useful of course, but you only recognize it once you know it!

Algorithmic panic - take 2

Shame, shame, shame.

The Blog Writing Rule

Writing at least one blog entry per month. That should be feasible, no? Unfortunately, last Friday, I was already more than 10 days late, so I had to do something!

And, when pressing the "Post" button, I had this awkward feeling that something was not really right with my "solution". Not really on the "speed" side because I knew that my list accesses were sloppy, but even on the "correctness" side. But I was past the date, damn it!

Fortunately, someone's watching

Ok, enough justifications. I was in fact very pleased to get so many insightful comments, be it on my blog, or on reddit. It's a pleasure to see that so many people care for The Truth on the Internet ( ;-) )

So I tried my best to review things up and to derive the maximum number of lessons out of theses errors. I'll start with that, then I'll copy the revised version of the algorithm.

Lessons learned

1. Check your definitions. The median definition is very straightforward,... if you have an odd number of elements. For an even number, you can elaborate several definitions, such as that one. For the code below, I choose to take define the median as a pair of values, which is both "middle" elements if the list is even. Pairs of values are then conveyed across all computations, so that it's always possible to eventually take the mean value if you feel like it


2. Don't trust test generation. I promise my tests were green when I posted! However, the generation parameters were very bad: the number of tests was too small, the number and range of elements in the lists also. The solution here was tested on larger samples (up to 8000 tests), with larger lists (of same size or different sizes)


3. Check your invariants. In this kind of "Divide and Conquer" approach, I should make sure that I have indeed the same problem, on a smaller scale, after dividing it. The code below shows that it is possible to take the same number of elements on each list, so that the median of the resulting lists is unchanged


4. Know your data structures. I knew that the most basic list was not necessarily efficient but I must say that I didn't check anything. And, in fact, when Jorge pointed out that List had an O(n) access for size, I was quite surprised. The other interesting thing is that optimizing an operation can be totally feasible but quite tricky to implement (see the medianOfSortedListAndOneElement function below). And of course, takeWhile(_) was pretty bad!

The algorithm

First of all, lets start by the definition of the median as a couple of values, which may be the same if the list has an odd number of elements:

case class ExtendedSeq(list: Seq[Int]) {
def median: (Int, Int) = {
if (list.isOdd)
(list(list.size / 2), list(list.size / 2))
else
(list(list.size / 2 - 1), list(list.size / 2))
}
def isOdd = list.size % 2 != 0
}

If elements access and size are O(1) as in the RandomAccessSeq collection, that should be ok in terms of complexity (thanks Jorge).

Here is the heart of the algorithm:

def median(l1: Seq[Int], l2:Seq[Int]): (Int, Int) = {
if (l1.isEmpty) l2.median
else if (l2.isEmpty) l1.median
else if (l1.size == 1) medianOfSortedListPlusOneElement(l2, l1(0))
else if (l2.size == 1) medianOfSortedListPlusOneElement(l1, l2(0))
else {
val (m1, m2) = (l1.median, l2.median)
m1.intersection(m2) match {
case Some(m) => m
case None => if (m1 < m2) median(trimLists(l1, l2))
else median(trimLists(l2, l1))
}
}
}

The code above starts by checking some special cases: one empty list, a list + one element, then 2 "regular" median values.

The first observation is that, if the median values (which are pairs of values) have a non-empty intersection, this intersection will be the median of both sorted lists.

This is the not-so readable code for checking if there is an overlap between 2 pairs:

case class ExtendedPair(pair1: (Int, Int)) {
def intersection(pair2: (Int, Int)) = {
if (pair2._1 <= pair1._1 && pair1._2 <= pair2._2) Some(pair1)
else if (pair1._1 <= pair2._1 && pair2._2 <= pair1._2) Some(pair2)
else if (pair1._1 <= pair2._1 && pair2._1 <= pair1._2 && pair1._2 <= pair2._2)
Some((pair2._1, pair1._2))
else if (pair2._1 <= pair1._1 && pair1._1 <= pair2._2 && pair2._2 <= pair1._2)
Some((pair1._1, pair2._2))
else None
}
}

Finally, if there is no overlap, we need to remove elements from both lists so that the median of the remaining 2 lists will be the same as the median of the original lists:

def trimLists(l1: Seq[int], l2: Seq[Int]) = {
val removableElements = if (l1.size < l2.size) l1.size / 2 else l2.size / 2
(l1.drop(removableElements), l2.take(l2.size - removableElements))
}

This time, I removed the same number of elements below the lower median and above the higher one (but still, a written mathematical proof would be more satisfying than "it's obvious that,...").

On the complexity side, I still have to check that I can find a data structure (a Seq) providing drop and take operations in O(1). From what I read of the implementation of RandomAccessSeq, that should be the case. RandomAccessSeq is using indexes to provide O(1) "slicing" operations as suggested by Mikko in the comments of the previous post.

Brutal force

And finally, I can't resist showing you the horror of finding the median of a sorted list and a supplementary element in O(1). There may be a subtle trick to do that (please say, if you find it!), but I used the brutal force:

def medianOfSortedListPlusOneElement(list: Seq[Int], element: Int) = {
if (list.size == 1){
if (element >= list(0))
(list(0), element)
else
(element, list(0))
}
else if (list.isOdd) {
if (element >= list(list.size / 2)) {
if (element <= list(list.size / 2 + 1))
(list(list.size / 2), element)
else
(list(list.size / 2), list(list.size / 2 + 1))
} else {
if (element >= list(list.size / 2 - 1))
(element, list(list.size / 2))
else
(list(list.size / 2 - 1), list(list.size / 2))
}
} else {
if (element >= list(list.size / 2 - 1)) {
if (element <= list(list.size / 2))
(element, element)
else
(list(list.size / 2), list(list.size / 2))
} else {
(list(list.size / 2 - 1), list(list.size / 2 - 1))
}
}
}

Not pretty, uh? Like the final approach for the 4-colors theorem.

Measuring the complexity

On my way to this post, I wanted to implement something which would store the test data and say if the curve looks like an O(log(n)), an O(nlog(n)),... In fact it doesn't seem so easy to do without any human validation and tweaking of the test data size until the complexity really shows. I'm thinking that using an utility to graph the performances and compare it with known complexities would be a good first step.

But maybe driving the test generation so as to get a good confidence on the possible complexity of the algorithm would be even better. I'll try to have a look at that,... Maybe for my next post?

Thanks for all the comments, I hope I didn't add more and more silliness to my previous errors.

PS: I liked a lot Tom's comments about the 2-keys indexing problem. Indeed, the problem was underspecified, in terms of access frequency, memory, etc,... And a log(n) solution would certainly be ok in my specific case.