The Uncarved Block

Haskell arrays are amazing

Functional languages are all about lists... but Haskell has incredible arrays

For certain classes of problems where speed is of the essence, table lookups can be a fantastic solution. However, many functional programmers eschew these approaches partly because they are thinking in terms of lists, and lookups in lists (as we all know) are a bit rubbish most of the time. However, Haskell has a fantastic array implementation that is flexible enough to put many imperative languages to shame (especially when combined with list comprehensions). Hopefully an example will illustrate what I mean. The actual array lookup I'm doing is into arrays of much larger size, but the principle is exactly the same.

Say you want to model a regular deck of 52 playing cards and want to write functions to convert them to and from integers. We want to proceed by rank (two to ace), then for each rank, by suit (clubs, diamonds, hearts, spades) so the two of clubs is going to be 0, the two of spades is going to be 3 and the ace of spades is going to be 51.

Let's start by defining the datatypes for Rank and Suit and a simple generic card type:

data Rank =
    Two
    | Three
    | Four
    | Five
    | Six
    | Seven
    | Eight
    | Nine
    | Ten
    | Jack
    | Queen
    | King
    | Ace
    deriving (Eq, Ord, Show, Read, Enum, Ix)

data Suit =
    Clubs
    | Diamonds
    | Hearts
    | Spades
    deriving (Eq, Ord, Show, Read, Enum, Ix)

data GenCard = GenCard Rank Suit
    deriving (Eq, Ord, Show, Read)

...so far so not very interesting. Now the two functions we want to define are:

genCardOfInt :: Int -> GenCard
genCardToInt :: GenCard -> Int

...and obviously, you could do:

genCardOfInt 0 = GenCard Two Clubs
genCardOfInt 1 = GenCard Two Diamonds

...etc all the way to...

genCardOfInt 51 = GenCard Ace Spades

...and then...

genCardToInt (GenCard Two Clubs) = 0

...and so on.

This would work, but it is extremely smelly. The thing that alerts us to the fact that this is fishy is that there is a lot of repetitive typing. The gods of Haskell frown on this and generally if you're doing a bunch of it, there's probably something wrong. And indeed there is. This approach leads to a linear search through cases every time until we find a match. If we care about performance, this will be horrible and imagine how bad it will get if we get a table with (say) 4 million entries? The actual problem I am interested in requires this size of table. Happily we can use an array-driven method to not only speed up our algorithm, but also to scrap all this tedious repitition. First, some background.

For genCardOfInt what we want to do is create an array of 52 cards and then use the int passed in to look up into this array. This will give us constant time access. So what we want is:

genCardOfInt x = lookup ! x
    where
        lookup = .... to be discussed

One possible solution is to do

        lookup = listArray (0,51) [GenCard Two Clubs, GenCard Two Diamonds etc etc etc]

... but that's almost as unhaskelly as what we had before! So, here's how we do it.

genCardOfInt x = lookup ! x
    where
        lookup = listArray (0,51) [GenCard r s|r<-enumFrom Two, s<-enumFrom Clubs]

Say whaaaa? Well, let's see what ghci has to say:

*Card> enumFrom Two
[Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten,Jack,Queen,King,Ace]
*Card> enumFrom Clubs
[Clubs,Diamonds,Hearts,Spades]

So because we said "deriving Enum" for our types, we get this ability for free. Mighty 'andy.

*Card> [GenCard r s|r<-enumFrom Two, s<-enumFrom Clubs]
[GenCard Two Clubs,GenCard Two Diamonds,GenCard Two Hearts,GenCard Two Spades,GenCard Three Clubs,GenCard Three Diamonds,GenCard Three Hearts,GenCard Three Spades,GenCard Four Clubs,GenCard Four Diamonds,GenCard Four Hearts,GenCard Four Spades,GenCard Five Clubs,GenCard Five Diamonds,GenCard Five Hearts,GenCard Five Spades,GenCard Six Clubs,GenCard Six Diamonds,GenCard Six Hearts,GenCard Six Spades,GenCard Seven Clubs,GenCard Seven Diamonds,GenCard Seven Hearts,GenCard Seven Spades,GenCard Eight Clubs,GenCard Eight Diamonds,GenCard Eight Hearts,GenCard Eight Spades,GenCard Nine Clubs,GenCard Nine Diamonds,GenCard Nine Hearts,GenCard Nine Spades,GenCard Ten Clubs,GenCard Ten Diamonds,GenCard Ten Hearts,GenCard Ten Spades,GenCard Jack Clubs,GenCard Jack Diamonds,GenCard Jack Hearts,GenCard Jack Spades,GenCard Queen Clubs,GenCard Queen Diamonds,GenCard Queen Hearts,GenCard Queen Spades,GenCard King Clubs,GenCard King Diamonds,GenCard King Hearts,GenCard King Spades,GenCard Ace Clubs,GenCard Ace Diamonds,GenCard Ace Hearts,GenCard Ace Spades]

