Sunday, August 8, 2010

Pydistutil Configuration

Something I rather hate is mixing up my development environment with my ordinary system environment. I like having things I installed by myself separated from system stuff (and I'm not the only one... think about the convention of /usr/local).

My development machine are essentially single user machines. I am the only user. Thus, global installations are simply an overkill. And I hate pervasive package management systems. For the MacOS, for example, I chose brew. I love it. It does the minimum to automate tedious stuff, but it leaves me complete freedom.

Of course, since I want my package manager to do the minimum necessary stuff, I don't want it to touch my python modules.
Don't touch my Python!
 Well... anyway. As a consequence I heavily rely on python own module manager, that is to say pip/easy_install. You can tell distutils the default base install directory. For example, on Linux, I have in my home directory this .pydistutil.cfg file:

% cat .pydistutils.cfg  
[install]
prefix = ~/.local
install_scripts = ~/bin
install_lib = ~/.local/lib/python$py_version_short/site-packages

which basically makes all the installation in my home-directory.
More information here and here. I use this for modules I plan to use in multiple projects.
Otherwise, virtualenv.

Saturday, August 7, 2010

Racket, can it be a Good Thing?

On June, 7th 2010 PLT Scheme changed name. Now it is called Racket.
I will not indulge in how bad does it sound (especially if you live in countries where racket is a social plague). I bet if some developer had his project dubbed "Rape" (Recursive Anamorphic Program Environment?) some people would argue. However, the authors do like it, and that's it.

PLT Scheme was a very good programming environment. I love scheme. I love Planet. There is a huge quantity of modules (which is something I really appreciate as a pythonista) and they are not only academic tools. Which is good. This comes with Racket as well.

Before the name change, PLT scheme was somewhat different from most scheme environments out there. First, it came batteries included. And as someone new to the platform I often found myself using functions which were srfi or even non standard. In both cases some alternative platforms I used did non include them. This was somewhat annoying...

Moreover the #scheme first line makes a "scheme" source non valid scheme. Thus this was not a minor issue to work with multiple environments. These stuffs is rather trivial to solve, simply I would have preferred to spend my time currying rather that understanding how each scheme environment extended r5rs with modules and make things r6rs compatible as well.

At least now it is Racket. It's based on scheme, but it's not scheme (which is something that, when done by Microsoft is regarded as a criminal offense). So I don't have to expect my racket programs work with scheme compilers.

Moreover, you can write in your CV:
IT skills: racket

... and that is not something to underestimate.

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...

Thursday, August 5, 2010

Reverse in continuation passing style

A simple reverse in continuation passing style...

(define cp-reverse
  (lambda (lst)
    (letrec ([cpr (lambda (lst k)
                 (cond 
                   ((null? lst) (k '()))
                   (else 
                    (cpr (cdr lst)
                         (lambda (v)
                           (cons (car lst) (k v)))))))])
      (cpr lst (lambda (v) v)))))

This is a rather classical example, indeed. Performance wise it's much slower than the reverse builtin in PLT-Scheme/Racket, still it's a lot faster than even more classical trivial implementations such as:


(define (sl-reverse lst)
  (cond
    ((null? lst) '())
    (else (append (sl-reverse (cdr lst)) (list (car lst))))))

In fact the example relies on TCO modulo cons.

In order to test the whole thing in racket I used this snippet:

(define (print-times min max step)
  (for ([m (in-range min max step)])
    (let ([lst (for/list ([i (in-range m)]) i)])
      (time (reverse lst))
      (time (cp-reverse lst))
      (time (sl-reverse lst)))))


> (print-times 5000 10000 1000)
cpu time: 0 real time: 10 gc time: 0
cpu time: 5 real time: 61 gc time: 0
cpu time: 1320 real time: 2088 gc time: 1099
cpu time: 1 real time: 0 gc time: 0
cpu time: 2 real time: 2 gc time: 0
cpu time: 463 real time: 493 gc time: 185
cpu time: 0 real time: 0 gc time: 0
cpu time: 2 real time: 21 gc time: 0
cpu time: 888 real time: 926 gc time: 505
cpu time: 0 real time: 0 gc time: 0
cpu time: 2 real time: 2 gc time: 0
cpu time: 2729 real time: 2766 gc time: 2235
cpu time: 0 real time: 0 gc time: 0
cpu time: 2 real time: 2 gc time: 0
cpu time: 2157 real time: 2227 gc time: 1522

Thursday, July 8, 2010

How much do you cost me?

Introduction

Study of algorithms and data structures is a major field in computer science. Efficiently representing and processing data has been of paramount importance since the early days of computing. Many problems are simply not accessible without proper study, since the computational complexity used to overwhelm older machines.

Today, we have powerful machines. Tomorrow, we will have more powerful and “more parallel” machines. One day, our computers will surpass the computing power of our own brain. However, no matters how fast the computer it is, there are always problems which are too hard, too expensive and, in short, intractable.

Essentially, the cost of the algorithm depends from the size of its input. For example, searching for an element in a list intuitively takes time proportional to the number of elements: after all, we don’t know where it is, we have to start somewhere and the element could be at the other end of the list.

In this case we say that the worst case is exactly A∙N, if A is the cost of comparing two elements and N is the length of the list. The average case is “only” A∙(N/2) since “in average” the element will be found after examining only half of the elements. Intuitively, the probability the element is the last one, is equal to the probability it is the first one, and so one. Computing average case cost is usually harder than worst case cost, because involves probability and random variables.

Landau notation and estimates

The computational cost of algorithms is usually represented in a unified notation (called Landau notation, from the name of the famous mathematician who introduced the notation in the field of calculus). Essentially, we don’t care about the multiplicative constants in the costs (e.g., A, “the cost of comparing two elements”) as long as they are, well… constants. We usually say that an looking for an element in a list costs O(N), which means that the time is proportional to the size of the list itself.

Most of the times, computing the cost of an algorithm is rather straightforward. Considering a strictly imperative setting, you essentially estimate the number of times you run the body of a loop. For example, let us consider
the following code:

M = {}
n = 4

for i in xrange(0, n):
    for j in xrange(0, i):
        M[i,j] = M[j,i] = i + j

print M

It is rather easy to see that if A is the cost of M[i,j] = M[j,i] = i + j, the whole thing costs (using presumed Gauss's formula):



That is to say it's O(N^2). A (correct) but coarse argument is that the internal loop runs i times, and i is at most n. So in the worse case it is n. The upper bound is used as an estimate n*n (again O(n*n)).

Things are not always that easy: sometimes the worse case is completely uninteresting. E.g., the most widely used simplex algorithm[0] has an exponential worse case cost. Exponential means that it essentially sucks. 2^1000 is a huge[1] number and 1000 is not an extremely huge input size. In pratice, the simplex algorithm is used, since the average complexity is polynomial. There is still active research on the estimation of the complexity of that algorithm.

The Master Theorem

A large class of algorithms are called "divide et impera"; a typical example is mergesort[2]. Here a trivial python implementation is given:

def merge_sort(L):
    if len(L) <= 1:
        return L
    else:
        middle = len(L) / 2
        left = merge_sort(L[:middle])
        right = merge_sort(L[middle:])
        return merge(left, right)

def merge(left, right):
    i = j = 0
    max_i = len(left)
    max_j = len(right)
    sorted_ = []
    while i < max_i and j < max_j:
        if left[i] < right[j]:
            sorted_.append(left[i])
            i += 1
        else:
            sorted_.append(right[j])
            j += 1
    if i < max_i:
        sorted_.extend(left[i:])
    else:
        sorted_.extend(right[j:])
    return sorted_
Of course, no Python user would use such functions. They use recursion instead of iteration, they are not particularly efficient, the "Timsort"[3] implemented in the Python library is far more efficient both algorithmically and code-wise. Algorithm books tell us that merge-sort has worse and average case O(N log N). As with most divide et impera algorithms, the result can be proved with ad hoc demostrations. However, I found extremely useful using the so-called master theorem[4]. This is essentially a swiss-knife result that gives us the complexity of an algorithm given the way the data is "divided" and "recomposed" (to rule, of course!). The complexity of the algorithm on an input of size n is T(n). f(n) is the cost done outside the recursive calls (e.g., the time needed to decompose and put together the data); a is the number of subproblems in the recursion and n/b is the size of the input of the subproblems. For example, in merge-sort a = 2 and b = 2. f(n) is proportional to the size of the input vector. The master theorem states that if:


The big theta notation essentially says that the specified bound is asymptotically both a lower and an upper bound. The big omega means that the specified bound is a lower bound.

Notes

[0] http://en.wikipedia.org/wiki/Simplex_algorithm
[1] 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
[2] http://en.wikipedia.org/wiki/Merge_sort [3] http://bugs.python.org/file4451/timsort.txt
[4] Introduction to Algorithms; Cormen, Leiserson, Rivest, and Stein; MIT Press

More on Word

Some pieces of software are unbelievable: as soon as you start noticing how nice are certain features, you discover that some other parts are so incredibly broken that you would want to put the whole thing in the thrash.

As is customary, I received the template to format my article, which was an IEEE standard. Nice. I opened it in Word 2007 for Mac (and the template should be compatible with that) and I removed the contents of the template (as advised on the template itself) and filled with my own text. Shit happens.

For no good reason, the styles built in the template were imported utterly broken. For example the base font size was plainly wrong, some styles did not use the “keep together with next paragraph” thing and so on. Essentially, I had to manually check the styles I used with ones in the unmodified template in order to fix things. And I was able to, however, it took time.

Then I sent the file to my advisor (as a word file, since perhaps he would like to modify some parts). And he told me I was not following the standards, for example because the image labels did not start with “Figure x.”; I found this quite strange, since I was pleasantly surprised by the fact that the text was auto-inserted choosing the style (I think as a kind of numerated list style). I checked my file and the figure x. was still there. So I opened the document with word for Windows. And the Figure x. was missing.

Apparently Office 2008 for Mac and Office 2007 for Windows are not compatible at all. The question is: how many things are going to be fucked up? If I got a lot of nice features which I can’t use because otherwise I can’t share files with colleagues, then the whole point of using Word instead of Latex (and especially the “revisions” feature, which is very useful) is lost.

As a nice side note, I discovered the bibliography IEEE style I downloaded to “auto-manage” the bibliography is utterly broken. First, the references like [4] have italic numbers and the standard I have to follow has not… perhaps there are multiple standards, but then… how to discover the right file?

Second the bibliography itself is formatted as a huge double column table, which is completely unusable in 2-column mode. Other styles (such as the ones built-in with Word) do not format the references, so essentially applying the “references” style to the unformatted references simply solves the problem. But not with the style I have to use. And I don’t understand why word does not come with such widely used styles and I have to download broken third party implementations. The solution was to format the references by hand. Besides, Papers “import into Word” feature also gave me some problems, as every single link contained the pdf url, even if that is not standard. I don’t understand if it’s Papers fault, Word fault or the style fault. Well… doing things by hand, I solved the problem.

I still don’t understand the part “Word is user friendly”.

Thursday, June 24, 2010

Portable pleasure

This spring, I got convinced I needed a netbook. Essentially, I'm a satisfied Apple user (mostly). However, I needed a portable machine to work with. Work means coding. No code, no work. Moreover, I was curious of the new Windows 7 operating system and I also wanted to check Linux progresses as a desktop OS. I don't consider myself an Apple zealot: no Apple netbook, does not mean no netbook at all. I think Virtual Machines are a great tool. However, I don't really relate with virtual os's. I don't customize them. I just use them in order to perform some given tasks (e.g., testing). I thought that using Linux and Windows as primary operating systems when I only had the netbook was a great way to try them. I was pretty sure I could work with Linux (and I was confident I could overcome every installation issue that would have faced). I bought an Asus EEE 1005PE with 14h of autonomy declared. Other specs were completely uninteresting to me. CPU, RAM. I didn't give a fuck. If I need power, I've got a Mac Pro. If I need portable power, I have a Mac Book Pro. The EEE has plenty of HD space, and comes already partitioned (which means installing Linux is easy and you retain the EEE auto-reinstal clean windows capabilities). Should I need more ram, I could buy it. Ubuntu 10.4 is a pleasant surprise: I have never seen such a nice DE in a non Apple OS. The desktop is nice, the apps work great. NetworkManager does not suck anymore, some Ubuntu specific apps are very nice. It is rather pleasurable to work with (even if I don't use the netbook remix variant). Of course, some apps are completely lacking. However, I need it to code, write docs, email, browse... chat. Good. Windows was both a nice surprise and a disappointment. Compared to other MS OS's it is a nice surprise. It is fast, small, easy to use. Compared to what I read about Windows 7 it is a delusion: I think that Ubuntu is way nicer. However, I need some MS applications and that's it. With a little bit of work I installed everything I needed and I have got a nice working system. Unfortunately, Windows seems to slow down progressively and I don't understand why (I stopped installing stuff, I don't think that could be the cause). Hardware support is excellent both with Windows and with Ubuntu. Battery life lasts longer with Windows. I think about 12 hours of effective light work, wireless on. Ubuntu is probably around 9, which is sufficient for my goals.

Location:Via Montebello,Parma,Italy

Monday, June 21, 2010

The power of the power of the power...

Though I really don’t particularly like writing about hardware in my blog, I provide some updates on my current hardware situation. Some years ago, I wrote a short series of blog posts about my iMac. That was a first series machine and I criticized lack of support for scientific software. This is partly related to GMP not being too much collaborative (as it appears they did not have intel macs, bac in the day).

Of course those problems have been entirely solved quite soon. I suppose that by the end of the year they were no more and the faster intel chips were just the right choice for everyone. A couple of years after those events I bought a MacBook Pro 15” laptop and that became my main work machine. Among the things I really liked best about the intel processor (and all the lovely virtualization software packets out there) was that I could run Linux in a virtual machine. Since I did not stop working on a certain “almost-linux-only” project, that was really good (and saved me from ugly partitioning). Although I had a Parallels license (which I used on the iMac), I preferred to use the simpler VirtualBox. Back in the day, Parallels did not support Linux exceptionally well and I could just live with the free VB. And that is the machine I used (along with the iMac) for my M.Sc. dissertation.

4 years after i first turned it on, my iMac blew off. Well, not literally. I suppose I can fix it, however, I chose to buy a new Mac. This time it is a MacPro. Essentially, I hated that if I thrashed my iMac (because of other damages) I could not save the monitor. So I chose a monitor-less computer. Essentially the thing is powerful. It is extremely evident when encoding my CDs into mp3’s. I don’t do much video editing.

Nonetheless, Parallels is completely capable of sucking all the resources. I had to stop using it, since the MP became useless. I suppose it may be related to having only 3 GB, even though the iMac had only 2. If I will need a virtualization software, I suppose I’ll go with VMWare or VirtualBox. But… that is not the last machine I bought.

Glad I chose Blogger

This new command line thing is awesome.

Probably I would have regretted if I chose another blog service.
This also makes me think writing the crappy little blogging application for Linux should be pretty easy. :)

Sunday, June 13, 2010

Prolog Homoiconicity

Although when talking about homoiconicity the first language that springs into mind is Lisp (or a Lisp dialect), Prolog is a completely homoiconic language. In Prolog everything is a term. Variables, atoms, numbers are terms. And everything else is a compound term: e.g., foo(a, b) is a term. Lists are the main data structure of Prolog, and a list is a term as well. In fact a Prolog list is i) the empty list []; ii) a term '.'(a, b), where a is a term and b is a list or a variable. We say that a compound term foo(x, y) has functor foo and arity 2. Clauses (the prolog equivalent of "functions") are represented as terms whose main functor is :-, their head is a prolog term and their body is a prolog term as well. Moreover, the read predicate reads from file (or standard input) some prolog predicates. E.g., this simple Prolog program reads a Prolog program into memory (and the program is represented as a list of prolog terms):
read_all([Predicate|T]) :-
	read(Predicate),
	Predicate \== end_of_file, !,
	read_all(T).
