Showing posts with label sorting. Show all posts
Showing posts with label sorting. Show all posts

Wednesday, December 15, 2010

Introducing Page Rank Maths

In my previous post, I started introducing Google PageRank. Now I'm starting to dig into PageRank maths. PageRank definition is not very "computation oriented".

r(Pi) = ∑ Pj ∈ Li  r(Pj)

|Pj|
(2)
There is not a clear temporal ordering of pages or any other sensible starting point. Hopefully, an iterative method may converge to the result. So we subscript the rank with a number indicating the step and express the (i)-th step in terms of the (i−1)-step:

rk(Pi) = ∑ Pj ∈ Li  rk−1(Pj)

|Pj|
(3)
Of course, we are not happy with the whole "may-converge" thing and would rather have some precise results. The first step is using some convenient notation: formula 2.2 is easy to grasp but awkward to manipulate.
The idea here is express the equations in matrix form. Let L = (eij) be the adjacency matrix of our graph. In other words:

eij =
1    if there is a link from i to j
0    otherwise
(4)
Then we consider the matrix H = ―D−1·A, where ―D = max(I, D) and D is the diagonal matrix of the degrees of the nodes. In Python we can compute H with:
import numpy as np

def h_matrix(adjacency_matrix):
    degrees = np.sum(adjacency_matrix, axis=1)
    degrees[degrees==0] = 1
    return np.dot(np.diag(1./degrees), adjacency_matrix)


# ...

The H j-th row can be interpreted as the distribution of probability for jumping from node j to each other node. If our graph is the web, H has billions of rows (and columns) but in average each page has about 10 out links. H is a very sparse matrix, indeed.
Let π be the rank-vector. Then we can reformulate the problem as:
πT = πT·H
(5)
Basically we are looking for the eigenvector associated with the eigenvalue λ1 = 1. The idea here is to use the power method. The power method is an iterative method that computes the eigenpair (λ1, π) of a matrix A. A must satisfy some basic hypotheses for the power method to work. Let {λ1,…,λn} be the set of eigenvalues of A, then it must be true that |λ1| > |λ2| ≥ … ≥ |λm|. As a consequence, λ1 ∈ \mathbbR.
Then it is possible to compute (λ1, π) with:

yn=Axn
      
xn + 1 = yn

cn
(6)
Then xn→ x and cn→λ1.
Here cn is a scaling factor; e.g. λ1, ||Anx0||, the component of yn of maximal magnitude. Now, the fact that we suggested λ1 as a possible scaling factor is not a mistake: it is proved that the actual Google matrix has λ1.
Moreover x0 ∉ R(A−λ1I), where R(M)={v|v=Ax}. However, in practice "most" starting vectors x0 are good.
If H were a usable Google Matrix, things would be very nice. H is essentially the matrix of a random walk over the graph. We could just see the whole thing as a Markov process. More on this can be found in [3].
The analogy with Markov processes here can give us some insight on the problems with the H matrix. If we interpret nodes as states, we can see that a closed class of state (possibly a singleton) would be a serious issue. And any node with no outgoing links constitutes such a class.
In the Markov chain interpretation, if such a node is reachable, then we would just remain there. Such nodes would be rank sinks. In such a setting, every web site with no outgoing links would acquire a rank far superior to its actual importance. Using Markov chains terminology, the H matrix is substochastic.
Another problem is related to cycles. E.g., consider the matrix
C =
0
1
1
0

Then: C2n = I and C2n+1 = C. As a consequence we cannot expect to converge (notice, here we have eigenvalues with same magnitude, so we know that the power method is not going to work).
Luckily enough, these problems are rather easy to fix modifying the matrix. More on this, later!

References

[1]
M. Franceschet. PageRank: Standing on the shoulders of giants. arxiv.org.
[2]
M.E.J. Newman. The structure and function of complex networks. Arxiv preprint cond-mat/0303516, 2003.
[3]
M.E.J. Newman. Networks: An Introduction. Oxford University Press, 2010.
[4]
L Page and S Brin. The pagerank citation ranking: Bringing order to the web. 1998.



Sunday, November 7, 2010

Java Collections sorting: copy, copy & copy in every sense!

This morning I was hacking into OpenJDK some implementation details. Curiosity, perhaps.