So that bit is a list comprehension. Basically it's a simple way of defining a list and avoids all the tedium and boilerplate. And "listArray" takes some dimensions and a list and returns an array. But the best is yet to come, and ghci hints at it now:

*Card> :t listArray (0,51) [GenCard r s|r<-enumFrom Two, s<-enumFrom Clubs]
listArray (0,51) [GenCard r s|r<-enumFrom Two, s<-enumFrom Clubs]
  :: (Num t, Ix t) => Array t GenCard

So the type of that "lookup" variable is an Array t GenCard. The "t" is some numeric type that we can use to index into our array, and the GenCards (as we know) are what's in the array. So the type of the array index is polymorphic.

This means we can do:

genCardToInt :: GenCard -> Int
genCardToInt (GenCard r s) = lookup ! (r,s)
    where
        lookup = listArray ((Two,Clubs),(Ace,Spades)) [x|x<-[0..51]]

Gosh! Let's try that expression in ghci:

*Card>  listArray ((Two,Clubs),(Ace,Spades)) [x|x<-[0..51]]
array ((Two,Clubs),(Ace,Spades)) [((Two,Clubs),0),((Two,Diamonds),1),((Two,Hearts),2),((Two,Spades),3),((Three,Clubs),4),((Three,Diamonds),5),((Three,Hearts),6),((Three,Spades),7),((Four,Clubs),8),((Four,Diamonds),9),((Four,Hearts),10),((Four,Spades),11),((Five,Clubs),12),((Five,Diamonds),13),((Five,Hearts),14),((Five,Spades),15),((Six,Clubs),16),((Six,Diamonds),17),((Six,Hearts),18),((Six,Spades),19),((Seven,Clubs),20),((Seven,Diamonds),21),((Seven,Hearts),22),((Seven,Spades),23),((Eight,Clubs),24),((Eight,Diamonds),25),((Eight,Hearts),26),((Eight,Spades),27),((Nine,Clubs),28),((Nine,Diamonds),29),((Nine,Hearts),30),((Nine,Spades),31),((Ten,Clubs),32),((Ten,Diamonds),33),((Ten,Hearts),34),((Ten,Spades),35),((Jack,Clubs),36),((Jack,Diamonds),37),((Jack,Hearts),38),((Jack,Spades),39),((Queen,Clubs),40),((Queen,Diamonds),41),((Queen,Hearts),42),((Queen,Spades),43),((King,Clubs),44),((King,Diamonds),45),((King,Hearts),46),((King,Spades),47),((Ace,Clubs),48),((Ace,Diamonds),49),((Ace,Hearts),50),((Ace,Spades),51)]

So we have an array which is indexed by a pair of (Rank, Suit), and the values are numbers from 0 to 51. And no boilerplate. If we wanted to test this, of course, we could write some quickcheck properties that verify this implementation against the naive one I gave before. This sort of model-based testing is going to be essential to verify the complex but fast implementation I have for my real problem against the obvious but insanely verbose and slow implementation that I can put in my test.

Aren't Haskell arrays amazing though? I was genuinely stunned when I realised I could do table-driven methods so elegantly. There's a nice tutorial to Haskell arrays in the "gentle introduction"

Using a Linux laptop in text mode

Linux has everything you need to use your laptop even if you choose not to use the pointy-clicky interfaces

So having recently got a netbook, I have been learning how to do everything on it. It's obviously all easy if you use something like xfce or gnome, or some netbook-specific spin like the Ubuntu netbook remix, but I like to use xmonad, so it helps to learn how to do all the laptoppy things in text mode.