read_all([]).

read_file(File, L) :-
	see(File), read_all(L), seen.
If we call read_file to read "read_file.pl" (spaces added by me):
?- read_file('read_all.pl', L).
L = [(read_all([_A|_B]):-read(_A),_A\==end_of_file,!,read_all(_B)),
     read_all([]),
     (read_file(_C,_D):-see(_C),read_all(_D),seen)] ? y
The display predicate some prolog interpreters provide gives us more insight on internal representation of terms. We call it on just one rule in order not to clutter the output with list representations:
?- read_file('read_all.pl', [H|_]), display(H).
:-(read_all(.(_932,_952)),,(read(_932),,(\==(_932,end_of_file),,(!,read_all(_952)))))
And here on the full program (notice... lists are "consed" with the '.' functor):
?- read_file('read_all.pl', L), display(L).
.(:-(read_all(.(_900,_920)),,(read(_900),,(\==(_900,end_of_file),,(!,read_all(_920))))),.(read_all([]),.(:-(read_file(_1410,_1430),,(see(_1410),,(read_all(_1430),seen))),[])))
Indeed, infix operators have been placed in prefix form like every other predicate (indeed, they are not special, just a parsing trick). Another "useful" predicate (not standard) is portray_clause/1, in order to do some "pretty-printing":
?- read_file('read_all.pl', [H|_]), portray_clause(H).
read_all([A|B]) :-
        read(A),
        A\==end_of_file, !,
        read_all(B).

The assert family of predicates can “store” into the Prolog database facts (this is very common, in fact, memoization is often implemented this way) and rules (this is done less frequently – and there may be a slight efficiency penality [CHECK]). Thus we can manipulate terms (indeed that is what Prolog is about) and load them as programs. Notice that I showed the read predicate only to show that “standard” Prolog programs are represented that way, but it is not necessary at all to use programs saved in files.

4 ?- assert((foo(X) :- bar(X), \+ baz(X))).
true.

5 ?- assert(bar(1)), assert(bar(2)), assert(baz(2)).
true.

6 ?- foo(X).
X = 1 .