After the announcement that Java is going to use timsort for generic sorting (which is a very good piece of news).
By the way, this is an actual implementation. I'm also wondering why these things are mostly ignored in Algorithm courses, but this is another story.

So... I checked the actual implementation. I had to download the OpenJDK source, as apparently on OS X IntelliJ does not auto-show implementation of base classes when I "jump to implementation". Anyway... the implementation of sorting is basically in Arrays.java.

I want to point out the funny comment relative to sort1 (which is the private function which actually does the job for non floating point stuff):
    /*
     * The code for each of the seven primitive types is largely identical.
     * C'est la vie.
     */


which really saved the day... Yes... it is almost cut and paste for the primitive types (char, int, byte, etc). I was amazed... it is 2010, we have generics but, as primitive types are not objects we have to actually write it 7 times. As a side note, code is recursive quicksort variant.

Ok, nice. You have implementation constraints and you have to cope with that. Yes... C++ templates perfectly solve this crap. But then, there is another very funny discovery (well, in fact I have known this for quite a few time). But that is not the only (technically) unnecessary copy which happen(ed|s) when dealing with sorting in Java; this is the sort function in Collections.java

public static <T extends Comparable<? super T>> void sort(List<T> list) {
Object[] a = list.toArray();
Arrays.sort(a);
ListIterator<T> i = list.listIterator();
for (int j=0; j<a.length; j++) {
i.next();
i.set((T)a[j]);
}
}


Wow! It is clear what it does? For every collection you sort, a copy is made. To be honest, the documentation states this clearly:


This implementation dumps the specified list into an array, sorts
the array, and iterates over the list resetting each element
from the corresponding position in the array. This avoids the
n^2 log(n) performance that would result from attempting
to sort a linked list in place.


Nice! Very nice... this means that the default sort function may be unnecessarily slow. Unless you use arrays (but using arrays in place of collections is considered a bad practice). Moreover, Collection is an interface: its implementation may be such that a full copy into main memory is not feasible. BTW, I should double check because the comments in Arrays.sort state that the algorithm is a modified quicksort, not a modified mergesort.

Oh... nice! Tomorrow I am discussing other methods in Arrays and Collections.

Sunday, October 17, 2010

Anagrams: fast python algorithm

So in this post I showed how the trivial algorithm to compute anagrams is stupid (and consequently slow) beyond hope. On the other hand I promised a nice algorithm to do that. In fact, it is just a matter of pre-processing the data. The idea is to compute a key which is unique to all the words which have the very same letters and then group all the words with the same key. One of the simplest keys could be a bag of the letters (if your language supports bags -- or multi-sets --). Python sets have unique keys, so they are not a good candidate. However a "sorted" variant of the word has the right properties. Every word with the same letters (which is the definition of palindrome) yields the same string if we sort its letters. A good structure is a dictionary with this key as key and a list of strings as value. When we find a string with a previously found key, we add it to the list of strings, otherwise we generate a new list. This is the very popular setdefault idiom:
d.setdefault(key, []).append(w)
The key for a word of length n is computed in O(n log n), insertion in the dictionary is O(1) (get the key is O(1) in average and append to the list is O(1)). So the algorithm uses rather efficient elementary operations. At this point "general" questions on the processed wordlist can be answered in O(M) if M is the number of the words (e.g., longest anagram, word with most anagrams). Simply getting the anagrams of a given word is an O(n log n) operation. We are not using an excessive amount of memory. We basically have an hash-table and some "pointers" to the original words. And now the code:
class Anagrams(object):
    def __init__(self, dictionary):
        self.dictionary = dictionary
        self.indexed_words = {}
        self.preprocess()

    def make_key(self, word):
        sorted_word = list(word)
        sorted_word.sort()
        return tuple(sorted_word)

    def preprocess(self):
        for word in self.dictionary:
            key = self.make_key(word)
            self.indexed_words.setdefault(key, []).append(word)

    def anagram(self, word):
        key = self.make_key(word)
        return tuple(self.indexed_words[key])

    def anagrams(self):
        return [anagrams for key, anagrams in self.indexed_words.iteritems()
                if len(anagrams) > 1]