Here's a list of things it's useful to know:

nm-applet - finds a network and connects to it automagically
pm-suspend - sends your laptop into a light sleep
pm-hibernate - sends your laptop into a heavy sleep where it uses no battery

I've got a new netbook

So finally I give in to temptation...

This week I bought a new netbook (an Acer Aspire One Pro P531). I was determined to reward someone who sells a Linux netbook but in spite of my best efforts noone was selling the sort of configuration I wanted (big SSD drive, 2GB memory) without Windows XP. So I bought one and installed Fedora on it myself.

It's a fantasticlittle device and everything on it just works with Fedora. Apart from the ethernet card which seems not to want to dhcp (although that might be some incompatibility with my home hub/router/adsl thingy).

I'm using it with Xmonad to do random bits of Haskell tinkering when I take the tube to work in the morning. The keyboard is a bit fiddly, but manageable, and the small form factor and light weight are terrific. In this respect it reminds me a little of my old sony vaio which was a sort of expensive predecessor to a netbook that I got second hand.

One thing I would have done differently if I had done more research is to get the other sort of Atom processor. Ones with an "N" or a "Z" at the front are based on the i686 architecture, whereas ones without that are x86_64 based. Since I have an x86_64 desktop pc at home it would have been slightly more convenient not to have to download two seperate fedora versions. The upside of this one is I think it has slightly better battery life. In fact, battery life seems amazing.

n-Bit Gray codes in Haskell

A first step in what will become a combinatorics library

I have been playing around with Gray's reflected binary code (aka Gray codes) and similar things a bit. Before I reveal why I'm doing this lets just dive in and write some code. Gray's algorithm is described well here. The code which follows is in haskell, because it's a really fantastic language and I'm playing around with it at the moment. For scala fans, don't worry. I haven't abandoned scala, this is a parallel effort.

So for starters we need a datatype for representing these things. This is how you define an algebraic datatype in haskell. In what follows, lines beginning "--" are single-line comments

-- | 'Bit' is the datatype for representing a bit in a gray code.  
data Bit = Zero | One deriving Show

Alright. So we have a type "Bit" with two constructors Zero and One and a "deriving Show" which means haskell figures out how to turn it into a string. This is useful when you're in ghci (the interactive haskell environment) debugging.

-- prepend a given item onto each of a list of lists (probably something to do this in the prelude)
prepend :: a -> [[a]] -> [[a]]
prepend t xs = map (t:) xs

A teeny helper function. Given a list of lists and a thing it sticks the thing on the front of each list in the outer list. This would append the thing on the end of each list:

append :: a -> [[a]] -> [[a]]
append t xs = map (++[t]) xs

Note I'm writing the type signatures explicitly but there's absolutely no problem if you leave them off. So let's generate our Gray codes:

-- | 'gray' generates the gray code sequence of length 'n'
gray :: Int -> [[Bit]]
gray 1 = [ [Zero], [One] ]
gray n = prepend Zero (gray (n-1)) ++ prepend One (descGray (n-1))

-- | 'descgray' generates the reversed gray code sequence of length 'n'
descGray :: Int -> [[Bit]]
descGray 1 = [ [One], [Zero] ]
descGray n = prepend One (gray (n-1)) ++ prepend Zero (descGray (n-1))

So we get an ascending and a descending one for free. Since the descending one is just the ascending one in reverse why (you might say) don't I just define descGray as descGray = reverse.gray ? Indeed, that may be a reasonable thing to do. I'm doing it this way to try to preserve as much laziness as possible, and (although my haskell-fu is still very weak at the moment) I think that if you reverse a list you pretty much have to evaluate each thing in the list. If you read the paper you'll see that this is Gray's (naive) algorithm and there has been an astonishing amount of research in this area leading to more efficient algorithms. I'll give those a crack at some point.

Why am I doing this? You'll see. This is at the heart of building a really cool combinatorics library. I needed something that could enumerate all combinations and permutations of various generic distribution-type things. There are similar but more recent orderings that are comparable to gray codes which I'm also looking into. They'll all be presented here in due course.

Which language would you use?

It depends