As a matter of fact, if we extend the corpus of words, the required time increases. In fact, using the dictionary from here it took 4 whole seconds. Which is nonetheless a lot less time than some million of years. I did not use that dictionary because it contains lots of misspells which can have sense when computing Levenshtein distances, but are quite a PITA if we compute anagrams (we don't want non existing words to popup as solutions).

Technorati Tags: , , , ,

Monday, October 11, 2010

Design of Sorting in Java and Python (pt. 2)

In this post I presented a pattern which arises very quickly when doing comparisons among objects. The pattern is essentially:


def name_cmp(fst, snd):
 return cmp(f(fst), f(snd))

In fact I see huge repetition in this kind of code, because we write code in that shape for many different "f"s (that is to say, getter or more complicated computation, but it boils down to that pattern).

For example, let us define a Person object in Java and also define comparators to sort a list of persons according to both the number and the phone number.

public class Person {
    private String name;
    private String phone;

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

    public String getPhone() {
        return phone;
    }

    public String getName() {
        return name;
    }

    static final private PhoneComparator PHONE_COMPARATOR = 
     new PhoneComparator();
    static final private NameComparator NAME_COMPARATOR = 
     new NameComparator();
    
    static private class PhoneComparator implements Comparator<Person> {
        public int compare(Person o1, Person o2) {
            return o1.getPhone().compareTo(o2.getPhone());
        }
    }

    static private class NameComparator implements Comparator<Person> {
        public int compare(Person o1, Person o2) {
            return o1.getName().compareTo(o2.getName());
        }
    }

    static Comparator<Person> ageComparator() { return PHONE_COMPARATOR; }
    static Comparator<Person> nameComparator() { return NAME_COMPARATOR; }
    
    // ide generated toString, equals etc. here\
};

I see huge repetition in this kind of code. But even in Python, defining all the compatators is boring and repetive. Of course in Python something like:

def cmp_buider(f):
 def cmp_(fst, snd):
  return cmp(f(fst), f(snd))
 return cmp_

would come quite natural. For the ones who do not like functionally flavoured code, this is the Java variant (which is longer and uglier).

public interface Key<K extends Comparable<K>, E> {
    public K index(E obj);
}

public class KeyToComparator {
    public static <K extends Comparable<K>, T> Comparator<T> build(final Key<K, T> key) {
        return new Comparator<T>() {
            public int compare(T o1, T o2) {
                return key.index(o1).compareTo(key.index(o2));
            }
        };
    }
}

and then we use it like:

static public void sortWithKey(List<Person> lstPerson) {
    Collections.sort(lstPerson, KeyToComparator.build(
            new Key<String, Person>() {
                public String index(Person obj) {
                    return obj.getName();
                }
            }));
}

I hope I'm not the only one feeling that this stuff is ugly as hell. However, it is quite general and powerful. And it lets you avoid code repetition. At the very least, the code to write a key is way shorter than the one to write a Comparator. Of course, if methods where first class objects...

Sunday, October 10, 2010

Parametrized Unit Tests (in Java)

While "creating" by hand parametrized test functionality is not extremely hard in Python (and in fact there are multiple people who did that, plus the "blessed" implementation which is going to be included in unittest2 some time soon, I hope), doing the same thing in Java is a nightmare.

Luckily enough, JUnit 4 offers the functionality out of the box.

In fact, I have to admit that it seems cleaner than Python alternative. Let us start with a Java quicksort (this time ugly, but not buggy, as far as my tests are correct).

package sorting.quicksort;

import java.util.List;

public class QuickSorter<E extends Comparable<E>> {
    List<E> lst;

    public synchronized void sort(List<E> lst) {
        this.lst = lst;
        sort(0, lst.size()-1);
    }

    private void sort(final int start, final int end) {
        if(start >= end) {
            return;
        }

        int pivot = partition(start, end);

        sort(start, pivot-1);
        sort(pivot+1, end);
    }

    private int partition(int start, int end) {
        int pivot = choosePivot(start, end);
        swap(pivot, start);

        int leftProcessed  = start, rightProcessed = end + 1;
        E pivotElement = lst.get(start);

        while(true) {
            do {
                leftProcessed++;
            } while(leftProcessed <= end &&
                lst.get(leftProcessed).compareTo(pivotElement) < 0);
            do {
                rightProcessed--;
            } while(lst.get(rightProcessed).compareTo(pivotElement) > 0);
            if(leftProcessed > rightProcessed) {
                break;
            }
            swap(leftProcessed, rightProcessed);
        }

        swap(start, rightProcessed);
        return rightProcessed;
    }

    private void swap(final int i, final int j) {
        E tmp = lst.get(i);
        lst.set(i, lst.get(j));
        lst.set(j, tmp);
    }

    private int choosePivot(final int left, final int right) {
        return left + (right - left)/2;
    }
}

The test case is:
package sorting.quicksort;

import org.junit.Assert;
import org.junit.Before;
import org.junit.Test;
import org.junit.runner.RunWith;
import org.junit.runners.Parameterized;

import java.util.*;

@RunWith(value=Parameterized.class)
public class QuickSorterTest<E extends Comparable<E>> {
    private List<E> lst;
    private List<E> copyOfLst;
    private QuickSorter<E> sorter;

    public QuickSorterTest(List<E> lst) {
        this.lst = lst;
    }

    @Parameterized.Parameters
    public static Collection<Object[]> getParameters() {
        return Arrays.asList(
                new Object[]{ Arrays.asList(1, 2, 3, 4) },
                new Object[]{ Arrays.asList(1, 2, 4, 3) },
                new Object[]{ Arrays.asList(1, 3, 2, 4) },
                new Object[]{ Arrays.asList(1, 3, 2, 3) },
                new Object[]{ Arrays.asList(4, 2, 3, 1) },
                new Object[]{ Arrays.asList(1, 4, 3, 2) },
                new Object[]{ Arrays.asList(3, 2, 1, 4) },
                new Object[]{ Arrays.asList(3, 2, 2, 4) },
                new Object[]{ Arrays.asList(3, 1, 1, 3) },
                new Object[]{ Arrays.asList(1, 1, 1, 1) },
                new Object[]{ Arrays.asList(1, 2, 1, 1) },
                new Object[]{ Arrays.asList(1, 2, 3) },
                new Object[]{ Arrays.asList(1, 3, 2) },
                new Object[]{ Arrays.asList(2, 1, 3) },
                new Object[]{ Arrays.asList(2, 3, 1) },
                new Object[]{ Arrays.asList(3, 1, 2) },
                new Object[]{ Arrays.asList(3, 2, 1) },
                new Object[]{ Arrays.asList(3, 2, 3) },
                new Object[]{ Arrays.asList(3, 1, 1) },
                new Object[]{ Arrays.asList() },
                new Object[]{ Arrays.asList(1) },
                new Object[]{ Arrays.asList("hello", "cruel", "world") });
    }

    @Before
    public void setUp() {
        lst = new ArrayList<E>(lst);
        copyOfLst = new ArrayList<E>(lst);
        sorter = new QuickSorter<E>();
    }

    @Test
    public void testSort() throws Exception {
        sorter.sort(lst);
        Collections.sort(copyOfLst);
        Assert.assertEquals(copyOfLst, lst);
    }
}

Now, the idea is quite clear. The code to be written is ugly. The parametrized test case is annotated with @RunWith and the Parametrized.class. The test case must have a method annotated with @Parameterized.Parameters. This method return a collection of array of objects.

The idea is that for each element returned by that method a new test case object is constructed using the array of objects as actual arguments. In fact, the design is not extremely different from what we saw in Python.

In Python we had a generator yielding the tests. Here we have a static method which does quite the same job. Of course, here a Collection is returned. And Java syntax here is quite cumbersome. Luckily enough there is some flexibility with types so that I do not need to make different test case classes just because of that.

In the example we sort both arrays of integers and arrays of strings.

The setUp method and the test method do roughly the same thing they do in Python. And also the Java version uses a constructor (in fact, I'm led to believe that it is simply unavoidable).

Friday, October 8, 2010

Parametrized Unit Testing (in Python)

In my last post I briefly discussed parametrized tests. Now it is time to see them in practice.

Right now, in Python there is no way to do parametrized testing. Well, actually there are many ways. Python is an extremely flexible and high level languange and adding that features is a matter of less than a hundred lines of code (or somewhat more, perhaps, if you want an even richer semantics). Just the bare bone functionality may take just a dozen lines of code.


For this reason, many people implemented their own version (see the linked thread). Plus there are nose and py.test variants. Actually, I have used nose to perform that kind of tests.


Nose allows tests to be simple functions (and uses python decorators to add setUp and tearDown functions when necessary).This is a (buggy) quicksort implementation in Python:

def quicksort(lst, select_pivot=lambda s, e: s + (s-e)/2):
    def partition(start, end):
        pivot = select_pivot(start, end)
        lst[pivot], lst[start] = lst[start], lst[pivot]
        left = start
        right = end + 1
        pivot_element = lst[pivot]
        
        while 1:
            left += 1
            right -= 1
            while (left <= end) and (lst[left] < pivot_element):
                left += 1
            while lst[right] > pivot_element:
                right -= 1
            if left > right:
                break
            lst[left], lst[right] = lst[right], lst[left]
        lst[start], lst[right] = lst[right], lst[start]
        return right
    
    def sort(start, end):
        if start >= end:
            return
        pivot = partition(start, end)
        sort(start, pivot-1)
        sort(pivot+1, start)
        
    sort(0, len(lst)-1)

And this is the test:

import qsort
from nose import tools

fixture_lsts = [
    [1, 2, 3, 4], [1, 2, 4, 3], [1, 3, 2, 4], [1, 3, 2, 3],
    [4, 2, 3, 1], [1, 4, 3, 2], [3, 2, 1, 4], [3, 2, 2, 4],
    [3, 1, 1, 3], [1, 1, 1, 1], [1, 2, 1, 1], [1, 2, 3],
    [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1],
    [3, 2, 3], [3, 1, 1], [], [1]
]

def check_sorted(lst):
    sorted_lst = sorted(lst)
    lst = lst[:] # do not modify lists
    qsort.quicksort(lst)
    tools.assert_equals(lst, sorted_lst)
    
def test_quicksort():
    for lst in fixture_lsts:
        yield check_sorted, lst
        


This is what I see from WingIDE:

Buggy quicksort test output


In fact, the test function is allowed to be a generator. In this case, it is assumed to yield tuples where the first argument is a callable and the others are its parameters. This way it is very easy to do parametrized testing.

Please notice (the documentation points that out) that if you attach setup and teardown stuff to the test function then they will be run just once, before and after the whole test function and not before and after every parametrized test. This is what I usually would like, in fact. In order to have this behavior, it is sufficient to add the functions to the yielded callable.

For example, I see that in the proposed "check" function we are doing stuff which should be done in a setup method. Unfortunately one of my greatest dissatisfactions with nose (and with_setup in particular) is that it completely forces you to use a stateful model.

From the setup function to the executed test function the only way to communicate is through the global environment. This is half bad (I've not yet decided if I prefer the overkill of TestCase classes -- where classes are just used for grouping purposes without "true" object behavior from the user point of view -- or using module variables -- which seems very pythonic as modules are objects but somewhat encourages the use of global variables).

This works quite well if you model all your tests with this idea in mind (and basically it's not different from xUnit model, just no classes). Unfortunately when you use generator based testing the test function receives the part of the SUT to test as a parameter, and the parameter is not accessible from the setup function. The easiest solution (the one which does not involve writing some custom test loader or at least a new piece of code at nose level), is useing classes.

In fact, the yielded callable is allowed to be a callable object. In this case, it is possible to define a setup function "as if" it were a TestCase (still, it is not a TestCase.

The resulting code is:

class CheckSorted(object):
    def __init__(self, lst):
        self.lst = lst

    def setup(self):
        self.sorted_lst = sorted(self.lst)
        self.lst = self.lst[:]

    def __call__(self):
        qsort.quicksort(self.lst)
        tools.assert_equals(self.lst, self.sorted_lst)


def test_quicksort2():
    for lst in fixture_lsts:
        yield CheckSorted(lst),

Please notice the , (comma) at the end of the yield statement. It is necessaty as a tuple must be yielded. Of course code could be made more explicit with

yield (CheckSorted(lst), )

I don't actually like this solution very much, since it is somewhat confusing what the costructor should do: usually TestCases do not have a constructor. Setup methods just set them up. However, in this case, part of the job is performed by the constructor (and there is no other way -- unless we use a global or something like that, but I believe the cure is worse than the illness). Anyway, this is what we have.

Friday, August 6, 2010

Quicksort in Continuation Passing Style

Quicksort is ugly. I hate quicksort. Still everyone praises quicksort...
... and Hoare did worse things. Such as Hoare logic. I had to manually verify some short programs using Hoare logic and sometimes I still wake up screaming at night after dreaming pages of Hoare triples.

Well, this is a problem of mine. Back to quicksort... among the elementary sorting algorithms it is one of the more difficult to understand. But it is rather efficient despite its worse case N^2 running time. Moreover, quicksort can be implemented in constant space (well, almost... you still have double recursive calls which means log n, if you consider the stack).

That is why (the idea that quicksort modifies an array in place) some believe functional implementations of quicksort are not quicksort. And I'm not going to argue. I don't care, in fact. Here my concern are the double recursive calls. Normally you could not tail call optimize such a thing (pseudo code from wikipedia):

procedure quicksort(array, left, right)
     if right > left
         select a pivot index //(e.g. pivotIndex := left+(right-left)/2)
         pivotNewIndex := partition(array, left, right, pivotIndex)
         quicksort(array, left, pivotNewIndex - 1)
         quicksort(array, pivotNewIndex + 1, right)

Luckily for us, we can use continuation passing style to trick this into tail call optimizable. Code in continuation passing style is always tail recursive.

Append:
(define cp-append 
  (lambda (lst1 lst2 k)
    (cond ((null? lst1) (k lst2))
          (else (cp-append (cdr lst1) lst2 
                           (lambda (rest)
                             (k (cons (car lst1) rest))))))))


The continuation function for the partition function shall take two parameters.
The first one "accumulates" less than items, the second one greater than items.

(define cp-partition
  (lambda (lst p? k)
    (letrec ([cpp 
              (lambda (lst k)
                (cond
                  ((null? lst) (k '() '()))
                  ((p? (car lst)) 
                   (cpp (cdr lst)
                        (lambda (p-true p-false)
                          (k (cons (car lst) p-true)
                             p-false))))
                  (else
                   (cpp (cdr lst)
                        (lambda (p-true p-false)
                          (k p-true
                             (cons (car lst) p-false)))))))])
      (cpp lst k))))

And eventually, his majesty the quicksort:
(define quicksort
  (lambda  (lst less?)
    (letrec 
        ([qs
          (lambda (lst k)
            (cond
              ((null? lst) (k '()))
              (else
               (let ([pivot (car lst)]
                     [rest (cdr lst)])
                 (cp-partition 
                  rest
                  (lambda (x) (less? x pivot))
                  (lambda (less-than greater-than)
                    (qs greater-than
                        (lambda (sorted-gt)
                          (qs less-than
                              (lambda (sorted-lt)
                                (cp-append
                                 sorted-lt
                                 (cons pivot sorted-gt) k)))))
      (qs lst (lambda (v) v)))))


Although the testing based on random data is a Very Bad Idea, this was the easiest way to rough up my code.

(define (random-list max length)
  (letrec ([rl (lambda (length)
                 (cond ((= length 0) '())
                       (else (cons (random max)
                                   (rl (- length 1))))))])
    (rl length)))

(define (test-quicksort)
  (let ([tests-per-length 50]
        [random-top-integer 1000]
        [list-lengths '(0 1 10 11 25 40 64 128)])
    (for-each 
     (lambda (length)
       (do ([i tests-per-length (- i 1)])
         ((zero? i))
         (let* ([the-list (random-list random-top-integer length)]
                [lt-sorted-list (sort the-list <)]
                [lt-qsorted-list (quicksort the-list <)]
                [gt-sorted-list (sort the-list >)]
                [gt-qsorted-list (quicksort the-list >)])
           (when (not (equal? lt-sorted-list lt-qsorted-list))
             (printf "~a[<]:~n~a~n~a~n~n" the-list
                     lt-sorted-list lt-qsorted-list))
           (when (not (equal? gt-sorted-list gt-qsorted-list))
             (printf "~a[>]:~n~a~n~a~n~n" the-list
                     gt-sorted-list gt-qsorted-list)))))
     list-lengths)))

I tried to avoid racket specific constructs, but I'm not sure that one or two slipped in the testing code.

BTW, here a version with named-let:
(define quicksort
  (lambda  (lst less?) 
    (let qs ([lst lst] [k (lambda (v) v)])
      (cond 
        ((null? lst) (k '()))
        (else
         (let ([pivot (car lst)]
               [rest (cdr lst)])
           (cp-partition 
            rest 
            (lambda (x) (less? x pivot))
            (lambda (less-than greater-than)
              (qs greater-than
                  (lambda (sorted-gt)
                    (qs less-than
                        (lambda (sorted-lt)
                          (cp-append 
                           sorted-lt
                           (cons pivot sorted-gt) k)))))))))))))

Seems I already wrote on Quicksort...

Sunday, April 23, 2006

Quicksort (C and Python)

The quick-sort algorithm

Quick-sort is a sorting algorithm developed by Hoare that in average does O(n log n). The O() notation means that if the list length is n, the number of operations used to sort the list is in the order of magnitudo of n log n (in this case).

To be more formal, we shall say that there exist two constants c1 and c2 such that the number of operations m isc1 * n * log n <= m <= c2 * n * log n. Notice that c1 and c2 are not specified in the O notation
We stop with math here. If you want more precise informations about the O notation read the wikipedia.

Now I'm going to briefly explain the quick-sort algorithm, and then we will use that algorithm to exploit the expressive power of some programming languages. Keep in mind that more expressive does not necessarily mean faster to run. It means faster to write. But being faster to write means that we can try more efficient algorithms, thus improving performance in a more effective way.
Expressive languages are usually high-level. In fact we can also profile the code, track the bottleneck and rewrite only a couple of functions in C or in ASM. This lead to fast programs with fast development cycles.

Quick-sort strategy

Quick-sort uses a divide and conquer approach. We divide a list in two sublists using a pivot element. All the elements less than the pivot come before it, those greater after. Then we recursively apply quick-sort to the two sublists.
Notice that the pivot is chosen arbitrarily. Some implementations chose the first element of the list. Some chose it randomly, some use other criteria.
Moreover there are implementations of quick-sort that sort the list in place (that is to say memory used is bounded), some need extra storage. Guess which are the more efficient.

On wikipedia you can find more informations about the quick-sort algorithm. You can find also some implementations, and a link to a page full of implementations in different languages.

I almost did not write code. I took code from the wikipedia or around the web and I'm going to comment it. This is the "new" work I did. Comments, explanations, not code. I also wrote a couple of implementations, however I don't remember which ones.

Python

2.4.2 (#1, Mar 23 2006, 22:00:44)  
[GCC 4.0.1 (Apple Computer, Inc. build 5250)] 

Now lets examine a Python high level implementation.

def qsort(L): 
    if L == []: return [] 
    return qsort([x for x in L[1:] if x< L[0]]) + L[0:1] + \ 
            qsort([x for x in L[1:] if x>=L[0]]) 

This is the exact description of the quick-sort algorithm, if you know what Python list comprehension are. When I write [x for x in L if p(x)], Python computer the list of elements of the L container that make p(x) true.

L[a:b] is the sublist of L[L[a], L[a+1], ..., L[b-1], L[b]]. L[a:]==L[a:len(L)], L[:b]==L[0:b]

[x for x in L[1:] if x< L[0]] is the sublist of the elements of L a part from the first one (that is the pivot) such that x is less than the pivot (L[0]).

[x for x in L[1:] if x>=L[0]] is the sublist of the elements of L a part from the first one (that is the pivot) such that x is greater or equal than the pivot (L[0]).

Then we pass those sublists to qsort that recursively sorts them, and eventually we return a append all those sublists and return them.

This implementation, although beautiful to read and understand is horrible from a performance point of view. It creates, concatenates and destroys a large number of lists. Moreover it is purely recursive. For large lists, we can even exhaust the stack space.

A version of quick-sort in Python that has no need for more memory is this:

def partition(array, begin, end, cmp): 
    while begin < end: 
        if cmp(array[begin], array[end]): 
            (array[begin], array[end]) = (array[end], array[begin]) 
            break 
        end -= 1 
    while begin < end: 
        if cmp(array[begin], array[end]): 
            (array[begin], array[end]) = (array[end], array[begin]) 
        break 
        begin += 1 
    return begin 
    
    def sort(array, cmp=lambda x, y: x > y, begin=None, end=None): 
        if begin is None: begin = 0 
        if end   is None: end   = len(array) 
        if begin < end: 
            i = partition(array, begin, end-1, cmp) 
        sort(array, cmp, begin, i) 
        sort(array, cmp, i+1, end) 

[NOTE: code has been reformatted without running it again: it is likely
to be wrong!]

If you want to run this code, remember that in Python indentation matters. If you just copy and paste the code, you may end with unusable code.
Anyway, I benchmarked all of them. "Elegant quicksort" is the first implementation. "Standard quicksort" is the latter. "Bultin sort" is standard python sort algorithm (that is not a quicksort, but a dedicated sorting algorithm).

Elegant qsort: [21.20, 19.75, 20.08, 20.46]
Builtin qsort: [1.61, 2.29, 2.29, 2.29]
Standard qsort: [31.15, 30.60, 30.68, 30.38]

Probably due to heavy optimization performed in list comprehension the "elegant" algorithm is faster than the standard one. [from the future: I hardly believe that]. However both of them are more than ten times slower than the python built-in sort algorithm.
We have learnt:

  1. High level languages are great to express algorithms in a readable way (and we can use them to test them and see how they work and all)
  2. You should not implement "basic" algorithms in high level languages: use their built-ins or write them in C (most built-ins are written in C).

C

Now it is time for a couple of pure C quick-sort:

void qsort1(int a[], int lo, int hi ){ 
    int h, l, p, t; 
    if (lo < hi) { 
        l = lo; 
        h = hi; 
        p = a[hi]; 
        do { 
            while ((l < h) && (a[l] <= p))  
                l = l+1; 
            while ((h > l) && (a[h] >= p)) 
                h = h-1; 
            if (l < h) { 
                t = a[l]; 
                a[l] = a[h]; 
                a[h] = t; 
            } 
        } while (l < h); 
        t = a[l]; 
        a[l] = a[hi]; 
        a[hi] = t; 
        qsort( a, lo, l-1 ); 
        qsort( a, l+1, hi ); 
    } 
} 

Another implementation is:

void quicksort(int x[], int first, int last) { 
    int pivIndex = 0; 
    if(first < last) { 
        pivIndex = partition(x,first, last); 
        quicksort(x,first,(pivIndex-1)); 
        quicksort(x,(pivIndex+1),last); 
    } 
} 

int partition(int y[], int f, int l) { 
    int up,down,temp; 
    int piv = y[f]; 
    up = f; 
    down = l; 
    goto partLS; 
    do {  
        temp = y[up]; 
        y[up] = y[down]; 
        y[down] = temp; 
        partLS: 
        while (y[up] <= piv && up < l) { 
         up++; 
        } 
        while (y[down] > piv  && down > f ) { 
            down--; 
        } 
    } while (down > up); 
    y[f] = y[down]; 
    y[down] = piv; 
    return down; 
} 

C also have a builtin quick-sort, that is much more advanced than this ones, since it allows to specify a comparison function, thus working with all kinds of data types.

Benchmarking these was hard. They too fast to be benchmarked with a list of a million elements and the standard posix clock function. However if you want to run the bench yourself, I could not use MacOS X specific bench functions.
So I used a list of 100000000 elements.
This is the (ugly) code I used for benchmarking. The test is not really interesting. Of course since C integers are machine integers it is dramatically faster than other solutions. It could have been more interesting if we sorted other data structures, where the gap could be reduced.

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include "qsort1.h" 
#include "qsort2.h" 
#define LEN 100000000 

int intcmp(const void *a, const void *b){ 
    return *(int*)a - *(int*)b; 
} 

int main(){ 
    size_t index; 
    int *l = malloc(sizeof(int)*LEN); 
    clock_t last, current; 
    srand(time(NULL)); 
    for(index = 0; index < LEN; ++index){ 
        l[index] = rand(); 
    } 
    
    last = clock(); 
    qsort(l,LEN, sizeof(int), intcmp); 
    current = clock(); 
    printf("C qsort: %f\n", (current-last)/(float)CLOCKS_PER_SEC); 
    
    for(index = 0; index < LEN; ++index){ 
        l[index] = rand(); 
    } 
    last = clock(); 
    qsort1(l,0, LEN-1); 
    current = clock(); 
    printf("int qsort: %f\n", (current-last)/(float)CLOCKS_PER_SEC); 
    for(index = 0; index < LEN; ++index){ 
        l[index] = rand(); 
    } 
    last = clock(); 
    quicksort(l,0, LEN-1); 
    current = clock(); 
    printf("qsort 2: %f\n", (current-last)/(float)CLOCKS_PER_SEC); 
    free(l); 
    return 0; 
}