I got a mail a few days ago about how I wrote a couple of articles with code for an options pricer in ocaml but all my latest articles were about scala. So which language would I use today if I were to write an options pricer. The answer of course is "it depends".

Now in and of themselves both ocaml and scala are fine choices for writing just about anything. But there are tradeoffs. If I was in charge of development of a brand new pricing and risk infrastructure for a big bank that had to be able to price everything from a stock to a digital multi-asset multi-barrier callable bermudan range-accrual thingummybobber that was going to be worked on by a thousand people ranging from the lowliest intern to the most brilliant genius then I have no hesitation in saying I would use scala.

In fact a friend of mine who is in charge of the development of risk and pricing systems at a major wall st firm told me that if he was building this infrastructure far a bank from scratch today he would use scala. He's the guy who persuaded me to try scala in the first place as a matter of fact.

The reasons this would be a fine choice should be obvious- it's very simple to write serious software in scala, it shouldn't take anyone of reasonable ability much time at all to learn, the syntax is not overly burdensome and tedious, performance is adequate or better for most things, the concurrency paradigm is tractable by normal human beings and there is fantastic library support because you can just use java stuff. The extensible syntax would help for various things and the cross-platform support is always a nice thing to have if you want to have number-crunching on a Linux compute farm and desktop apps on Windows or whatever.

On the other hand if I was setting up a small quant trading shop/hedge fund or doing it for my own benefit then the choices are much more varied. I might use scala still (it would still be an excellent choice), I might do it in ocaml (or even Haskell in fact), particularly if I was going to be all or most of the programming myself or I had access to recruit a decent pipeline of smart functional-programming aware people.

The benefits of doing it in ocaml (or Haskell) would be that I would probably have a more mentally-stimulating time doing it (this is can be an important motivation also if you have a super-bright team), and would probably end up with something more aesthetically pleasing from a pure comp-sci point of view.

On the other hand I would certainly have more frustrations (eg Why has no-one noticed that you can only do one request through the ocaml curl library because there is a memory scribble? But I digress). I wouldn't really want to lead a group of 100 guys and have to keep teaching haskell monad combinators or whatever every time a new person joined. And maintaining/code reviewing etc could become excruciating when you were dealing with people of average ability less one or two standard deviations.

So horses for courses. Ultimately writing good software always requires thought, discipline and some skill. The right language fits the problem domain and suits the group and organization. Good programmers can learn to write good software in any language.

Website Changes

An occasionally-updated summary of stuff I've changed

So when I had a day off work this week I decided to fix a few things on this site that have been bugging me for a while:

most of the content is dynamic, yet the sitemap (used by search engines to figure out how to slurp your site) was not
I never bothered to put proper meta description tags on pages
now that there were quite a few articles things were getting lost in painful navigation
The site didn't do gzip encoding
I've finally added comments thanks to intense debate

So I've had a go at fixing them. I now generate my sitemaps automatically using the code I wrote to generate the atom and rss feeds, transparently gzip stuff if your client can accept that, and there is a new "treeview" page on the side to make things easier to get to. Oh, and the pages have proper description tags which should make them more useful on search engines.

Hope this all helps.

Adding comments

So I've thought about doing this for some time but never had the spare time. Finally I add comments.

I have added comments for the first time. I don't expect there will be a great deal of traffic and at the moment just to get started I have moderation on. If things seem to be resonably calm and the wingnuts stay away, I'll turn moderation off.

I'm using the comment system from intense debate, which has been fantastically easy to set up and use so far. Let me know what you think.

Thank goodness for greasemonkey

Google recently added an extremely annoying feature to google reader. Thankfully you can decide for yourself how their website works.

Every morning I drink my coffee and read various news sites using google reader, so I was horrified recently to see a big yellow smiley popping up on a bunch of posts. This is some sort of social-networking-inspired feature which allows people to say which posts they "like", then anyone else reading that post gets a smiley to say "6 people liked this post" or whatever.

Perhaps this is just the ticket if you're desperately insecure about whether or not your taste in news is in tune with the mainstream of the internets, but as someone who does not give a flying... err monkeys about what other people think (by and large) I hate this feature. In fact, clicking through to the google forum on this feature I find that it turns out there are a lot of people who do not "like" reader's new "like" feature.

Thanks to the miracle of greasemonkey I can decide for myself what google's site looks (and works) like in my browser, and in fact Jason Fager has already posted a greasemonkey script that removes the offending smileys.

helper classes for scala

My first opensource scala package

So I've released my helpers code via sbaz (the package is uncarved-helpers), or you can just get a tarball. Scaladoc is available here, although I warn you that my patience with hacking the stylesheet only goes so far and I don't yet know how to make scaladoc that is not foul to look at. To build it, and run the tests, do...

ant build run_tests

...after unpacking the tarball somewhere

The motivation behind this package is that I wanted to do some http stuff and wanted to make my client code be as simple as possible while still being polite to servers. So that meant supporting conditional get using Etags, Last-Modified and also supporting gzip encoding. I do all this by using the apache commons httpclient-4.x and httpcore-4.x libraries and wrapping them all up in a class that's convenient and simple to use. Here's a taster:

 import com.uncarved.helpers.http._

 val helper = new BasicClient()
 //Get a webpage as a string (if you look at the apache http log4j messages you can see it
 //does conditional get and has transparent gzip support)
 val str = helper.get("http://www.theflautadors.org/")

 //Get some XML (with request parameters supplied)
 val params = List(("tag"->""), ("limit"->"5"))
 val xml = helper.getXML(Request(RequestType.GET, "http://www.uncarved.com/index.py/rss1.1.xml", params))
 val items = xml \\ "item"

Not Memoization

Sometimes you need a better algorithm, rather than papering over the cracks with optimisation

Sometimes memoization is a useful technique for optimisation, but when you see it explained,it is often demonstrated using little explanatory code or algorithms that show the technique itself well, but are not the best way of achieving their supposed purpose. So you may see a small Fibonacci function for example, that demonstrates memoization well, but is a bad way of calculating Fibonacci numbers.

The same is true of recursion. You will often see something like this factorial function used to demonstrate recursion:

def fact(n: Int) : BigInt = {
    assert(n>=0, "factorial is only defined over non-negative integers")
    if(n==0) then 1
    else fact(n-1) * n
}

...and the article will then go on to demonstrate how to get the tail call optimisation to work so you don't run out of stack space. Or alternatively, will use this as an example of how you memoize a recursive function. The fact is that this is just a really bad way to write a factorial function. Here's one that is not bad:

def fact(n: Int) = {
    assert(n>=0, "factorial is only defined over non-negative integers")
    (BigInt(1) /: (1 to n)) (_*_)
}

Note that we don't need to rely on any special optimisations and this solution runs in constant space and in O(n) time. All we are doing is folding over the range from 1 to the n and multiplying as we go. To put it another way, we are going from the bottom up. Amusingly when we go from the bottom up, we get a function application that looks like butt-cheeks (I can't be the only person to have noticed that).

Now it doesn't take a genius to figure out how we use a memoization-like technique to optimise this if we want to. (I know that it can't because I thought of how to optimise it and I'm no genius)

object Factorial {
    import scala.collection.mutable.ArrayBuffer

    private val cache : ArrayBuffer[BigInt] = new ArrayBuffer()
    cache += 1
    cache += 1

    def apply(n: Int) = {
        assert(n>=0, "Factorial is only defined over non-negative integers")
        if(n>cache.length-1) {
            for(idx<-cache.length to n)
                cache += cache(idx-1) * idx
        }

        cache(n)
    }
}

So we make an array and store the intermediate results in it so that if we are called later for a value that we have already computed, we can simply return it from the array. Because we are using an array, the lookup happens in constant time with little or no constant overhead, rather than using a map or hash as is usually required for general memoization. We can do this because we know the space of inputs is not sparse - we need to compute every intermediate value anyway so we may as well cache them.

This is not memoization because what we are caching is not (just) the result of the function invocation but all partial results. As such, our optimisation has more in common with techniques used in dynamic programming.

The point is that writing the best algorithm is almost always the best optimisation you can do, and if you're jumping through hoops to optimise something, it may be because there's a better way to do it. In the words of Nadia Boulanger, "Never strain to avoid the obvious".

P.S. If you're computing Fibonacci numbers recursively, you are living in a state of sin. Please use the closed form solution. Thanks.